text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.NonlinearEquations (ModiaMath/NonlinearEquations/_module.jl)
#

"""
    y = solveOneNonlinearEquation(f, u_min, u_max; u_nominal=1.0, tolerance=100.0*eps())

This function determines the solution of one non-linear algebraic equation `y=f(u)`
in one unknown `u` in a reliable and efficient way. It is one of the best
numerical algorithms for this purpose. As input, the nonlinear function `f(u)`
has to be given, as well as an interval `u_min`, `u_max` that contains the solution,
that is `f(u_min)` and `f(u_max)` must have a different sign.
The function computes iteratively smaller intervals
in which a sign change is present and terminates when the following test
is fulfilled (`u_nominal` is the nominal value of `u`, so the order of magnitude of `u`
around the solution and `tolerance` is the relative tolerance;
for example, tolerance=1e-10 means that the solution is correct up to 10 digits):

```
absTol = 0.1*u_nominal*tolerance
if abs(length_of_u_interval) <= max(tolerance*abs(u), absTol) || f(u) == 0.0
    # root u found (interval is small enough)
```

The interval reduction is performed
using inverse quadratic interpolation (interpolating with a quadratic polynomial
through the last three points and computing the zero). If this fails, bisection
is used, which always reduces the interval by a factor of two. The inverse quadratic
interpolation method has superlinear convergence. This is roughly the same convergence
rate as a globally convergent Newton method, but without the need to compute derivatives
of the non-linear function. The solver function is a direct mapping of the Algol 60 procedure
`zero` to Julia, from:

- Brent R.P. (1973): *Algorithms for Minimization without derivatives.*
  Prentice Hall, 1973, pp. 58-59.\\
  Download: [http://wwwmaths.anu.edu.au/~brent/pd/rpb011i.pdf](http://wwwmaths.anu.edu.au/~brent/pd/rpb011i.pdf)\\
  Errata and new print: [http://wwwmaths.anu.edu.au/~brent/pub/pub011.html](http://wwwmaths.anu.edu.au/~brent/pub/pub011.html)

# Examples
```julia
import ModiaMath: solveOneNonlinearEquation

fun(u; w=1.0) = 3*u - sin(w*u) - 1
u_zero = solveOneNonlinearEquation(u->fun(u; w=3.0), 0.0, 5.0)
println("residue = ", fun(u_zero; w=3.0))
```

# Remarks

- The interface was made such that it is identical to function
  [Modelica.Math.Nonlinear.solveOneNonlinearEquation](https://doc.modelica.org/Modelica%203.2.3/Resources/helpDymola/Modelica_Math_Nonlinear.html#Modelica.Math.Nonlinear.solveOneNonlinearEquation)
  in order that automatic translation of Modelica to Modia is simplified.
  However, the termination condition was changed:
  The original Brent algorithm uses an absolute tolerance for termination
  (abs(length_of_interval) <= 2*eps()*abs(u) + tolerance || f(u) == 0.0), whereas
  this was changed here to take a relative tolerance into account and use a similar
  definition of tolerances as used by Modia and Modelica,
  because easier to understand by a user (with relative tolerance and nominal value, instead of
  relative and absolute tolerances).

- Newer algorithms for the problem are presented in
  *Alefeld, Potra, Shi (1995):* [Algorithm 748: enclosing zeros of continuous functions](https://dl.acm.org/citation.cfm?id=210111).
  Here, an inverse cubic interpolation is used instead of an inverse quadratic interpolation
  as in case of the Brent algorithm. The numerical experiments in this article with
  15 test problems show that Brents algorithm needs about 5% more function evaluations
  for all the 15 test problems as the newly presented algorithm 4.2
  (in some cases Brents algorithm needs slightly less and in other cases
  slightly more function evaluations).

- Julia package [Roots.jl](https://github.com/JuliaMath/Roots.jl) provides various algorithms
  to solve the problem above but it seems that Brents algorithm is not yet included
  (as of Nov. 2018).
```
"""
function solveOneNonlinearEquation(f::Function, u_min, u_max; u_nominal=1.0, tolerance=100.0*eps())
  @assert(u_nominal > 0.0)
  @assert(tolerance > 0.0)
  @assert(u_max > u_min)
  u_min1::Float64     = convert(Float64, u_min)
  u_max1::Float64     = convert(Float64, u_max)
  u_nominal1::Float64 = convert(Float64, u_nominal)
  tolerance1::Float64 = convert(Float64, tolerance)
  absTol::Float64     = 0.1*u_nominal*tolerance1
  u::Float64 = 0.0
  a::Float64 = u_min1   # Current best minimum interval value
  b::Float64 = u_max1   # Current best maximum interval value
  c::Float64 = 0.0     # Intermediate point a <= c <= b"
  d::Float64 = 0.0
  e::Float64 = 0.0     # "b - a"
  m::Float64 = 0.0
  s::Float64 = 0.0
  p::Float64 = 0.0
  q::Float64 = 0.0
  r::Float64 = 0.0
  tol::Float64 = 0.0
  fa::Float64 = 0.0    # = f(a)
  fb::Float64 = 0.0    # = f(b)
  fc::Float64 = 0.0
  found::Bool = false

  # Check that f(u_min1) and f(u_max1) have different sign
  fa = f(u_min1)
  fb = f(u_max1)
  fc = fb
  if fa > 0.0 && fb > 0.0 || fa < 0.0 && fb < 0.0
    error("\nThe arguments u_min and u_max provided in the function call\n",
          "    solveOneNonlinearEquation(f,u_min,u_max)\n",
          "do not bracket the root of the single non-linear equation 0=f(u):\n",
          "  u_min = ", string(u_min1), "\n",
          "  u_max = ", string(u_max1), "\n",
          "  fa = f(u_min) = ", string(fa), "\n",
          "  fb = f(u_max) = ", string(fb), "\n",
          "fa and fb must have opposite sign which is not the case here.\n")
  end

  # Initialize variables
  c  = a
  fc = fa
  e  = b - a
  d  = e

  # Search
  while !found
    if abs(fc) < abs(fb)
      a  = b
      b  = c
      c  = a
      fa = fb
      fb = fc
      fc = fa
    end

    # Original formulation in Brents algorithm:
    # tol = 2*eps()*abs(b) + tolerance1
    tol = max(tolerance1*abs(b), absTol)
    m = (c - b)/2

    if abs(m) <= tol || fb == 0.0
      # root found (interval is small enough)
      found = true
      u = b
    else
      # Determine if a bisection is needed
      if abs(e) < tol || abs(fa) <= abs(fb)
        e = m
        d = e
      else
        s = fb/fa
        if a == c
          # linear interpolation
          p = 2*m*s
          q = 1 - s
        else
          # inverse quadratic interpolation
          q = fa/fc
          r = fb/fc
          p = s*(2*m*q*(q - r) - (b - a)*(r - 1))
          q = (q - 1)*(r - 1)*(s - 1)
        end

        if p > 0
          q = -q
        else
          p = -p
        end

        s = e
        e = d
        if 2*p < 3*m*q - abs(tol*q) && p < abs(0.5*s*q)
          # interpolation successful
          d = p/q
        else
          # use bi-section
          e = m
          d = e
        end
      end

      # Best guess value is defined as "a"
      a  = b
      fa = fb
      b  = b + (abs(d) > tol ? d : (m > 0 ? tol : -tol))
      fb = f(b)

      if fb > 0 && fc > 0 || fb < 0 && fc < 0
        # initialize variables
        c  = a
        fc = fa
        e  = b - a
        d  = e
      end
    end
  end

  return u
end


text:
arrayof(::Type{S}, d::NTuple{N, Any}) where {S, N} = similar(defaultarray(S, N), d)

struct StructArrayInitializer{F, G}
    unwrap::F
    arrayof::G
end
StructArrayInitializer(unwrap = t->false) = StructArrayInitializer(unwrap, arrayof)

const default_initializer = StructArrayInitializer()

function (s::StructArrayInitializer)(S, d)
    ai = ArrayInitializer(s.unwrap, s.arrayof)
    buildfromschema(typ -> ai(typ, d), S)
end

struct ArrayInitializer{F, G}
    unwrap::F
    arrayof::G
end
ArrayInitializer(unwrap = t->false) = ArrayInitializer(unwrap, arrayof)

(s::ArrayInitializer)(S, d) = s.unwrap(S) ? buildfromschema(typ -> s(typ, d), S) : s.arrayof(S, d)

_axes(itr) = _axes(itr, Base.IteratorSize(itr))
_axes(itr, ::Base.SizeUnknown) = nothing
_axes(itr, ::Base.HasLength) = (Base.OneTo(length(itr)),)
_axes(itr, ::Base.HasShape) = axes(itr)

"""
    collect_structarray(itr; initializer = default_initializer)

Collects `itr` into a `StructArray`. The user can optionally pass a `initializer`, that is to say
a function `(S, d) -> v` that associates to a type and a size an array of eltype `S`
and size `d`. By default `initializer` returns a `StructArray` of `Array` but custom array types
may be used.
"""
function collect_structarray(itr; initializer = default_initializer)
    ax = _axes(itr)
    elem = iterate(itr)
    _collect_structarray(itr, elem, ax; initializer = initializer)
end

function _collect_structarray(itr::T, ::Nothing, ax; initializer = default_initializer) where {T}
    S = Core.Compiler.return_type(first, Tuple{T})
    return initializer(S, something(ax, (Base.OneTo(0),)))
end

function _collect_structarray(itr, elem, ax; initializer = default_initializer)
    el, st = elem
    S = typeof(el)
    dest = initializer(S, something(ax, (Base.OneTo(1),)))
    offs = first(LinearIndices(dest))
    @inbounds dest[offs] = el
    return _collect_structarray!(dest, itr, st, ax)
end

function _collect_structarray!(dest, itr, st, ax)
    offs = first(LinearIndices(dest)) + 1
    return collect_to_structarray!(dest, itr, offs, st)
end

_collect_structarray!(dest, itr, st, ::Nothing) =
    grow_to_structarray!(dest, itr, iterate(itr, st))

function collect_to_structarray!(dest::AbstractArray, itr, offs, st)
    # collect to dest array, checking the type of each result. if a result does not
    # match, widen_from_instance the result type and re-dispatch.
    i = offs
    while true
        elem = iterate(itr, st)
        elem === nothing && break
        el, st = elem
        if iscompatible(el, dest)
            @inbounds dest[i] = el
            i += 1
        else
            new = widen_from_instance(dest, i, el)
            @inbounds new[i] = el
            return collect_to_structarray!(new, itr, i+1, st)
        end
    end
    return dest
end

function grow_to_structarray!(dest::AbstractArray, itr, elem = iterate(itr))
    # collect to dest array, checking the type of each result. if a result does not
    # match, widen_from_instance the result type and re-dispatch.
    i = length(dest)+1
    while elem !== nothing
        el, st = elem
        if iscompatible(el, dest)
            push!(dest, el)
            elem = iterate(itr, st)
            i += 1
        else
            new = widen_from_instance(dest, i, el)
            push!(new, el)
            return grow_to_structarray!(new, itr, iterate(itr, st))
        end
    end
    return dest
end

# Widen `dest` to contain `el` and copy until index `i-1`
widen_from_instance(dest::AbstractArray, i, el::T) where {T} = widen_from_type(dest, i, T)
# Widen `dest` to contain elements of type `T` and copy until index `i-1`
function widen_from_type(dest::AbstractArray{S}, i, ::Type{T}) where {S, T}
    U = _promote_typejoin(S, T)
    return _widenstructarray(dest, i, U)
end

function _widenstructarray(dest::StructArray, i, ::Type{T}) where {T}
    sch = hasfields(T) ? staticschema(T) : nothing
    sch !== nothing && fieldnames(sch) == propertynames(dest) || return _widenarray(dest, i, T)
    types = ntuple(x -> fieldtype(sch, x), fieldcount(sch))
    cols = Tuple(fieldarrays(dest))
    newcols = map((a, b) -> _widenstructarray(a, i, b), cols, types)
    return StructArray{T}(newcols)
end

_widenstructarray(dest::AbstractArray, i, ::Type{T}) where {T} = _widenarray(dest, i, T)

_widenarray(dest::AbstractArray{T}, i, ::Type{T}) where {T} = dest
function _widenarray(dest::AbstractArray, i, ::Type{T}) where T
    new = similar(dest, T, length(dest))
    copyto!(new, firstindex(new), dest, firstindex(dest), i-1)
    new
end

"""
    dest = StructArrays.append!!(dest, itr)

Try to append `itr` into a vector `dest`, widening the element type of `dest` if
it cannot hold the elements of `itr`. That is to say,
```julia
vcat(dest, StructVector(itr)) == append!!(dest, itr)
```
holds. Note that the `dest` argument may or may not be the same object as the
returned value.

The state of `dest` is unpredictable after `append!!` is called (e.g., it may
contain some, none or all the elements from `itr`).
"""
append!!(dest::AbstractVector, itr) =
    _append!!(dest, itr, Base.IteratorSize(itr))

function _append!!(dest::AbstractVector, itr, ::Union{Base.HasShape, Base.HasLength})
    n = length(itr)  # itr may be stateful so do this first
    fr = iterate(itr)
    fr === nothing && return dest
    el, st = fr
    i = lastindex(dest) + 1
    new = iscompatible(el, dest) ? dest : widen_from_instance(dest, i, el)
    resize!(new, length(dest) + n)
    @inbounds new[i] = el
    return collect_to_structarray!(new, itr, i + 1, st)
end

_append!!(dest::AbstractVector, itr, ::Base.SizeUnknown) =
    grow_to_structarray!(dest, itr)

# Optimized version when element collection is an `AbstractVector`
# This only works for julia 1.3 or greater, which has `append!` for `AbstractVector`
@static if VERSION ≥ v"1.3.0"
    function append!!(dest::V, v::AbstractVector{T}) where {V<:AbstractVector, T}
        new = iscompatible(T, V) ? dest : widen_from_type(dest, length(dest) + 1, T)
        return append!(new, v)
    end
end

ext: jl


text:
"""
    Masked{K<:Kernel} <: Kernel

A wrapper for kernels so that they are only applied along certain dimensions.

This is similar to the `active_dims` kernel attribute in the python GPy package and to the
`covMask` function in the matlab gpml package.

The implementation is very simple: any function of the kernel that takes an `X::Matrix`
input is delegated to the wrapped kernel along with a view of `X` that only includes the
active dimensions.
"""
struct Masked{K<:Kernel, DIM} <: Kernel
    kernel::K
    active_dims::SVector{DIM, Int}
end
num_dims(masked::Masked{K, DIM}) where {K, DIM} = DIM
function Masked(kern::Kernel, active_dims::AbstractVector)
    dim = length(active_dims)
    K = typeof(kern)
    return Masked{K,dim}(kern, SVector{dim, Int}(active_dims))
end
xview(masked::Masked, X::AbstractMatrix) = view(X, masked.active_dims, :)

################
## KernelData ##
################
struct MaskedData{M1<:AbstractMatrix,M2<:AbstractMatrix,KD<:KernelData} <: KernelData
    X1view::M1
    X2view::M2
    wrappeddata::KD
end
function KernelData(masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix)
    X1view = xview(masked,X1)
    if X1 === X2
        X2view = X1view
    else
        X2view = xview(masked,X2)
    end
	wrappeddata = KernelData(masked.kernel, X1view, X2view)
    return MaskedData(X1view, X2view, wrappeddata)
end

@inline @inbounds function cov_ij(masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::MaskedData, i::Int, j::Int, dim::Int)
    return cov_ij(masked.kernel, data.X1view, data.X2view, data.wrappeddata, i, j, num_dims(masked))
end
@inline @inbounds function cov_ij(masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, dim::Int)
    return cov_ij(masked.kernel, xview(masked,X1), xview(masked,X2), data, i, j, num_dims(masked))
end

@inline function dKij_dθp(masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::MaskedData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(masked.kernel, data.X1view, data.X2view, data.wrappeddata, i, j, p, num_dims(masked))
end
@inline function dKij_dθp(masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(masked.kernel, xview(masked,X1), xview(masked,X2), data, i, j, p, num_dims(masked))
end
@inline @inbounds function dKij_dθ!(dK::AbstractVector, masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::KernelData, i::Int, j::Int, dim::Int, npars::Int)
    X1view, X2view = xview(masked,X1), xview(masked,X2)
    return dKij_dθ!(dK, masked.kernel, X1view, X2view, data, i, j, num_dims(masked), npars)
end
@inline @inbounds function dKij_dθ!(dK::AbstractVector, masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::MaskedData, i::Int, j::Int, dim::Int, npars::Int)
    return dKij_dθ!(dK, masked.kernel, data.X1view, data.X2view, data.wrappeddata, i, j, num_dims(masked), npars)
end

function kernel_data_key(masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix)
    k = kernel_data_key(masked.kernel, xview(masked,X1), xview(masked,X2))
    return @sprintf("%s_active=%s", k, masked.active_dims)
end

get_params(masked::Masked) = get_params(masked.kernel)
get_param_names(masked::Masked) = get_param_names(masked.kernel)
num_params(masked::Masked) = num_params(masked.kernel)
set_params!(masked::Masked, hyp) = set_params!(masked.kernel, hyp)

function cov(masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::MaskedData)
    return cov(masked.kernel, data.X1view, data.X2view, data.wrappeddata)
end
function cov!(s::AbstractMatrix, masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::MaskedData)
    return cov!(s, masked.kernel, data.X1view, data.X2view, data.wrappeddata)
end
function grad_slice!(dK::AbstractMatrix, masked::Masked, X1::AbstractMatrix, X2::AbstractMatrix, data::MaskedData, iparam::Int)
    return grad_slice!(dK, masked.kernel, data.X1view, data.X2view, data.wrappeddata, iparam)
end

function cov(masked::Masked, x1::AbstractVector, x2::AbstractVector)
    return cov(masked.kernel, view(x1,masked.active_dims), view(x2,masked.active_dims))
end

# priors
get_priors(masked::Masked) = get_priors(masked.kernel)
set_priors!(masked::Masked, priors) = set_priors!(masked.kernel, priors)


text:
ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

using FastAI
using ReTest
FastAI.runtests([ReTest.fail, ReTest.not(ReTest.pass)])

module FastAITests

using InlineTest

using ..FastAI
import ..FastAI: Block, Encoding, encodedblock, decodedblock, encode, decode,
    testencoding, test_method_show, checkblock
using ..FastAI.Tabular: EncodedTableRow
using Flux.Optimise: Optimiser, ADAM, apply!
import Makie

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"
include("testdata.jl")

include("encodingapi.jl")
include("fasterai.jl")
include("training.jl")

include("makie.jl")

end

FastAITests.runtests([ReTest.fail, ReTest.not(ReTest.pass)])


text:
module MultiFieldFESpacesWithLinearConstraintsTests

using Gridap.Algebra
using Gridap.Arrays
using Gridap.Fields
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.TensorValues
using Gridap.MultiField
using Gridap.CellData
using Gridap.ReferenceFEs
using Test

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet_1",[1,2,5])
add_tag_from_tags!(labels,"neumann_1",[6,7,8])
add_tag_from_tags!(labels,"dirichlet_2",[1,3,7])
add_tag_from_tags!(labels,"neumann_2",[5,6,8])

Ω = Triangulation(model)
Γ1 = BoundaryTriangulation(model,tags="neumann_1")
Γ2 = BoundaryTriangulation(model,tags="neumann_2")
Λ = SkeletonTriangulation(model)

dΩ = Measure(Ω,2)
dΓ1 = Measure(Γ1,2)
dΓ2 = Measure(Γ2,2)
dΛ = Measure(Λ,2)

n_Γ1 = get_normal_vector(Γ1)
n_Γ2 = get_normal_vector(Γ2)

reffe = ReferenceFE(lagrangian,Float64,1)
V1 = FESpace(model,reffe,conformity=:H1,dirichlet_tags="dirichlet_1")
V2 = FESpace(model,reffe,conformity=:H1,dirichlet_tags="dirichlet_2")

sDOF_to_dof = [1,5,-2]
sDOF_to_dofs = Table([[-1,4],[4,6],[-1,-3]])
sDOF_to_coeffs = Table([[0.5,0.5],[0.5,0.5],[0.5,0.5]])
V1c = FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V1)

sDOF_to_dof = [4,-2]
sDOF_to_dofs = Table([[2,6],[-1,-3]])
sDOF_to_coeffs = Table([[0.5,0.5],[0.5,0.5]])
V2c = FESpaceWithLinearConstraints(sDOF_to_dof, sDOF_to_dofs, sDOF_to_coeffs, V2)

u1(x) = x[1] + 2*x[2]
u2(x) = 4*x[1] + x[2]
f1(x) = -Δ(u1)(x) + u2(x)
#f2(x) = -Δ(u2)(x)

U1c = TrialFESpace(V1c,u1)
U2c = TrialFESpace(V2c,u2)

U = MultiFieldFESpace([U1c,U2c])
V = MultiFieldFESpace([V1c,V2c])
test_fe_space(U)
test_fe_space(V)
@test has_constraints(U)
@test has_constraints(V)

a((u1,u2),(v1,v2)) =
  ∫( ∇(v1)⋅∇(u1) + ∇(v2)⋅∇(u2) + v1*u2 )*dΩ +
  ∫( jump(v1)*jump(u2) + jump(v2)*jump(u2) )*dΛ

b((v1,v2)) =
  ∫( v1*f1 )*dΩ +
  ∫( v1*(n_Γ1⋅∇(u1)) )*dΓ1 +
  ∫( v2*(n_Γ2⋅∇(u2)) )*dΓ2 +
  ∫( jump(v2) )*dΛ

op = AffineFEOperator(a,b,U,V)
uh = solve(op)

uh1, uh2 = uh
e1 = u1 - uh1
e2 = u2 - uh2

#using Gridap.Visualization
#v1 = FEFunction(V1,rand(num_free_dofs(V1)))
#v2 = FEFunction(V2,rand(num_free_dofs(V2)))
#v1c = FEFunction(V1c,rand(num_free_dofs(V1c)))
#v2c = FEFunction(V2c,rand(num_free_dofs(V2c)))
#writevtk(Ω,"trian",nsubcells=10,cellfields=["v1"=>v1,"v2"=>v2,"v1c"=>v1c,"v2c"=>v2c])

#using Gridap.Visualization
#writevtk(Ω,"trian",nsubcells=10,cellfields=["uh1"=>uh1,"uh2"=>uh2,"e1"=>e1,"e2"=>e2])

e1_l2 = sqrt(sum(∫(e1*e1)*dΩ))
e2_l2 = sqrt(sum(∫(e2*e2)*dΩ))

tol = 1.0e-9
@test e1_l2 < tol
@test e2_l2 < tol

end # module


text:
using Documenter, JuliaGrid

makedocs(
    sitename = "JuliaGrid",
    modules = [JuliaGrid],
    clean = false,
    doctest = false,
    format = Documenter.HTML(assets=["assets/style.css"]),
    pages = [
        "Home" => "index.md",
        "Input Data" => "man/input.md",
        "Power Flow" => "man/flow.md",
        "Optimal Power Flow" => "man/optimalflow.md",
        "State Estimation" => "man/estimation.md",
        "Measurement Generator" => "man/generator.md",
        "Theoretical Background" => [
            "Network Equations" => "man/tbbranch.md",
            "Power Flow" => "man/tbflow.md",
            "Optimal Power Flow" => "man/tboptimal.md",
            "State Estimation" => "man/tbestimate.md",
            "State Estimation Routines" => "man/tbestimaterest.md",
        ],
    ]
)

deploydocs(
    repo = "github.com/mcosovic/JuliaGrid.jl.git",
    target = "build",
)

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


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

"""
    Sphere(center, radius)

A sphere with `center` and `radius`.
"""
struct Sphere{Dim,T} <: Primitive{Dim,T}
  center::Point{Dim,T}
  radius::T
end

Sphere(center::Tuple, radius) = Sphere(Point(center), radius)

"""
    Sphere(p1, p2, p3)

A 2D sphere passing through points `p1`, `p2` and `p3`.
"""
function Sphere(p1::Point{2}, p2::Point{2}, p3::Point{2})
  x1, y1 = p2 - p1
  x2, y2 = p3 - p2
  c1 = centroid(Segment(p1, p2))
  c2 = centroid(Segment(p2, p3))
  l1 = Line(c1, c1 + Vec(y1, -x1))
  l2 = Line(c2, c2 + Vec(y2, -x2))
  center = l1 ∩ l2
  radius = norm(center - p2)
  Sphere(center, radius)
end

Sphere(p1::Tuple, p2::Tuple, p3::Tuple) = Sphere(Point(p1), Point(p2), Point(p3))

paramdim(::Type{<:Sphere{Dim}}) where {Dim} = Dim - 1

isconvex(::Type{<:Sphere}) = false

center(s::Sphere) = s.center
radius(s::Sphere) = s.radius

# https://en.wikipedia.org/wiki/N-sphere#Volume_and_surface_area
function measure(s::Sphere{Dim}) where {Dim}
  r, n = s.radius, Dim
  2π^(n/2)*r^(n-1) / gamma(n/2)
end

function Base.in(p::Point, s::Sphere)
  x = coordinates(p)
  c = coordinates(s.center)
  r = s.radius
  sum(abs2, x - c) ≈ r^2
end


text:
#https://en.wikipedia.org/wiki/Closure_%28computer_programming%29

function start(x)
    if x < 10
        x = 10 #setting conditional environment
    end
    function increment(y)
        return x + y
    end
    return increment
end

function h()
    f = start(15)
    g = start(1)

    println("15 + 4 = ",f(4))
    println("15 + 20 = ", f(20))
    println("10 + 4 = ", g(4))
    println("10 + 24 = ", g(24))
end

precompile(h)

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
struct ChangedAsset
    url ::String
    modified ::Int
    is_css ::Bool
end
mutable struct StateReload
    hash::Union{String, Nothing}
    hard::Bool
    changed_assets::Vector{ChangedAsset}
    task ::Union{Nothing, Task}
    StateReload(hash) = new(hash, false, ChangedAsset[], nothing)
end

mutable struct StateCache
    resources::ApplicationResources
    index_string ::String
    dependencies_json ::String
    need_recache ::Bool
    StateCache(app, registry) = new(_cache_tuple(app, registry)..., false)
end

_dep_clientside_func(func::ClientsideFunction) = func
_dep_clientside_func(func) = nothing
function _dependencies_json(app::DashApp)
    id_prop_named(p::IdProp) = (id = p[1], property = p[2])
    result = map(values(app.callbacks)) do dep
        (inputs = id_prop_named.(dep.id.input),
        state = id_prop_named.(dep.id.state),
        output = output_string(dep.id),
        clientside_function = _dep_clientside_func(dep.func)
        )
    end
    return JSON2.write(result)
end

function _cache_tuple(app::DashApp, registry::ResourcesRegistry)
    app_resources = ApplicationResources(app, registry)
    index_string::String = index_page(app, app_resources)
    dependencies_json = _dependencies_json(app)
    return (app_resources, index_string, dependencies_json)
end

struct HandlerState
    app::DashApp
    registry::ResourcesRegistry
    cache::StateCache
    reload::StateReload
    HandlerState(app, registry = main_registry()) = new(app, registry, StateCache(app, registry), make_reload_state(app))
end

make_reload_state(app::DashApp) = get_devsetting(app, :hot_reload) ? StateReload(generate_hash()) : StateReload(nothing)

get_cache(state::HandlerState) = state.cache

function rebuild_cache!(state::HandlerState)
    cache = get_cache(state)
    (cache.resources, cache.index_string, cache.dependencies_json) = _cache_tuple(state.app, state.registry)
    cache.need_recache = false
end


text:
include("./typedefs.jl")
include("./conversions.jl")

"""
    setglobals(;isOK, extant=scriptexists, excpn = nothing)

Set the globals of the `LabView0mqJl` module. By default do not change the value of `scriptexists`.
Use this function from the top-level scripts, e.g. as executing from LabVIEW.

# Examples
```julia-repl

julia> setglobals(;isOK=true, extant=true);

```

"""
function setglobals(; isOK, extant=scriptexists, excpn = nothing)
    global scriptOK = isOK
    global scriptexists = extant
    global scriptexcep = excpn
    return nothing
end

"""
    get_script_path(p)

Accept a Julia script (with user defined functions) as full path for user functions of just
file name for examples delivered with this package. If the file not found, return the path
to the empty file "dummy.jl". If `p` is empty string, return path to "Examples-UserFn.jl",
which in turn includes some example scripts. Set the global `scriptexists` to `true`
if everything OK.

# Examples
```julia-repl

julia> include(get_script_path(raw"C:\\Users\\eben\\Desktop\\Julia\\my_functions.jl"))
julia> server_0mq4lv((;foo=foo, bar, baz))
```
"""
function get_script_path(p="")
    global scriptexists
    srcdir = @__DIR__
    dummy = joinpath(srcdir, "dummy.jl")
    _, suffix = splitext(p)
    if p == ""
        scriptexists = true
        p1 = joinpath(srcdir, "Examples-UserFn.jl")
        return p1
    elseif suffix != ".jl"
        scriptexists = false
        return dummy
    elseif isfile(p)
        scriptexists = true
        return realpath(p)
    end

    p1 = joinpath(srcdir, p) # used for example scripts residing in the package

    if isfile(p1)
        scriptexists = true
        return p1
    else
        scriptexists = false
        return dummy
    end
end

"""
    get_LVlib_path()

A helper function: return path to the LabVIEW library `LV-Julia_ZMQ-Roundtrip_lib.lvlib`

# Examples
```julia-repl

julia> using LabView0mqJl
julia> get_LVlib_path()
C:/_LabView_projects/ZMQ/LabView0mqJl.jl/src/LabVIEW
```
"""
function get_LVlib_path()
    srcdir = @__DIR__
    return joinpath(dirname(srcdir), "LabVIEW")
end

function StackFrame_to_NamedTuple(fm)
    #= StackFrame, see:
    https://docs.julialang.org/en/v1/base/stacktraces/

    func::Symbol
       The name of the function containing the execution context.
    linfo::Union{Core.MethodInstance, CodeInfo, Nothing}
       The MethodInstance containing the execution context (if it could be found).
    file::Symbol
       The path to the file containing the execution context.
    line::Int
       The line number in the file containing the execution context.
    from_c::Bool
       True if the code is from C.
    inlined::Bool
       True if the code is from an inlined frame.
    pointer::UInt64
    =#
    return (
        frame_summary = string(fm),
        func = fm.func,
        linfo = string(fm.linfo),
        file = fm.file,
        line = fm.line,
        from_c = fm.from_c,
        inlined = fm.inlined,
        pointer = fm.pointer,
    )
end

"""
    build_err_info(;err, errcode, source, stack_trace, excep)

Combine data for LV-style error cluster (the first three arguments), and detailed
exception information.
"""
function build_err_info(;
    err::Bool = false,
    errcode::Int = 0,
    source::String = "",
    stack_trace = [],
    excep = "",
)
    if err
        # default values on error
        if errcode == 0
            errcode = 5235805
        end
    end

    if !isempty(stack_trace)
        stack_trace = [StackFrame_to_NamedTuple(e) for e in stack_trace]
    end

    return (;
        status = err,
        code = errcode,
        source = source,
        detailed_info = (exception_or_info = string(excep), trace = stack_trace),
    )
end

"""
    puttogether(;y, err, opt_header, returncode)

Combine data to be sent to LabVIEW into a byte vector.
# Arguments
- `y`: data returned by the user function. Can be Dict, NamedTuple or other data type
    convertible to pairs with  `key::T where T<: Union{Symbol, AbstractString}`
- other arguments are self-explaining
"""
function puttogether(;
    y = Dict{Symbol,Any}(),
    err = build_err_info(),
    opt_header::ByteArr = UInt8[],
    returncode::Int = 0,
)

    returncode = UInt8(returncode)
    y = Dict{Symbol,Any}(pairs(y)) # y could be Dict or named tuple

    @assert haskey(err, :status) &
            haskey(err, :code) &
            haskey(err, :source) &
            haskey(err, :detailed_info)

    if haskey(y, :bigarrs)
        arrs = pop!(y, :bigarrs)
        bin_data, bindata_descr = nums2bin(arrs)
        push!(y, :bindata_descr => bindata_descr)
    elseif haskey(y, :bin_data)
        bin_data = y[:bin_data]
    else
        bin_data = UInt8[]
    end
    push!(y, :errorinfo => err)
    jsonstring = ByteArr(JSON3.write(y))

    o_h_lng = int2bytar(length(opt_header))
    bin_lng = int2bytar(length(bin_data))
    js_lng = int2bytar(length(jsonstring))

    r = vcat(returncode, PROTOC_V, o_h_lng, bin_lng, opt_header, bin_data, jsonstring)
    return r
end

function parse_cmnd(b)
    c = b[1]
    prot_v = b[2]
    prot_OK = prot_v <= PROTOC_V
    if c == UInt8('p')
        command = :ping
    elseif c == UInt8('s')
        command = :stop
    elseif c == UInt8('c')
        command = :callfun
    else
        command = :undef
    end
    return (; command, prot_OK, prot_v)
end

# # # # # # #
"""
    parse_REQ(b)

Parse data received from LabVIEW, first pass. If there are any binary data, these will be
further parsed downstream.
# Arguments
- `b::UInt8[]`: data received
"""
function parse_REQ(b)
    opt_header = fun2call = json_data = json_dict = args = bin_data = bytearr_lng = nothing
    o_h_lng_start = 3
    bin_lng_start = o_h_lng_start + 4
    o_h_lng = bytar2int(b[o_h_lng_start:3+o_h_lng_start])
    bin_lng = bytar2int(b[bin_lng_start:3+bin_lng_start])

    o_h_start = bin_lng_start + 4
    if o_h_lng > 0
        opt_header = b[o_h_start:o_h_start+o_h_lng-1]
    end

    bin_start = o_h_start + o_h_lng
    if bin_lng > 0
        bin_data = b[bin_start:bin_start+bin_lng-1]
    end

    json_start = bin_start + bin_lng
    json_data = String(b[json_start:end])
    json_dict = Dict(JSON3.read(json_data))
    fun2call = Symbol(pop!(json_dict, :fun2call))
    args = json_dict # the rest
    if haskey(args, :bindata_descr)
        bindata_descr = pop!(args, :bindata_descr)
        numarrs = bin2nums(bin_data = bin_data, bindata_descr = bindata_descr) #TODO add semicolon before kwargs?
        args = merge(args, numarrs)
    elseif !isnothing(bin_data)
        push!(args, :bin_data => bin_data)
    end
    return (; opt_header, fun2call, args)
end

function version_this_pkg()
    v = PkgVersion.Version(LabView0mqJl)
    return (; major=v.major, minor=v.minor, patch=v.patch)
end


text:
# standard binomial
function pmf(x::TI, n::TI, p::TR; partial::Tuple{TI,TI} = (0, 0)) where {TI<:Integer,TR<:Real}

    if !(0 <= partial[1] <= partial[2] <= n)
        throw(DomainError((partial..., n), @sprintf "0 <= xpartial=%i <= npartial=%i <= n=%i violated" partial[1] partial[2] n))
    end
    deltax = x - partial[1]
    deltan = n - partial[2]
    if deltax < 0; return 0.0 end
    if deltan < deltax; return 0.0 end
    if x > n
        throw(DomainError((x, n), @sprintf "x=%i <= n=%i violated" x n))
    end
    return (gamma(deltan + 1)/gamma(deltax + 1)/gamma(deltan - deltax + 1)) *
        p^deltax * (1 - p)^(deltan - deltax)
end

function pmf(x::TI, n::TI, p::TP; partial::Tuple{TI,TI} = (0, 0)) where {TI<:Integer,TP<:Prior}
    expectation(
        p -> pmf(x, n, p; partial = partial),
        update(p, partial[1], partial[2])
    )
end

function cdf(x::TR1, n::TI, p::TR2; partial::Tuple{TI,TI} = (0, 0)) where {TI<:Integer,TR1,TR2<:Real}

    if x < 0; return 0.0 end
    if x > n; return 1.0 end
    if x == n; return 1.0 end
    if !(0 <= partial[1] <= partial[2] <= n)
        throw(DomainError((partial..., n), @sprintf "0 <= xpartial=%i <= npartial=%i <= n=%i violated" partial[1] partial[2] n))
    end
    deltax = x - partial[1]
    deltan = n - partial[2]
    if deltax < 0; return 0.0 end # is x < xpartial, so is x-1 .. 0
    if deltax > deltan; return 1.0 end
    return beta_inc(deltan - deltax, deltax + 1, 1 - p, p)[1]
end

function cdf(x::TR, n::TI, p::TP; partial::Tuple{TI,TI} = (0, 0)) where {TI<:Integer,TR<:Real,TP<:bad.Prior}
    expectation(
        p -> cdf(x, n, p; partial = partial),
        update(p, partial[1], partial[2])
    )
end


text:
using Documenter, MovingWeightedLeastSquares

makedocs(
  modules = [MovingWeightedLeastSquares],
  format = :html,
  sitename ="MovingWeightedLeastSquares.jl",
  authors = "Tung Anh Vu",
  pages = Any[
    "Home" => "index.md",
    "Constructors" => "constructors.md",
    "Approximation" => "approximation.md"
  ]
)

deploydocs(
  repo = "github.com/vutunganh/MovingWeightedLeastSquares.jl.git",
  target = "build",
  julia = "0.6",
  osname = "linux",
  deps = nothing,
  make = nothing
)


text:
# Compute extension velocity field for Fisher-Stefan
# Alex Tam, 09/11/2021

"Compute velocity extension by orthogonal extrapolation"
function extend_velocity(D, dΩ, U, ϕ, par, dx, dy)
    dτ = 0.5*(dx/5 + dy/5) # "Time" step for extrapolation PDE
    vi = interface_velocity(dΩ, U, ϕ, par, dx, dy) # Find speed at interface using Stefan condition
    V = velocity_init(D, dΩ, par, vi) # Obtain "initial" condition for V(x,y,t) using speed at interface
    v = build_vector(V, D) # Create vector using data at interior grid points
    # Solve PDE to extrapolate V outwards from interface
    prob = ODEProblem((du, u, p, t) -> vel_outward_rhs!(du, u, p, t, D, dΩ, vi, par, ϕ, dx, dy), v, (0, par.V_Iterations*dτ))
    sol = solve(prob, Tsit5(), reltol = 1e-3, abstol = 1e-6, saveat = par.V_Iterations*dτ)
    v = sol[:,end] # Update vector of V
    # # Solve PDE to extrapolate V inwards from interface
    prob = ODEProblem((du, u, p, t) -> vel_inward_rhs!(du, u, p, t, D, dΩ, vi, par, ϕ, dx, dy), v, (0, par.V_Iterations*dτ))
    sol = solve(prob, Tsit5(), reltol = 1e-3, abstol = 1e-6, saveat = par.V_Iterations*dτ)
    # Reshape solution to matrix
    V = build_v_matrix(sol[:,end], par, D)
    return V
end

"Apply Stefan condition to find speed at interface points"
function interface_velocity(dΩ, U, ϕ, par, dx, dy)
    vi = Vector{Float64}() # Pre-allocate empty vector of interface speeds
    for ip in dΩ # Loop over interface points
        # Compute du/dx, dϕ/dx, du/dy and dϕ/dy
        if ip.m.xInd != ip.p.xInd # If interface lies in x-direction
            # du/dx
            if (ip.m.Ω == true) && (ip.m.dΩ == false) # Left point lies inside Ω and not close to interface
                if ϕ[ip.m.xInd-1,ip.m.yInd] < 0 # Check adjacent grid point
                    um = U[ip.m.xInd-1,ip.m.yInd] # Set to grid value if inside Ω
                else
                    um = par.u_f # Set to interface value if outisde Ω
                end
                ux = (1/dx)*(ip.m.θxp*um/(ip.m.θxm*(ip.m.θxm+ip.m.θxp)) - (ip.m.θxm+ip.m.θxp)*U[ip.m.xInd,ip.m.yInd]/(ip.m.θxm*ip.m.θxp) + par.u_f*(2*ip.m.θxp+ip.m.θxm)/(ip.m.θxp*(ip.m.θxp+ip.m.θxm))) # Biased stencil
            elseif (ip.m.Ω == true) && (ip.m.dΩ == true) # Left point lies inside Ω and close to interface
                ϕm = ϕ[ip.m.xInd-1,ip.m.yInd]
                ϕmm = ϕ[ip.m.xInd-2,ip.m.yInd]
                if (ϕm < 0) && (ϕmm < 0) # Check two adjacent grid points
                    ux = 1/(2*dx)*(3*U[ip.m.xInd,ip.m.yInd] - 4*U[ip.m.xInd-1,ip.m.yInd] + U[ip.m.xInd-2,ip.m.yInd]) # Second-order
                elseif (ϕm < 0) && (ϕmm >= 0)
                    θ = ϕm/(ϕm - ϕmm)
                    if θ < par.θb # If opposite ghost node is close to interface
                        ux = 0.0
                    else
                        ux = (1/dx)*(par.u_f/(θ*(1+θ)) - (U[ip.m.xInd-1,ip.m.yInd]*(1+θ))/θ + (par.u_f*(2+θ))/(1+θ)) # Second-order difference
                    end
                else
                    ux = 0.0
                end
            elseif (ip.p.Ω == true) && (ip.p.dΩ == false) # Right point lies inside Ω and not close to interface
                if ϕ[ip.p.xInd+1,ip.p.yInd] < 0 # Check adjacent grid point
                    up = U[ip.p.xInd+1,ip.p.yInd] # Set to grid value if inside Ω
                else
                    up = par.u_f # Set to grid value if outside Ω
                end
                ux = (1/dx)*(-par.u_f*(2*ip.p.θxm + ip.p.θxp)/(ip.p.θxm*(ip.p.θxm+ip.p.θxp)) + (ip.p.θxm+ip.p.θxp)*U[ip.p.xInd,ip.p.yInd]/(ip.p.θxm*ip.p.θxp) - ip.p.θxm*up/(ip.p.θxp*(ip.p.θxm+ip.p.θxp))) # Biased stencil
            else # Right point lies inside Ω and close to interface
                ϕp = ϕ[ip.p.xInd+1,ip.p.yInd]
                ϕpp = ϕ[ip.p.xInd+2,ip.p.yInd]
                if (ϕp < 0) && (ϕpp < 0) # Check two adjacent grid points
                    ux = 1/(2*dx)*(-U[ip.p.xInd+2,ip.p.yInd] + 4*U[ip.p.xInd+1,ip.p.yInd] - 3*U[ip.p.xInd,ip.p.yInd]) # Second-order
                elseif (ϕp < 0) && (ϕpp >= 0)
                    θ = ϕp/(ϕp - ϕpp)
                    if θ < par.θb # If opposite ghost node is close to interface
                        ux = 0.0
                    else
                        ux = (1/dx)*(-(par.u_f*(2+θ))/(1+θ) + (U[ip.m.xInd+1,ip.m.yInd]*(1+θ))/θ - par.u_f/(θ*(1+θ))) # Second-order difference
                    end
                else
                    ux = 0.0
                end
            end
            # dϕ/dx
            ϕxl = (ϕ[ip.m.xInd+1, ip.m.yInd]-ϕ[ip.m.xInd-1, ip.m.yInd])/(2*dx)
            ϕxr = (ϕ[ip.p.xInd+1, ip.p.yInd]-ϕ[ip.p.xInd-1, ip.p.yInd])/(2*dx)
            ϕx = ip.m.θxp*ϕxr + (1-ip.m.θxp)*ϕxl
            # dϕ/dy
            ϕp = ip.m.θxp*ϕ[ip.p.xInd, ip.p.yInd+1] + (1-ip.m.θxp)*ϕ[ip.m.xInd, ip.m.yInd+1]
            ϕm = ip.m.θxp*ϕ[ip.p.xInd, ip.p.yInd-1] + (1-ip.m.θxp)*ϕ[ip.m.xInd, ip.m.yInd-1]
            ϕy = (ϕp-ϕm)/(2*dy) # Use ghost point method
            # du/dy
            if ϕp < 0 # Positive ghost point inside Ω
                up =  ip.m.θxp*U[ip.p.xInd, ip.p.yInd+1] + (1-ip.m.θxp)*U[ip.m.xInd, ip.m.yInd+1]
                ϕpp =  ip.m.θxp*ϕ[ip.p.xInd, ip.p.yInd+2] + (1-ip.m.θxp)*ϕ[ip.m.xInd, ip.m.yInd+2] # Second ghost point
                if ϕpp < 0
                    upp =  ip.m.θxp*U[ip.p.xInd, ip.p.yInd+2] + (1-ip.m.θxp)*U[ip.m.xInd, ip.m.yInd+2]
                    uy = (-upp + 4*up - 3*par.u_f)/(2*dy) # Second-order difference
                else
                    θ = ϕp/(ϕp - ϕpp) # Distance between middle ghost point and interface
                    if θ < par.θb
                        uy = 0.0 # Set to zero if ghost point close to interface
                    else
                        uy = ((-(2+θ)*par.u_f)/(1+θ) + (1+θ)*up/θ - par.u_f/(θ*(1+θ)))/dy # Second-order difference
                    end
                end
            elseif ϕm < 0 # Negative ghost point inside Ω
                um =  ip.m.θxp*U[ip.p.xInd, ip.p.yInd-1] + (1-ip.m.θxp)*U[ip.m.xInd, ip.m.yInd-1]
                ϕmm = ip.m.θxp*ϕ[ip.p.xInd, ip.p.yInd-2] + (1-ip.m.θxp)*ϕ[ip.m.xInd, ip.m.yInd-2]
                if ϕmm < 0
                    umm = ip.m.θxp*U[ip.p.xInd, ip.p.yInd-2] + (1-ip.m.θxp)*U[ip.m.xInd, ip.m.yInd-2]
                    uy = (3par.u_f - 4*um + umm)/(2*dy) # Second-order difference
                else
                    θ = ϕm/(ϕm - ϕmm) # Distance between middle ghost point and interface
                    if θ < par.θb
                        uy = 0.0 # Set to zero if ghost point close to interface
                    else
                        uy = ((2+θ)*par.u_f/(1+θ) - (1+θ)*um/θ + par.u_f/(θ*(θ+1)))/dy # Second-order difference
                    end
                end
            else # Both ghost points outside Ω
                uy = 0.0
            end
        else # If interface lies in y-direction
            # du/dx and du/dy
            if (ip.m.Ω == true) && (ip.m.dΩ == false) # Bottom point lies inside Ω and not close to interface
                if ϕ[ip.m.xInd,ip.m.yInd-1] < 0 # Check adjacent grid point
                    ub = U[ip.m.xInd,ip.m.yInd-1] # Set to grid value if inside Ω
                else
                    ub = par.u_f # Set to interface value if outside Ω
                end
                uy = (1/dy)*(ip.m.θyp*ub/(ip.m.θym*(ip.m.θym+ip.m.θyp)) - (ip.m.θym+ip.m.θyp)*U[ip.m.xInd,ip.m.yInd]/(ip.m.θym*ip.m.θyp) + par.u_f*(2*ip.m.θyp+ip.m.θym)/(ip.m.θyp*(ip.m.θyp+ip.m.θym)) ) # Biased stencil
            elseif (ip.m.Ω == true) && (ip.m.dΩ == true) # Bottom point lies inside Ω and close to interface
                ϕm = ϕ[ip.m.xInd,ip.m.yInd-1]
                ϕmm = ϕ[ip.m.xInd,ip.m.yInd-2]
                if (ϕm < 0) && (ϕmm < 0) # Check two adjacent points
                    uy = 1/(2*dy)*(3*U[ip.m.xInd,ip.m.yInd] - 4*U[ip.m.xInd,ip.m.yInd-1] + U[ip.m.xInd,ip.m.yInd-2]) # Second-order
                elseif (ϕm < 0) && (ϕmm >= 0)
                    θ = ϕm/(ϕm - ϕmm)
                    if θ < par.θb # If opposite ghost node is close to interface
                        uy = 0.0
                    else
                        uy = (1/dy)*(par.u_f/(θ*(1+θ)) - (U[ip.m.xInd,ip.m.yInd-1]*(1+θ))/θ + (par.u_f*(2+θ))/(1+θ)) # Second-order difference
                    end
                else
                    uy = 0.0
                end
            elseif (ip.p.Ω == true) && (ip.p.dΩ == false) # Top point lies inside Ω and not close to interface
                if ϕ[ip.p.xInd,ip.p.yInd+1] < 0 # Check adjacent grid point
                    ut = U[ip.p.xInd,ip.p.yInd+1] # Set to grid value if inside Ω
                else
                    ut = par.u_f # Set to interface value if outside Ω
                end
                uy = (1/dy)*(-par.u_f*(2*ip.p.θym + ip.p.θyp)/(ip.p.θym*(ip.p.θym+ip.p.θyp)) + (ip.p.θym+ip.p.θyp)*U[ip.p.xInd, ip.p.yInd]/(ip.p.θym*ip.p.θyp) - ip.p.θym*ut/(ip.p.θyp*(ip.p.θym+ip.p.θyp)) ) # Biased stencil
            else
                ϕp = ϕ[ip.m.xInd,ip.m.yInd+1]
                ϕpp = ϕ[ip.m.xInd,ip.m.yInd+2]
                if (ϕp < 0) && (ϕpp < 0) # Check two adjacent nodes
                    uy = 1/(2*dy)*(-U[ip.p.xInd,ip.p.yInd+2] + 4*U[ip.p.xInd,ip.p.yInd+1] - 3*U[ip.p.xInd,ip.p.yInd]) # Second-order
                elseif (ϕp < 0) && (ϕpp >= 0)
                    θ = ϕp/(ϕp - ϕpp)
                    if θ < par.θb # If opposite ghost node is close to interface
                        uy = 0.0
                    else
                        uy = (1/dy)*(-(par.u_f*(2+θ))/(1+θ) + (U[ip.m.xInd,ip.m.yInd+1]*(1+θ))/θ - par.u_f/(θ*(1+θ))) # Second-order difference
                    end
                else
                    uy = 0.0
                end
            end
            # dϕ/dx
            ϕp = ip.m.θyp*ϕ[ip.p.xInd+1, ip.p.yInd] + (1-ip.m.θyp)*ϕ[ip.m.xInd+1, ip.m.yInd]
            ϕm = ip.m.θyp*ϕ[ip.p.xInd-1, ip.p.yInd] + (1-ip.m.θyp)*ϕ[ip.m.xInd-1, ip.m.yInd]
            ϕx = (ϕp-ϕm)/(2*dx) # Use ghost point method
            # dϕ/dy
            ϕyb = (ϕ[ip.m.xInd, ip.m.yInd+1]-ϕ[ip.m.xInd, ip.m.yInd-1])/(2*dy)
            ϕyt = (ϕ[ip.p.xInd, ip.p.yInd+1]-ϕ[ip.p.xInd, ip.p.yInd-1])/(2*dy)
            ϕy = ip.m.θyp*ϕyt + (1-ip.m.θyp)*ϕyb
            # du/dx
            if ϕp < 0 # Positive ghost point inside Ω
                up = ip.m.θyp*U[ip.p.xInd+1, ip.p.yInd] + (1-ip.m.θyp)*U[ip.m.xInd+1, ip.m.yInd]
                ϕpp = ip.m.θyp*ϕ[ip.p.xInd+2, ip.p.yInd] + (1-ip.m.θyp)*ϕ[ip.m.xInd+2, ip.m.yInd]
                if ϕpp < 0
                    upp = ip.m.θyp*U[ip.m.xInd+2, ip.p.yInd] + (1-ip.m.θyp)*U[ip.m.xInd+2, ip.m.yInd]
                    ux = (-upp + 4*up -3*par.u_f)/(2*dx) # Second-order
                else
                    θ = ϕp/(ϕp - ϕpp) # Distance between middle ghost point and interface
                    if θ < par.θb
                        ux = 0.0 # Set to zero if ghost point close to interface
                    else
                        ux = ((-(2+θ)*par.u_f)/(1+θ) + (1+θ)*up/θ - par.u_f/(θ*(1+θ)))/dx # Second-order
                    end
                end
            elseif ϕm < 0 # Negative ghost point inside Ω
                um = ip.m.θyp*U[ip.m.xInd-1, ip.p.yInd] + (1-ip.m.θyp)*U[ip.m.xInd-1, ip.m.yInd]
                ϕmm = ip.m.θyp*ϕ[ip.p.xInd-2, ip.p.yInd] + (1-ip.m.θyp)*ϕ[ip.m.xInd-2, ip.m.yInd]
                if ϕmm < 0
                    umm = ip.m.θyp*U[ip.p.xInd-2, ip.p.yInd] + (1-ip.m.θyp)*U[ip.m.xInd-2, ip.m.yInd]
                    ux = (3*par.u_f - 4um + umm)/(2*dx) # Second-order
                else
                    θ = ϕm/(ϕm - ϕmm) # Distance between middle ghost point and interface
                    if θ < par.θb
                        ux = 0.0 # Set to zero if ghost point close to interface
                    else
                        ux = ((2+θ)*par.u_f/(1+θ) - (1+θ)*um/θ + par.u_f/(θ*(θ+1)))/dx # Second-order
                    end
                end
            else # Both ghost points outside Ω
                ux = 0.0
            end
        end
        push!(vi, -par.κ*(ux*ϕx + uy*ϕy)) # Append speed at current point to vector
    end
    return vi
end

"Obtain initial condition for V(x,y,t) using speed at interface"
function velocity_init(D, dΩ, par, vi)
    V = zeros(par.Nx, par.Ny) # Pre-allocate matrix of V(x,y,t)
    # Manually overwrite speed close to interface
    for gp in D # Loop over interior grid points
        if (gp.dΩ == true) # If grid point is close to interface
            # Locate closest interface point to current grid point
            if (gp.θxp < gp.θxm) && (gp.θxp < gp.θyp) && (gp.θxp < gp.θym) # Closest point in positive x-direction
                xm = gp.xInd; xp = gp.xInd+1; ym = gp.yInd; yp = gp.yInd
            elseif (gp.θxm < gp.θxp) && (gp.θxm < gp.θyp) && (gp.θxm < gp.θym) # Closest point in negative x-direction
                xm = gp.xInd-1; xp = gp.xInd; ym = gp.yInd; yp = gp.yInd
            elseif (gp.θyp < gp.θxm) && (gp.θyp < gp.θxp) && (gp.θyp < gp.θym) # Closest point in positive y-direction
                xm = gp.xInd; xp = gp.xInd; ym = gp.yInd; yp = gp.yInd+1
            else # Closest point in negative y-direction
                xm = gp.xInd; xp = gp.xInd; ym = gp.yInd-1; yp = gp.yInd
            end
            # Assign V to interface speed at closest point
            for i = 1:length(dΩ)
                if (xm == dΩ[i].m.xInd) && (ym == dΩ[i].m.yInd) && (xp == dΩ[i].p.xInd) && (yp == dΩ[i].p.yInd)
                    V[gp.xInd, gp.yInd] = vi[i]
                end
            end
        end
    end
    return V
end

"Obtain velocity at a specified interface point"
function obtain_interface_velocity(xm, ym, xp, yp, dΩ, vi)
    for i = 1:length(dΩ) # Loop over candidate points
        if (xm == dΩ[i].m.xInd) && (ym == dΩ[i].m.yInd) && (xp == dΩ[i].p.xInd) && (yp == dΩ[i].p.yInd) # If all indices are correct
            return vi[i]
        end
    end
end

"Construct right-hand vector of outward velocity extension PDE, for use in DifferentialEquations.jl"
function vel_outward_rhs!(du, u, p, t, D, dΩ, vi, par, ϕ, dx, dy)
    V = build_v_matrix(u, par, D) # Generate matrix
    for i = 1:length(D) # Loop over interior grid points
        if (D[i].Ω == true) || (D[i].dΩ == true) # Grid points in Ω or close to interface
            du[i] = 0.0 # Set RHS to zero
        else
            # Compute advection coefficients
            ϕx = 1/(2*dx)*(ϕ[D[i].xInd+1, D[i].yInd]-ϕ[D[i].xInd-1, D[i].yInd]) # dϕ/dx
            ϕy = 1/(2*dy)*(ϕ[D[i].xInd, D[i].yInd+1]-ϕ[D[i].xInd, D[i].yInd-1]) # dϕ/dy
            n = sqrt(ϕx^2 + ϕy^2) # Normalisation
            if n == 0 # Prevent division by zero
                a = 0; b = 0
            else
                a = ϕx/n # x-component
                b = ϕy/n # y-component
            end
            # Discretise dV/dx (first-order upwind)
            if a < 0 # Forward difference
                if D[i].θxp == 1.0 # Regular discretisation
                    Vx = (V[D[i].xInd+1, D[i].yInd] - V[D[i].xInd, D[i].yInd])/dx
                else # Irregular discretisation
                    v = obtain_interface_velocity(D[i].xInd, D[i].yInd, D[i].xInd+1, D[i].yInd, dΩ, vi)
                    Vx = (v - V[D[i].xInd, D[i].yInd])/(D[i].θxp*dx)
                end
            else # Backward difference
                if D[i].θxm == 1.0 # Regular discretisation
                    Vx = (V[D[i].xInd, D[i].yInd] - V[D[i].xInd-1, D[i].yInd])/dx
                else # Irregular discretisation
                    v = obtain_interface_velocity(D[i].xInd-1, D[i].yInd, D[i].xInd, D[i].yInd, dΩ, vi)
                    Vx = (V[D[i].xInd, D[i].yInd]-v)/(D[i].θxm*dx)
                end
            end
            # Discretise dV/dy (first-order upwind)
            if b < 0 # Forward difference
                if D[i].θyp == 1.0 # Regular discretisation
                    Vy = (V[D[i].xInd, D[i].yInd+1] - V[D[i].xInd, D[i].yInd])/dy
                else # Irregular discretisation
                    v = obtain_interface_velocity(D[i].xInd, D[i].yInd, D[i].xInd, D[i].yInd+1, dΩ, vi)
                    Vy = (v - V[D[i].xInd, D[i].yInd])/(D[i].θyp*dy)
                end
            else # Backward difference
                if D[i].θym == 1.0
                    Vy = (V[D[i].xInd, D[i].yInd] - V[D[i].xInd, D[i].yInd-1])/dy
                else
                    # Obtain interface velocity
                    v = obtain_interface_velocity(D[i].xInd, D[i].yInd-1, D[i].xInd, D[i].yInd, dΩ, vi)
                    Vy = (V[D[i].xInd, D[i].yInd]-v)/(D[i].θym*dy)
                end
            end
            du[i] = -a*Vx - b*Vy # Update RHS term
        end
    end
end

"Construct right-hand vector of inward velocity extension PDE, for use in DifferentialEquations.jl"
function vel_inward_rhs!(du, u, p, t, D, dΩ, vi, par, ϕ, dx, dy)
    V = build_v_matrix(u, par, D) # Generate matrix
    for i = 1:length(D) # Loop over interior grid points
        if (D[i].Ω == false) || (D[i].dΩ == true) # Grid points outside Ω or close to interface
            du[i] = 0.0 # Set RHS to zero
        else
            # Compute advection coefficients
            ϕx = 1/(2*dx)*(ϕ[D[i].xInd+1, D[i].yInd]-ϕ[D[i].xInd-1, D[i].yInd]) # dϕ/dx
            ϕy = 1/(2*dy)*(ϕ[D[i].xInd, D[i].yInd+1]-ϕ[D[i].xInd, D[i].yInd-1]) # dϕ/dy
            n = sqrt(ϕx^2 + ϕy^2) # Normalisation
            if n == 0 # Prevent division by zero
                a = 0; b = 0
            else
                a = -ϕx/n # x-component
                b = -ϕy/n # y-component
            end
            # Discretise dV/dx (first-order upwind)
            if a < 0 # Forward difference
                if D[i].θxp == 1.0 # Regular discretisation
                    Vx = (V[D[i].xInd+1, D[i].yInd] - V[D[i].xInd, D[i].yInd])/dx
                else # Irregular discretisation
                    v = obtain_interface_velocity(D[i].xInd, D[i].yInd, D[i].xInd+1, D[i].yInd, dΩ, vi)
                    Vx = (v - V[D[i].xInd, D[i].yInd])/(D[i].θxp*dx)
                end
            else # Backward difference
                if D[i].θxm == 1.0 # Regular discretisation
                    Vx = (V[D[i].xInd, D[i].yInd] - V[D[i].xInd-1, D[i].yInd])/dx
                else # Irregular discretisation
                    v = obtain_interface_velocity(D[i].xInd-1, D[i].yInd, D[i].xInd, D[i].yInd, dΩ, vi)
                    Vx = (V[D[i].xInd, D[i].yInd]-v)/(D[i].θxm*dx)
                end
            end
            # Discretise dV/dy (first-order upwind)
            if b < 0 # Forward difference
                if D[i].θyp == 1.0 # Regular discretisation
                    Vy = (V[D[i].xInd, D[i].yInd+1] - V[D[i].xInd, D[i].yInd])/dy
                else # Irregular discretisation
                    v = obtain_interface_velocity(D[i].xInd, D[i].yInd, D[i].xInd, D[i].yInd+1, dΩ, vi)
                    Vy = (v - V[D[i].xInd, D[i].yInd])/(D[i].θyp*dy)
                end
            else # Backward difference
                if D[i].θym == 1.0
                    Vy = (V[D[i].xInd, D[i].yInd] - V[D[i].xInd, D[i].yInd-1])/dy
                else
                    # Obtain interface velocity
                    v = obtain_interface_velocity(D[i].xInd, D[i].yInd-1, D[i].xInd, D[i].yInd, dΩ, vi)
                    Vy = (V[D[i].xInd, D[i].yInd]-v)/(D[i].θym*dy)
                end
            end
            du[i] = -a*Vx - b*Vy # Update RHS term
        end
    end
end


text:
@testset "Vec" begin
    for T ∈ Floats
        a = Vec{3,T}(1,2,3)
        b = Vec{3,T}(2,4,6)
        @test a ⊙ b ≈ [2, 8, 18]
        @test b ⊘ a ≈ [2, 2, 2]
        @test inv(a)*a ≈ 1
        @test norm²(a) ≈ 1 + 4 + 9
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
checkDir()

@testset "File loading with compensation" begin
    _, data = loadFCS("test-compensation.fcs")
    @test isapprox(
        data[:, 7:18],
        [
            3332.2309256243007 8088.362793640517 538.3720631079816 144.5766025806829 -5.026619216289584 656.7371177910801 34824.02845241197 1939.2085663518994 268.10446790549787 7164.229365371252 1487.6587206861805 -103.62798381469148
            466.6162342252232 6100.635084662845 471.12011441761274 4270.515980761523 237.1142731684298 154.2816233797197 8063.704826604625 2363.291912539149 250.09636142954616 -5.275245824878045 546.3837818701431 7.180492945860638
            4988.506790656656 5836.121341276924 305.25138438318146 -78.2911490240366 -130.7712323646924 256.75533709823753 -28.969688160222205 2463.006063022074 -107.24661567025781 212.81623260036247 50.45721112838109 946.1078887373078
            2474.009284729641 3662.411462672774 499.30026183606464 217.49388411442334 273.06602846319055 102.19931337770262 195.20362476100473 2593.704979451791 1597.8609042177852 -99.00019479802567 326.5676563660152 2481.6467014782347
            3737.3002300540215 92.63719998587527 103.7345736486239 34.18875829231229 -23.820984710851445 71.25989624438593 492.24250044483324 934.9753295720772 -32.31068183338777 114.67958090123807 239.36797617353315 23.052151803822486
        ],
    )
end


text:
"""
ULMFiT - testing file
"""

# Test the language model on the validation set
function test_lm(lm, data_gen)
    model_layers = Chain(
        lm.embedding_layer,
        lm.lstm_layer1,
        lm.lstm_layer2,
        lm.lstm_layer3,
        x -> lm.embedding_layer(x, true),
        softmax
    )
    testmode!(model_layers)
    sum_l, l_vect = 0, []
    for iter=1:num_of_iters
        x, y = take!(data_gen)
        h = broadcast(w -> model_layers[1](indices([w], lm.vocab, "_unk_")), x)
        h = model_layers[2:end].(h)
        y = broadcast(x -> Flux.onehotbatch([x], lm.vocab, "_unk_"), y)
        Flux.reset!(model_layers)
        l = sum(crossentropy.(Tracker.data.(h), y))
        sum_l += l
        push!(l_vect, l)
    end
    return sum_l/num_of_iters, l_vect
end

# Sampling
function sampling(starting_text::String, lm::LanguageModel=LanguageModel())
    model_layers = Chain(
        lm.lstm_layer1,
        lm.lstm_layer2,
        lm.lstm_layer3,
    )
    testmode!(model_layers)
    tokens = tokenize(starting_text)
    word_indices = map(x -> indices([x], lm.vocab, "_unk_"), tokens)
    embeddings = lm.embedding_layer.(word_indices)
    h = (model_layers.(embeddings))[end]
    probabilities = softmax(lm.embedding_layer(h, true))
    prediction = lm.vocab[findall(isequal(maximum(probabilities)), probabilities)[1]]
    println("SAMPLING...")
    print(prediction, ' ')
    while true
        h = lm.embedding_layer(indices([prediction], lm.vocab, "_unk_"))
        h = model_layers(h)
        probabilities = softmax(lm.embedding_layer(h, true))
        prediction = lm.vocab[findall(isequal(maximum(probabilities)), probabilities)[1]]
        print(prediction, ' ')
        prediction == "_pad_" && break
    end
end

# Test Classifier
function test_classifier(classifier::TextClassifier, data_gen)
    testmode!(classifier)
    sum_l, l_vect = 0, []
    for iter=1:num_of_iters
        x = take!(data_gen)
        y = take!(data_gen)
        h = broadcast(x -> classifier.rnn_layers[1](indices([x], lm.vocab, "_unk_")), x)
        h = classifier.rnn_layers[2:end].(h)
        h = classifier.linear_layers(h)
        Flux.reset!(classifier)
        l = crossentropy(Tracker.data(h), y)
        sum_l += l
        push!(l_vect, l)
    end
    return sum_l/num_of_iters, l_vect
end


text:
"""
convert to AutomotiveDrivingModels.Scene
"""
function state2scene(mdp::PedMDP, s::PedMDPState, ped_type::VehicleDef = mdp.ped_type)
    scene = Scene()
    push!(scene, Vehicle(s.ego, mdp.ego_type, EGO_ID))
    push!(scene, Vehicle(s.ped, mdp.ped_type, PED_ID))
    return scene
end

function animate_states(mdp::PedMDP, states::Vector{PedMDPState}, actions::Vector{PedMDPAction};
                        overlays=SceneOverlay[IDOverlay()],
                        cam=StaticCamera(VecE2(0, -5.), 17.))
    duration = length(states)*mdp.ΔT
    fps = Int(1/mdp.ΔT)
    function render_states(t, dt)
        frame_index = Int(floor(t/dt)) + 1
        scene = state2scene(mdp, states[frame_index])
        return AutoViz.render(scene,
                mdp.env,
                cat(1, overlays, TextOverlay(text = ["Acc: $(actions[frame_index].acc)"],
                                            font_size=20,
                                            pos=VecE2(0.,8.),
                                            incameraframe=true)),
                cam=cam,
                car_colors=get_colors(scene))
    end
    return duration, fps, render_states
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using Plots, GraphRecipes
import Plots: plot
export plot, spectrum

function plot(ϕ::FourierField{S, T}; ylims = ()) where {S, T <: FourierMetaData}
    dims = length(ϕ.metadata.grid.grid)
    if dims == 1
        x = ϕ.metadata.grid.grid[1][:]
        dd = ϕ.metadata.transform.backward * ϕ.data
        plot(x, real.(dd))
        plot!(xlabel = "x")
        plot!(ylabel = ϕ.metadata.name)
        plot!(legend = false, ylims = ylims)
    elseif dims == 2
        x = ϕ.metadata.grid.grid[1][:]
        y = ϕ.metadata.grid.grid[2][:]
        dd = ϕ.metadata.transform.backward * ϕ.data
        contourf(x, y, real.(dd)', linewidth = 0, color = :thermometer)
        contourf!(xlabel = "x")
        contourf!(ylabel = "y")
        contourf!(title =  ϕ.metadata.name)
    else
        print("Plotting is not supported for fields ")
        print("with dimensions greater ≥ 3")
    end
end

plot(ϕ::Field{S, T}; kwargs...) where {S <: FourierField, T} = plot(ϕ.data, kwargs)

function spectrum(ϕ::FourierField)
    dims = length(ϕ.metadata.grid.grid)
    if dims == 1
        f = ϕ.data
        g = log10.(abs.(f))[1:div(length(f),2)+1]
        ymax = maximum(g) * 1.1
        ymin = maximum(g)  - 16
        wvi = collect(1:div(length(f),2)+1) .- 1
        p1 = scatter(wvi, g, label = ϕ.metadata.name)
        scatter!(xlabel = "wavenumber index")
        scatter!(ylabel = "log10(spectral amplitude)")
        scatter!(ylims = (ymin, ymax))
        return p1
    elseif dims == 2
        g = log10.(abs.(ϕ.data))
        cmax = maximum(g) * 1.1
        cmin = maximum(g)  - 16
        clims = (cmin, cmax)
        p1 = heatmap(g, clims = clims, color = :thermometer)
        scatter!(xlabel = "x, wavenumber index")
        scatter!(ylabel = "y, wavenumber index")
        scatter!(title = "log10(spectral amplitude)")
        return p1
    end
    return nothing
end

plot(a::AbstractExpression) = plot(to_expr(a))
plot(a::Equation) = plot(to_expr(a))

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
export decompose, ms_decompose, cst_rkf, eps_rkf, weights, normlz

#------------------------------------------------------------------------
eps_rkf = eps::Float64 -> function (S)
  i :: Int = 1;
  while i<= length(S) && S[i]/S[1] > eps
    i+= 1;
  end
  i-1;
end

cst_rkf = r::Int64 -> function (S) return r end

#------------------------------------------------------------------------
# Decomposition of the pencil of matrices
function decompose(H::Vector{Matrix{C}}, lambda::Vector, rkf::Function) where C
    n = length(H)

    H0 = sum(H[i]*lambda[i] for i in 1:length(lambda))

    U, S, V = svd(H0)       # H0= U*diag(S)*V'
    r = rkf(S)

    Sr  = S[1:r]
    Sri = diagm([one(C)/S[i] for i in 1:r])

    M = Matrix{C}[]
    for i in 1:length(H)
    	push!(M, Sri*(U[:,1:r]')*H[i]*(V[:,1:r]))
    end

    if r > 1
        Xi, E, DiagInfo = diagonalization(M)
    else
        Xi = fill(zero(C),n,r)
        for i in 1:n
            Xi[i,1] = M[i][1,1]
        end
        E  = fill(1.0,1,1)
        DiagInfo = Dict{String,Any}( "case" => "1x1" )
    end

    Uxi = (U[:,1:r].*Sr')*E
    Vxi = (E\ V[:,1:r]')

    Info = Dict{String,Any}( "diagonalization" => DiagInfo)
    return Xi, Uxi, Vxi, Info
end

#------------------------------------------------------------------------
"""
```
decompose(σ :: Series{C,M}, rkf :: Function)
```
Decompose the series ``σ`` as a weighted sum of exponentials.
Return ``ω``, ``Ξ`` where
 - ``ω`` is the vector of weights,
 - ``Ξ`` is the matrix of frequency points, stored per row.
The list of monomials of degree ``\\leq {d-1 \\over 2}`` are used to construct
the Hankel matrix, where ``d`` is the maximal degree of the moments in ``σ``.

The optional argument `rkf` is the rank function used to determine the numerical rank from the vector S of singular values. Its default value `eps_rkf(1.e-6)` determines the rank as the first i s.t. S[i+1]/S[i]< 1.e-6 where S is the vector of singular values.

If the rank function cst_rkf(r) is used, the SVD is truncated at rank r.
"""
function decompose(sigma::Series{R,M},
                   rkf::Function = eps_rkf(1.e-6),
                   weps::Float64=1.e-5) where {R, M}

    d = maxdegree(sigma)
    X = variables(sigma)

    d0 = div(d-1,2); d1 = d-1-d0
    B0 = monomials(X, seq(0:d0))
    B1 = monomials(X, seq(0:d1))

    H = Matrix{R}[hankel(sigma, B0, B1)]
    for x in X
        push!(H, hankel(sigma, B0, [b*x for b in B1]))
    end

    lambda = [1.0]
    Xi, Uxi, Vxi = decompose(H, lambda,  rkf)

    n, r = size(Xi)

    w = fill(one(eltype(Xi)),r)

    for i in 1:r
        w[i] = Xi[1,i]
        Xi[:,i]/= Xi[1,i]
        w[i]*= Uxi[1,i]*Vxi[i,1]
    end

    # remove weights below threshold weps
    I = Bool[abs(w[i])>weps for i in 1:length(w)]
    return w[I], Xi[2:end,I]

end

#------------------------------------------------------------------------
function ms_decompose(sigma::Series{R,M},
                      rkf::Function = eps_rkf(1.e-6),
                      weps::Float64=1.e-5) where {R, M}
    return decompose(sigma,rkf,weps)
end
#------------------------------------------------------------------------
function normlz(M::AbstractMatrix,i)
    diagm(0 => [1/M[i,j] for j in 1:size(M,1)])*M
end


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# llplot() is our workhorse plotting function
function llplot()
    global gnuplot_state
    global gaston_config

    # select current figure
    c = findfigure(gnuplot_state.current)
    if c == 0
        println("No current figure")
        return
    end
    fig = gnuplot_state.figs[c]
    config = fig.conf

    # if figure has no data, stop here
    if isempty(fig.isempty)
        return
    end

    # Reset gnuplot settable options.
    gnuplot_send("\nreset\n")

    # Build terminal setup string and send it to gnuplot
    gnuplot_send(termstring())

    # Datafile filename. This is where we store the coordinates to plot.
    # This file is then read by gnuplot to do the actual plotting. One file
    # per figure handle is used; this avoids polutting /tmp with too many files.
    filename = joinpath(tempdir(),
                        "gaston-$(gaston_config.tmpprefix)-$(fig.handle)")
    f = open(filename,"w")

    # Send appropriate coordinates and data to gnuplot, depending on
    # whether we are doing 2-d, 3-d or image plots.

    # 2-d plot: Z is empty or plostyle is {,rgb}image
    if isempty(fig.curves[1].Z) ||
        fig.curves[1].conf.plotstyle == "image" ||
        fig.curves[1].conf.plotstyle == "rgbimage"
        # create data file
        for i in fig.curves
            ps = i.conf.plotstyle
            if ps == "errorbars" || ps == "errorlines"
                if isempty(i.E.yhigh)
                    # ydelta (single error coordinate)
                    writedlm(f,[i.x i.y i.E.ylow],' ')
                else
                    # ylow, yhigh (double error coordinate)
                    writedlm(f,[i.x i.y i.E.ylow i.E.yhigh],' ')
                end
            elseif ps == "financebars"
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(length(i.x),5)
                # output matrix
                for col = 1:length(i.x)
                    tmparr[col,1] = i.x[col]
                    tmparr[col,2] = i.F.open[col]
                    tmparr[col,3] = i.F.low[col]
                    tmparr[col,4] = i.F.high[col]
                    tmparr[col,5] = i.F.close[col]
                end
                writedlm(f,tmparr,' ')
            elseif ps == "image"
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(length(i.x)*length(i.y),3)
                tmparr_row_index = 1  # index into tmparr row
                # output matrix
                for row = 1:length(i.y)
                    x = length(i.x)
                    for col = 1:length(i.x)
                        tmparr[tmparr_row_index,1] = i.x[col]
                        tmparr[tmparr_row_index,2] = i.y[row]
                        tmparr[tmparr_row_index,3] = i.Z[row,col]
                        tmparr_row_index = tmparr_row_index+1
                        x = x-1
                    end
                end
                writedlm(f,tmparr,' ')
            elseif ps == "rgbimage"
                # data is written to tmparr, which is then written to disk
                tmparr = zeros(length(i.x)*length(i.y), 5)
                tmparr_row_index = 1
                # output matrix
                for col = 1:length(i.x)
                    y = length(i.y)
                    for row = 1:length(i.y)
                        tmparr[tmparr_row_index,1] = i.x[col]
                        tmparr[tmparr_row_index,2] = i.y[row]
                        tmparr[tmparr_row_index,3] = i.Z[y,col,1]
                        tmparr[tmparr_row_index,4] = i.Z[y,col,2]
                        tmparr[tmparr_row_index,5] = i.Z[y,col,3]
                        tmparr_row_index = tmparr_row_index+1
                        y = y-1
                    end
                end
                writedlm(f,tmparr,' ')
            else
                writedlm(f,[i.x i.y],' ')
            end
            write(f,"\n\n")
        end
        close(f)
        # send figure configuration to gnuplot
        gnuplot_send_fig_config(config)
        # Send user command to gnuplot
        !isempty(fig.gpcom) && gnuplot_send(fig.gpcom)
        # send plot command to gnuplot
        gnuplot_send(linestr(fig.curves, "plot", filename))

    # 3-d plot: Z is not empty and plotstyle is not {,rgb}image
    elseif !isempty(fig.curves[1].Z) &&
            fig.curves[1].conf.plotstyle != "image" &&
            fig.curves[1].conf.plotstyle != "rgbimage"
        # create data file
        for i in fig.curves
            # data is written to tmparr, which is then written to disk
            tmparr = zeros(1, 3)
            tmparr_row_index = 1
            for row in 1:length(i.x)
                for col in 1:length(i.y)
                    tmparr[1,1] = i.x[row]
                    tmparr[1,2] = i.y[col]
                    tmparr[1,3] = i.Z[row,col]
                    writedlm(f,tmparr,' ')
                end
                write(f,"\n")
            end
            write(f,"\n\n")
        end
        close(f)
        # send figure configuration to gnuplot
        gnuplot_send_fig_config(config)
        # Send user command to gnuplot
        !isempty(fig.gpcom) && gnuplot_send(fig.gpcom)
        # send command to gnuplot
        gnuplot_send(linestr(fig.curves, "splot",filename))
    end

    # Make sure gnuplot is done; if terminal is text, read data
    # reset error handling
    err = ""
    gnuplot_state.gp_lasterror = err
    gnuplot_state.gp_error = false

    attempt_stderr = 20
    attempt_stdout = 100
    sleep_interval = 0.05
    sleep_increment = 1.2

    gnuplot_send("printerr \"GastonDone\"\n")
    sleep(sleep_interval)

    # wait for stderr channel to be ready
    si = sleep_interval
    count = 0
    while !isready(ChanStdErr)
        sleep(si)
        si = sleep_increment * si
        count = count + 1
        count > attempt_stderr &&
        error("Gnuplot is taking too long to respond.")
    end

    # read all data in channel
    si = sleep_interval
    while isready(ChanStdErr)
        err = err * take!(ChanStdErr)
        sleep(si)
        si = sleep_increment * si
    end

    # check for errors while plotting
    if err != "GastonDone\n"
        gnuplot_state.gp_lasterror = err
        gnuplot_state.gp_error = true
        @warn("Gnuplot returned an error message:\n  $err)")
    end

    # if there was no error and text terminal, read all data from stdout
    if err == "GastonDone\n"
        if (gaston_config.terminal ∈ term_text)
            # wait for stdout to be ready
            si = sleep_interval
            count = 0
            while !isready(ChanStdOut)
                sleep(si)
                si = sleep_increment * si
                count = count + 1
                count > attempt_stdout &&
                error("Gnuplot is taking too long to respond.")
            end
            svgdata = ""
            si = sleep_interval
            while isready(ChanStdOut)
                svgdata = svgdata * take!(ChanStdOut)
                sleep(si)
                si = sleep_increment * si
            end
            fig.svg = svgdata
        end
    end

    return nothing

end


text:
# Support files containing all functions for Kalman Filtering and MHTs

function InitKF(d1::Detection,d2::Detection)
    # d1 - first detection of object
    # d2 - second detection of object (assumed to be associated)

    # Z matrix
    Δt = 1                       # 1 Frame
    TL = zeros(2,2)              # Top Left
    TR = Matrix(I, 2, 2)         # Top Right
    BL = (-1/Δt)*Matrix(I, 2, 2) # Bottom Left
    BR = (1/Δt)*Matrix(I, 2, 2)  # Bottom Right
    Z₀ = [TL TR; BL BR]          # Make Z Matrix

    # Centroid Measurments
    y_prev = [d1.Centroid[1]; d1.Centroid[2]]
    y₀     = [d2.Centroid[1]; d2.Centroid[2]]

    # Initialized State x₀
    x₀ = Z₀*[y_prev ; y₀]

    # Centroid Uncertainties
    R_prev = [d1.CentroidUncertainty[1]; d1.CentroidUncertainty[2]].*Matrix(I, 2, 2)
    R₀     = [d2.CentroidUncertainty[1]; d2.CentroidUncertainty[2]].*Matrix(I, 2, 2)
    R      = [R_prev TL; TL R₀]

    # Initialize Covariance P₀
    P₀ = Z₀*R*Z₀'

    # Initialize residuals r₀
    r₀ = zeros(2,1)

    # Create Track Object
    T = Track([d1,d2],SMatrix{4,1}(x₀[:]), SMatrix{4,4}(P₀[:]), SMatrix{2,1}(r₀[:]))

    return TrackHistory([T])
end # end InitKF

# Init Kalman Filter with KF method
function KF(d1::Detection,d2::Detection)
    return InitKF(d1,d2)
end # end KF

function TimeUpdate(x::Track, tk::Float64, T::TuningParams)

    σw² = T.ProcessNoise

    tk₋ = x.Measurements[end].Time

    # State Transition Matix
    Φ(tk,tk₋) = [Matrix(I, 2, 2) (tk-tk₋)*Matrix(I, 2, 2); zeros(2,2) Matrix(I, 2, 2)]

    # Discrete-time Equivalent of B
    Γ(tk,tk₋) = [(exp(tk-tk₋)-1)*Matrix(I, 2, 2) ; (tk-tk₋)*Matrix(I, 2, 2)]

    # Process Noise
    function Qk(tk,tk₋,σw²)#,W)
        Δt = tk - tk₋
        #σw² = 0.1 #var(W)
        Q  = σw²*Matrix(I, 2, 2)
        return [(Q*Δt^3)/3 (Q*Δt^2)/2; (Q*Δt^2)/2 Q*Δt]
    end

    xk₋⁺ = x.State[end] # Adding [end] here probs breaks KF!
    Pk₋⁺ = x.Covariance[end]
    uk₋  = zeros(2,1) # No input
    #W    = ModTrack[tk][2][1][5]
    xk⁻  = Φ(tk,tk₋)*xk₋⁺ + Γ(tk,tk-1)*uk₋
    Pk⁻  = Φ(tk,tk₋)*Pk₋⁺*(Φ(tk,tk-1))' + Qk(tk,tk-1,σw²)#,W)

    return SMatrix{4,1}(xk⁻[:]), SMatrix{4,4}(Pk⁻[:])
end # End TimeUpdate

function MeasurementUpdate(xk⁻::SMatrix{4,1}, Pk⁻::SMatrix{4,4}, yk::SMatrix{2,1}, Rk::SMatrix{2,2})
    # Residuals, Kalman Gain
    H   = [Matrix(I, 2, 2) zeros(2,2)]
    rk⁻ = yk - H*xk⁻
    Kk  = Pk⁻*H'*inv(H*Pk⁻*H' + Rk) # Changed from Notes

    # Measurement Update
    xk⁺ = xk⁻ + Kk*rk⁻ # Changed from Notes
    rk⁺ = yk - H*xk⁺
    Pk⁺ = (Matrix(I, 4, 4) - Kk*H)*Pk⁻*(Matrix(I, 4,4)-Kk*H)' + Kk*Rk*Kk'

    return SMatrix{4,1}(xk⁺[:]),SMatrix{4,4}(Pk⁺[:]),SMatrix{2,1}(rk⁺[:])

end # End MeasurementUpdate

# Linear Kalman Filter for Image Space # CURRENTLY UNTESTED AFTER CHANGES
function KF!(T::TrackHistory,d::Detection,TP::TuningParams)

    # Time Update
    xk⁻, Pk⁻ = TimeUpdate(T.Tracks[end], d.Time, TP)

    # Read Next Observation
    yk = SMatrix{2,1}(d.Centroid[1], d.Centroid[2])
    R = [d.CentroidUncertainty[1]; d.CentroidUncertainty[2]].*Matrix(I, 2, 2)
    Rk = SMatrix{2,2}(R[:])

    # Measurement Update
    xk⁺,Pk⁺, rk⁺ = MeasurementUpdate(xk⁻, Pk⁻, yk, Rk)

    # Measurement History
    detections = [T.Tracks[end].Measurements;d]

    # Update Tracks
    push!(T.Tracks,Track(detections,SMatrix{4,1}(xk⁺[:]), SMatrix{4,4}(Pk⁺[:]), SMatrix{2,1}(rk⁺[:])))

    return T
end # end KF!

# Make a Dummy Detection (for missed detections)
function dummydetect(Time::Float64;Pd::Float64=1.)
    return Detection(0, Time, Pd, SVector(0.,0.),SVector(0.,0.),SVector(0.),SVector(0.),SVector(0., 0., 0., 0.))
end

# Birth Model
function BirthModel(d::Detection;n_dummys::Int64=0)
    # Starlink Birth Model - velocity
    tanFOV = (11/50)    # Camera half height is 11mm and 60mm focal length
    r = 6378.1+500      # Km radius of starlink
    v = sqrt(398600/r)  # km/s velocity (circular orbit)
    x = r*tanFOV        # half FOV in km
    pixperkm = 1080/x   # Pixels per KM
    vpix = v*pixperkm   # velocity in pixels/sec
    vpixframe = vpix/30 # velcity pixels per frame (30Hz framerate)

    state = vcat(d.Centroid[:], 0.0, 0.0)
    covar = Matrix(1.0I, 4, 4)
    covar[1,1] = d.CentroidUncertainty[1]
    covar[2,2] = d.CentroidUncertainty[2]
    covar[3,3] = sqrt(5) # This is hard coded for now
    covar[4,4] = sqrt(5) # This is hard coded for now
    r = [0. 0.]

    # Add a bunch of dummy detects for prior missed frames when birthing new track
    detectarray = [d]
    states = [SMatrix{4,1}(state[:])]
    covars = [SMatrix{4,4}(covar[:])]
    resids = [SMatrix{2,1}(r[:])]
    if n_dummys>0
        for ii = 1:n_dummys
            newdummy = dummydetect(d.Time-ii)
            detectarray = vcat(newdummy,detectarray)
            states = vcat([SMatrix{4,1}([0.;0.;0.;0.])],states)
            covars = vcat([SMatrix{4,4}([0. 0. 0. 0.; 0. 0. 0. 0.; 0. 0. 0. 0.; 0. 0. 0. 0.])],covars)
            resids = vcat([SMatrix{2,1}([0.;0.])],resids)
        end
    end

    return Track(detectarray,states,covars,resids)
end

# TrackCost without a PrevCost input calculates the cost of the current
# time step using all the available history in the track and assignments
function TrackCost(track::Track, assignments::Array{Int64,1})

    S = length(assignments)
    Cost = 0.
    λν = 0.2 # I dont know what this should be
    λc = 1080*1920*1e-5 # Hard code for now but should pass this in
    λex = λν + λc
    for ii in 1:S
        #Pd = clamp(track.Measurements[end].Pd,0.001,.999)
        Pd = .999
        Σ = track.Covariance[ii][1:2,1:2]
        r = track.Residuals[ii]
        if assignments[ii] == 0
            δ = 0
            Term1 = (δ-1)*log(1-Pd) +100
            Term2 = 0
            Term3 = 0
        else
            δ = 1
            Term1 = (δ-1)*log(1-Pd)
            Term2 = -δ*log(Pd/(λex*sqrt(norm(2*π*Σ))))
            Term3 = δ*(1/2*r'*inv(Σ)*r)[end] # Do this in the IF becuase zeros'*inv(zeros)*zeros fails for dummy detections
        end
        c = Term1 + Term2 + Term3
        Cost += c
    end
    return Cost
end # End TrackCost

# TrackCost with a PrevCost input calculates updates the cost assuming the
# previous cost is a running sum of all previous cost
function TrackCost(track::Track, assignments::Array{Int64,1}, PrevCost::Float64)

    λν = 0. # I dont know what this should be but since im cheating on gamma anyways its probably fine for now
    λc = 1080*1920*1e-5 # Hard code for now but should pass this in
    λex = λν + λc
    #Pd = clamp(track.Measurements[end].Pd,0.001,.999)
    Pd = .999
    Σ = track.Covariance[end][1:2,1:2]
    r = track.Residuals[end]
    if assignments[end] == 0
        δ = 0
        Term1 = (δ-1)*log(1-Pd) + 100
        Term2 = 0
        Term3 = 0
    else
        δ = 1
        Term1 = (δ-1)*log(1-Pd)
        Term2 = -δ*log(Pd/(λex*sqrt(norm(2*π*Σ))))
        Term3 = δ*(1/2*r'*inv(Σ)*r)[end] # Do this in the IF becuase zeros'*inv(zeros)*zeros fails for dummy detections
    end
    #show δ, Term1, Term2, Term3, assignments
    c = Term1 + Term2 + Term3
    Cost = c + PrevCost

    return Cost
end # End TrackCost

# Init TOMHT on first frame with only Detections
function TOMHT(Detections::Array{Detection,1})
    # Number of Detections
    nd = length(Detections)

    # Create a dummy track for each detection and assign it to an AssociationMatrix
    tracks = Array{Track,1}(undef,nd)
    AssMat = Array{Int,2}(undef,nd,1)
    cost = zeros(Float64,size(AssMat))
    for ii in 1:nd
        tracks[ii] = BirthModel(Detections[ii])
        AssMat[ii] = Detections[ii].ID
        cost[ii]   = TrackCost(tracks[ii], [AssMat[ii]], 0.0)
    end

    return tracks, AssociationMatrix(AssMat), HypothesisCost(cost), Catalog(Array{Int64,1}(undef,0),Array{Track,1}(undef,0))
end # end TOMHT

# Track Oriented Multiple Hypothesis Tracker (TOMHT)
function TOMHT!(Detections::Array{Detection,1},Tracks::Array{Track,1},AssMat::AssociationMatrix,Cost::HypothesisCost, catalog::Catalog, T::TuningParams, FrameNum::Int64)

    # Check AssMat and Tracks are the same size
    assmat  = copy(AssMat.TOMHT)
    CostMat = copy(Cost.CostMatrix)
    nTracks = length(Tracks)
    nAssMat = size(assmat,1)
    ndetect = length(Detections)

    # Add a new empty column to AssMat and CostMat that needs to be filled
    TrackHypoth = copy(Tracks) # Dont want to return a reference - or do i?
    newcol      = Array{Int,2}(undef,nAssMat,1)
    newassmat   = hcat(assmat,newcol)
    newcostmat  = hcat(CostMat,newcol)
    newcat      = copy(catalog.Tracks)
    frameends = copy(catalog.EndFrames)

    # Create Hypotheses for each track-detection pair in each cluster
    clusters = Array{Cluster,1}(undef,nTracks)
    rowoffset = 0 # Row Offset Value (add to this when a new row (hypothesis) is inserted for a track)
    for ii in 1:nTracks

        # Time Update to get predicted state at current frame
        tk = Detections[1].Time
        xk⁻, Pk⁻ = TimeUpdate(Tracks[ii], tk, T)

        # Measurement Gating - make clusters of detections near each track (in measurment space)
        H = [Matrix(I, 2, 2) zeros(2,2)]
        P = H*Pk⁻*H'
        MD = zeros(Float64,size(Detections))
        for ii in 1:length(Detections)
            d = Detections[ii]
            R = [d.CentroidUncertainty[1]; d.CentroidUncertainty[2]].*Matrix(I, 2, 2)
            S = P + R
            rm = reshape(d.Centroid,2,1) - H*xk⁻
            MD[ii] = (rm'*inv(S)*rm)[end] # Mahalanobis Distance
        end
        boolindx = vec(MD.<5^2) # 2 sigma mahalanobis distance
        clusters[ii] = Cluster(Tracks[ii], Detections[boolindx])

        # Missed Detection Hypothesis - will always be the first hypothesis
        dummy = dummydetect(Detections[1].Time)#,Pd=Tracks[ii].Measurements[end].Pd)
        TrackHypoth[ii+rowoffset] = Track(vcat(Tracks[ii].Measurements,dummy),vcat(Tracks[ii].State,[SMatrix{4,1}(xk⁻[:])]), vcat(Tracks[ii].Covariance,[SMatrix{4,4}(Pk⁻[:])]), vcat(Tracks[ii].Residuals,[SMatrix{2,1}([0.;0.])]))
        newassmat[ii+rowoffset,end] = dummy.ID # Missed Detection
        newcostmat[ii+rowoffset,end] = TrackCost(TrackHypoth[ii+rowoffset], newassmat[ii+rowoffset,:], CostMat[ii,end])

        # Assign each detection in cluster as possible association hypothesis
        if !isempty(clusters[ii].Detections)
            for d in clusters[ii].Detections
                # Assignment
                newrow    = copy(newassmat[ii+rowoffset,:]) # Make a copy of the track (row)
                newrow[end] = d.ID                          # Assign detection to the track
                frontass  = newassmat[1:ii+rowoffset,:]
                backass   = newassmat[ii+rowoffset+1:end,:]
                newassmat = vcat(frontass, newrow', backass)

                # Track Hypothesis Creation
                xk⁻ = TrackHypoth[ii+rowoffset].State[end] # Remember current state has been time updated
                Pk⁻ = TrackHypoth[ii+rowoffset].Covariance[end]
                yk  = SMatrix{2,1}(d.Centroid[1], d.Centroid[2])
                R   = [d.CentroidUncertainty[1]; d.CentroidUncertainty[2]].*Matrix(I, 2, 2)
                Rk  = SMatrix{2,2}(R[:])
                xk⁺,Pk⁺,rk⁺ = MeasurementUpdate(xk⁻, Pk⁻, yk, Rk)
                newTrack    = Track(vcat(Tracks[ii].Measurements[1:end],d),vcat(Tracks[ii].State[1:end],[SMatrix{4,1}(xk⁺[:])]), vcat(Tracks[ii].Covariance[1:end],[SMatrix{4,4}(Pk⁺[:])]), vcat(Tracks[ii].Residuals[1:end],[SMatrix{2,1}(rk⁺)]))
                fronttrack  = TrackHypoth[1:ii+rowoffset]
                backtrack   = TrackHypoth[ii+rowoffset+1:end]
                TrackHypoth = vcat(fronttrack, newTrack, backtrack)

                # Cost Calculation
                newcostrow = copy(newcostmat[ii+rowoffset,:])
                newcostrow[end] = TrackCost(TrackHypoth[ii+rowoffset], newrow, CostMat[ii,end])
                frontcost  = newcostmat[1:ii+rowoffset,:]
                backcost   = newcostmat[ii+rowoffset+1:end,:]
                newcostmat = vcat(frontcost, newcostrow', backcost)

                # Increment Counter
                rowoffset += 1
            end
        end
    end

    # Lastly, detections could be a new track
    n_prevframe = T.ScanningWindow
    for ii in 1:ndetect
        TrackHypoth = vcat(TrackHypoth, BirthModel(Detections[ii],n_dummys=n_prevframe))
        newrow = zeros(Int64,size(newassmat[1,:]))
        newrow[end] = Detections[ii].ID
        newassmat = vcat(newassmat,newrow')
        newcostrow = zeros(Float64,size(newcostmat[1,:]))
        TrackMultiplyer = max(1.,FrameNum-T.ScanningWindow)
        newcostrow[1] = TrackMultiplyer*TrackCost(TrackHypoth[end], [newrow[1]], 0.0)
        for jj in 2:length(newcostrow)-1
            newcostrow[jj] = newcostrow[jj-1] + newcostrow[1]
        end
        newcostrow[end] = TrackCost(TrackHypoth[end], [newrow[end]], newcostrow[end-1])
        newcostmat = vcat(newcostmat,newcostrow')
    end

    # Only save info for 5 frames
    if size(newassmat,2) > T.ScanningWindow
        # Cut First Column
        newassmat  = newassmat[:,2:end]
        newcostmat = newcostmat[:,2:end]
        # A rows of all dummy tracks means track is gone (ex: [1 0 0 0 0 0] -> [0 0 0 0 0] -> Write to catalog)
        rmidx       = findall(x -> x==0 , sum(newassmat,dims=2)[:])
        newassmat   = newassmat[setdiff(1:end, rmidx), :]
        newcostmat  = newcostmat[setdiff(1:end, rmidx), :]
        # Put Cut Tracks into the catalog
        newcat      = vcat(catalog.Tracks, TrackHypoth[rmidx])
        frames      = Vector{Int64}(undef,length(rmidx))
        frames     .= FrameNum
        frameends   = vcat(catalog.EndFrames,frames)
        TrackHypoth = TrackHypoth[setdiff(1:end, rmidx), :][:]
    end

return TrackHypoth, AssociationMatrix(newassmat), HypothesisCost(newcostmat), Catalog(frameends,newcat)

end # end TOMHT!

export KF, KF!, TOMHT, TOMHT!


text:
function lpc2c!(c::AbstractVector, a::AbstractVector, loggain::Bool=false)
    @assert length(a) == length(c)
    m = length(a) - 1
    fill!(c, zero(eltype(c)))
    if loggain
        c[1] = a[1]
    else
        c[1] = log(a[1])
    end
    c[2] = -a[2]
    for k=3:m+1
        d = 0.0
        upperlimit = (k >= m+1) ? m+1 : k
        for i=2:upperlimit
            d += (i-1) *  c[i] * a[k-i+1]
        end
        c[k] = -d / (k-1) - a[k]
    end

    c
end

lpc2c(a::AbstractVector, loggain::Bool=false) = lpc2c!(similar(a), a, loggain)

function lpc2c(state::SpectralParamState{LinearPredictionCoef})
    assert_not_ready_to_filt(state)

    def = paramdef(state)
    newdef = LinearCepstrum(param_order(def))
    data = rawdata(state)
    SpectralParamState(newdef, lpc2c(data, has_loggain(state)),
                       true, gain_normalized(state))
end


text:
using QuantExQASM
using PicoQuant

using Compat.Test
using TestSetExtensions
using PyCall

include("test_utils.jl")

@testset "All the tests" begin
    @includetests ARGS
end


text:
using Vulkan
using Dictionaries
using Rhyolite
using Test
using Accessors

resource(filename) = joinpath(@__DIR__, "resources", filename)

@testset "Rhyolite.jl" begin
    @testset "Initialization" begin
        config = dictionary([QUEUE_GRAPHICS_BIT | QUEUE_COMPUTE_BIT => 1])
        instance, device = init(; queue_config = config, device_extensions = ["VK_KHR_synchronization2"])
        disp = QueueDispatch(device, config)
        pools = ThreadedCommandPool(device, disp, dictionary(1:4 .=> collect(keys(config))))
        ev = Event(device)
        buffer = first(unwrap(allocate_command_buffers(device, CommandBufferAllocateInfo(CommandPool(pools), COMMAND_BUFFER_LEVEL_PRIMARY, 1))))
        @record buffer begin
            cmd_set_event(ev, PIPELINE_STAGE_TOP_OF_PIPE_BIT)
        end
        queue = submit(disp, QUEUE_GRAPHICS_BIT, [SubmitInfo2KHR([], [CommandBufferSubmitInfoKHR(buffer, 0)], [])])
        @test unwrap(queue_wait_idle(queue)) == SUCCESS
    end

    @testset "Shaders" begin
        include("shaders.jl")
    end

    @testset "Pipelines" begin
        include("pipelines.jl")
    end

    @testset "Rendering" begin
        @testset "Headless rendering" begin
            include("headless.jl")
        end
    end
end


text:
## extend RegionTrees

"""
    getnthcell(c::Cell, n)

Get the cell at the n-th level (at the depth `n`) from the parents of a cell.
If `n` is larger than the depth of `c`, return nothing
(it means, that the asked depth is below the current level).
Also return nothing, if the asked level is below the minimum level, currently set to 1.
"""
function getnthcell(c::Cell, n)
    #n < 1 && throw(BoundsError(c, n))
    n < 1 && return nothing
    maxd = c.data.depth
    n == maxd && return c
    while true
        c = parent(c)
        c === nothing && return nothing
        c.data.depth == n && return c
    end
    return nothing
end

"""
    struct RANSACCloud{A<:AbstractArray, B<:AbstractArray, C<:AbstractArray}

A struct to wrap a point cloud. Stores the vertices, the normals,
    the octree of the vertices,
    the subsets (as an array of arrays of vertice indexes),
    an array to indicate if a point is part of an already extracted primitive,
    the size of the point cloud (number of vertices),
    the weight of each octree level (populated during the construction of the octree,
    by copying the given element),
    an array to store the sum of the score of already
    extracted primitives for each octree level.
"""
struct RANSACCloud{A<:AbstractArray, B<:AbstractArray, C<:AbstractArray}
    vertices::A
    normals::B
    octree::Cell
    subsets::Vector{Vector{Int}}
    isenabled::BitArray{1}
    size::Int
    levelweight::C
    levelscore::C
    #=
    # Arguments
    - `vertices::AbstractArray`: an array of vertices.
    - `normals::AbstractArray`: an array of surface normals.
    - `subsets::Vector{Vector{Int}}`: an array consisting of index-arrays,
        that describe which point is in which subset.
    - `isenabled::BitArray{1}`: an array indicating if the given point is enabled
        (`true`) or already has been extracted (`false`).
    - `size::Int`: number of points in the cloud.
    - `levelweight::Vector{Float}`: the weight of every octree level.
    - `levelscore::Vector{Float}`: for every octree level,
        store the sum of the scores of primitives that have been extracted.
    =#
end

Base.show(io::IO, x::RANSACCloud{A,B,C}) where {A,B,C} =
    print(io, "RANSACCloud of size $(x.size)")

Base.show(io::IO, ::MIME"text/plain", x::RANSACCloud{A,B,C}) where {A,B,C} =
    print(io, "RANSACCloud of size $(x.size) & $(length(x.subsets)) subsets")

"""
    nomodRANSACCloud(vertices, normals, subsets)

Construct a `RANSACCloud` without touching the vertices, normals and subsets.
Other fields are computed.

# Arguments
- `vertices`: an array of vertices.
- `normals`: an array of surface normals.
- `subsets::Vector{Vector{Int}}`: a list of indexes for each subset.
"""
function nomodRANSACCloud(vertices, normals, subsets)
    octree = buildoctree(vertices)
    octree_d = octreedepth(octree)
    s = size(vertices, 1)
    levelscore = zeros(eltype(eltype(vertices)), (octree_d,))
    levelweight = fill!(similar(levelscore), 1/octree_d)
    RANSACCloud(vertices, normals, octree, subsets, trues(s), s, levelscore, levelweight)
end

"""
    RANSACCloud(vertices, normals, subsets; force_eltype::Union{Nothing,DataType}=nothing)

Construct a `RANSACCloud` with the given `subsets`.
Vertices and normals are converted to array of `SVector`s.

# Arguments
- `vertices`: an array of vertices.
- `normals`: an array of surface normals.
- `subsets::Vector{Vector{Int}}`: a list of indexes for each subset.
- `force_eltype::Union{Nothing,DataType}=nothing`:
    an element type can be forced for the normals and vertices
    (in practice `Float64` or `Float32`). If not specified,
    the element type of the passed `vertices` and `normals` will be used.
"""
function RANSACCloud(vertices, normals, subsets; force_eltype::Union{Nothing,DataType}=nothing)
    @assert size(vertices) == size(normals) "Every point must have a normal."
    v_type = force_eltype === nothing ? eltype(eltype(vertices)) : force_eltype
    n_type = force_eltype === nothing ? eltype(eltype(normals)) : force_eltype
    vc = [SVector{3,v_type}(v) for v in vertices]
    nc = [SVector{3,n_type}(v) for v in normals]
    nomodRANSACCloud(vc, nc, subsets)
end

"""
    RANSACCloud(vertices, normals, numofsubsets::Int)

Construct a `RANSACCloud` with `numofsubsets` number of random subsets.
Vertices and normals are converted to array of `SVector`s.

# Arguments
- `vertices`: an array of vertices.
- `normals`: an array of surface normals.
- `numofsubsets::Int`: number of subsets.
- `force_eltype::Union{Nothing,DataType}=nothing`:
    an element type can be forced for the normals and vertices
    (in practice `Float64` or `Float32`). If not specified,
    the element type of the passed `vertices` and `normals` will be used.
"""
function RANSACCloud(vertices, normals, numofsubsets::Int; force_eltype::Union{Nothing,DataType}=nothing)
    @assert numofsubsets > 0 "At least 1 subset please!"
    # subset length
    function makesubset(l, n)
        ssl = fld(l, n)
        alls = randperm(l)
        subsets = [ alls[(i-1)*ssl+1:i*ssl] for i in 1:n-1]
        push!(subsets, alls[(n-1)*ssl+1:end])
        subsets
    end
    subs = makesubset(length(vertices), numofsubsets)
    return RANSACCloud(vertices, normals, subs; force_eltype=force_eltype)
end

"""
    struct OctreeNode{B<:AbstractArray}

Store the data attached to one cell of the octree.
Stores the indices of the points that are inside of the cell,
and its own depth (depth of the root is 1).
"""
struct OctreeNode{B<:AbstractArray, I<:Integer}
    incellpoints::B
    depth::I
end

Base.show(io::IO, x::OctreeNode{A}) where {A} =
    print(io, "OctreeNode: $(length(x.incellpoints)) ps, $(x.depth) d")

Base.show(io::IO, ::MIME"text/plain", x::OctreeNode{A}) where {A} =
    print(io, """OctreeNode{$A}\n with $(length(x.incellpoints)) points, at depth $(x.depth)""")

struct OctreeRefinery{A<:AbstractArray} <: AbstractRefinery
    count::Int
    vertices::A
end

function needs_refinement(r::OctreeRefinery, cell)
    length(cell.data.incellpoints) > r.count
end

function refine_data(r::OctreeRefinery, cell::Cell, indices)
    # new boundary
    boundary = child_boundary(cell, indices)
    # points that are in the current cell
    points = @view r.vertices[cell.data.incellpoints]
    # points that are in the new cell
    bolarr = map(x -> iswithinrectangle(boundary, x), points)
    # new cell inherits ++ of the depth
    d = cell.data.depth + 1
    OctreeNode(cell.data.incellpoints[bolarr], d)
end

"""
    iswithinrectangle(rect::HyperRectangle, p)

Decide if p (3 dimensional point) is within rect.

Works only with those prisms which axes are aligned with the coordinate frame.
"Bottom/left" gives `false` and "top/right" gives `true`.
"""
function iswithinrectangle(rect::HyperRectangle, p)
    vs = vertices(rect)
    vmin = vs[1,1,1]
    vmax = vs[2,2,2]
    for i in 1:3
        vmin[i] < p[i] || return false
        vmax[i] >= p[i] || return false
    end
    return true
end

function updatelevelweight(pc, x = 9//10)
    P = pc.levelweight
    σ = pc.levelscore
    w = sum( σ[i]/P[i] for i in eachindex(P) )
    for i in eachindex(P)
        pc.levelweight[i] = x*σ[i]/(w*P[i]) + (1-x)/length(P)
    end
end

"""
    octreedepth(pc::RANSACCloud)

Compute the depth of the octree stored in the cloud.
"""
function octreedepth(pc::RANSACCloud)
    return octreedepth(pc.octree)
end

"""
    octreedepth(cell::Cell)

Iterate through the leafes and search the maximum depth (`cell.data.depth`).
"""
function octreedepth(cell::Cell)
    maxdepth = cell.data.depth
    alll = allleaves(cell)
    for a in alll
        if a.data.depth > maxdepth
            maxdepth = a.data.depth
        end
    end
    return maxdepth
end

"""
    buildoctree(vertices)

Build octree of the vertices.
"""
function buildoctree(vertices)
    l = size(vertices, 1)
    minV, maxV = findAABB(vertices)
    root=Cell(SVector{3}(minV), SVector{3}(maxV), OctreeNode(collect(1:l), 1))
    r = OctreeRefinery(8, vertices)
    adaptivesampling!(root, r)
    return root
end


text:
ENV["GKSwstype"]="100" #to enable saving figure in server

#this file is to calculate prediction accuracy for NN-LMM
using DataFrames,CSV,Random,DelimitedFiles,Statistics
using Plots

Random.seed!(123)
method="omics_gblup_1200"  #NN-GBLUP
main_path="/group/qtlchenggrp/tianjing/omics_fix_all/full_omics/"
cd(main_path*method)

nRep=20  #data1,...,data20
chain_length=5000
final_accuracy_all=ones(nRep,2)*999.0

for i in 1:nRep
    #read bv
    bv=vec(readdlm("/group/qtlchenggrp/tianjing/omics_data/andrew1055/data$i/tbv_m.txt"));

    #read results for repi
    if ispath("data$i/results/MCMC_samples_EBV_NonLinear.txt")
        mcmc=readdlm("data$i/results/MCMC_samples_EBV_NonLinear.txt",',')[2:end,:]
        if size(mcmc,1)==1000 #1000 iterations were saved
            mcmc_mean=zeros(size(mcmc))
            for i in 1:1000
                mcmc_mean[i,:]=mean(mcmc[1:i,:],dims=1)
            end
            accuracy_datai=zeros(1000)
            for i in 1:1000
                accuracy_datai[i]=cor(mcmc_mean[i,:][956:1055],bv[956:1055]) #testing individual
            end

            #save results
            open(main_path*"accuracy_$method/accuracy.$method.data$i.txt", "w") do io
                writedlm(io, accuracy_datai)
            end
            final_accuracy=round(accuracy_datai[end],digits=3)
            final_accuracy_all[i,:]=[i,accuracy_datai[end]]

            #plot results
            myfig=plot(collect(chain_length/1000:chain_length/1000:chain_length),accuracy_datai,xlabel="iteration",ylabel="accuracy",title="full_omics.data$i.accu=$final_accuracy")
            savefig(myfig,main_path*"accuracy_$method/accuracy.$method.data$i.png")
        end
    end
end
open(main_path*"accuracy_$method/accuracy_final.$method.txt", "w") do io
    writedlm(io, final_accuracy_all)
end

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## dummy stub for https://github.com/JuliaBinaryWrappers/LibCURL_jll.jl

module LibCURL_jll

using Libdl

const PATH_list = String[]
const LIBPATH_list = String[]

export libcurl

# These get calculated in __init__()
libcurl_handle = C_NULL
libcurl_path = ""

if Sys.iswindows()
    const libnghttp2 = "libnghttp2-14.dll"
    const libcurl = "libcurl-4.dll"
elseif Sys.isapple()
    const libnghttp2 = "libnghttp2.14.dylib"
    const libcurl = "libcurl.4.dylib"
else
    const libnghttp2 = "libnghttp2.so"
    const libcurl = "libcurl.so"
end

function __init__()
    dlopen(libnghttp2)
    global libcurl_handle = dlopen(libcurl)
    global libcurl_path = dlpath(libcurl_handle)
end

end  # module LibCURL_jll

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia