text:
using Test
using LinearAlgebra
using KernelFunctions
using SpecialFunctions

x = rand()*2; v1 = rand(3); v2 = rand(3); id = IdentityTransform()
@testset "Kappa functions of kernels" begin
    @testset "Constant" begin
        @testset "ZeroKernel" begin
            k = ZeroKernel()
            @test eltype(k) == Any
            @test kappa(k,2.0) == 0.0
        end
        @testset "WhiteKernel" begin
            k = WhiteKernel()
            @test eltype(k) == Any
            @test kappa(k,1.0) == 1.0
            @test kappa(k,0.0) == 0.0
            @test EyeKernel == WhiteKernel
        end
        @testset "ConstantKernel" begin
            c = 2.0
            k = ConstantKernel(c=c)
            @test eltype(k) == Any
            @test kappa(k,1.0) == c
            @test kappa(k,0.5) == c
        end
    end
    @testset "Cosine" begin
        k = CosineKernel()
        @test eltype(k) == Any
        @test kappa(k, 1.0) ≈ -1.0 atol=1e-5
        @test kappa(k, 2.0) ≈ 1.0 atol=1e-5
        @test kappa(k, 1.5) ≈ 0.0 atol=1e-5
        @test kappa(k,x) ≈ cospi(x) atol=1e-5
        @test k(v1, v2) ≈ cospi(sqrt(sum(abs2.(v1-v2)))) atol=1e-5
    end
    @testset "Exponential" begin
        @testset "SqExponentialKernel" begin
            k = SqExponentialKernel()
            @test kappa(k,x) ≈ exp(-x)
            @test k(v1,v2) ≈ exp(-norm(v1-v2)^2)
            @test kappa(SqExponentialKernel(),x) == kappa(k,x)
        end
        @testset "ExponentialKernel" begin
            k = ExponentialKernel()
            @test kappa(k,x) ≈ exp(-x)
            @test k(v1,v2) ≈ exp(-norm(v1-v2))
            @test kappa(ExponentialKernel(),x) == kappa(k,x)
        end
        @testset "GammaExponentialKernel" begin
            γ = 2.0
            k = GammaExponentialKernel(γ=γ)
            @test kappa(k,x) ≈ exp(-(x)^(γ))
            @test k(v1,v2) ≈ exp(-norm(v1-v2)^(2γ))
            @test kappa(GammaExponentialKernel(),x) == kappa(k,x)
            @test GammaExponentialKernel(gamma=γ).γ == [γ]
            #Coherence :
            @test KernelFunctions._kernel(GammaExponentialKernel(γ=1.0),v1,v2) ≈ KernelFunctions._kernel(SqExponentialKernel(),v1,v2)
            @test KernelFunctions._kernel(GammaExponentialKernel(γ=0.5),v1,v2) ≈ KernelFunctions._kernel(ExponentialKernel(),v1,v2)
        end
    end
    @testset "Exponentiated" begin
        @testset "ExponentiatedKernel" begin
            k = ExponentiatedKernel()
            @test kappa(k,x) ≈ exp(x)
            @test kappa(k,-x) ≈ exp(-x)
            @test k(v1,v2) ≈ exp(dot(v1,v2))
        end
    end
    @testset "Matern" begin
        @testset "MaternKernel" begin
            ν = 2.0
            k = MaternKernel(ν=ν)
            matern(x,ν) = 2^(1-ν)/gamma(ν)*(sqrt(2ν)*x)^ν*besselk(ν,sqrt(2ν)*x)
            @test MaternKernel(nu=ν).ν == [ν]
            @test kappa(k,x) ≈ matern(x,ν)
            @test kappa(k,0.0) == 1.0
            @test kappa(MaternKernel(ν=ν),x) == kappa(k,x)
        end
        @testset "Matern32Kernel" begin
            k = Matern32Kernel()
            @test kappa(k,x) ≈ (1+sqrt(3)*x)exp(-sqrt(3)*x)
            @test k(v1,v2) ≈ (1+sqrt(3)*norm(v1-v2))exp(-sqrt(3)*norm(v1-v2))
            @test kappa(Matern32Kernel(),x) == kappa(k,x)
        end
        @testset "Matern52Kernel" begin
            k = Matern52Kernel()
            @test kappa(k,x) ≈ (1+sqrt(5)*x+5/3*x^2)exp(-sqrt(5)*x)
            @test k(v1,v2) ≈ (1+sqrt(5)*norm(v1-v2)+5/3*norm(v1-v2)^2)exp(-sqrt(5)*norm(v1-v2))
            @test kappa(Matern52Kernel(),x) == kappa(k,x)
        end
        @testset "Coherence Materns" begin
            @test kappa(MaternKernel(ν=0.5),x) ≈ kappa(ExponentialKernel(),x)
            @test kappa(MaternKernel(ν=1.5),x) ≈ kappa(Matern32Kernel(),x)
            @test kappa(MaternKernel(ν=2.5),x) ≈ kappa(Matern52Kernel(),x)
        end
    end
    @testset "Polynomial" begin
        c = randn();
        @testset "LinearKernel" begin
            k = LinearKernel()
            @test kappa(k,x) ≈ x
            @test k(v1,v2) ≈ dot(v1,v2)
            @test kappa(LinearKernel(),x) == kappa(k,x)
        end
        @testset "PolynomialKernel" begin
            k = PolynomialKernel()
            @test kappa(k,x) ≈ x^2
            @test k(v1,v2) ≈ dot(v1,v2)^2
            @test kappa(PolynomialKernel(),x) == kappa(k,x)
            #Coherence test
            @test kappa(PolynomialKernel(d=1.0,c=c),x) ≈ kappa(LinearKernel(c=c),x)
        end
    end
    @testset "RationalQuadratic" begin
        @testset "RationalQuadraticKernel" begin
            α = 2.0
            k = RationalQuadraticKernel(α=α)
            @test RationalQuadraticKernel(alpha=α).α == [α]
            @test kappa(k,x) ≈ (1.0+x/2.0)^-2
            @test k(v1,v2) ≈ (1.0+norm(v1-v2)^2/2.0)^-2
            @test kappa(RationalQuadraticKernel(α=α),x) == kappa(k,x)
        end
        @testset "GammaRationalQuadraticKernel" begin
            k = GammaRationalQuadraticKernel()
            @test kappa(k,x) ≈ (1.0+x^2.0/2.0)^-2
            @test k(v1,v2) ≈ (1.0+norm(v1-v2)^4.0/2.0)^-2
            @test kappa(GammaRationalQuadraticKernel(),x) == kappa(k,x)
            a = 1.0 + rand()
            @test GammaRationalQuadraticKernel(alpha=a).α == [a]
            #Coherence test
            @test kappa(GammaRationalQuadraticKernel(α=a,γ=1.0),x) ≈ kappa(RationalQuadraticKernel(α=a),x)
        end
    end
    @testset "Transformed/Scaled Kernel" begin
        s = rand()
        v = rand(3)
        k = SqExponentialKernel()
        kt = TransformedKernel(k,ScaleTransform(s))
        ktard = TransformedKernel(k,ARDTransform(v))
        ks = ScaledKernel(k,s)
        @test kappa(kt,v1,v2) == kappa(transform(k,ScaleTransform(s)),v1,v2)
        @test kappa(kt,v1,v2) == kappa(transform(k,s),v1,v2)
        @test kappa(kt,v1,v2) ≈ kappa(k,s*v1,s*v2) atol=1e-5
        @test kappa(ktard,v1,v2) ≈ kappa(transform(k,ARDTransform(v)),v1,v2) atol=1e-5
        @test kappa(ktard,v1,v2) == kappa(transform(k,v),v1,v2)
        @test kappa(ktard,v1,v2) == kappa(k,v.*v1,v.*v2)
        @test KernelFunctions.metric(kt) == KernelFunctions.metric(k)
        @test kappa(ks,x) == s*kappa(k,x)
        @test kappa(ks,x) == kappa(s*k,x)
    end
    @testset "KernelCombinations" begin
        k1 = LinearKernel()
        k2 = SqExponentialKernel()
        k3 = RationalQuadraticKernel()
        X = rand(2,2)
        @testset "KernelSum" begin
            w = [2.0,0.5]
            k = KernelSum([k1,k2],w)
            ks1 = 2.0*k1
            ks2 = 0.5*k2
            @test length(k) == 2
            @test kappa(k,v1,v2) == kappa(2.0*k1+0.5*k2,v1,v2)
            @test kappa(k+k3,v1,v2) ≈ kappa(k3+k,v1,v2)
            @test kappa(k1+k2,v1,v2) == kappa(KernelSum([k1,k2]),v1,v2)
            @test kappa(k+ks1,v1,v2) ≈ kappa(ks1+k,v1,v2)
            @test kappa(k+k,v1,v2) == kappa(KernelSum([k1,k2,k1,k2],vcat(w,w)),v1,v2)
        end
        @testset "KernelProduct" begin
            k = KernelProduct([k1,k2])
            @test length(k) == 2
            @test kappa(k,v1,v2) == kappa(k1*k2,v1,v2)
            @test kappa(k*k,v1,v2) ≈ kappa(k,v1,v2)^2
            @test kappa(k*k3,v1,v2) ≈ kappa(k3*k,v1,v2)
        end
    end
end


text:
n = [1,2,3,4]
k = ["a","b","c"]
path = Path(n, k)

# path to arcs
function test_path_to_arcs(path::Path)
    res = arcs(path)[1]
    for a in arcs(path)[2:end]
        @test res.v == a.u
        res = a
    end
    return nothing
end

test_path_to_arcs(path)

# path to graph
function test_path_to_graphs(path)
    G = Graph(path)
    for a in arcs(G)
        @test G[a,:flow] == 1
        @test G[a.u,a.v,a.key,:flow] == 1
    end
    return nothing
end

test_path_to_graphs(path)

@test nodelist(path) == path.nodelist
@test keylist(path) == path.keylist

G = Graph(6)
add_arcs_from!(G, [(1, 2), (2, 3), (3, 5), (3, 6), (4, 6), (6, 4), (2, 4)])
@test nodelist(depth_first_search(G, 1, 2)) == [1,2]
@test nodelist(depth_first_search(G, 1, 3)) == [1,2,3]
@test nodelist(depth_first_search(G, 1, 4)) == [1,2,4]
@test nodelist(depth_first_search(G, 1, 5)) == [1,2,3,5]
@test nodelist(depth_first_search(G, 1, 6)) == [1,2,4,6] ||
nodelist(depth_first_search(G, 1, 6)) == [1,2,3,6]

@test isempty(depth_first_search(G, 2, 1))
@test isempty(depth_first_search(G, 2, 1))
@test isempty(depth_first_search(G, 3, 1))
@test isempty(depth_first_search(G, 4, 1))
@test isempty(depth_first_search(G, 5, 1))
@test isempty(depth_first_search(G, 6, 1))

for a in arcs(G)
    G[a,:flow] = 1
end

@test nodelist(dijkstra_shortest_path(G, 1, 2)) == [1,2]
@test nodelist(dijkstra_shortest_path(G, 1, 3)) == [1,2,3]
@test nodelist(dijkstra_shortest_path(G, 1, 4)) == [1,2,4]
@test nodelist(dijkstra_shortest_path(G, 1, 5)) == [1,2,3,5]
@test nodelist(dijkstra_shortest_path(G, 1, 6)) == [1,2,4,6] ||
nodelist(dijkstra_shortest_path(G, 1, 6)) == [1,2,3,6]

@test isempty(dijkstra_shortest_path(G, 2, 1))
@test isempty(dijkstra_shortest_path(G, 2, 1))
@test isempty(dijkstra_shortest_path(G, 3, 1))
@test isempty(dijkstra_shortest_path(G, 4, 1))
@test isempty(dijkstra_shortest_path(G, 5, 1))
@test isempty(dijkstra_shortest_path(G, 6, 1))

@test dijkstra_shortest_path(G, 1)[1][6] == 3

ext: jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayPathRulePropertiesFormat <: SwaggerModel
    paths::Any # spec type: Union{ Nothing, Vector{String} } # spec name: paths
    backendAddressPool::Any # spec type: Union{ Nothing, SubResource } # spec name: backendAddressPool
    backendHttpSettings::Any # spec type: Union{ Nothing, SubResource } # spec name: backendHttpSettings
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState

    function ApplicationGatewayPathRulePropertiesFormat(;paths=nothing, backendAddressPool=nothing, backendHttpSettings=nothing, provisioningState=nothing)
        o = new()
        validate_property(ApplicationGatewayPathRulePropertiesFormat, Symbol("paths"), paths)
        setfield!(o, Symbol("paths"), paths)
        validate_property(ApplicationGatewayPathRulePropertiesFormat, Symbol("backendAddressPool"), backendAddressPool)
        setfield!(o, Symbol("backendAddressPool"), backendAddressPool)
        validate_property(ApplicationGatewayPathRulePropertiesFormat, Symbol("backendHttpSettings"), backendHttpSettings)
        setfield!(o, Symbol("backendHttpSettings"), backendHttpSettings)
        validate_property(ApplicationGatewayPathRulePropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ApplicationGatewayPathRulePropertiesFormat

const _property_map_ApplicationGatewayPathRulePropertiesFormat = Dict{Symbol,Symbol}(Symbol("paths")=>Symbol("paths"), Symbol("backendAddressPool")=>Symbol("backendAddressPool"), Symbol("backendHttpSettings")=>Symbol("backendHttpSettings"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ApplicationGatewayPathRulePropertiesFormat = Dict{Symbol,String}(Symbol("paths")=>"Vector{String}", Symbol("backendAddressPool")=>"SubResource", Symbol("backendHttpSettings")=>"SubResource", Symbol("provisioningState")=>"String")
Base.propertynames(::Type{ ApplicationGatewayPathRulePropertiesFormat }) = collect(keys(_property_map_ApplicationGatewayPathRulePropertiesFormat))
Swagger.property_type(::Type{ ApplicationGatewayPathRulePropertiesFormat }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ApplicationGatewayPathRulePropertiesFormat[name]))}
Swagger.field_name(::Type{ ApplicationGatewayPathRulePropertiesFormat }, property_name::Symbol) =  _property_map_ApplicationGatewayPathRulePropertiesFormat[property_name]

function check_required(o::ApplicationGatewayPathRulePropertiesFormat)
    true
end

function validate_property(::Type{ ApplicationGatewayPathRulePropertiesFormat }, name::Symbol, val)
end


text:
# We define a simple Gaussian model with unknown mean and variance.

# \toc

# ## Model

using Turing

@model function gdemo(x, y)
  s ~ InverseGamma(2, 3)
  m ~ Normal(0, sqrt(s))
  x ~ Normal(m, sqrt(s))
  y ~ Normal(m, sqrt(s))
end;

# ## Output

# and run the sampler:

chains = sample(gdemo(1.5, 2), NUTS(0.65), 1000)

# \defaultoutput{}


text:
"""
Created on Mon Apr 13 2020

@author: Yoann Pradat

Code supporting presentation on NMF.
"""

include("../../../src/VNMF.jl")
using .VNMF
using Distributions
using Random
include("../util/plot_cost.jl")

gr() # set the backend to GR

#### # 1. SIMULATE DATA
#### # ################################################################################################################

#### simulate matrix for testing NMF
F = 10
N = 25
K = 5
rng = MersenneTwister(123)

#### truncated-normal factors with noise
μ = 2
σ = √2
dist_W = truncated(Normal(μ, σ), 0, Inf)
dist_H = truncated(Normal(μ, σ), 0, Inf)

dist_E = Normal(μ, 0.2*σ)
dist_E_type = "Normal"

# α = μ
# θ = 1
# dist_E = Gamma(α, θ)
# dist_E_type = "Gamma"

W = rand(rng, dist_W, F, K)
H = rand(rng, dist_H, K, N)
E = rand(rng, dist_E, F, N)

V = W * H + E

fo = "../../plot/simulation/nmf_noise"
mkpath(fo)

#### global params
global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0,1),
    rank          = K,
    max_iter      = 20_000,
    stopping_crit = :none,
    verbose       = true,
)

#### # 2. FEVOTTE-IDIER NMF ALGORITHMS WITH KNOWN RANK
#### # ################################################################################################################

list_res = Array{NMFResults{Float64}, 1}()
list_lab = Array{String, 1}()
list_ls  = Array{Symbol, 1}()
list_col = Array{String, 1}()

local_params = FIParams(
    scale_W = true,
    α_H     = 0,
    α_W     = 0,
)

####
#### 2.1 β divergences, FI algorithms
####

for (β, col) in zip([2, 1, 0], ["black", "red", "green"])
    local_params.β = β

    run_mm = true
    run_me = true
    run_h  = true

    if β == 2
        run_h = false
    elseif β == 1
        run_h  = false
        run_me = false
    end

    if run_mm
        local_params.alg = :mm

        push!(list_res, nmf_FI(V, global_params, local_params))
        push!(list_lab, raw"$\mathrm{MM} \quad \beta =" * "$(β)" * raw"$")
        push!(list_ls, :dash)
        push!(list_col, col)
    end

    if run_me
        local_params.alg = :me

        push!(list_res, nmf_FI(V, global_params, local_params))
        push!(list_lab, raw"$\mathrm{ME} \quad \beta =" * "$(β)" * raw"$")
        push!(list_ls, :solid)
        push!(list_col, col)
    end

    if run_h
        local_params.alg = :h

        push!(list_res, nmf_FI(V, global_params, local_params))
        push!(list_lab, raw"$\mathrm{H} \quad \beta =" * "$(β)" * raw"$")
        push!(list_ls, :dot)
        push!(list_col, col)
    end

end

####
#### 2.2 PLOT
####

fn = joinpath(fo, "nmf_FI_noise_$(dist_E_type).pdf")
plot_cost(list_res, list_lab, list_ls, list_col, fn, relative=true)

#### # 3. MU NMF ALGORITHMS WITH KNOWN RANK
#### # ################################################################################################################

list_res = Array{NMFResults{Float64}, 1}()
list_lab = Array{String, 1}()
list_ls  = Array{Symbol, 1}()
list_col = Array{String, 1}()

####
#### 2.1 β divergences
####

for (β, col) in zip([2, 1, 0], ["black", "red", "green"])
    local_params = MUParams(
        β       = β,
        scale_W = true,
        div     = :β,
        alg     = :mu
    )

    push!(list_res, nmf_MU(V, global_params, local_params))
    push!(list_lab, raw"$\mathrm{MU} \quad \beta =" * "$(β)" * raw"$")
    push!(list_ls, :solid)
    push!(list_col, col)
end

####
#### 2.2 α divergences
####

for (α, col) in zip([2, 1, 0], ["black", "red", "green"])
    local_params = MUParams(
        α       = α,
        scale_W = true,
        div     = :α,
        alg     = :mu
    )

    push!(list_res, nmf_MU(V, global_params, local_params))
    push!(list_lab, raw"$\mathrm{MU} \quad \alpha =" * "$(α)" * raw"$")
    push!(list_ls, :dash)
    push!(list_col, col)
end

####
#### 2.3 PLOT
####

fn = joinpath(fo, "nmf_MU_noise_$(dist_E_type).pdf")
plot_cost(list_res, list_lab, list_ls, list_col, fn, relative=true)

#### # Cichocki NMF algorithms
#### #################################################################################################################

#params = CParams(
#    α             = 2,
#    max_iter      = convert(Int, 1e4),
#    stopping_crit = :none,
#    stopping_tol  = 1e-6,
#    α_H           = 0,
#    α_W           = 0,
#    l₁ratio_H     = 1,
#    l₁ratio_W     = 1,
#    μ_H           = 0,
#    μ_W           = 0,
#    α_sparse_H    = 0,
#    α_sparse_W    = 0,
#    ω             = 1,
#    verbose       = true,
#)
#
#
#res_α = nmf_α(
#    V,
#    W_init,
#    H_init,
#    params
#)


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _precess{T<:AbstractFloat}(ra::T, dec::T, equinox1::T, equinox2::T,
                                    FK4::Bool, radians::Bool)
    if radians
        ra_rad  = ra
        dec_rad = dec
    else
        # Convert to radians.
        ra_rad  = deg2rad(ra)
        dec_rad = deg2rad(dec)
    end
    a = cos(dec_rad)
    x = [a*cos(ra_rad), a*sin(ra_rad), sin(dec_rad)]
    r = premat(equinox1, equinox2, FK4=FK4)
    x2 = r*x
    ra_rad  = atan2(x2[2], x2[1])
    dec_rad = asin(x2[3])
    if radians
        return cirrange(ra_rad, 2.0*pi), dec_rad
    else
        return cirrange(rad2deg(ra_rad)), rad2deg(dec_rad)
    end
end

"""
    precess(ra, dec, equinox1, equinox2[, FK4=true, radians=true]) -> prec_ra, prec_dec

### Purpose ###

Precess coordinates from `equinox1` to `equinox2`.

### Explanation ###

The default `(ra, dec)` system is FK5 based on epoch J2000.0 but FK4 based on
B1950.0 is available via the `FK4` boolean keyword.

### Arguments ###

* `ra`: input right ascension, scalar or vector, in degrees, unless the
  `radians` keyword is set to `true`
* `dec`: input declination, scalar or vector, in degrees, unless the `radians`
  keyword is set to `true`
* `equinox1`: original equinox of coordinates, numeric scalar.
* `equinox2`: equinox of precessed coordinates.
* `FK4` (optional boolean keyword): if this keyword is set to `true`, the FK4
  (B1950.0) system precession angles are used to compute the precession matrix.
  When it is `false`, the default, use FK5 (J2000.0) precession angles.
* `radians` (optional boolean keyword): if this keyword is set to `true`, then
  the input and output right ascension and declination vectors are in radians
  rather than degrees.

### Output ###

The 2-tuple `(ra, dec)` of coordinates modified by precession.

### Example ###

The Pole Star has J2000.0 coordinates (2h, 31m, 46.3s, 89d 15' 50.6"); compute
its coordinates at J1985.0

``` julia
ra, dec = ten(2,31,46.3)*15, ten(89,15,50.6)
# => (37.94291666666666,89.26405555555556)
adstring(precess(ra, dec, 2000, 1985), precision=1)
# => " 02 16 22.73  +89 11 47.3"
```

Precess the B1950 coordinates of Eps Ind (RA = 21h 59m,33.053s, DEC = (-56d,
59', 33.053") to equinox B1975.

``` julia
ra, dec = ten(21, 59, 33.053)*15, ten(-56, 59, 33.053)
# => (329.88772083333333,-56.992514722222225)
adstring(precess(ra, dec, 1950, 1975, FK4=true), precision=1)
# => " 22 01 15.46  -56 52 18.7"
```

### Method ###

Algorithm from "Computational Spherical Astronomy" by Taff (1983),
p. 24. (FK4). FK5 constants from "Explanatory Supplement To The Astronomical
Almanac" 1992, page 104 Table 3.211.1
(https://archive.org/details/131123ExplanatorySupplementAstronomicalAlmanac).

### Notes ###

Accuracy of precession decreases for declination values near 90 degrees.
`precess` should not be used more than 2.5 centuries from 2000 on the FK5 system
(1950.0 on the FK4 system).  If you need better accuracy, use `bprecess` or
`jprecess` as needed.

Code of this function is based on IDL Astronomy User's Library.
"""
precess(ra::Real, dec::Real, equinox1::Real, equinox2::Real;
        FK4::Bool=false, radians::Bool=false) =
            _precess(promote(float(ra), float(dec),
                             float(equinox1), float(equinox2))...,
                     FK4, radians)

precess(radec::Tuple{Real, Real}, equinox1::Real, equinox2::Real;
        FK4::Bool=false,
        radians::Bool=false) =
            precess(radec[1], radec[2], equinox1, equinox2,
                    FK4=FK4, radians=radians)

function precess{R<:Real, D<:Real}(ra::AbstractArray{R},
                                   dec::AbstractArray{D}, equinox1::Real,
                                   equinox2::Real; FK4::Bool=false,
                                   radians::Bool=false)
    @assert length(ra) == length(dec) "ra and dec arrays should be of the same length"
    typera = typeof(float(one(R)))
    ra_out  = similar(ra,  typera)
    dec_out = similar(dec, typera)
    for i in eachindex(ra)
        ra_out[i], dec_out[i] = precess(ra[i], dec[i], equinox1, equinox2,
                                        FK4=FK4, radians=radians)
    end
    return ra_out, dec_out
end


text:
###########
# Generic methods for obtaining Jacobi matrices with one of the a, b and c parameters lowered by 1.
# passing symbols :a, :b or :c into lowindex determines which parameter is lowered
######
function initialα_gen(P::SemiclassicalJacobi, lowindex::Symbol)
    t = P.t; a = P.a; b = P.b; c = P.c;
    A,B,_ = recurrencecoefficients(P)
    # this makes use of Q_1 == A[1]x + B[1]
    if lowindex == :a
        return -(A[1]*sum(SemiclassicalJacobiWeight(t, a, b, c)) + B[1]*sum(SemiclassicalJacobiWeight(t, a-1, b, c)))/(sum(SemiclassicalJacobiWeight(t, a-1, b, c)))
    elseif lowindex == :b
        return -(A[1]*sum(SemiclassicalJacobiWeight(t, a+1, b-1, c)) + B[1]*sum(SemiclassicalJacobiWeight(t, a, b-1, c)))/(sum(SemiclassicalJacobiWeight(t, a, b-1, c)))
    else # lowindex == :c
        return -(A[1]*sum(SemiclassicalJacobiWeight(t, a+1, b, c-1)) + B[1]*sum(SemiclassicalJacobiWeight(t, a, b, c-1)))/(sum(SemiclassicalJacobiWeight(t, a, b, c-1)))
    end
end

mutable struct LoweringCoefficients{T} <: AbstractCachedVector{T}
    P::SemiclassicalJacobi{T}
    lowind::Symbol # options for lowering are :a, :b and :c
    data::Vector{T}
    datasize::Tuple{Int}
    LoweringCoefficients{T}(P::SemiclassicalJacobi{T}, lowindex::Symbol) where T = new{T}(P, lowindex, [initialα_gen(P, lowindex)], (1,))
end
LoweringCoefficients(P::SemiclassicalJacobi{T}, lowindex::Symbol) where T = LoweringCoefficients{T}(P, lowindex)
size(::LoweringCoefficients) = (ℵ₀,)

function getindex(α::LoweringCoefficients, I::UnitRange{Int})
    resizedata!(α, maximum(I))
    return getindex(α.data,I)
end

function getindex(α::LoweringCoefficients, I::Int)
    resizedata!(α, I)
    return getindex(α.data,I)
end

function resizedata!(α::LoweringCoefficients, nm::Int)
    νμ = length(α.data)
    if nm > νμ
        olddata = copy(α.data)
        nmax = maximum(nm)
        α.data = similar(olddata,nmax)
        α.data[1:νμ] = olddata[1:νμ]
        inds = νμ:nmax
        cache_filldata!(α, inds)
        α.datasize = (nmax,)
    end
    return α
end
cache_filldata!(α::LoweringCoefficients, inds::UnitRange{Int}) =  αgenfillerbackwards!(α.data, 100, α.P, α.lowind, inds)

# fill in coefficients via Miller recurrence
function αgenfillerbackwards!(α::Vector{T}, addscale::Int, P::SemiclassicalJacobi{T}, lowindex::Symbol, inds::UnitRange{Int}) where T
    maxI::Int = maximum(inds)
    minI::Int = minimum(inds)
    oldval = α[minI]
    n = maxI
    k = one(T)
    k2 = zero(T)
    if lowindex == :a
        _,B,C = recurrencecoefficients(P)
        while abs(k2-k) > 1e-10
            k2, k = k, one(T)
            n += addscale
            @inbounds for j = n:-1:maxI
                k = -C[j+1]/(k+B[j+1])
            end
        end
        α[maxI] = k
        @inbounds for j = maxI:-1:minI+1
            α[j-1] = -C[j]/(α[j]+B[j])
        end
    elseif lowindex == :b
        A,B,C = recurrencecoefficients(P)
        while abs(k2-k) > 1e-10
            k2, k = k, one(T)
            n += addscale
            @inbounds for j = n:-1:maxI
                k = -C[j+1]/(k+A[j+1]+B[j+1])
            end
        end
        α[maxI] = k
        @inbounds for j = maxI:-1:minI+1
            α[j-1] = -C[j]/(α[j]+A[j]+B[j])
        end
    else # lowindex == :c
        A,B,C = recurrencecoefficients(P)
        while abs(k2-k) > 1e-10
            k2, k = k, one(T)
            n += addscale
            @inbounds for j = n:-1:maxI
                k = -C[j+1]/(k+A[j+1]*P.t+B[j+1])
            end
        end
        α[maxI] = k
        @inbounds for j = maxI:-1:minI+1
            α[j-1] = -C[j]/(α[j]+A[j]*P.t+B[j])
        end
    end
    α[minI:maxI] = ((oldval)/α[inds[1]]).*α[minI:maxI]
    α
end

function symlowered_jacobimatrix(Q::SemiclassicalJacobi, lowindex::Symbol)
    bands = SemiclassicalLoweredJacobiBands(Q,lowindex)
    return SymTridiagonal(bands[1,:],bands[2,:])
end

mutable struct SemiclassicalLoweredJacobiBands{T} <: AbstractCachedMatrix{T}
    P::SemiclassicalJacobi{T}
    data::Array{T} # to avoid redundant re-computations, the bands are stored as a (∞,2)-array
    αcfs::LoweringCoefficients{T}
    datasize::Tuple{Int,Int}
end
size(r::SemiclassicalLoweredJacobiBands) = (2,ℵ₀)

function SemiclassicalLoweredJacobiBands(P::SemiclassicalJacobi{T}, lowindex::Symbol) where T
    cachedα = LoweringCoefficients(P,lowindex)
    a, b, c, t, αcfs = P.a, P.b, P.c, P.t, cachedα[1]
    C, A = P.X[2,1], P.X[1,1]
    if lowindex == :a
        SemiclassicalLoweredJacobiBands{T}(P, [A-C*αcfs, C/(sqrt(sum(SemiclassicalJacobiWeight(t,a-1,b,c)))/sqrt((((αcfs^2-2*αcfs*A/C+A^2/C^2)*sum(SemiclassicalJacobiWeight(t, a-1, b, c))+(2*αcfs/C-2*A/C^2)*sum(SemiclassicalJacobiWeight(t, a, b, c))))+(1/C^2)*(sum(SemiclassicalJacobiWeight(t, a+1, b, c)))))], cachedα, (2,1))
    elseif lowindex == :b
        SemiclassicalLoweredJacobiBands{T}(P, [A-C*αcfs, C/(sqrt(sum(SemiclassicalJacobiWeight(t,a,b-1,c)))/sqrt((((αcfs^2-2*αcfs*A/C+A^2/C^2)*sum(SemiclassicalJacobiWeight(t, a, b-1, c))+(2*αcfs/C-2*A/C^2)*sum(SemiclassicalJacobiWeight(t, a+1, b-1, c))))+(1/C^2)*(sum(SemiclassicalJacobiWeight(t, a+2, b-1, c)))))], cachedα, (2,1))
    elseif lowindex == :c
        SemiclassicalLoweredJacobiBands{T}(P, [A-C*αcfs, C/(sqrt(sum(SemiclassicalJacobiWeight(t,a,b,c-1)))/sqrt((((αcfs^2-2*αcfs*A/C+A^2/C^2)*sum(SemiclassicalJacobiWeight(t, a, b, c-1))+(2*αcfs/C-2*A/C^2)*sum(SemiclassicalJacobiWeight(t, a+1, b, c-1))))+(1/C^2)*(sum(SemiclassicalJacobiWeight(t, a+2, b, c-1)))))], cachedα, (2,1))
    end
end

function getindex(r::SemiclassicalLoweredJacobiBands, I::Vararg{Int,2})
    r.αcfs[maximum(I)+1] # expanding the cache of coefficients here prevents redundant recomputation later
    resizedata!(r,maximum(I))
    getindex(r.data,I[1],I[2])
end
function getindex(r::SemiclassicalLoweredJacobiBands, I::Int, J::UnitRange{Int})
    r.αcfs[maximum(J)+1] # expanding the cache of coefficients here prevents redundant recomputation later
    resizedata!(r,maximum(J))
    getindex(r.data,I,J)
end
function resizedata!(r::SemiclassicalLoweredJacobiBands, nm::Int)
    νμ = length(r.data[1,:])
    if nm > νμ
        olddata = copy(r.data)
        r.data = similar(olddata,2,nm)
        r.data[1:2,1:νμ] = olddata
        inds = νμ:nm
        cache_filldata!(r, inds)
        r.datasize = (2,nm)
    end
    return r
end
cache_filldata!(r::SemiclassicalLoweredJacobiBands, inds::UnitRange{Int}) =  loweringjacobibandfill!(r.data, r.αcfs, r.P, inds)

function loweringjacobibandfill!(data::Array, αcfs::LoweringCoefficients, P::SemiclassicalJacobi, inds::UnitRange{Int})
    lowI = minimum(inds):maximum(inds)-1
    shiftI = minimum(inds)+1:maximum(inds)
    C,A = subdiagonaldata(P.X), diagonaldata(P.X)
    data[1,shiftI] = -C[shiftI].*αcfs[shiftI]+C[lowI].*αcfs[lowI]+A[shiftI]
    data[2,shiftI] = sqrt.(C[lowI].*C[shiftI].*αcfs[shiftI]./αcfs[lowI])
    data
end

###########
# Methods for the special case of computing SemiclassicalJacobi(t,0,0,-1)
######
# As this can be built from SemiclassicalJacobi(t,0,0,0) which is just shifted and normalized Legendre(), we have more explicit methods at our disposal due to explicitly known coefficients for the Legendre bases.

###
#   α coefficients implementation, cached, direct and via recurrences
###
# inital value n=0 for α_{n,n-1}(t) coefficients
initialα(t) = t-2/(log1p(t)-log(t-1))

# compute n-th coefficient from direct evaluation formula
αdirect(n, t) = 2*n*_₂F₁general2((1+n)/2,(n+2)/2,(2*n+3)/2,1/t^2)/(t*2*(1+2*n)*_₂F₁general2(n/2,(n+1)/2,(2*n+1)/2,1/t^2))
# this version takes a pre-existing vector v and fills in the missing data guided by indices in inds using explicit formula
function αdirect!(α, t, inds)
    @inbounds for n in inds
        α[n] = 2*n*_₂F₁general2((1+n)/2,(2+n)/2,(2*n+3)/2,1/t^2)/(t*2*(1+2*n)*_₂F₁general2(n/2,(n+1)/2,(2*n+1)/2,1/t^2))
    end
end

# takes a previously computed vector of α_{n,n-1}(t) that has been increased in size and fills in the missing data guided by indices in inds via forward recurrence
function αcoefficients!(α, t, inds)
    @inbounds for n in inds
       α[n] = (t*(2*n-1)/n-(n-1)/(n*α[n-1]))
    end
end

# takes a previously computed vector of α_{n,n-1}(t) that has been increased in size and fills in the missing data guided by indices in inds, using a final condition computation and backward recurrence
function backαcoeff!(α, t, inds)
    α[end] = αdirect(length(α),t)
    @inbounds for n in reverse(inds)[1:end-1]
       α[n-1] = (n-1)/(t*(2*n-1)-n*α[n])
    end
end

# cached implementation using stable back recurrence to fill data
mutable struct neg1c_αcfs{T} <: AbstractCachedVector{T}
    t::T
    data::Vector{T}
    datasize::Tuple{Int}
    neg1c_αcfs{T}(t::T) where T = new{T}(t, [initialα(2*t-1),αdirect(2,2*t-1)], (2,))
end
neg1c_αcfs(t::T) where T = neg1c_αcfs{T}(t)
size(α::neg1c_αcfs) = (ℵ₀,)
cache_filldata!(α::neg1c_αcfs, inds) = backαcoeff!(α.data, 2*α.t-1, inds)

function resizedata!(α::neg1c_αcfs, nm)
    νμ = length(α.data)
    if nm > νμ
        olddata = copy(α.data)
        α.data = similar(olddata,maximum(nm))
        α.data[1:νμ] = olddata[1:νμ]
        inds = Array(νμ:maximum(nm))
        cache_filldata!(α, inds)
        α.datasize = (nm,)
    end
    α
end

###
#   cached implementation of normalization constants
###
mutable struct neg1c_normconstant{T} <: AbstractCachedVector{T}
    t::T
    data::Vector{T}
    datasize::Tuple{Int}
    neg1c_normconstant{T}(t::T) where T = new{T}(t, neg1c_normconstinitial(t, 10), (10,))
end
neg1c_normconstant(t::T) where T = neg1c_normconstant{T}(t)
size(B::neg1c_normconstant) = (ℵ₀,)
cache_filldata!(B::neg1c_normconstant, inds) = neg1c_normconstextension!(B.data, inds, B.t)

function resizedata!(B::neg1c_normconstant, nm)
    νμ = length(B.data)
    if nm > νμ
        olddata = copy(B.data)
        B.data = similar(olddata,maximum(nm))
        B.data[1:νμ] = olddata[1:νμ]
        inds = Array(νμ-1:maximum(nm))
        cache_filldata!(B, inds)
        B.datasize = (nm,)
    end
    B
end

function neg1c_normconstinitial(t::T, N::Integer) where T
    # generate α coefficients for OPs via recurrence
    α = zeros(T,N-1)
    α[1] = initialα(2*t-1)
    backαcoeff!(α,2*t-1,(2:N-1))
    # normalization constants
    B = [one(T)]
    append!(B,sqrt(2*acoth(2*t-1))*sqrt.((1:N-1)./(2 .*α)))
    return B
end

function neg1c_normconstextension!(B::Vector, inds, t::T) where T
    n = maximum(inds)
    m = minimum(inds)-1
    # generate missing α coefficients
    α = zeros(T, n)
    backαcoeff!(α,2*t-1,(m:n))
    # normalization constants
    norm = sqrt(2*acoth(2*t-1))*sqrt.((m:n)./(2 .*α[m:n]))
    B[m+1:n] = norm[1:end-1]
    B
end

###
#   bidiagonal operator which converts from OPs wrt (t-x)^-1 to shifted Legendre
###
function neg1c_tolegendre(t::T) where T
    norm = neg1c_normconstant(t)
    α = neg1c_αcfs(t)
    supD = norm.*Vcat(zeros(T,1),-α)
    D = norm
    return Bidiagonal(D,supD[2:end],:U)
end

###
#   explicit polynomial evaluations
###
# Evaluate the n-th normalized OP wrt 1/(t-x) at point x
function evalQn(n::Integer, x, t::T) where T
    # this version recomputes α based on t
    t <= 1 && throw(ArgumentError("t must be greater than 1."))
    n == 0 && return one(T)
    α = αdirect(n,2*t-1)
    return sqrt(2*acoth(2*t-1))*(-α*jacobip(n-1,0,0,2*x-1)+jacobip(n,0,0,2*x-1))*sqrt(n/(2*α))
end
# Evaluate the n-th non-normalized OP wrt 1/(t-x) at point x
function evalϕn(n::Integer, x, t::T) where T
    # this version recomputes α based on t
    t <= 1 && throw(ArgumentError("t must be greater than 1."))
    n == 0 && return one(T)
    return -αdirect(n,2*t-1)*jacobip(n-1,0,0,2*x-1)+jacobip(n,0,0,2*x-1)
end


text:
using UnPack

const SPECIES_NAMES = ["taxon", "species"]
const FLOAT_TYPES = [Real, Union{Missing, <:Real}]
const TAXON_NAME = "taxon"

function process_column(col::AbstractVector{<:Real})::Vector{Float64}
    n = length(col)
    return [convert(Float64, x) for x in col]
end

function process_column(col::AbstractVector{Union{Missing, T}})::Vector{Float64} where T <: Real
    n = length(col)
    v = zeros(n)
    for i = 1:n
        if ismissing(col[i])
            v[i] = NaN
        else
            v[i] = convert(Float64, col[i])
        end
    end
    return v
end

function parse_extension(path::AbstractString)
    last_dot = findlast(isequal('.'), path)
    ext = path[(last_dot + 1):end]
    return ext
end

function parse_continuous_data(datapath::AbstractString)::DataFrame
    ext = parse_extension(datapath)
    if ext in keys(FILE_EXTENSIONS)
        f = FILE_EXTENSIONS[ext]
        return f(datapath)
    end
end

function parse_continuous_data_csv(datapath::AbstractString)
    df = DataFrame(CSV.File(datapath))
    processed_df = parse_continuous_df(df)

    return processed_df
end

const TAXA = "TAXA"
const CHARACTERS = "CHARACTERS"
const DATATYPE = "DATATYPE"

function contains_taxon(line::AbstractString)
    line = strip(line)
    return !isempty(line) && line != ";"
end

function parse_taxon_values(s::AbstractString)
    comps = split(strip(s))
    taxon = comps[1]
    values = missing_parse.(comps[2:end])
    return (taxon = taxon, values = values)
end

function missing_parse(s::AbstractString)
    if s in ["?"]
        return NaN
    else
        return parse(Float64, s)
    end
end

# function parse_taxa_nexus(lines::Vector{<:AbstractString}, ind::Int)
#     @assert strip(lines[ind]) == nexus_begin(TAXA)

# function parse_continuous_data_nexus(datapath::AbstractString; title::String = "")
#     lines = readlines(data_path)
#     char_start = -1
#     char_begin = nexus_begin(CHARACTERS)
#     taxa_begin = nexus_begin(TAXA)

#     i = 1
#     while i <= length(lines)
#         line = strip(lines[i])
#         if line == char_begin
#             @unpack title, parse_nexus_character(lines, i)
#         else
#             i+= 1
#         end
# end

const NEXUS_COMMENT = r"\[.*?\]"
const NEXUS_BLOCK = r"\s*BEGIN\s+(.*?);\s*(.*?)\s*END;"s
const NEXUS_TAXA = r"\s*DIMENSIONS\s+NTAX\s*=\s*(\d+?)\s*;.*TAXLABELS\s*(.*?)\s*;"s
const NEXUS_TITLE = r"\s*TITLE\s+(.*?)\s*;"
const NEXUS_DIM = r"DIMENSIONS\s+NCHAR\s*=\s*(\d+)\s*;"
const NEXUS_CHARLABELS = r"\s+CHARSTATELABELS\s+(.*?);"s
const NEXUS_MATRIX = r"\s+MATRIX\s+(.*?)\s+;"s

function nexus_format_regex(s::AbstractString)
    space_or_other = "(?:\\s+|\\s+.+\\s+)"
    r = "FORMAT" * space_or_other * s * "\\s*=\\s*(.*?)(?:\\s+|;)"
    return Regex(r)
end

function parse_characters_nexus(s::AbstractString, n_taxa::Int)
    title = ""
    title_match = match(NEXUS_TITLE, s)
    if !isnothing(title_match)
        title = title_match[1]
    end

    dim = parse(Int, match(NEXUS_DIM, s)[1])

    datatype = match(nexus_format_regex(DATATYPE), s)[1]

    labels_string = match(NEXUS_CHARLABELS, s)[1]
    labels_and_inds = strip.(split(labels_string, ','))
    @assert length(labels_and_inds) == dim
    labels = fill("", dim)
    for i = 1:dim
        split_label = split(labels_and_inds[i])
        @assert parse(Int, split_label[1]) == i
        label = strip(join(split_label[2:end], '_'), ['\'','"'])
        labels[i] = label
    end

    taxa = fill("", n_taxa)
    data = fill(NaN, n_taxa, dim)
    data_lines = split(match(NEXUS_MATRIX, s)[1], '\n')
    @assert length(data_lines) == n_taxa

    for i = 1:n_taxa
        @unpack taxon, values = parse_taxon_values(data_lines[i])
        taxa[i] = taxon
        data[i, :] .= values
    end

    df = [DataFrame(taxon = taxa) DataFrame(data, labels)]

    return (title = title, datatype = datatype, data = df)
end

function split_by_pattern(s::AbstractString, r::Regex)
    matches = RegexMatch[]
    between = String[]

    m = match(r, s)
    start = 1
    while !isnothing(m)
        between_inds = start:(m.offset - 1)
        push!(between, s[between_inds])
        push!(matches, m)

        s = @view s[(m.offset + length(m.match)):end]

        m = match(r, s)
    end
    push!(between, s)
    return (matches = matches, between = between)
end

function remove_comments_nexus(s::AbstractString)
    @unpack between = split_by_pattern(s, NEXUS_COMMENT)
    return join(between)
end

function parse_continuous_data_nexus(path::String)
    s = read(path, String)
    s = remove_comments_nexus(s)

    @unpack matches = split_by_pattern(s, NEXUS_BLOCK)

    taxon_ind = findall(x -> x[1] == TAXA, matches)
    @assert length(taxon_ind) == 1
    taxa = parse_nexus_taxa(matches[taxon_ind[1]][2])
    n_taxa = length(taxa)

    character_blocks = findall(x -> x[1] == CHARACTERS, matches)

    df = DataFrame(taxon = taxa)

    for block in @view matches[character_blocks]
        content = block[2]
        @unpack datatype, title, data = parse_characters_nexus(content, n_taxa)
        @assert df.taxon == data.taxon
        df = [df data[!, 2:end]]
    end

    return df
end

function parse_nexus_taxa(s::AbstractString)
    m = match(NEXUS_TAXA, s)
    n_taxa = parse(Int, m[1])
    taxa = split(m[2])
    if length(taxa) != n_taxa
        error("could not parse nexus taxa")
    end
    return taxa
end

# function parse_first_block(s::AbstractString)

# function parse_taxa_block(lines::Vector{<:AbstractString})

function parse_continuous_df(df::DataFrame)

    check_valid(df)
    colnames = names(df)

    taxa = df[!, taxon_inds[1]]
    data_inds = Int[]

    processed_df = DataFrame()
    processed_df[TAXON_NAME] = taxa

    nms = names(df)

    for i = 1:size(df, 2)
        col = df[!, i]
        for type in FLOAT_TYPES
            if eltype(col) <: type
                processed_df[nms[i]] = process_column(col)
            end
        end
    end

    return processed_df
end

function check_valid(df::DataFrame)
    colnames = names(df)
    taxon_inds = findall(x -> lowercase(x) in SPECIES_NAMES, colnames)
    if length(taxon_inds) > 1
        error("cannot determine taxon column. multiple matches.")
    elseif length(taxon_inds) == 0
        error("not a valid trait data set. cannot find list of taxa.")
    end
    return colnames[taxon_inds[1]]
end

const FILE_EXTENSIONS = Dict{String, Function}(
                            "csv" => parse_continuous_data_csv,
                            "nex" => parse_continuous_data_nexus
                            )


text:
net = init_network([2, 2, 2])

@assert abs(net.weights[1][1]) > 0.0001

# reset "random" weights so we get predictable results
net.weights = Array{Float64}[
  [
    -0.42700000000000005 0.7829999999999999;
    0.8089999999999999 0.02299999999999991;
    -0.06999999999999995 0.24;
  ],
  [
    -0.7010000000000001 -0.783;
    -0.917 0.4099999999999999;
    -0.834 0.9550000000000001;
  ]]

@assert train(net, [0.15,0.7],[0.1,0.9]) == 0.026566707596555807

# review data
@assert net.weights == Array{Float64}[
[
  -0.42731782254634215 0.7832271291136619;
  0.8075168281170702 0.02405993586375571;
  -0.07211881697561397 0.24151419409107971;
],
[
  -0.7017660587853684 -0.7754389164772143;
  -0.9177482539018814 0.4173853473864355;
  -0.8352632670143251 0.9674685828154825;
]]

@assert net.activation_nodes == Array{Float64}[
  [0.15,0.7,1.0],
  [0.6064108194715588,0.5923165042674361,1.0],
  [0.14157761100099062,0.6732739153591676]
]

@assert net.deltas == Array{Float64}[
  [-0.00847526790245606 0.006056776364318863 0.0],
  [-0.005053068057300355 0.049874331261929916]
]

@assert net.last_changes == Array{Float64}[
  [
    -0.0012712901853684088 0.0009085164546478294;
    -0.005932687531719241 0.004239743455023204;
    -0.00847526790245606 0.006056776364318863;
  ],
  [
    -0.0030642351414730666 0.030244334091142907;
    -0.002993015607525591 0.02954138954574243;
    -0.005053068057300355 0.049874331261929916;
  ]
]


text:
using Base.Threads
using DataFrames
using Discretizers
using Eolas
using LightGraphs
using LightGraphs.SimpleGraphs
using Random
using XLSX
using Plots
using MLBase

normalize(xs) = xs ./ sum(xs)

function filtrate!(frames, col)
    for (i, frame) in enumerate(frames[2:end])
        frames[1 + i] = filter(r -> all(g -> r[col] in g[:,col], frames[1:i]), frame)
    end
    frames
end

function kos()
    excel = XLSX.openxlsx("kos.xlsx")
    sheets = XLSX.sheetnames(excel)[2:4]
    frames = map(s -> DataFrame(XLSX.gettable(excel[s])...), sheets)
    sort!.(dropmissing!.(select!.(frames, :KO), disallowmissing=true))
end

function bootstrap(df, col, n)
    gf = DataFrame(Symbol(string(col) * "0") => df[:,col])
    for i in 1:n
        gf[!,Symbol(string(col) * string(i))] = sort!(rand(df[:,col], length(df[:,col])))
    end
    gf
end

function genus()
    excel = XLSX.openxlsx("genus.xlsx")
    sheets = XLSX.sheetnames(excel)[2:4]
    map(sheets) do sheet
        column = filter(!ismissing, excel[sheet]["J"])
        name, values = Symbol(column[1]), column[2:end]
        df = dropmissing!(DataFrame(name => values), disallowmissing=true)
        df[!,name] = strip.(df[:,name])
        sort!(df)
    end
end

discretizer(dfs, col) = CategoricalDiscretizer(sort!(unique!(vcat(dfs...)))[:,col])
function encode!(disc::CategoricalDiscretizer, df::DataFrame, col::Symbol)
    df[!,Symbol(string(col) * "code")] = encode(disc, df[:,col])
    df
end

function discretize(dfs, col)
    disc = discretizer(dfs, col)
    map(df -> encode!(disc, df, col), dfs), disc
end

struct Space
    n::Int
    v::Int
    Space(n) = n < 0 ? error("invalid dimension") : new(n, 2^n - 1)
end

Base.iterate(space::Space) = trues(space.n), (0, trues(space.n))

function Base.iterate(space::Space, state)
    n, s = state
    if n >= space.v - 1
        nothing
    else
        for i in 1:space.n
            if s[i]
                s[i] = zero(s[i])
                s[1:i-1] .= one(s[i])
                break
            end
        end
        copy(s), (n + 1, s)
    end
end

Base.length(space::Space) = space.v
Base.eltype(space::Space) = BitArray{1}

venn(xs, col) = venn(map(x -> x[:,col], xs))
function venn(xs::AbstractVector{T}) where T
    v = Dict{BitArray{1}, T}()
    for idx in Space(length(xs))
        v[idx] = setdiff(intersect(xs[idx]...), xs[.!idx]...)
    end
    v
end

function assignkos(genera::AbstractVector{DataFrame}, kos::AbstractVector{DataFrame})
    @assert length(genera) == length(kos)
    gvenn = venn(genera, :genuscode)
    kvenn = venn(kos, :KOcode)
    merge((assignkos(gvenn[idx], kvenn[idx]) for idx in keys(gvenn))...)
end

function assignkos(genera::AbstractVector{S}, kos::AbstractVector{T}) where {S, T}
    assignment = Dict{S, Set{T}}()
    for g in genera
        assignment[g] = Set(rand(kos, 10))
    end
    for k in kos
        push!(assignment[rand(genera)], k)
    end
    assignment
end

function group!(depends, genera, n, m)
    N = length(genera)
    for _ in 1:n
        U = rand(2:m)
        group = genera[randperm(N)[1:U]]
        for i in 1:U, j in 1:U
            if i != j
                add_edge!(depends, group[i], group[j])
            end
        end
    end
    depends
end

function group(genera::AbstractVector{DataFrame}, n, m)
    M = length(genera)
    N = length(unique!(vcat(genera...)).genuscode)
    gvenn = venn(genera, :genuscode)
    T = gvenn |> valtype |> eltype

    G = SimpleGraph(N)
    for i in 1:M
        a = Set{T}()
        for j in 1:M
            for s in Space(M)
                if s[i] && !s[j]
                    union!(a, Set(gvenn[s]))
                end
            end
        end
        group!(G, collect(a), n, m)
    end
    G
end

sample(G, n) = sort!(randperm(nv(G))[1:n])

redux(G, prev) = filter(i -> neighbors(G, i) ⊆ prev, prev)

function subsample(G, prev, n)
    next = redux(G, prev)
    N = length(next)
    if N > n
        sort!(next[randperm(N)][1:n])
    else
        while N < n
            next = unique!(append!(next, sample(G, n - N)))
            N = length(next)
        end
        sort!(next)
    end
    #  N = length(prev)
    #  if N > n
    #      next = sort!(prev[randperm(N)[1:n]])
    #      redux(G, next)
    #  else
    #      prev
    #  end
end

function occurances(G, samplings)
    N = length(samplings)
    occurs = Array{Int}(undef, N, nv(G))
    for j in 1:nv(G), i in 1:N
        occurs[i, j] = (j in samplings[i]) ? 2 : 1
    end
    occurs
end

cbw(G, m, dfs::AbstractVector{DataFrame}) = cbw(G, m, nrow.(dfs))

function cbw(G, m, ns)
    oc = Array{Int}(undef, length(ns), m, nv(G))
    @threads for i in 1:m
        samplings = [sample(G, ns[1])]
        for n in ns[2:end]
            push!(samplings, subsample(G, samplings[end], n))
        end
        oc[:,i,:] .= occurances(G, samplings)
    end
    oc
end

function Eolas.mutualinfo(xs::AbstractMatrix, ys::AbstractMatrix)
    @assert size(xs) == size(ys)
    m = MutualInfo(2, 2)
    for k in 1:size(xs,2)
        @views observe!(m, xs[1:end-1, k], ys[2:end, k])
    end
    estimate(m)
end

function Eolas.transferentropy(xs::AbstractMatrix, ys::AbstractMatrix)
    @assert size(xs) == size(ys)
    te = TransferEntropy(2, 2, 1)
    for k in 1:size(xs, 2)
        @views observe!(te, xs[:, k], ys[:, k])
    end
    estimate(te)
end

function sigtest(measure, xs, ys; nperms=100)
    xsperm = copy(xs)
    gt = measure(xs, ys)
    c = 1
    for i in 1:nperms
        @views for k in 1:size(xsperm, 2)
            shuffle!(xsperm[:, k])
        end
        c += (measure(xsperm, ys) ≥ gt)
    end
    gt, c / (nperms + 1)
end

function lmi(data; nperms=0)
    mi = Array{Float64}(undef, size(data,3), size(data,3))
    @threads for i in 1:size(data, 3)
        for j in i:size(data, 3)
            a, b = if !iszero(nperms)
                a, p = @views sigtest(mutualinfo, data[:, :, i], data[:, :, j]; nperms=nperms)
                b, p = @views sigtest(mutualinfo, data[:, :, j], data[:, :, i]; nperms=nperms)
                a = p < 0.05 ? a : zero(a)
                b = p < 0.05 ? b : zero(b)
                a, b
            else
                a = @views mutualinfo(data[:, :, i], data[:, :, j])
                b = @views mutualinfo(data[:, :, j], data[:, :, i])
                a, b
            end
            mi[i, j] = mi[j, i] = 0.5 * (a + b)
        end
    end
    mi
end

function nte(data; nperms=0)
    te = Array{Float64}(undef, size(data,3), size(data,3))
    @threads for i in 1:size(data, 3)
        for j in i:size(data, 3)
            a, b = if !iszero(nperms)
                a, p = @views sigtest(transferentropy, data[:, :, i], data[:, :, j]; nperms=nperms)
                b, p = @views sigtest(transferentropy, data[:, :, j], data[:, :, i]; nperms=nperms)
                a = p < 0.05 ? a : zero(a)
                b = p < 0.05 ? b : zero(b)
                a, b
            else
                a = @views transferentropy(data[:, :, i], data[:, :, j])
                b = @views transferentropy(data[:, :, j], data[:, :, i])
                a, b
            end
            te[i, j] = te[j, i] = 0.5 * (a + b)
        end
    end
    te
end

function neighboring(G)
    ns = falses(nv(G), nv(G))
    for i in 1:nv(G)
        for j in neighbors(G, i)
            ns[i, j] = true
        end
    end
    ns
end

function main()
    genera = genus()
    discretize(genera, :genus)
    G = group(genera, 20, 5)
    gt = vec(neighboring(G))
    @time te = vec(nte(cbw(G, 12, genera); nperms=100))
    idx = te .!= zero(eltype(te))
    gt = gt[idx]
    te = te[idx]
    p = rocplot(gt, te)
    for _ in 1:10
        G = group(genera, 20, 5)
        gt = vec(neighboring(G))
        @time te = vec(nte(cbw(G, 12, genera); nperms=100))
        idx = te .!= zero(eltype(te))
        gt = gt[idx]
        te = te[idx]
        rocplot!(p, gt, te)
    end
    savefig(p, "roc.svg")
    p
end

function rocplot!(p, gt, score)
    c = roc(gt, score, 100)
    fpr, tpr = false_positive_rate.(c), true_positive_rate.(c)
    plot!(p, fpr, tpr; linecolor=:black)
end

function rocplot(gt, score)
    p = plot(title="Receiver Operating Characteristic", xlabel="FPR", ylabel="TPR",
             xlim=(0,1), ylim=(0,1), aspect_ratio=1, size=(600,600), legend=nothing)
    plot!(p, [0, 1], [0, 1]; linecolor=:gray, line=:dash)
    rocplot!(p, gt, score)
end

max_stars_repo_name: ELIFE-ASU/dglmoore

max_issues_repo_name: ELIFE-ASU/dglmoore

max_forks_repo_name: ELIFE-ASU/dglmoore


text:
"""
Financial blotter type: keeps track of filled transactions
for one specific instrument (in metadata).
"""
typealias Blotter{T<:Float64} TimeArray{T,2,FinancialSeries.Ticker}

"Standard blotter columns: order quantity and fill price."
const blottercolnames = ["Qty", "FillPrice"]

"Placeholder ticker"
const blotterticker = FinancialSeries.Ticker("ticker")

"Blotter object without any filled transactions"
emptyblotter() = TimeSeries.TimeArray([DateTime(1795,10,31)], [0.0 0.0], blottercolnames, blotterticker)

# """
# Add a single transaction at the end (time-wise) of the Blotter-object.
# Modifies `fb`.
# """
# function push!(fb::Blotter, tstamp::DateTime, qty::Float64, fillprice::Float64)
#   fb.colnames == blottercolnames || error("Unexpected blotter column names")

#   ### this does not modify the object but creates another one
#   ### same with TimeSeries merge

#   return TimeArray([fb.timestamp, tstamp],
#                    [fb.values, [qty fillprice]],
#                    fb.colnames,
#                    fb.meta)
# end

# push!(fb::Blotter, tstamp::DateTime, qty::Int64, fillprice::Float64) =
#   push!(fb, tstamp, float(qty), fillprice)

### TODO use merge, append!
# this probably blows up if you try to insert vs add the end
#add!(b::Blotter, entry::Blotter) = Blotter(vcat(b.timestamp, entry.timestamp), vcat(b.values, entry.values), blottercolnames, blotterticker)

# """
# Backtesting blotter for instrument series `fts` based on `signal`.
# CAUTION: not tested.
# ...
# Returns `Blotter`-object.
# """
# function Blotter(signal::TimeArray{Bool,1}, fts::FinancialTimeSeries{Float64,2};
#                  quantity = 100, price = "open", slippage = .1)
#   op, cl = fts["Open"], fts["Close"]
#     tt = lag(signal)              # 495 row of bools, the next day
#     t  = discretesignal(tt)  # 78 rows of first true and first false, as floats though

#     price == "open" ?
#     fills = op[datetolastsecond([t.timestamp])] :
#     price == "close" ?
#     fills = cl[datetolastsecond([t.timestamp])] :
#     error("only open and close prices supported for naive backtests")

#     notsigned  = TimeArray(t.timestamp, quantity * ones(length(t)), ["Qty"])
#     exitdates  = findwhen(t.==0)
#     qty        = quantity * ones(length(t.timestamp))

#     for i in 1:length(notsigned)
#         for j in 1:length(exitdates)
#             if notsigned[i].timestamp[1] == exitdates[j]
#                 qty[i] = qty[i] * -1
#             end
#         end
#     end

#    datetimes = datetolastsecond(t.timestamp)

#    Blotter(datetimes, [qty fills.values], blottercolnames, fts.instrument.ticker)
# end


text:
precompile(Tuple{typeof(PlayGround.include), String})
precompile(Tuple{typeof(PlayGround.World.include), String})
precompile(Tuple{typeof(Base._iterate), Base.Dict{Int32, Symbol}, Int64})
precompile(Tuple{Type{PlayGround.World.ResourceName}, Int32})
precompile(Tuple{typeof(RuntimeGeneratedFunctions.init), Module})
precompile(Tuple{typeof(Base.lock), RuntimeGeneratedFunctions.var"#7#8"{Module}, Base.Threads.SpinLock})
precompile(Tuple{typeof(PlayGround.World.eval), Expr})
precompile(Tuple{Type{PlayGround.World.OperatorType}, Int32})
precompile(Tuple{Type{PlayGround.World.Operator}, Symbol, PlayGround.World.OperatorType})
precompile(Tuple{typeof(Base.vect), PlayGround.World.Operator, Vararg{PlayGround.World.Operator, N} where N})
precompile(Tuple{typeof(Base.vect), Array{PlayGround.World.Operator, 1}, Vararg{Array{PlayGround.World.Operator, 1}, N} where N})
precompile(Tuple{typeof(Base.allocatedinline), Type{PlayGround.World.Operator}})
precompile(Tuple{typeof(Base.reduce), typeof(Base.vcat), Array{Array{PlayGround.World.Operator, 1}, 1}})
precompile(Tuple{RuntimeGeneratedFunctions.var"#@RuntimeGeneratedFunction", LineNumberNode, Module, Any})
precompile(Tuple{typeof(RuntimeGeneratedFunctions._RGF_constructor_code), Expr, Expr})
precompile(Tuple{typeof(GLFW_jll.__init__)})
precompile(Tuple{Core.var"#Type##kw", NamedTuple{(:libc, :compiler_abi), Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.FreeBSD}, Symbol})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.FreeBSD, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.FreeBSD})
precompile(Tuple{typeof(Libdl.dlopen), String, UInt32})
precompile(Tuple{typeof(JLLWrappers.get_julia_libpaths)})
precompile(Tuple{typeof(GLFW._ErrorCallbackWrapper), Int32, Base.Cstring})
precompile(Tuple{typeof(GLFW.__init__)})
precompile(Tuple{typeof(CImGui_jll.__init__)})
precompile(Tuple{typeof(Libdl.dlopen), String})
precompile(Tuple{typeof(Base.filter!), Base.var"#64#65"{typeof(Base.isempty)}, Array{String, 1}})
precompile(Tuple{typeof(Base.join), Array{String, 1}, Char})
precompile(Tuple{typeof(CImGui.GLFWBackend.__init__)})
precompile(Tuple{typeof(CImPlot_jll.__init__)})
precompile(Tuple{typeof(Base.foreach), CImPlot_jll.var"#7#9", Tuple{Array{String, 1}}})
precompile(Tuple{typeof(Base.foreach), CImPlot_jll.var"#8#10", Tuple{Array{String, 1}}})
precompile(Tuple{typeof(GLFW.WindowHint), UInt32, Int64})
precompile(Tuple{typeof(GLFW.WindowHint), UInt32, UInt32})
precompile(Tuple{typeof(Base.join), Tuple{Symbol, Symbol, Symbol}, String})
precompile(Tuple{Base.var"##sprint#355", Nothing, Int64, typeof(Base.sprint), Function, Tuple{Symbol, Symbol, Symbol}, Vararg{Any, N} where N})
precompile(Tuple{typeof(Base.join), Base.GenericIOBuffer{Array{UInt8, 1}}, Tuple{Symbol, Symbol, Symbol}, String})
precompile(Tuple{typeof(GLFW.SetErrorCallback), Function})
precompile(Tuple{Type{PlayGround.World.GlyphRange}, UInt16, UInt16})
precompile(Tuple{typeof(CImGui.LibCImGui.ImFontConfig_ImFontConfig)})
precompile(Tuple{typeof(Base.tail), Tuple{String, Expr, Expr}})
precompile(Tuple{typeof(Printf._printf), String, Expr, String, Tuple{Expr, Expr}})
precompile(Tuple{typeof(Printf.gen_e), String, Int64, Int64, Char})
precompile(Tuple{typeof(Base.tail), Tuple{String, Symbol, Symbol}})
precompile(Tuple{typeof(Printf._printf), String, Expr, String, Tuple{Symbol, Symbol}})
precompile(Tuple{typeof(Printf._printf), String, Expr, String, Tuple{}})
precompile(Tuple{Type{PlayGround.World.StateFloatIndices}, Int32})
precompile(Tuple{Type{PlayGround.World.StateIntIndices}, Int32})
precompile(Tuple{Base.var"#@ccallable", LineNumberNode, Module, Any})
precompile(Tuple{typeof(Base.expand_ccallable), Nothing, Expr})
precompile(Tuple{Type{Base.Generator{I, F} where F where I}, Base.var"#283#284", Array{Any, 1}})
precompile(Tuple{typeof(Base.collect_similar), Array{Any, 1}, Base.Generator{Array{Any, 1}, Base.var"#283#284"}})
precompile(Tuple{typeof(Base._ccallable), Type{T} where T, Type{T} where T})
precompile(Tuple{typeof(GLFW.Terminate)})


text:
using POMDPs

#POMDPs.add("POMDPToolbox")
#POMDPs.add("POMDPXFiles")


text:
module Manipulation

include("types.jl")
include("arrays.jl")
include("tables.jl")
include("piecewise.jl")

export DataGrid
export merge_grids, merge_knots
export unique_knots
export piecewise_function
export merge_table_slave_to_master_sigmoid

end  # module Manipulation


text:
module CLagrangianFESpaces

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery
using Gridap.DOFBases: _length
using Gridap.ConformingFESpaces: _CellField
using Gridap.ConformingFESpaces: CellEqClass
using Gridap.BoundaryGrids: _setup_tags

import Gridap: num_free_dofs
import Gridap: num_diri_dofs
import Gridap: diri_tags
import Gridap: apply_constraints
import Gridap: apply_constraints_rows
import Gridap: apply_constraints_cols
import Gridap: celldofids
import Gridap: interpolate_values
import Gridap: interpolate_diri_values
import Gridap: CellField
import Gridap: CellBasis
import Gridap: Triangulation

export CLagrangianFESpace

struct CLagrangianFESpace{D,Z,T,S,B} <: FESpace{D,Z,T}
  grid::Grid{D,Z}
  node_to_label::Vector{Int}
  tag_to_labels::Vector{Vector{Int}}
  diritags::Vector{Int}
  dirimasks::Vector{B}
  nfreedofs::Int
  ndiridofs::Int
  node_and_comp_to_dof::Vector{S}
  cell_to_dofs::CellVector{Int}
  reffe::LagrangianRefFE{Z,T}
  cellbasis::CellBasis{Z,T}
  dirinodes::Vector{Int}
end

function CLagrangianFESpace(
  ::Type{T},
  grid::Grid{D,Z},
  node_to_label::AbstractVector{<:Integer},
  tag_to_labels::Vector{Vector{Int}},
  diritags::Vector{Int},
  dirimasks::Vector{B}) where {D,Z,T,B}

  node_and_comp_to_dof, nfreedofs, ndiridofs, dirinodes =
    _setup_node_and_comp_to_dof(
    T, node_to_label, tag_to_labels, diritags, dirimasks)

  cell_to_nodes = cells(grid)

  reffe = _setup_reffe(T,celltypes(grid),cellorders(grid))
  lnode_and_comp_to_ldof = reffe.dofbasis.node_and_comp_to_dof

  cell_to_dofs = _setup_cell_to_dofs(
    cell_to_nodes,
    node_and_comp_to_dof,
    lnode_and_comp_to_ldof)

  cellbasis = _setup_cellbasis(reffe,grid)

  S = _S(T)

  CLagrangianFESpace{D,Z,T,S,B}(
    grid,
    node_to_label,
    tag_to_labels,
    diritags,
    dirimasks,
    nfreedofs,
    ndiridofs,
    node_and_comp_to_dof,
    cell_to_dofs,
    reffe,
    cellbasis,
    dirinodes)

end

function CLagrangianFESpace(
  ::Type{T},
  model::DiscreteModel,
  order::Integer,
  diritags=Int[],
  dirimasks=nothing) where T
  CLagrangianFESpace(T,model,FaceLabels(model),order,diritags,dirimasks)
end

function CLagrangianFESpace(
  ::Type{T},
  model::DiscreteModel,
  facelabels::FaceLabels,
  order::Integer,
  diritags=Int[],
  dirimasks=nothing) where T

  _diri_tags = _setup_tags(facelabels,diritags)
  _diri_masks = _setup_masks(T,_diri_tags,dirimasks)

  grid, node_to_label, tag_to_labels = _setup_grid(model,facelabels,order)

  CLagrangianFESpace(
    T,grid,node_to_label,tag_to_labels,_diri_tags,_diri_masks)

end

num_free_dofs(fesp::CLagrangianFESpace) = fesp.nfreedofs

num_diri_dofs(fesp::CLagrangianFESpace) = fesp.ndiridofs

diri_tags(fesp::CLagrangianFESpace) = fesp.diritags

apply_constraints(
  ::CLagrangianFESpace, cellvec::CellVector, cellids::CellNumber) = cellvec

apply_constraints_rows(
  ::CLagrangianFESpace, cellmat::CellMatrix, cellids::CellNumber) = cellmat

apply_constraints_cols(
  ::CLagrangianFESpace, cellmat::CellMatrix, cellids::CellNumber) = cellmat

celldofids(fesp::CLagrangianFESpace) = fesp.cell_to_dofs

function interpolate_values(fesp::CLagrangianFESpace,fun::Function)
  _interpolate_values(fesp,fun)
end

function interpolate_values(fesp::CLagrangianFESpace{D,Z,T},val::T) where {D,Z,T}
  _interpolate_values(fesp,val)
end

function _interpolate_values(fesp,fun)
  zh = zero(fesp)
  fdof_to_val = free_dofs(zh)
  ddof_to_val = diri_dofs(zh)
  _fill_interpolated_vals!(fdof_to_val,ddof_to_val,fesp,fun)
  (fdof_to_val, ddof_to_val)
end

function interpolate_diri_values(fesp::CLagrangianFESpace, funs::Vector{<:Function})
  _interpolate_diri_values(fesp, funs)
end

function interpolate_diri_values(fesp::CLagrangianFESpace{D,Z,T},vals::Vector{T}) where {D,Z,T}
  _interpolate_diri_values(fesp, vals)
end

function _interpolate_diri_values(fesp, funs)

  T = value_type(fesp)
  E = eltype(T)
  ndiri = num_diri_dofs(fesp)
  ddof_to_val = zeros(E,ndiri)
  node_to_coords = points(fesp.grid)
  node_and_comp_to_dof = fesp.node_and_comp_to_dof
  node_to_label = fesp.node_to_label
  tag_to_labels = fesp.tag_to_labels
  dirinodes = fesp.dirinodes

  for (fun,diritag,dirimask) in zip(funs,fesp.diritags,fesp.dirimasks)
    _fill_diri_values!(
      ddof_to_val,
      fun,
      diritag,
      dirimask,
      dirinodes,
      node_to_label,
      tag_to_labels,
      node_to_coords,
      node_and_comp_to_dof)
  end

  ddof_to_val
end

function CellField(
  fesp::CLagrangianFESpace{D,Z,T},
  free_dofs::AbstractVector{E},
  diri_dofs::AbstractVector{E}) where {D,Z,T,E}

  _CellField(fesp, free_dofs, diri_dofs, T, E)
end

CellBasis(fesp::CLagrangianFESpace{D,Z,T}) where {D,Z,T} = fesp.cellbasis

Triangulation(fesp::CLagrangianFESpace) = Triangulation(fesp.grid)

# Helpers

function _setup_masks(T,diritags,dirimasks)
  if dirimasks === nothing
    ncomps = _length(T)
    ntags = length(diritags)
    t = tuple(fill(true,ncomps)...)
    return fill(t,ntags)
  else
    @assert length(diritags) == length(dirimasks)
    return dirimasks
  end
end

function _setup_node_and_comp_to_dof(
  ::Type{T}, node_to_label, tag_to_labels, diritags, dirimasks) where T

  S = _S(T)
  nnodes = length(node_to_label)
  node_and_comp_to_dof = zeros(S,nnodes)
  dirinodes = Int[]

  nf, nd = _fill_node_and_comp_to_dof!(
    node_and_comp_to_dof, dirinodes, node_to_label, tag_to_labels, diritags, dirimasks)

  (node_and_comp_to_dof, nf, nd, dirinodes)

end

function _fill_node_and_comp_to_dof!(
  node_and_comp_to_dof::Vector{S},
  dirinodes,
  node_to_label,
  tag_to_labels,
  diritags,
  dirimasks) where S

  nnodes = length(node_to_label)
  fdof = 1
  ddof = -1

  for node in 1:nnodes
    label = node_to_label[node]

    comp_to_dof, fdof, ddof, isdirinode =
      _compute_comp_to_dof(S,label,tag_to_labels,diritags,dirimasks,fdof,ddof)

    if isdirinode
      push!(dirinodes,node)
    end

    node_and_comp_to_dof[node] = comp_to_dof

  end

  (fdof-1,-(ddof+1))

end

function _compute_comp_to_dof(
  ::Type{S}, label, tag_to_labels, diritags, dirimasks, fdof, ddof) where S <: MultiValue

  comp_to_dof = zero(mutable(S))
  comps = eachindex(comp_to_dof)

  isdirinode = false
  for comp in comps
    dof, fdof, ddof = _compute_dof(label,tag_to_labels,diritags,dirimasks,comp,fdof,ddof)
    comp_to_dof[comp] = dof
    if dof < 0
      isdirinode = true
    end
  end

  (comp_to_dof, fdof, ddof, isdirinode)

end

function _compute_comp_to_dof(
  ::Type{S}, label, tag_to_labels, diritags, dirimasks, fdof, ddof) where S <: Real

  comp = 1
  dof, fdof, ddof =  _compute_dof(
    label,tag_to_labels,diritags,dirimasks,comp,fdof,ddof)

  isdirinode = dof < 0

  (dof, fdof, ddof, isdirinode)

end

function _compute_dof(label,tag_to_labels,diritags,dirimasks,comp,fdof,ddof)

  isfree = true
  for (i,diritag) in enumerate(diritags)

    if label in tag_to_labels[diritag]
      comp_to_mask = dirimasks[i]
      mask = comp_to_mask[comp]
      if mask
        dof = ddof
        ddof -= 1
        isfree = false
        break
      end
    end
  end

  if isfree
    dof = fdof
    fdof += 1
  end

  (dof,fdof,ddof)

end

_S(::Type{<:Real}) = Int

_S(::Type{MultiValue{S,T,N,L}}) where {S,T,N,L} = MultiValue{S,Int,N,L}

function _setup_reffe(ct,co)
  @notimplemented
end

function _setup_reffe(
  ::Type{T},
  ct::ConstantCellValue{<:NTuple{Z}},
  co::ConstantCellValue) where {T,Z}
  t = ct.value
  o = co.value
  p = Polytope(t)
  LagrangianRefFE(T,p,fill(o,Z))
end

function _setup_cellbasis(reffe,grid)
  lshapes = shfbasis(reffe)
  nc = ncells(grid)
  shb = ConstantCellValue(lshapes, nc)
  trian = Triangulation(grid)
  phi = CellGeomap(trian)
  basis = attachgeomap(shb,phi)
  basis
end

function _setup_cell_to_dofs(
  cell_to_nodes,
  node_and_comp_to_dof::Vector{S},
  lnode_and_comp_to_ldof) where S

  ncomp = _length(S)
  ncells = length(cell_to_nodes)
  nlnodes, ncomp = size(lnode_and_comp_to_ldof)
  nldofs = nlnodes * ncomp

  ndata = ncells*nldofs
  I = eltype(S)

  cell_to_dofs_data = zeros(I,ndata)

  _fill_cell_to_dofs!(
    cell_to_dofs_data,
    cell_to_nodes,
    node_and_comp_to_dof,
    lnode_and_comp_to_ldof)

  CellVectorFromDataAndStride(cell_to_dofs_data,nldofs)

end

function  _fill_cell_to_dofs!(
  cell_to_dofs_data,
  cell_to_nodes,
  node_and_comp_to_dof,
  lnode_and_comp_to_ldof)

  k = 0

  nlnodes, ncomps = size(lnode_and_comp_to_ldof)
  nldofs = nlnodes * ncomps

  for nodes in cell_to_nodes

    for (lnode,node) in enumerate(nodes)
      for comp in 1:ncomps
        ldof = lnode_and_comp_to_ldof[lnode,comp]
        dof = node_and_comp_to_dof[node][comp]
        cell_to_dofs_data[k+ldof] = dof
      end
    end

    k += nldofs

  end

end

function _fill_interpolated_vals!(fdof_to_val,ddof_to_val,fesp,fun)

  x = points(fesp.grid)
  fx = _eval_bc(fun,x)
  T = value_type(fesp)
  @assert eltype(fx) == T
  z = zero(T)
  comps = eachindex(z)

  _fill_interpolated_vals_kernel!(
    fdof_to_val,
    ddof_to_val,
    fesp.node_and_comp_to_dof,
    fx,
    comps)

end

_eval_bc(fun::Function,x) = fun.(x)

_eval_bc(val,x) = fill(val,length(x))

function  _fill_interpolated_vals_kernel!(
  fdof_to_val,
  ddof_to_val,
  node_and_comp_to_dof,
  node_and_comp_to_val,
  comps)

  nnodes = length(node_and_comp_to_dof)
  @assert nnodes == length(node_and_comp_to_val)

  for node in 1:nnodes
    comp_to_dof = node_and_comp_to_dof[node]
    comp_to_val = node_and_comp_to_val[node]
    for comp in comps
      dof = comp_to_dof[comp]
      val = comp_to_val[comp]
      if dof > 0
        fdof_to_val[dof] = val
      elseif dof < 0
        ddof_to_val[-dof] = val
      else
        @unreachable
      end
    end
  end

end

function _fill_diri_values!(
  ddof_to_val,
  fun,
  diritag,
  dirimask,
  dirinodes,
  node_to_label,
  tag_to_labels,
  node_to_coords,
  node_and_comp_to_dof)

  for node in dirinodes

    label = node_to_label[node]
    if label in tag_to_labels[diritag]
      x = node_to_coords[node]
      fx = _eval(fun,x)
      comp_to_dof = node_and_comp_to_dof[node]
      for comp in eachindex(fx)
        mask = dirimask[comp]
        if mask
          val = fx[comp]
          dof = comp_to_dof[comp]
          ddof = -dof
          ddof_to_val[ddof] = val
        end
      end
    end

  end

end

@inline _eval(fun::Function,x) = fun(x)

@inline _eval(val,x) = val

function _same_order(order,co)
  @notimplemented
end

function _same_order(order,co::ConstantCellValue)
  order == co.value
end

function _setup_grid(model,facelabels,order)

  _grid = Grid(model)

  if _same_order(order,cellorders(_grid))
    @notimplementedif order !=1
    grid = _grid
    node_dim = 0
    node_to_label = labels_on_dim(facelabels,node_dim)
  else
    grid, node_to_label = grid_from_model_and_order(model,order)
  end

  tag_to_labels = facelabels.tag_to_labels

  (grid, node_to_label, tag_to_labels)

end

function grid_from_model_and_order(
  model::DiscreteModel{D}, order::Integer) where D

  node_to_coords, dim_to_nface_to_nodes, node_to_label =
    _setup_nodes_to_coords(model,D,order)

  graph = GridGraph(model)

  dim_to_cell_to_nface = [connections(graph,D,i) for i in 0:D]

  offset = length.(dim_to_nface_to_nodes)
  for i in 2:length(offset)
    offset[i] += offset[i-1]
  end
  offset = tuple(offset[1:(end-1)]...)

  cell_to_faces = local_append(offset, dim_to_cell_to_nface...)

  face_to_nodes = append(dim_to_nface_to_nodes...)

  grid = Grid(model,D)
  reffe = _setup_reffe(order,celltypes(grid))

  cell_to_nodes = CellEqClass(cell_to_faces, face_to_nodes, reffe)

  cell_to_nodes_data, cell_to_nodes_ptrs = compress(cell_to_nodes)

  cell_to_extrusion = celltypes(grid)

  n = ncells(grid)
  cell_to_order = ConstantCellValue(order,n)

  grid = UnstructuredGrid(
    node_to_coords,
    cell_to_nodes_data,
    cell_to_nodes_ptrs,
    cell_to_extrusion,
    cell_to_order)

  (grid, node_to_label)

end

function _setup_nodes_to_coords(model,D,order)

  node_to_coords = Point{D,Float64}[]
  node_to_label = Int[]

  facelabels = FaceLabels(model)

  vertex_to_label = labels_on_dim(facelabels,0)

  grid = Grid(model,0)
  _push_nodes!(node_to_coords,points(grid))
  _push_nodes!(node_to_label,vertex_to_label)
  node_to_nodes = cells(grid)

  dim_to_nface_to_nodes = []
  push!(dim_to_nface_to_nodes,node_to_nodes)

  for d in 1:D

    grid = Grid(model,d)
    trian = Triangulation(grid)
    reffe = _setup_reffe(order,celltypes(grid))

    ilnodes = reffe.nfacenodes[end]
    linode_to_qcoords = reffe.dofbasis.nodes[ilnodes]
    n = ncells(grid)
    q = ConstantCellValue(linode_to_qcoords,n)
    phi = CellGeomap(trian)
    nface_ilnode_to_coords = evaluate(phi,q)

    nface_to_nodes_data = Int[]
    nface_to_nodes_ptrs = zeros(Int,n+1)

    nface_to_label = labels_on_dim(facelabels,d)

    _push_nface_nodes!(
      node_to_coords,
      node_to_label,
      nface_to_nodes_data,
      nface_to_nodes_ptrs,
      nface_ilnode_to_coords,
      nface_to_label)

    nface_to_nodes = CellVectorFromDataAndPtrs(
      nface_to_nodes_data, nface_to_nodes_ptrs)

    push!(dim_to_nface_to_nodes,nface_to_nodes)

  end

  (node_to_coords, dim_to_nface_to_nodes, node_to_label)

end

function _setup_reffe(order,ct::ConstantCellValue)
  e = ct.value
  p = Polytope(e)
  D = length(e)
  LagrangianRefFE(Float64,p,fill(order,D))
end

function _push_nodes!(node_to_coords,coords)
  for x in coords
    push!(node_to_coords,x)
  end
end

function _push_nface_nodes!(
  node_to_coords,
  node_to_label,
  nface_to_nodes_data,
  nface_to_nodes_ptrs,
  nface_ilnode_to_coords,
  nface_to_label)

  node = length(node_to_coords)+1
  for (nface,ilnode_to_coords) in enumerate(nface_ilnode_to_coords)

    label = nface_to_label[nface]
    nface_to_nodes_ptrs[nface+1] = length(ilnode_to_coords)

    for coords in ilnode_to_coords
      push!(node_to_coords,coords)
      push!(node_to_label,label)
      push!(nface_to_nodes_data,node)
      node += 1
    end
  end

  length_to_ptrs!(nface_to_nodes_ptrs)

end

end # module


text:
using Test, Alexya

include("widgets.jl")
include("canvas.jl")
include("grid.jl")
include("others.jl")
include("controls.jl")


text:
module HuffmanCode

using BufferedStreams, FileIO

include("Node.jl")

export ZIP

const DICT = 8

global HuffmanCode = Dict()

function zipFile(src::String, dir::String)
    bis = BufferedInputStream(open(src, "r"))
    b = read(bis, UInt8)
    close(bis)
    # 进行压缩
    byteZip = HuffmanZip(b)
    # 输出文件
    bos = BufferedOutputStream(open(dir, "w"))
    write(bos, byteZip)
    flush(bos)
    close(bos)
end

function huffmanZip(bytes::UInt8)
    # 获取节点
    nodes = getNodes(bytes)
    # 创建树
    tree = createHuffmanTree(nodes)
    # 获取编码表
    huffmanCodes = getCodes(tree)
    # 返回数据和编码表
    return zip(bytes, huffmanCodes)
end

function getCodes(tree::Node)
    if tree === nothing
        error("你搁这整活呐?")
    end
    getCodes(tree.left, "0")
    getCodes(tree.right, "1")
    return huffmanCodes
end

function getCodes(node::Node, code::String)
    str = Array{String, 1}
    push!(str, code)
    if node.data === nothing
        getCodes(node.left, "0")
        getCodes(node.right, "1")
    else
        push!(huffmanCodes, node.data, str)
    end
end

end

max_stars_repo_name: Colcactus/Huffman

max_issues_repo_name: Colcactus/Huffman

max_forks_repo_name: Colcactus/Huffman


text:
@testset "Selection lower bound tests" begin
    @test selectionLowerBound([1, 3], [3, 1], [2, 2], 2, getKnapsackConstraints([1, 2], 1), 0.5, KnapsackProblemDescriptor(2)) == 3.5

    @test selectionLowerBound([2 3; 1 4], [2 2; 3 1], [5 4; 3 4], 6, getAssignmentConstraints(2), 0.5, AssignmentProblemDescriptor(2)) == 11
end


text:
import Distributions.partype
import PDMats.dim
import Base.Random.rand
doc"""
    LKJ(η, p)

The *LKJ correlation matrix* with shape parameter \eta.

```julia
LKJ(eta, p)        # LKJ distribution with shape parameter \eta and dimension p

dim(d)            # Get the dimensions of LKJ correlation matrix
size(d)           # Get the size of correlation matrix as a tuple
params(d)         # Get the parameters, i.e. (mu, sig)

rand(d, n)          # Draw random samples from LKJ distribution, by default n=1
meanlogdet(d)       # Expected log determinant
meanlog(d, h)
entropy(d)
```

"""

immutable LKJ{T<:Real} <: ContinuousMatrixDistribution
    η::T
    p::Int

    function LKJ(η, p)
        η > zero(η) || error("η should be greater than zero.")
        new(η, p)
    end
end

LKJ(η::Real, p::Integer) = LKJ{typeof(η)}(η, p)
LKJ(η::Integer, p::Integer) = LKJ{Float64}(Float64(η), p)

@inline partype{T<:Real}(d::LKJ{T}) = T

"""
Draw n random samples with replacement from the LKJ distribution (n = 1 by default).
"""
function rand(d::LKJ, n::Integer)
    η = d.η
    p = d.p
    idxvec = [η + (p - i)/2. for i in 2:p] # generate Beta distribution indices
    L = zeros(p, p)
    for j in 1:length(idxvec)
        L[j+1:end, j] = 2 * rand(Distributions.Beta(idxvec[j], idxvec[j]), p - j) - 1
    end
    L = _cpc_to_corr(L, p)
    S = PDMats.PDMat(Base.LinAlg.Cholesky(full(L), :L)) # consistent with Stan Lib
    rand(S.mat, n)
end

rand(d::LKJ) = rand(d::LKJ, 1)

function _cpc_to_corr(L::AbstractMatrix, p::Real)  # consistent with Stan Lib
    #p = dim(d)
    LL = eye(p)
    for j in 1:p
        for i in j+1:p
            sumsqs = 0 # faster than dot(vec(L[i, 1:j]), vec(L[i, 1:j]))
            for ij in 1:j
                sumsqs += LL[i, ij]^2
            end
            LL[i, j] = L[i, j] * sqrt(1 - sumsqs)
        end
        sumsqs = 0
        for ij in 1:j-1
            sumsqs += LL[j, ij]^2
        end
        LL[j, j] = sqrt(1 - sumsqs)
    end
    LowerTriangular(LL)
end

dim(d::LKJ) = d.p
size(d::LKJ) = (d.p, d.p)
params(d::LKJ) = (d.η, d.p)

"""
E[log |A|] where
A ~ LKJ(eta) of dimension p
"""
function meanlogdet(d::LKJ)
    η = d.η
    p = d.p
    idxvec = [η + (p - i)/2. for i in 2:p]
    ldsum = 0
    for i in 1:length(idxvec)
        ldsum += (digamma(idxvec[i]) - digamma(2 * idxvec[i])) * (p - i)
    end
    ldsum
end

"""
E_q[log p(x)] where
p(x) = LKJ(h, p)  # p = matrix dimension
q(x) = LKJ(eta, p)
"""
function meanlog(d::LKJ, h::Real)
    η = d.η
    p = d.p
    betaq = [η + (p - i)/2. for i in 2:p]
    betap = [h + (p - i)/2. for i in 2:p]
    Elogx = 0
    for i in 1:length(betaq) # compute from Dirichlet mean log
        Elogx += (p - i) * (2(betap[i] - 1) .* (digamma(betaq[i]) - digamma(2betaq[i])) - (2lgamma(betap[i]) - lgamma(2betap[i])))
    end
    Elogx
end

"""
Compute entropy for LKJ distribution by Dirichelet entropy
"""
function entropy(d::LKJ)
    η = d.η
    p = d.p
    idxvec = [η + (p - i)/2. for i in 2:p]
    H = 0
    for i in 1:length(idxvec) # compute from Dirichlet entropy
        H += (p - i) * (2*lgamma(idxvec[i]) - lgamma(2idxvec[i]) + 2(idxvec[i] - 1) .* digamma(2idxvec[i]) - 2 * (idxvec[i] - 1) .* digamma(idxvec[i]))
    end
    H
end


text:
########################################
## File Name: state_transition_observation_models.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: State Transition and Observation Models for SACBP
########################################

import ForwardDiff: gradient, jacobian, JacobianConfig, Chunk
using LinearAlgebra
using Random

##### State Transition Models ####
abstract type TransModel end
abstract type StateTransModel <: TransModel end

struct TransModel_Pos <: StateTransModel
# Linear state transition model for 2D position.
end
function trans(model::TransModel_Pos,x::PhysPos{T},u::PosControl{<:Real},dt::Real)::PhysPos{T} where T <: Real
# Euler Approximation.
    if x.dim != u.dim
        error(ArgumentError("PhysState and Control have inconsistent dimension values."));
    elseif x.t != u.t
        error(ArgumentError("PhysState and Control have inconsistent time parameters."));
    else
        pos = x.pos + u.vel*dt;
        t = round(x.t + dt, digits=5);
        return PhysPos(t,pos,x.dim)
    end
end
function trans(model::TransModel_Pos,x::PhysPos{T},u::PosControl{<:Real},dt::Real,
               Q::Matrix{Float64},rng::MersenneTwister)::PhysPos{T} where T <: Real
# Euler Approximation with stochastic disturbance.
    if x.dim != u.dim
        error(ArgumentError("PhysState and Control have inconsistent dimension values."));
    elseif x.t != u.t
        error(ArgumentError("PhysState and Control have inconsistent time parameters."));
    else
        pos = x.pos + u.vel*dt + rand(rng,MvNormal(zeros(x.dim),Q*dt));
        t = round(x.t + dt, digits=5);
        return PhysPos(t,pos,x.dim)
    end
end
trans_jacobi(model::TransModel_Pos,x::PhysPos,u::PosControl{<:Real}) = diagm(zeros(x.dim)); # Jacobian of dynamics.
trans_jacobi_euler(model::TransModel_Pos,x::PhysPos,u::PosControl{<:Real},dt::Real) = Matrix(1.0I, x.dim, x.dim) + trans_jacobi(model,x,u)*dt; # Jacobian of euler state transition for discrete-time EKF update..
trans_u_coeff(model::TransModel_Pos,x::PhysPos) = Matrix(1.0I, x.dim, x.dim); # Liner control coefficient for computing optimal control value.

struct TransModel_Manipulate2D <: StateTransModel
# Non-linear state transition model for planar object transport.
end
function trans(model::TransModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real},dt::Real)::PhysManipulate2D{T} where T <: Real
# Euler Approximation.
    if x.t != u.t
        error(ArgumentError("PhysState and Control have inconsistent time parameters."))
    else
        pos = x.pos + x.vel*dt; # position of the center of mass.
        θ   = x.θ + x.ω*dt;     # angular velocity of the center of mass
        vel = x.vel - x.μ/x.m*x.vel*dt + 1/x.m*[u.fx,u.fy]*dt;
        ω   = x.ω - 1/x.J*dot(x.r,[sin(x.θ),cos(x.θ)])*u.fx*dt + 1/x.J*dot(x.r,[cos(x.θ),-sin(x.θ)])*u.fy*dt + 1/x.J*u.tr*dt;
        m = x.m;
        J = x.J;
        r = x.r;
        μ = x.μ;
        t = round(x.t + dt, digits=5);
        return PhysManipulate2D(t,pos,θ,vel,ω,m,J,r,μ,2)
    end
end
function trans(model::TransModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real},dt::Real,
               Q::Matrix{Float64},rng::MersenneTwister)::PhysManipulate2D{T} where T <: Real
# Euler Approximation with stochastic disturbance.
    if x.t != u.t
        error(ArgumentError("PhysState and Control have inconsistent time parameters."))
    else
        vect = rand(rng,MvNormal(zeros(11),Q*dt));
        vect[1:2] += x.pos + x.vel*dt; # position of the center of mass.
        vect[3]   += x.θ + x.ω*dt;     # angular velocity of the center of mass
        vect[4:5] += x.vel - x.μ/x.m*x.vel*dt + 1/x.m*[u.fx,u.fy]*dt;
        vect[6]   += x.ω - 1/x.J*dot(x.r,[sin(x.θ),cos(x.θ)])*u.fx*dt + 1/x.J*dot(x.r,[cos(x.θ),-sin(x.θ)])*u.fy*dt + 1/x.J*u.tr*dt;
        vect[7]   += x.m;
        vect[8]   += x.J;
        vect[9:10] += x.r;
        vect[11]  += x.μ;
        t = round(x.t + dt, digits=5);
        return PhysManipulate2D(t,vect)
    end
end
function trans_jacobi(model::TransModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real}) where T <: Real
    if x.t != u.t
        error(ArgumentError("PhysState and Control have inconsistent time parameters."))
    else
        Jacobian = zeros(T,11,11);
        Jacobian[1,4] = 1.;
        Jacobian[2,5] = 1.;
        Jacobian[3,6] = 1.;
        Jacobian[4,4] = -x.μ/x.m;
        Jacobian[4,7] = x.μ/x.m^2*x.vel[1] - 1/x.m^2*u.fx;
        Jacobian[4,11]= -1/x.m*x.vel[1];
        Jacobian[5,5] = -x.μ/x.m;
        Jacobian[5,7] = x.μ/x.m^2*x.vel[2] - 1/x.m^2*u.fy;
        Jacobian[5,11] = -1/x.m*x.vel[2];
        Jacobian[6,3] = -1/x.J*dot(x.r,[cos(x.θ),-sin(x.θ)])*u.fx + 1/x.J*dot(x.r,[-sin(x.θ),-cos(x.θ)])*u.fy;
        Jacobian[6,8] = 1/x.J^2*dot(x.r,[sin(x.θ),cos(x.θ)])*u.fx - 1/x.J^2*dot(x.r,[cos(x.θ),-sin(x.θ)])*u.fy - 1/x.J^2*u.tr;
        Jacobian[6,9] = -1/x.J*sin(x.θ)*u.fx + 1/x.J*cos(x.θ)*u.fy;
        Jacobian[6,10] = -1/x.J*cos(x.θ)*u.fx - 1/x.J*sin(x.θ)*u.fy;
        return Jacobian
    end
end
function trans_jacobi_auto(model::TransModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real}) where T <: Real
    function f(xVec)
        x = PhysManipulate2D(x.t,xVec)
        return vec(trans(model,x,u,1.));
    end
    return ForwardDiff.jacobian(f,vec(x)) - Matrix(1.0I, 11, 11)
end
function trans_jacobi_euler(model::TransModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real},dt::Real) where T <: Real
    return Matrix{T}(1.0I, 11, 11) + trans_jacobi(model,x,u)*dt;
end
function trans_u_coeff(model::TransModel_Manipulate2D,x::PhysManipulate2D{T}) where T <: Real
    H = zeros(T,11,3);
    H[4,1] = 1/x.m;
    H[5,2] = 1/x.m;
    H[6,1] = -1/x.J*dot(x.r,[sin(x.θ),cos(x.θ)]);
    H[6,2] = 1/x.J*dot(x.r,[cos(x.θ),-sin(x.θ)]);
    H[6,3] = 1/x.J;
    return H
end

#### Observation Models ####
abstract type ObserveModel end

struct ObserveModel_Range <: ObserveModel
# Range-only measurements from relative positions.
end
function observe(model::ObserveModel_Range,p_robot::PhysPos,q_pos::Vector{<:Real})
    # Noiseless Observation
    return norm(p_robot.pos - q_pos)
end
function observe(model::ObserveModel_Range,p_robot::PhysPos,q_pos::Vector{<:Real},
                 R::Matrix{Float64},rng::MersenneTwister)
    # Observation with non-additive Gaussian noise.
    return norm(p_robot.pos - q_pos - rand(rng,MvNormal(zeros(p_robot.dim),R)))
end
function observe_jacobi(model::ObserveModel_Range,p_robot::PhysPos{T},q_pos::Vector{<:Real}) where T <: Real
    return transpose(q_pos - p_robot.pos)./(norm(q_pos - p_robot.pos) + 1e-15); # Prevent singularity.
end
function observe_jacobi_auto(model::ObserveModel_Range,p_robot::PhysPos{T},q_pos::Vector{<:Real}) where T <: Real
    return ForwardDiff.gradient(q_pos -> observe(model,p_robot,q_pos), q_pos)'
end

struct ObserveModel_Manipulate2D <: ObserveModel
end
function observe(model::ObserveModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real})::Vector{T} where T <: Real
    if x.t != u.t
        error(ArgumentError("PhysState and Control have inconsistent time parameters."))
    else
        px_robot = x.pos[1] + x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ);
        py_robot = x.pos[2] + x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ);
        θ = x.θ;
        vx_robot = x.vel[1] - x.ω*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ));
        vy_robot = x.vel[2] + x.ω*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ));
        ω = x.ω;
        α_robot = 1/x.J*u.tr - 1/x.J*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ))*u.fx + 1/x.J*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ))*u.fy;
        ax_robot = 1/x.m*(u.fx - x.μ*x.vel[1]) - α_robot*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ)) - x.ω^2*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ));
        ay_robot = 1/x.m*(u.fy - x.μ*x.vel[2]) + α_robot*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ)) - x.ω^2*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ));
        return [px_robot,py_robot,θ,vx_robot,vy_robot,ω,ax_robot,ay_robot,α_robot]
    end
end
function observe(model::ObserveModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real},
                 R::Matrix{Float64},rng::MersenneTwister)::Vector{T} where T <: Real
    obs = observe(model,x,u);
    obs =  obs + rand(rng,MvNormal(zeros(length(obs)),R));
    return obs
end
function observe_jacobi(model::ObserveModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D) where T <: Real
    if x.t != u.t
        error(ArgumentError("PhysState and Control have inconsistent time parameters."))
    else
        Jacobian = zeros(T,9,11);
        Jacobian[1,1] = 1.;
        Jacobian[1,3] = x.r[1]*(-sin(x.θ)) - x.r[2]*cos(x.θ);
        Jacobian[1,9] = cos(x.θ);
        Jacobian[1,10] = -sin(x.θ);
        Jacobian[2,2] = 1.;
        Jacobian[2,3] = x.r[1]*cos(x.θ) + x.r[2]*(-sin(x.θ));
        Jacobian[2,9] = sin(x.θ);
        Jacobian[2,10] = cos(x.θ);
        Jacobian[3,3] = 1.;
        Jacobian[4,3] = -x.ω*x.r[1]*cos(x.θ) - x.ω*x.r[2]*(-sin(x.θ));
        Jacobian[4,4] = 1.;
        Jacobian[4,6] = -x.r[1]*sin(x.θ) - x.r[2]*cos(x.θ);
        Jacobian[4,9] = -x.ω*sin(x.θ);
        Jacobian[4,10] = -x.ω*cos(x.θ);
        Jacobian[5,3] = x.ω*x.r[1]*(-sin(x.θ)) - x.ω*x.r[2]*cos(x.θ);
        Jacobian[5,5] = 1.;
        Jacobian[5,6] = x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ);
        Jacobian[5,9] = x.ω*cos(x.θ);
        Jacobian[5,10] = -x.ω*sin(x.θ);
        Jacobian[6,6] = 1.;
        Jacobian[9,3] = -1/x.J*(x.r[1]*cos(x.θ) + x.r[2]*(-sin(x.θ)))*u.fx + 1/x.J*(x.r[1]*(-sin(x.θ)) - x.r[2]*cos(x.θ))*u.fy;
        Jacobian[9,8] = -1/x.J^2*u.tr + 1/x.J^2*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ))*u.fx - 1/x.J^2*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ))*u.fy;
        Jacobian[9,9] = -1/x.J*sin(x.θ)*u.fx + 1/x.J*cos(x.θ)*u.fy;
        Jacobian[9,10] = -1/x.J*cos(x.θ)*u.fx - 1/x.J*sin(x.θ)*u.fy;
        α_robot = 1/x.J*u.tr - 1/x.J*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ))*u.fx + 1/x.J*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ))*u.fy;
        Jacobian[7,3] = -Jacobian[9,3]*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ)) - α_robot*(x.r[1]*cos(x.θ) + x.r[2]*(-sin(x.θ))) - x.ω^2*(x.r[1]*(-sin(x.θ)) - x.r[2]*cos(x.θ));
        Jacobian[7,4] = -x.μ/x.m;
        Jacobian[7,6] = -2*x.ω*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ));
        Jacobian[7,7] = -1/x.m^2*(u.fx - x.μ*x.vel[1]);
        Jacobian[7,8] = -Jacobian[9,8]*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ));
        Jacobian[7,9] = -Jacobian[9,9]*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ)) - α_robot*sin(x.θ) - x.ω^2*cos(x.θ);
        Jacobian[7,10] = -Jacobian[9,10]*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ)) - α_robot*cos(x.θ) + x.ω^2*sin(x.θ);
        Jacobian[7,11] = -1/x.m*x.vel[1];
        Jacobian[8,3] = Jacobian[9,3]*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ)) + α_robot*(x.r[1]*(-sin(x.θ)) - x.r[2]*cos(x.θ)) - x.ω^2*(x.r[1]*cos(x.θ) + x.r[2]*(-sin(x.θ)));
        Jacobian[8,5] = -x.μ/x.m;
        Jacobian[8,6] = -2*x.ω*(x.r[1]*sin(x.θ) + x.r[2]*cos(x.θ));
        Jacobian[8,7] = -1/x.m^2*(u.fy - x.μ*x.vel[2]);
        Jacobian[8,8] = Jacobian[9,8]*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ));
        Jacobian[8,9] = Jacobian[9,9]*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ)) + α_robot*cos(x.θ) - x.ω^2*sin(x.θ);
        Jacobian[8,10] = Jacobian[9,10]*(x.r[1]*cos(x.θ) - x.r[2]*sin(x.θ)) + α_robot*(-sin(x.θ)) - x.ω^2*cos(x.θ);
        Jacobian[8,11] = -1/x.m*x.vel[2];
        return Jacobian
    end
end
function observe_jacobi_auto(model::ObserveModel_Manipulate2D,x::PhysManipulate2D{T},u::MControl2D{<:Real}) where T <: Real
    function f(xVec)
        x = PhysManipulate2D(x.t,xVec);
        return vec(observe(model,x,u));
    end
    return ForwardDiff.jacobian(f,vec(x))
end


text:
using Random
using Dates
using Printf
using MPI

mutable struct MonteCarloStatistics
    sweeps::Int
    attemptedLocalUpdates::Int
    acceptedLocalUpdates::Int
    attemptedReplicaExchanges::Int
    acceptedReplicaExchanges::Int
    initializationTime::Float64

    MonteCarloStatistics() = new(0, 0, 0, 0, 0, time())
end

mutable struct MonteCarlo{T<:Lattice}
    lattice::T

    beta::Float64
    thermalizationSweeps::Int
    measurementSweeps::Int
    measurementRate::Int
    replicaExchangeRate::Int
    reportInterval::Int
    checkpointInterval::Int

    rng::MersenneTwister
    seed::UInt
    sweep::Int

    observables::Observables
end

function MonteCarlo(
    lattice::T,
    beta::Float64,
    thermalizationSweeps::Int,
    measurementSweeps::Int;
    measurementRate::Int = 1,
    replicaExchangeRate::Int = 10,
    reportInterval::Int = round(Int, 0.05 * (thermalizationSweeps + measurementSweeps)),
    checkpointInterval::Int = 3600,
    rng::MersenneTwister = copy(Random.GLOBAL_RNG),
    seed::UInt = rand(Random.RandomDevice(),UInt)
    ) where T<:Lattice

    mc = MonteCarlo(deepcopy(lattice), beta, thermalizationSweeps, measurementSweeps, measurementRate, replicaExchangeRate, reportInterval, checkpointInterval, rng, seed, 0, Observables(lattice))
    Random.seed!(mc.rng, mc.seed)

    return mc
end

function run!(mc::MonteCarlo{T}; outfile::Union{String,Nothing}=nothing) where T<:Lattice
    #init MPI
    rank = 0
    commSize = 1
    allBetas = zeros(0)
    enableMPI = false
    if MPI.Initialized()
        commSize = MPI.Comm_size(MPI.COMM_WORLD)
        rank = MPI.Comm_rank(MPI.COMM_WORLD)
        if commSize > 1
            allBetas = zeros(commSize)
            allBetas[rank + 1] = mc.beta
            MPI.Allgather!(UBuffer(allBetas, 1), MPI.COMM_WORLD)
            enableMPI = true
            rank == 0 && @printf("MPI detected. Enabling replica exchanges across %d simulations.\n", commSize)
        end
    end

    #init IO
    enableOutput = typeof(outfile) != Nothing
    if enableOutput
        enableMPI && (outfile *= "." * string(rank))
        isfile(outfile) && error("File ", outfile, " already exists. Terminating.")
    end

    #init spin configuration
    if mc.sweep == 0
        for i in 1:length(mc.lattice)
            setSpin!(mc.lattice, i, uniformOnSphere(mc.rng))
        end
    end

    #init Monte Carlo run
    totalSweeps = mc.thermalizationSweeps + mc.measurementSweeps
    partnerSpinConfiguration = deepcopy(mc.lattice.spins)
    energy = getEnergy(mc.lattice)

    #launch Monte Carlo run
    lastCheckpointTime = time()
    statistics = MonteCarloStatistics()
    rank == 0 && @printf("Simulation started on %s.\n\n", Dates.format(Dates.now(), "dd u yyyy HH:MM:SS"))

    while mc.sweep < totalSweeps
        #perform local sweep
        for i in 1:length(mc.lattice)
            #select random spin
            site = rand(mc.rng, 1:length(mc.lattice))

            #propose new spin configuration
            newSpinState = uniformOnSphere(mc.rng)
            energyDifference = getEnergyDifference(mc.lattice, site, newSpinState)

            #check acceptance of new configuration
            statistics.attemptedLocalUpdates += 1
            p = exp(-mc.beta * energyDifference)
            if (rand(mc.rng) < min(1.0, p))
                setSpin!(mc.lattice, site, newSpinState)
                energy += energyDifference
                statistics.acceptedLocalUpdates += 1
            end
        end
        statistics.sweeps += 1

        #perform replica exchange
        if enableMPI && mc.sweep % mc.replicaExchangeRate == 0
            #determine MPI rank to exchagne configuration with
            if iseven(mc.sweep ÷ mc.replicaExchangeRate)
                partnerRank = iseven(rank) ? rank + 1 : rank - 1
            else
                partnerRank = iseven(rank) ? rank - 1 : rank + 1
            end

            if partnerRank >= 0 && partnerRank < commSize
                #obtain energy of new configuration
                partnerEnergy = MPISendrecvFloat(energy, partnerRank, MPI.COMM_WORLD)

                #check acceptance of new configuration
                statistics.attemptedReplicaExchanges += 1
                exchangeAccepted = false
                if iseven(rank)
                    p = exp(-(allBetas[rank + 1] - allBetas[partnerRank + 1]) * (partnerEnergy - energy))
                    exchangeAccepted = (rand(mc.rng) < min(1.0, p)) ? true : false
                    MPISendBool(exchangeAccepted, partnerRank, MPI.COMM_WORLD)
                else
                    exchangeAccepted = MPIRecvBool(partnerRank, MPI.COMM_WORLD)
                end
                if (exchangeAccepted)
                    energy = partnerEnergy
                    MPI.Sendrecv!(mc.lattice.spins, partnerRank, 0, partnerSpinConfiguration, partnerRank, 0, MPI.COMM_WORLD)
                    (mc.lattice.spins, partnerSpinConfiguration) = (partnerSpinConfiguration, mc.lattice.spins)
                    statistics.acceptedReplicaExchanges += 1
                end
            end
        end

        #perform measurement
        if mc.sweep >= mc.thermalizationSweeps
            if mc.sweep % mc.measurementRate == 0
                performMeasurements!(mc.observables, mc.lattice, energy)
            end
        end

        #increment sweep
        statistics.sweeps += 1
        mc.sweep += 1

        #runtime statistics
        t = time()
        if mc.sweep % mc.reportInterval == 0
            #collect statistics
            progress = 100.0 * mc.sweep / totalSweeps
            thermalized = (mc.sweep >= mc.thermalizationSweeps) ? "YES" : "NO"
            sweeprate = statistics.sweeps / (t - statistics.initializationTime)
            sweeptime = 1.0 / sweeprate
            eta = (totalSweeps - mc.sweep) / sweeprate

            localUpdateAcceptanceRate = 100.0 * statistics.acceptedLocalUpdates / statistics.attemptedLocalUpdates
            if enableMPI
                replicaExchangeAcceptanceRate = 100.0 * statistics.acceptedReplicaExchanges / statistics.attemptedReplicaExchanges
                allLocalAppectanceRate = zeros(commSize)
                allLocalAppectanceRate[rank + 1] = localUpdateAcceptanceRate
                MPI.Allgather!(UBuffer(allLocalAppectanceRate, 1), MPI.COMM_WORLD)
                allReplicaExchangeAcceptanceRate = zeros(commSize)
                allReplicaExchangeAcceptanceRate[rank + 1] = replicaExchangeAcceptanceRate
                MPI.Allgather!(UBuffer(allReplicaExchangeAcceptanceRate, 1), MPI.COMM_WORLD)
            end

            #print statistics
            if rank == 0
                str = ""
                str *= @sprintf("Sweep %d / %d (%.1f%%)", mc.sweep, totalSweeps, progress)
                str *= @sprintf("\t\tETA : %s\n", Dates.format(Dates.now() + Dates.Second(round(Int64,eta)), "dd u yyyy HH:MM:SS"))
                str *= @sprintf("\t\tthermalized : %s\n", thermalized)
                str *= @sprintf("\t\tsweep rate : %.1f sweeps/s\n", sweeprate)
                str *= @sprintf("\t\tsweep duration : %.3f ms\n", sweeptime * 1000)

                if enableMPI
                    for n in 1:commSize
                        str *= @sprintf("\t\tsimulation %d update acceptance rate: %.2f%%\n", n - 1, allLocalAppectanceRate[n])
                        str *= @sprintf("\t\tsimulation %d replica exchange acceptance rate : %.2f%%\n", n - 1, allReplicaExchangeAcceptanceRate[n])
                    end
                else
                    str *= @sprintf("\t\tupdate acceptance rate: %.2f%%\n", localUpdateAcceptanceRate)
                end
                str *= @sprintf("\n")
                print(str)
            end

            #reset statistics
            statistics = MonteCarloStatistics()
        end

        #write checkpoint
        if enableOutput
            checkpointPending = time() - lastCheckpointTime >= mc.checkpointInterval
            enableMPI && (checkpointPending = MPIBcastBool(checkpointPending, 0, MPI.COMM_WORLD))
            if checkpointPending
                writeMonteCarlo(outfile, mc)
                lastCheckpointTime = time()
                rank == 0 && @printf("Checkpoint written on %s.\n", Dates.format(Dates.now(), "dd u yyyy HH:MM:SS"))
            end
        end
    end

    #write final checkpoint
    if enableOutput
        writeMonteCarlo(outfile, mc)
        rank == 0 && @printf("Checkpoint written on %s.\n", Dates.format(Dates.now(), "dd u yyyy HH:MM:SS"))
    end

    #return
    rank == 0 && @printf("Simulation finished on %s.\n", Dates.format(Dates.now(), "dd u yyyy HH:MM:SS"))
    return nothing
end


text:
@recipe function f(pt::CartesianPoint)
    xguide --> "x"
    yguide --> "y"
    zguide --> "z"
    unitformat --> :slash
    if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
        aspect_ratio --> 1.0
    end
    @series begin
        seriestype --> :scatter
        [pt.x]*internal_length_unit, [pt.y]*internal_length_unit, [pt.z]*internal_length_unit
    end
end

@recipe function f(pt::CylindricalPoint)
    @series begin
        CartesianPoint(pt)
    end
end

@recipe function f(v::AbstractVector{<:CartesianPoint})
    xguide --> "x"
    yguide --> "y"
    zguide --> "z"
    unitformat --> :slash
    if occursin("GRBackend", string(typeof(plotattributes[:plot_object].backend)))
        aspect_ratio --> 1.0
    end
    @series begin
        seriestype --> :scatter
        [v[i].x for i in eachindex(v)]*internal_length_unit, [v[i].y for i in eachindex(v)]*internal_length_unit, [v[i].z for i in eachindex(v)]*internal_length_unit
    end
end

@recipe function f(v::AbstractVector{<:CylindricalPoint})
    @series begin
        map(x -> CartesianPoint(x), v)
    end
end


text:
using SYNJULIA
using Test


text:
module RankWrapper

using JavaCall

const AbstractStorelessUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.AbstractStorelessUnivariateStatistic
const AbstractUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.AbstractUnivariateStatistic
const AggregatableStatistic = @jimport org.hipparchus.stat.descriptive.AggregatableStatistic
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const DoubleConsumer = @jimport java.util.function.DoubleConsumer
const Iterable = @jimport java.lang.Iterable
const JString = @jimport java.lang.String
const KthSelector = @jimport org.hipparchus.util.KthSelector
const List = @jimport java.util.List
const Max = @jimport org.hipparchus.stat.descriptive.rank.Max
const Median = @jimport org.hipparchus.stat.descriptive.rank.Median
const Min = @jimport org.hipparchus.stat.descriptive.rank.Min
const NaNStrategy = @jimport org.hipparchus.stat.ranking.NaNStrategy
const Object = @jimport java.lang.Object
const PSquarePercentile = @jimport org.hipparchus.stat.descriptive.rank.PSquarePercentile
const PSquarePercentile_PSquareMarkers = @jimport org.hipparchus.stat.descriptive.rank.PSquarePercentile$PSquareMarkers
const Percentile = @jimport org.hipparchus.stat.descriptive.rank.Percentile
const Percentile_EstimationType = @jimport org.hipparchus.stat.descriptive.rank.Percentile$EstimationType
const PivotingStrategy = @jimport org.hipparchus.util.PivotingStrategy
const RandomGenerator = @jimport org.hipparchus.random.RandomGenerator
const RandomPercentile = @jimport org.hipparchus.stat.descriptive.rank.RandomPercentile
const StorelessUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.StorelessUnivariateStatistic
const UnivariateStatistic = @jimport org.hipparchus.stat.descriptive.UnivariateStatistic

export Max
export Median
export Min
export PSquarePercentile
export Percentile
export RandomPercentile
export accept
export aggregate
export and_then
export clear
export copy
export equals
export evaluate
export get_aggregate_n
export get_aggregate_quantile_rank
export get_aggregate_rank
export get_class
export get_data
export get_estimation_type
export get_kth_selector
export get_n
export get_na_n_strategy
export get_pivoting_strategy
export get_quantile
export get_quantile_rank
export get_rank
export get_result
export hash_code
export increment
export increment_all
export max_values_retained
export new_markers
export notify
export notify_all
export quantile
export reduce
export set_data
export set_quantile
export to_string
export wait
export with_estimation_type
export with_kth_selector
export with_na_n_strategy

include("max.jl")
include("median.jl")
include("min.jl")
include("p_square_percentile.jl")
include("percentile.jl")
include("random_percentile.jl")

end