text:
using Test

include("rational-numbers.jl")

@test RationalNumber <: Real
@test_throws ArgumentError RationalNumber(0, 0)

@testset "One- & Zero-elements" begin
    @test zero(RationalNumber{Int}) == RationalNumber(0, 1)
    @test one(RationalNumber{Int})  == RationalNumber(1, 1)
end

@testset "Arithmetic" begin
    @testset "Addition" begin
        @test RationalNumber( 1, 2) + RationalNumber( 2, 3) == RationalNumber( 7, 6)
        @test RationalNumber( 1, 2) + RationalNumber(-2, 3) == RationalNumber(-1, 6)
        @test RationalNumber(-1, 2) + RationalNumber(-2, 3) == RationalNumber(-7, 6)
        @test RationalNumber( 1, 2) + RationalNumber(-1, 2) == RationalNumber( 0, 1)
    end

    @testset "Subtraction" begin
        @test RationalNumber( 1, 2) - RationalNumber( 2, 3) == RationalNumber(-1, 6)
        @test RationalNumber( 1, 2) - RationalNumber(-2, 3) == RationalNumber( 7, 6)
        @test RationalNumber(-1, 2) - RationalNumber(-2, 3) == RationalNumber( 1, 6)
        @test RationalNumber( 1, 2) - RationalNumber( 1, 2) == RationalNumber( 0, 1)
    end

    @testset "Multiplication" begin
        @test RationalNumber( 1, 2) * RationalNumber( 2, 3) == RationalNumber( 1, 3)
        @test RationalNumber(-1, 2) * RationalNumber( 2, 3) == RationalNumber(-1, 3)
        @test RationalNumber(-1, 2) * RationalNumber(-2, 3) == RationalNumber( 1, 3)
        @test RationalNumber( 1, 2) * RationalNumber( 2, 1) == RationalNumber( 1, 1)
        @test RationalNumber( 1, 2) * RationalNumber( 1, 1) == RationalNumber( 1, 2)
        @test RationalNumber( 1, 2) * RationalNumber( 0, 1) == RationalNumber( 0, 1)
    end

    @testset "Exponentiation" begin
        @testset "Exponentiation of a rational number" begin
            @test RationalNumber( 1, 2)^3 == RationalNumber( 1, 8)
            @test RationalNumber(-1, 2)^3 == RationalNumber(-1, 8)
            @test RationalNumber( 0, 1)^5 == RationalNumber( 0, 1)
            @test RationalNumber( 1, 1)^4 == RationalNumber( 1, 1)
            @test RationalNumber( 1, 2)^0 == RationalNumber( 1, 1)
            @test RationalNumber(-1, 2)^0 == RationalNumber( 1, 1)
        end

        @testset "Exponentiation of a real number to a rational number" begin
            @test 8^RationalNumber( 4, 3) ≈ 15.999999999999998
            @test 9^RationalNumber(-1, 2) ≈ 0.3333333333333333
            @test 2^RationalNumber( 0, 1) ≈ 1.0
        end
    end

    @testset "Division" begin
        @test RationalNumber( 1, 2) / RationalNumber( 2, 3) == RationalNumber( 3, 4)
        @test RationalNumber( 1, 2) / RationalNumber(-2, 3) == RationalNumber(-3, 4)
        @test RationalNumber(-1, 2) / RationalNumber(-2, 3) == RationalNumber( 3, 4)
        @test RationalNumber( 1, 2) / RationalNumber( 1, 1) == RationalNumber( 1, 2)
    end
end

@testset "Absolute value" begin
    @test abs(RationalNumber( 1,  2)) == RationalNumber(1, 2)
    @test abs(RationalNumber(-1,  2)) == RationalNumber(1, 2)
    @test abs(RationalNumber( 0,  1)) == RationalNumber(0, 1)
    @test abs(RationalNumber( 1, -2)) == RationalNumber(1, 2)
end

@testset "Reduction to lowest terms" begin
    r = RationalNumber(2, 4)
    @test numerator(r)   == 1
    @test denominator(r) == 2

    r = RationalNumber(-4, 6)
    @test numerator(r)   == -2
    @test denominator(r) ==  3

    r = RationalNumber(3, -9)
    @test numerator(r)   == -1
    @test denominator(r) ==  3

    r = RationalNumber(0, 6)
    @test numerator(r)   == 0
    @test denominator(r) == 1

    r = RationalNumber(-14, 7)
    @test numerator(r)   == -2
    @test denominator(r) ==  1

    r = RationalNumber(13, 13)
    @test numerator(r)   == 1
    @test denominator(r) == 1

    r = RationalNumber(1, -1)
    @test numerator(r)   == -1
    @test denominator(r) ==  1
end

# TODO add to problem spec
# The following testset is based on the tests for rational numbers in Julia Base (MIT license)
# https://github.com/JuliaLang/julia/blob/52bafeb981bac548afd2264edb518d8d86944dca/test/rational.jl
# https://github.com/JuliaLang/julia/blob/52bafeb981bac548afd2264edb518d8d86944dca/LICENSE.md
@testset "Ordering" begin
    for a in -5:5, b in -5:5, c in -5:5
        a == b == 0 && continue

        r = RationalNumber(a, b)

        @test (r == c) == (a / b == c)
        @test (r != c) == (a / b != c)
        @test (r <= c) == (a / b <= c)
        @test (r <  c) == (a / b <  c)
        @test (r >= c) == (a / b >= c)
        @test (r >  c) == (a / b >  c)

        for d in -5:5
            c == d == 0 && continue

            s = RationalNumber(c, d)

            @test (r == s) == (a / b == c / d)
            @test (r != s) == (a / b != c / d)
            @test (r <= s) == (a / b <= c / d)
            @test (r <  s) == (a / b <  c / d)
            @test (r >= s) == (a / b >= c / d)
            @test (r >  s) == (a / b >  c / d)
        end
    end
end

@testset "Showing RationalNumbers" begin
    @test sprint(show, RationalNumber(23, 42)) == "23//42"
    @test sprint(show, RationalNumber(-2500, 5000)) == "-1//2"
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
abstract type AbstractIntervalInput end

"""
    struct ParamIntervalInput <: AbstractIntervalInput
        theta_init::Vector{Float64} # initial parameters vector
        theta_num::Int              # number of the parameter for identifiability analysis
        loss_func::Function         # loss function
        method::Symbol
        options::Any
    end

Structure storing input data for parameter identification
"""
struct ParamIntervalInput <: AbstractIntervalInput
    theta_init::Vector{Float64} # initial parameters vector
    theta_num::Int # number of the parameter for analysis
    loss_func::Function # loss function
    method::Symbol
    options::Dict{Symbol, Any}
end

struct PredictionIntervalInput <: AbstractIntervalInput
    theta_init::Vector{Float64}
    scan_func::Function
    loss_func::Function
    method::Symbol
    options::Dict{Symbol, Any}
end

"""
    struct ParamInterval
        input::ParamIntervalInput
        loss_init::Float64
        method::Symbol
        result::Tuple{EndPoint, EndPoint}
    end

Structure storing result of parameter identification
"""
struct ParamInterval
    input::AbstractIntervalInput
    loss_init::Float64
    method::Symbol
    result::Tuple{EndPoint, EndPoint}
end

"""
    function get_interval(
        theta_init::Vector{Float64},
        theta_num::Int,
        loss_func::Function,
        method::Symbol;

        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        scan_bounds::Tuple{Float64,Float64} = unscaling.(
            (-9.0, 9.0),
            scale[theta_num]
            ),
        kwargs...
    )
Computes confidence interval for single component `theta_num` of parameter vector.

## Return
`ParamInterval` structure storing all input data and estimated confidence interval.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values should not necessary be the optimum values of `loss_func` but `loss_func(theta_init)` should be lower than `loss_crit`.
- `theta_num`: index of vector component for identification: `theta_init(theta_num)`.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` for profile likelihood-based (PL) identification. Usually we use log-likelihood for PL analysis: ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.
- `method`: computational method to estimate confidence interval's endpoint. Currently the following methods are implemented: `:CICO_ONE_PASS`, `:LIN_EXTRAPOL`, `:QUADR_EXTRAPOL`.

## Keyword arguments
- `scale`: vector of scale transformations for each parameters' component. Possible values: `:direct` (`:lin`), `:log`, `:logit`. This option can speed up the optimization, especially for wide `theta_bounds`. The default value is `:direct` (no transformation) for all parameters.
- `scan_bounds`: scan bounds tuple for `theta_num` parameter. Should be within the `theta_bounds` for `theta_num` parameter. Default is `(-9.,9.)` for `:direct` scales and `(1e-9, 1e+9)` for `:log`.
- `kwargs...`: the additional arguments passed to [`get_endpoint`](@ref)
"""
function get_interval(
    theta_init::Vector{Float64},
    theta_num::Int,
    loss_func::Function,
    method::Symbol;

    scale::Vector{Symbol} = fill(:direct, length(theta_init)),
    scan_bounds::Tuple{Float64,Float64} = unscaling.(
        (-9.0, 9.0),
        scale[theta_num]
        ),
    kwargs... # other options for get_right_endpoint
)
    # both endpoints
    endpoints = [get_endpoint(
        theta_init,
        theta_num,
        loss_func,
        method,
        [:left,:right][i]; # method

        scale,
        scan_bound = scan_bounds[i],
        kwargs... # options for local fitter
        ) for i in 1:2]

    input = ParamIntervalInput(
        theta_init,
        theta_num,
        loss_func,
        method,
        Dict(:scale=>scale, :scan_bounds=>scan_bounds, kwargs...) # ? NamedTuple
        )

    ParamInterval(
        input,
        loss_func(theta_init),
        method,
        Tuple(endpoints)
        )
end

"""
    get_interval(
        theta_init::Vector{Float64},
        scan_func::Function,
        loss_func::Function,
        method::Symbol;

        loss_crit::Float64 = 0.0,
        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
            fill((-Inf, Inf), length(theta_init)),
            scale
            ),
        scan_bounds::Tuple{Float64,Float64} = unscaling.(
            (-9.0, 9.0),
            :direct
            ),
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 0.,
        local_alg::Symbol = :LN_NELDERMEAD,
        autodiff::Bool = true,
        kwargs...
        )
Computes confidence interval for function of parameters `scan_func`.

## Return
`ParamInterval` structure storing all input data and estimated confidence interval.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values should not necessary be the optimum values of `loss_func` but `loss_func(theta_init)` should be lower than `loss_crit`.
- `scan_func`: scan function of parameters.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` the profile of which is analyzed. Usually we use log-likelihood for profile analysis in form ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.
- `method`: computational method to estimate confidence interval's endpoint. Currently supports only `:CICO_ONE_PASS` method.

## Keyword arguments
- `scale`: vector of scale transformations for each parameters' component. Possible values: `:direct` (`:lin`), `:log`, `:logit`. This option can speed up the optimization, especially for wide `theta_bounds`. The default value is `:direct` (no transformation) for all parameters.
- `scan_bounds`: scan bounds tuple for `scan_func` values. Default is `(1e-9, 1e9)` .
- `kwargs...`: the additional arguments passed to [`get_endpoint`](@ref)
"""
function get_interval(
    theta_init::Vector{Float64},
    scan_func::Function,
    loss_func::Function,
    method::Symbol;

    scale::Vector{Symbol} = fill(:direct, length(theta_init)),
    scan_bounds::Tuple{Float64,Float64} = unscaling.(
        (-9.0, 9.0),
        scale[theta_num]
        ),
    kwargs... # other options for get_right_endpoint
)
    # both endpoints
    endpoints = [get_endpoint(
        theta_init,
        scan_func,
        loss_func,
        method,
        [:left,:right][i]; # method

        scale,
        scan_bound = scan_bounds[i],
        kwargs...
        ) for i in 1:2]

    input = PredictionIntervalInput(
        theta_init,
        scan_func,
        loss_func,
        method,
        Dict(:scale=>scale, :scan_bounds=>scan_bounds, kwargs...) # ? NamedTuple
        )

    ParamInterval(
        input,
        loss_func(theta_init),
        method,
        Tuple(endpoints)
        )
end

ext: jl


text:
using ProjAssistant
@quickactivate

@time begin
    import Chemostat_InSilico
    const Dyn = Chemostat_InSilico.Dynamic

    using Base.Threads
    using Plots
    import GR
    !isinteractive() && GR.inline("png")
end

## ------------------------------------------------------
# common globals
let
    sglob(Dyn;
        δ = 450 # the number of discrete bins
    )
end

## ------------------------------------------------------
@time let
    # input
    net2D = Dyn.ToyModel2D()

    freeids = [:z, :ug]
    @lglob Dyn δ

    @info("Building Vcell2D", nthreads())
    @show freeids δ
    fLb, fUb = Dyn.bounds(net2D, freeids)
    filter(net, v) = all(fLb .< v .< fUb) # avoid artifacts at extremes
    Vcell2D = Dyn.Space(net2D, freeids, δ; filter)
    @sglob Dyn Vcell2D net2D

    ps = Dyn.pol_plots(Vcell2D; bins = 10000)
    sfig(Dyn, ps, "Vcell2D.png")
end

## ------------------------------------------------------
@time let
    # input
    net3D = Dyn.ToyModel3D()

    freeids = [:z, :ug, :uo]
    @lglob Dyn δ

    @info("Building Vcell3D", nthreads())
    @show freeids δ
    fLb, fUb = Dyn.bounds(net3D, freeids)
    filter(net, v) = all(fLb .< v .< fUb) # avoid artifacts at extremes
    Vcell3D = Dyn.Space(net3D, freeids, δ; filter)
    @sglob Dyn Vcell3D net3D

    ps = Dyn.pol_plots(Vcell3D; bins = 10000)
    sfig(Dyn, ps, "Vcell3D.png")
end

## ------------------------------------------------------
# Test dimensionality
let
    net3D = Dyn.ToyModel3D()
    tol = 1e-3

    Dyn.fix!(net3D, :z, 0.3)
    Dyn.reducebox!(net3D)
    @assert any(net3D.ub .- net3D.lb .> tol)

    lb, ub = Dyn.bounds(net3D, :ug)
    Dyn.fix!(net3D, :ug, lb + (ub - lb) / 2)
    Dyn.reducebox!(net3D)
    @assert any(net3D.ub .- net3D.lb .> tol)

    lb, ub = Dyn.bounds(net3D, :uo)
    Dyn.fix!(net3D, :uo, lb + (ub - lb) / 2)
    Dyn.reducebox!(net3D)
    @assert all(net3D.ub .- net3D.lb .≈ 0.0)
    println(net3D)
end


text:
@testset "Ylm" begin

function explicit_shs(θ, φ)
   Y00 = 0.5 * sqrt(1/π)
   Y1m1 = 0.5 * sqrt(3/(2*π))*sin(θ)*exp(-im*φ)
   Y10 = 0.5 * sqrt(3/π)*cos(θ)
   Y11 = -0.5 * sqrt(3/(2*π))*sin(θ)*exp(im*φ)
   Y2m2 = 0.25 * sqrt(15/(2*π))*sin(θ)^2*exp(-2*im*φ)
   Y2m1 = 0.5 * sqrt(15/(2*π))*sin(θ)*cos(θ)*exp(-im*φ)
   Y20 = 0.25 * sqrt(5/π)*(3*cos(θ)^2 - 1)
   Y21 = -0.5 * sqrt(15/(2*π))*sin(θ)*cos(θ)*exp(im*φ)
   Y22 = 0.25 * sqrt(15/(2*π))*sin(θ)^2*exp(2*im*φ)
   Y3m3 = 1/8 * exp(-3 * im * φ) * sqrt(35/π) * sin(θ)^3
   Y3m2 = 1/4 * exp(-2 * im * φ) * sqrt(105/(2*π)) * cos(θ) * sin(θ)^2
   Y3m1 = 1/8 * exp(-im * φ) * sqrt(21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
   Y30 = 1/4 * sqrt(7/π) * (-3 * cos(θ) + 5 * cos(θ)^3)
   Y31 = -(1/8) * exp(im * φ) * sqrt(21/π) * (-1 + 5 * cos(θ)^2) * sin(θ)
   Y32 = 1/4 * exp(2 * im * φ) * sqrt(105/(2*π)) * cos(θ) * sin(θ)^2
   Y33 = -(1/8) * exp(3 * im * φ) * sqrt(35/π) * sin(θ)^3
   return [Y00, Y1m1, Y10, Y11, Y2m2, Y2m1, Y20, Y21, Y22,
           Y3m3, Y3m2, Y3m1, Y30, Y31, Y32, Y33]
end

@info("Test 1: check complex spherical harmonics against explicit expressions")
nsamples = 30
for n = 1:nsamples
   θ = rand() * π
   φ = (rand()-0.5) * 2*π
   Y = compute_y(3, cos(θ), φ)
   Yex = explicit_shs(θ, φ)
   print((@test Y ≈ Yex), " ")
end
println()

##
@info("Test 2: check faster implementation matches")
nsamples = 30
L = 30
for n = 1:nsamples
   θ = rand() * π
   φ = (rand()-0.5) * 2*π
   Y = compute_y(L, cos(θ), φ)
   r = 0.1+rand()
   R = SVector(r*sin(θ)*cos(φ), r*sin(θ)*sin(φ), r*cos(θ))
   Ynew = cYlm_from_cart(L, R)
   print((@test Y ≈ Ynew), " ")
end
println()

# ##  TODO: there is something wrong here! => test errors due to sqrt(neg numbers)
# @info("Test 3: Old vs New Timing")
# L = 15
# θ = rand() * π
# cos_θ = cos(θ)
# φ = (rand()-0.5) * 2*π
# r = 0.1+rand()
# R = SVector(r*sin(θ)*cos(φ), r*sin(θ)*sin(φ), r*cos(θ))
# x, y, z, s = SphericalHarmonics.compute_rxz(R)
# coeff = SphericalHarmonics.compute_coefficients(L)
# P = Vector{Float64}(undef, SphericalHarmonics.sizeP(L))
# SphericalHarmonics.compute_p!(L, z, coeff, P)
# Y = Vector{ComplexF64}(undef, SphericalHarmonics.sizeY(L))
# @info("Old Implementation:")
# @btime compute_y!($L, $cos_θ, $φ, $P, $Y)
# @info("New Implementation:")
# @btime cYlm_from_cart!($Y, $L, $r, $x, $z, $s, $P)
# @info("Experimental Real Spherical Harmonics:")
# Yr = Vector{Float64}(undef, SphericalHarmonics.sizeY(L))
# @btime SphericalHarmonics.rYlm_from_cart!($Y, $L, $r, $x, $z, $s, $P)

end # @testset


text:
module VariableTest

using Test
using ForneyLab

@testset "Variable" begin
    g = FactorGraph()

    # Variable call should construct a Variable that is assigned to graph
    var = Variable(id=:my_var)
    @test isa(var, Variable)
    @test length(g.variables) == 1
    @test ===(g.variables[:my_var], var)
end

@testset "associate!" begin
    FactorGraph()
    var = Variable()

    # Variable should be associated with a new edge
    node1 = Clamp(var, 0.0)
    iface1 = node1.interfaces[1]
    edge1 = iface1.edge
    @test length(var.edges) == 1 # Check Variable
    @test ===(var.edges[1], edge1)
    @test ===(edge1.variable, var) # Check Edge
    @test ===(edge1.a, iface1)
    @test edge1.b == nothing
    @test ===(iface1.edge, edge1) # Check Interface
    @test iface1.partner == nothing

    # Variable should be associated with an existing edge
    node2 = Clamp(var, 0.0)
    iface2 = node2.interfaces[1]
    @test length(var.edges) == 1 # Check Variable
    @test ===(var.edges[1], edge1)
    @test ===(edge1.variable, var) # Check edge
    @test ===(edge1.a, iface1)
    @test ===(edge1.b, iface2)
    @test ===(iface1.edge, edge1) # Check Interfaces
    @test ===(iface2.edge, edge1)
    @test ===(iface1.partner, iface2)
    @test ===(iface2.partner, iface1)

    # Variable should become equality constrained
    node3 = Clamp(var, 0.0)
    eq = node1.interfaces[1].partner.node
    @test isa(eq, Equality) # Check new Equality
    eq_iface1 = eq.interfaces[1]
    eq_iface2 = eq.interfaces[2]
    eq_iface3 = eq.interfaces[3]
    iface3 = node3.interfaces[1]
    edge2 = iface2.edge
    edge3 = iface3.edge
    @test length(var.edges) == 3 # Check Variable
    @test ===(var.edges[1], edge1)
    @test ===(var.edges[2], edge2)
    @test ===(var.edges[3], edge3)
    @test ===(edge1.variable, var) # Check Edges
    @test ===(edge2.variable, var)
    @test ===(edge3.variable, var)
    @test ===(edge1.a, iface1)
    @test ===(edge1.b, eq_iface1)
    @test ===(edge2.a, iface2)
    @test ===(edge2.b, eq_iface2)
    @test ===(edge3.a, iface3)
    @test ===(edge3.b, eq_iface3)
    @test ===(iface1.edge, edge1) # Check Interfaces
    @test ===(eq_iface1.edge, edge1)
    @test ===(iface1.partner, eq_iface1)
    @test ===(eq_iface1.partner, iface1)
    @test ===(iface2.edge, edge2)
    @test ===(eq_iface2.edge, edge2)
    @test ===(iface2.partner, eq_iface2)
    @test ===(eq_iface2.partner, iface2)
    @test ===(iface3.edge, edge3)
    @test ===(eq_iface3.edge, edge3)
    @test ===(iface3.partner, eq_iface3)
    @test ===(eq_iface3.partner, iface3)
end

@testset "@RV" begin
    g = FactorGraph()

    # @RV should construct a new variable
    x = constant(0.0)
    @RV y ~ GaussianMeanVariance(x, constant(1.0))
    @test length(g.variables) == 3 # including constants
    @test y.id == :y # automatically assigned id based on variable name in code
    @test haskey(g.variables, y.id)
    @test g.variables[y.id] === y

    # @RV ~ should reuse existing variable and handle keyword agruments
    y_old = y
    @RV y ~ GaussianMeanVariance(constant(0.0), constant(1.0); id=:g_node)
    @test length(g.variables) == 5 # including constants
    @test haskey(g.nodes, :g_node)
    @test y === y_old

    # @RV should handle array element assignments and explicit Variable ids
    g = FactorGraph()
    vars = Vector{Variable}(undef, 2)
    i = 1
    @RV [id=:v*i] vars[1] ~ GaussianMeanVariance(constant(0.0), constant(1.0); id=:tst1) # new Variable
    @test length(g.variables) == 3 # including constants
    @test vars[1].id == :v1
    @RV vars[1] ~ GaussianMeanVariance(constant(0.0), constant(1.0); id=:tst2) # existing Variable
    @test length(g.variables) == 5 # including constants
    @RV vars[2*i] ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @test vars[2*i].id == :vars_2
    varmatrix = Matrix{Variable}(undef,2,2)
    @RV varmatrix[1,2*i] ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @test varmatrix[1,2*i].id == :varmatrix_1_2
    vardict = Dict{Int,Variable}()
    @RV vardict[3*i+1] ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @test vardict[3*i+1].id == :vardict_4

    # @RV should work with '= syntax'
    g = FactorGraph()
    @RV x = constant(1.0) + constant(2.0)
    @test length(g.variables) == 3
    @test x.id == :x
    @RV [id=:my_y] y = x + constant(2.0)
    @test length(g.variables) == 5
    @test y.id == :my_y

    # @RV without node definition should create a new Variable
    g = FactorGraph()
    @RV x
    @test isa(x, Variable)
    @test g.variables[:x] === x
    @RV [id=:x_new] x
    @test g.variables[:x_new] === x
    @test length(g.variables) == 2

    # @RV without node definition should create a new Variable with a given index
    g = FactorGraph()
    xi = Vector{Variable}(undef, 1)
    t = 1
    @RV xi[t]
    @test isa(xi[1], Variable)
    @test g.variables[:xi_1] === xi[1]

    # @RV should throw an exception on incorrect usage
    @test_throws ErrorException @RV 1
    @test_throws ErrorException @RV [ 1 ]
    @test_throws ErrorException @RV [ 1 ] 1
    @test_throws ErrorException @RV [ 1 ] x
    @test_throws ErrorException @RV [ x ] x
    @test_throws ErrorException @RV [ x = :x_id ]
    @test_throws ErrorException @RV [ x = :x_id ] 1
end

end #module


text:
using FiniteDifferences: grad, jacobian, _jvp, jvp, j′vp, _j′vp, to_vec

@testset "grad" begin

    @testset "jvp(::$T)" for T in (Float64,)
        rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
        x, y = randn(rng, T, N), randn(rng, T, M)
        ẋ, ẏ = randn(rng, T, N), randn(rng, T, M)
        xy, ẋẏ = vcat(x, y), vcat(ẋ, ẏ)
        ż_manual = _jvp(fdm, (xy)->sum(sin, xy), xy, ẋẏ)[1]
        ż_auto = jvp(fdm, x->sum(sin, x[1]) + sum(sin, x[2]), ((x, y), (ẋ, ẏ)))
        ż_multi = jvp(fdm, (x, y)->sum(sin, x) + sum(sin, y), (x, ẋ), (y, ẏ))
        @test ż_manual ≈ ż_auto
        @test ż_manual ≈ ż_multi
    end

    @testset "grad(::$T)" for T in (Float64,)
        rng, fdm = MersenneTwister(123456), central_fdm(5, 1)
        x = randn(rng, T, 2)
        xc = copy(x)
        @test grad(fdm, x->sin(x[1]) + cos(x[2]), x)[1] ≈ [cos(x[1]), -sin(x[2])]
        @test xc == x
    end

    function check_jac_and_jvp_and_j′vp(fdm, f, ȳ::Array, x::Array, ẋ::Array, J_exact)
        xc = copy(x)

        # Validate inputs.
        @assert length(x) == length(ẋ)
        @assert length(ȳ) == length(f(x))

        # Check that the jacobian is as expected.
        J_fdm = jacobian(fdm, f, x)[1]
        @test size(J_fdm) == (length(ȳ), length(x))
        @test J_fdm ≈ J_exact
        @test J_fdm == jacobian(fdm, f, x)[1]

        # Check that the estimated jvp and j′vp are consistent with their definitions.
        @test _jvp(fdm, f, x, ẋ) ≈ J_exact * ẋ
        @test _j′vp(fdm, f, ȳ, x) ≈ transpose(J_exact) * ȳ

        # Check that no mutation occured that wasn't reverted.
        @test xc == x
    end

    @testset "jacobian / _jvp / _j′vp (::$T)" for T in (Float64,)
        rng, P, Q, fdm = MersenneTwister(123456), 3, 2, central_fdm(5, 1)
        ȳ = randn(rng, T, P)
        A = randn(rng, T, P, Q)
        x = randn(rng, T, Q)
        ẋ = randn(rng, T, Q)
        Ac = copy(A)

        check_jac_and_jvp_and_j′vp(fdm, x->A * x, ȳ, x, ẋ, A)
        @test Ac == A
        check_jac_and_jvp_and_j′vp(fdm, x->sin.(A * x), ȳ, x, ẋ, cos.(A * x) .* A)
        @test Ac == A

        # Prevent regression against https://github.com/JuliaDiff/FiniteDifferences.jl/issues/67
        J = first(jacobian(fdm, identity, x))
        @test J ≈ one(Matrix{T}(undef, size(J)))
    end

    @testset "multi vars jacobian/grad" begin
        rng, fdm = MersenneTwister(123456), central_fdm(5, 1)

        f1(x, y) = x * y + x
        f2(x, y) = sum(x * y + x)
        f3(x::Tuple) = sum(x[1]) + x[2]
        f4(d::Dict) = sum(d[:x]) + d[:y]

        @testset "jacobian" begin
            @testset "check multiple matrices" begin
                x, y = rand(rng, 3, 3), rand(rng, 3, 3)
                jac_xs = jacobian(fdm, f1, x, y)
                @test jac_xs[1] ≈ jacobian(fdm, x->f1(x, y), x)[1]
                @test jac_xs[2] ≈ jacobian(fdm, y->f1(x, y), y)[1]
            end

            @testset "check mixed scalar and matrices" begin
                x, y = rand(3, 3), 2.0
                jac_xs = jacobian(fdm, f1, x, y)
                @test jac_xs[1] ≈ jacobian(fdm, x->f1(x, y), x)[1]
                @test jac_xs[2] ≈ jacobian(fdm, y->f1(x, y), y)[1]
            end
        end

        @testset "grad" begin
            @testset "check multiple matrices" begin
                x, y = rand(rng, 3, 3), rand(rng, 3, 3)
                dxs = grad(fdm, f2, x, y)
                @test dxs[1] ≈ grad(fdm, x->f2(x, y), x)[1]
                @test dxs[2] ≈ grad(fdm, y->f2(x, y), y)[1]
            end

            @testset "check mixed scalar & matrices" begin
                x, y = rand(rng, 3, 3), 2.0
                dxs = grad(fdm, f2, x, y)
                @test dxs[1] ≈ grad(fdm, x->f2(x, y), x)[1]
                @test dxs[2] ≈ grad(fdm, y->f2(x, y), y)[1]
            end

            @testset "check tuple" begin
                x, y = rand(rng, 3, 3), 2.0
                dxs = grad(fdm, f3, (x, y))[1]
                @test dxs[1] ≈ grad(fdm, x->f3((x, y)), x)[1]
                @test dxs[2] ≈ grad(fdm, y->f3((x, y)), y)[1]
            end

            @testset "check dict" begin
                x, y = rand(rng, 3, 3), 2.0
                d = Dict(:x=>x, :y=>y)
                dxs = grad(fdm, f4, d)[1]
                @test dxs[:x] ≈ grad(fdm, x->f3((x, y)), x)[1]
                @test dxs[:y] ≈ grad(fdm, y->f3((x, y)), y)[1]
            end
        end
    end

    @testset "j′vp(::$T)" for T in (Float64,)
        rng, N, M, fdm = MersenneTwister(123456), 2, 3, central_fdm(5, 1)
        x, y = randn(rng, T, N), randn(rng, T, M)
        z̄ = randn(rng, T, N + M)
        xy = vcat(x, y)
        x̄ȳ_manual = j′vp(fdm, xy->sin.(xy), z̄, xy)[1]
        x̄ȳ_auto = j′vp(fdm, x->sin.(vcat(x[1], x[2])), z̄, (x, y))[1]
        x̄ȳ_multi = j′vp(fdm, (x, y)->sin.(vcat(x, y)), z̄, x, y)
        @test x̄ȳ_manual ≈ vcat(x̄ȳ_auto...)
        @test x̄ȳ_manual ≈ vcat(x̄ȳ_multi...)
    end

    # Tests for complex numbers, to prevent regressions against
    # https://github.com/JuliaDiff/FiniteDifferences.jl/pull/76 and
    # https://github.com/JuliaDiff/FiniteDifferences.jl/issues/75
    @testset "Complex correctness - $T" for T in (ComplexF64,)
        rng = MersenneTwister(123456)
        x = randn(rng, T)
        y = randn(rng, T)
        fdm = FiniteDifferences.Central(5, 1)

        # Addition.
        dx, dy = FiniteDifferences.jacobian(fdm, +, x, y)
        @test dx ≈ [1 0; 0 1]
        @test dy ≈ [1 0; 0 1]

        # Negation.
        dx, = FiniteDifferences.jacobian(fdm, -, x)
        @test dx ≈ [-1 0; 0 -1]

        # Multiplication.
        dx, dy = FiniteDifferences.jacobian(fdm, *, x, y)
        @test dx ≈ [real(y) -imag(y); imag(y) real(y)]
        @test dy ≈ [real(x) -imag(x); imag(x) real(x)]

        # Magnitude
        dx, = FiniteDifferences.grad(fdm, abs2, x)
        @test real(dx) ≈ 2 * real(x)
        @test imag(dx) ≈ 2 * imag(x)

        # Concatenate tangents / cotangents.
        ẋ = randn(rng, T)
        ẏ = randn(rng, T)
        ẋẏ = vcat(first(to_vec(ẋ)), first(to_vec(ẏ)))
        z̄ = randn(rng, T)
        z̄_vec = first(to_vec(z̄))

        # Check jvp against jacobian result.
        f(x, y) = 5 * sin(x) * cos(y)
        Jx, Jy = FiniteDifferences.jacobian(fdm, f, x, y)
        ż = jvp(fdm, f, (x, ẋ), (y, ẏ))
        ż_vec = [real(ż), imag(ż)]
        ż_manual = hcat(Jx, Jy) * ẋẏ
        @test ż_vec ≈ ż_manual

        # Check j′vp against jacobian result.
        x̄, ȳ = j′vp(fdm, f, z̄, x, y)
        @test [real(x̄), imag(x̄)] ≈ Jx'z̄_vec
        @test [real(ȳ), imag(ȳ)] ≈ Jy'z̄_vec
    end
end


text:
# Model
include_model("hopper")

# modify parameters to match implicit dynamics example
gravity = 9.81 # gravity
μ_world = 0.8 # coefficient of friction
μ_joint = 0.0

mb = 3.0 # body mass
ml = 0.3  # leg mass
Jb = 0.75 # body inertia
Jl = 0.075 # leg inertia

model = Hopper{Discrete, FixedTime}(n, m, d,
			   mb, ml, Jb, Jl,
			   μ_world, gravity,
			   qL, qU,
			   uL, uU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)

# Horizon
T = 5
# Time step
h = 0.05
tf = h * (T - 1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 1000.0 #model_ft.uU
_ul = zeros(model.m)
_ul[model.idx_u] .= -1000.0 #model_ft.uL
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
q1 = [0.0; 0.5; 0.0; 0.5]
# qM = [0.5; 0.5; 0.0; 0.5]
qT = q1#[1.0; 0.5; 0.0; 0.5]
q_ref = q1#[0.5; 0.75; 0.0; 0.25]
x_ref = [q_ref; q_ref]

xl, xu = state_bounds(model, T, [q1; q1], [q1; q1])

# Objective

# gait 1
obj_contact_penalty = PenaltyObjective(1.0e5, model.m)

obj = MultiObjective([obj_contact_penalty])

# Constraints
include_constraints(["contact"])
con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = [[q1; q1] for t = 1:T]
u0 = [[0.0; model.g * (model.mb + model.ml) * 0.5 * h; 0.0 * rand(model.m - model.nu)] for t = 1:T-1] # random controls

fd(model, x0[1], x0[1], u0[1], zeros(model.d), h, 1)

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
# include_snopt()

@time z̄, info = solve(prob, copy(z0),
	nlp = :ipopt,
	tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

x̄, ū = unpack(z̄, prob)
u_stand = ū

@show check_slack(z̄, prob)

using Plots
t = range(0, stop = h * (T - 1), length = T)
plot(t, hcat(ū..., ū[end])[1:2,:]', linetype=:steppost,
	xlabel="time (s)", ylabel = "control",
	label = ["angle" "length"],
	width = 2.0, legend = :top)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model,
    state_to_configuration(x̄),
	# state_to_configuration([[x̄[1] for i = 1:50]...,x̄..., [x̄[end] for i = 1:50]...]),
	Δt = h,
	scenario = :flip)

@save joinpath(pwd(), "examples/implicit_dynamics/examples/comparisons/hopper_stand.jld2") u_stand

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Base.Iterators: product

# --- Nodes and States ---

"""
    Node = Int

Primitive type for node. Alias for `Int`.
"""
const Node = Int

"""
    abstract type AbstractNode end

Node type for directed, acyclic graph.
"""
abstract type AbstractNode end

function validate_node(j::Node, I_j::Vector{Node})
    if !allunique(I_j)
        throw(DomainError("All information nodes should be unique."))
    end
    if !all(i < j for i in I_j)
        throw(DomainError("All nodes in the information set must be less than node j."))
    end
end

"""
    ChanceNode <: AbstractNode

Chance node type.

# Examples
```julia
c = ChanceNode(3, [1, 2])
```
"""
struct ChanceNode <: AbstractNode
    j::Node
    I_j::Vector{Node}
    function ChanceNode(j::Node, I_j::Vector{Node})
        validate_node(j, I_j)
        new(j, I_j)
    end
end

"""
    DecisionNode <: AbstractNode

Decision node type.

# Examples
```julia
d = DecisionNode(2, [1])
```
"""
struct DecisionNode <: AbstractNode
    j::Node
    I_j::Vector{Node}
    function DecisionNode(j::Node, I_j::Vector{Node})
        validate_node(j, I_j)
        new(j, I_j)
    end
end

"""
    ValueNode <: AbstractNode

Value node type.

# Examples
```julia
v = ValueNode(4, [1, 3])
```
"""
struct ValueNode <: AbstractNode
    j::Node
    I_j::Vector{Node}
    function ValueNode(j::Node, I_j::Vector{Node})
        validate_node(j, I_j)
        new(j, I_j)
    end
end

"""
    const State = Int

Primitive type for the number of states. Alias for `Int`.
"""
const State = Int

"""
    States <: AbstractArray{State, 1}

States type. Works like `Vector{State}`.

# Examples
```julia
S = States([2, 3, 2, 4])
```
"""
struct States <: AbstractArray{State, 1}
    vals::Vector{State}
    function States(vals::Vector{State})
        if !all(vals .≥ 1)
            throw(DomainError("All states must be ≥ 1."))
        end
        new(vals)
    end
end

Base.size(S::States) = size(S.vals)
Base.IndexStyle(::Type{<:States}) = IndexLinear()
Base.getindex(S::States, i::Int) = getindex(S.vals, i)
Base.length(S::States) = length(S.vals)
Base.eltype(S::States) = eltype(S.vals)

"""
    function States(states::Vector{Tuple{State, Vector{Node}}})

Construct states from vector of (state, nodes) tuples.

# Examples
```julia-repl
julia> S = States([(2, [1, 3]), (3, [2, 4, 5])])
States([2, 3, 2, 3, 3])
```
"""
function States(states::Vector{Tuple{State, Vector{Node}}})
    S_j = Vector{State}(undef, sum(length(j) for (_, j) in states))
    for (s, j) in states
        S_j[j] .= s
    end
    States(S_j)
end

"""
    function validate_influence_diagram(S::States, C::Vector{ChanceNode}, D::Vector{DecisionNode}, V::Vector{ValueNode})

Validate influence diagram.
"""
function validate_influence_diagram(S::States, C::Vector{ChanceNode}, D::Vector{DecisionNode}, V::Vector{ValueNode})
    n = length(C) + length(D)
    if length(S) != n
        throw(DomainError("Each change and decision node should have states."))
    end
    if Set(c.j for c in C) ∪ Set(d.j for d in D) != Set(1:n)
        throw(DomainError("Union of change and decision nodes should be {1,...,n}."))
    end
    if Set(v.j for v in V) != Set((n+1):(n+length(V)))
        throw(DomainError("Values nodes should be {n+1,...,n+|V|}."))
    end
    I_V = union((v.I_j for v in V)...)
    if !(I_V ⊆ Set(1:n))
        throw(DomainError("Each information set I(v) for value node v should be a subset of C∪D."))
    end
    # Check for redundant nodes.
    leaf_nodes = setdiff(1:n, (c.I_j for c in C)..., (d.I_j for d in D)...)
    for i in leaf_nodes
        if !(i∈I_V)
            @warn("Chance or decision node $i is redundant.")
        end
    end
    for v in V
        if isempty(v.I_j)
            @warn("Value node $(v.j) is redundant.")
        end
    end
end

# --- Paths ---

"""
    const Path{N} = NTuple{N, State} where N

Path type. Alias for `NTuple{N, State} where N`.
"""
const Path{N} = NTuple{N, State} where N

"""
    function paths(states::AbstractVector{State})

Iterate over paths in lexicographical order.

# Examples
```julia-repl
julia> states = States([2, 3])
julia> vec(collect(paths(states)))
[(1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3)]
```
"""
function paths(states::AbstractVector{State})
    product(UnitRange.(one(eltype(states)), states)...)
end

"""
    function paths(states::AbstractVector{State}, fixed::Dict{Node, State})

Iterate over paths with fixed states in lexicographical order.

# Examples
```julia-repl
julia> states = States([2, 3])
julia> vec(collect(paths(states, Dict(1=>2))))
[(2, 1), (2, 2), (2, 3)]
```
"""
function paths(states::AbstractVector{State}, fixed::Dict{Node, State})
    iters = collect(UnitRange.(one(eltype(states)), states))
    for (i, v) in fixed
        iters[i] = UnitRange(v, v)
    end
    product(iters...)
end

"""
    const ForbiddenPath = Tuple{Vector{Node}, Set{Path}}

ForbiddenPath type.

# Examples
```julia
ForbiddenPath[
    ([1, 2], Set([(1, 2)])),
    ([3, 4, 5], Set([(1, 2, 3), (3, 4, 5)]))
]
```
"""
const ForbiddenPath = Tuple{Vector{Node}, Set{Path}}

# --- Probabilities ---

"""
    struct Probabilities{N} <: AbstractArray{Float64, N}

Construct and validate stage probabilities.

# Examples
```julia-repl
julia> data = [0.5 0.5 ; 0.2 0.8]
julia> X = Probabilities(2, data)
julia> s = (1, 2)
julia> X(s)
0.5
```
"""
struct Probabilities{N} <: AbstractArray{Float64, N}
    j::Node
    data::Array{Float64, N}
    function Probabilities(j::Node, data::Array{Float64, N}) where N
        for i in CartesianIndices(size(data)[1:end-1])
            if !(sum(data[i, :]) ≈ 1)
                throw(DomainError("Probabilities should sum to one."))
            end
        end
        new{N}(j, data)
    end
end

Base.size(P::Probabilities) = size(P.data)
Base.IndexStyle(::Type{<:Probabilities}) = IndexLinear()
Base.getindex(P::Probabilities, i::Int) = getindex(P.data, i)
Base.getindex(P::Probabilities, I::Vararg{Int,N}) where N = getindex(P.data, I...)

(X::Probabilities)(s::Path) = X[s...]

# --- Path Probability ---

"""
    abstract type AbstractPathProbability end

Abstract path probability type.

# Examples
```julia
struct PathProbability <: AbstractPathProbability
    C::Vector{ChanceNode}
    # ...
end

(U::PathProbability)(s::Path) = ...
```
"""
abstract type AbstractPathProbability end

"""
    DefaultPathProbability <: AbstractPathProbability

Path probability.

# Examples
```julia
P = DefaultPathProbability(C, X)
s = (1, 2)
P(s)
```
"""
struct DefaultPathProbability <: AbstractPathProbability
    C::Vector{ChanceNode}
    X::Vector{Probabilities}
end

function (P::DefaultPathProbability)(s::Path)
    prod(X(s[[c.I_j; c.j]]) for (c, X) in zip(P.C, P.X))
end

# --- Consequences ---

"""
    Consequences{N} <: AbstractArray{Float64, N}

State utilities.

# Examples
```julia-repl
julia> vals = [1.0 -2.0; 3.0 4.0]
julia> Y = Consequences(3, vals)
julia> s = (1, 2)
julia> Y(s)
-2.0
```
"""
struct Consequences{N} <: AbstractArray{Float64, N}
    j::Node
    data::Array{Float64, N}
end

Base.size(Y::Consequences) = size(Y.data)
Base.IndexStyle(::Type{<:Consequences}) = IndexLinear()
Base.getindex(Y::Consequences, i::Int) = getindex(Y.data, i)
Base.getindex(Y::Consequences, I::Vararg{Int,N}) where N = getindex(Y.data, I...)

(Y::Consequences)(s::Path) = Y[s...]

# --- Path Utility ---

"""
    abstract type AbstractPathUtility end

Abstract path utility type.

# Examples
```julia
struct PathUtility <: AbstractPathUtility
    V::Vector{ValueNode}
    # ...
end

(U::PathUtility)(s::Path) = ...
```
"""
abstract type AbstractPathUtility end

"""
    DefaultPathUtility <: AbstractPathUtility

Default path utility.

# Examples
```julia
U = DefaultPathUtility(V, Y)
s = (1, 2)
U(s)
```
"""
struct DefaultPathUtility <: AbstractPathUtility
    V::Vector{ValueNode}
    Y::Vector{Consequences}
end

function (U::DefaultPathUtility)(s::Path)
    sum(Y(s[v.I_j]) for (v, Y) in zip(U.V, U.Y))
end

# --- Local Decision Strategy ---

"""
    LocalDecisionStrategy{N} <: AbstractArray{Int, N}

Local decision strategy type.

# Examples
```julia
Z = LocalDecisionStrategy(1, data)
Z(s_I)
```
"""
struct LocalDecisionStrategy{N} <: AbstractArray{Int, N}
    j::Node
    data::Array{Int, N}
    function LocalDecisionStrategy(j::Node, data::Array{Int, N}) where N
        if !all(0 ≤ x ≤ 1 for x in data)
            throw(DomainError("All values x must be 0 ≤ x ≤ 1."))
        end
        for s_I in CartesianIndices(size(data)[1:end-1])
            if !(sum(data[s_I, :]) == 1)
                throw(DomainError("Values should add to one."))
            end
        end
        new{N}(j, data)
    end
end

Base.size(Z::LocalDecisionStrategy) = size(Z.data)
Base.IndexStyle(::Type{<:LocalDecisionStrategy}) = IndexLinear()
Base.getindex(Z::LocalDecisionStrategy, i::Int) = getindex(Z.data, i)
Base.getindex(Z::LocalDecisionStrategy, I::Vararg{Int,N}) where N = getindex(Z.data, I...)

function (Z::LocalDecisionStrategy)(s_I::Path)::State
    findmax(Z[s_I..., :])[2]
end

# --- Decision Strategy ---

"""
    DecisionStrategy

Decision strategy type.
"""
struct DecisionStrategy
    D::Vector{DecisionNode}
    Z_j::Vector{LocalDecisionStrategy}
end


text:
include("testsystem1.jl")
srand(11)
Y, X = sample(20, 100, M)
kalmanrts(Y, M)


text:
# This file is a part of EmpiricalDistributions.jl, licensed under the MIT License (MIT).

"""
    UvBinnedDist{T<:Real} <: StatsBase.Distribution{Univariate,Continuous}

A binned univariate distribution, usually derived from a histogram.

Constructors:

```julia
    UvBinnedDist(h::StatsBase.Histogram{<:Real,1})
    UvBinnedDist{T<:Real}(h::StatsBase.Histogram{<:Real,1})
```

You can convert a `UvBinnedDist` back to a histogram via

```julia
convert(StatsBase.Histogram, dist::UvBinnedDist)
StatsBase.Histogram(dist::UvBinnedDist)
```
"""
struct UvBinnedDist{
    T <: Real,
    U <: Real,
    VT <: AbstractVector{T},
    VU <: AbstractVector{U}
} <: Distribution{Univariate,Continuous}
    _edge::VT
    _edge_cdf::VU
    _bin_pdf::VU
    _bin_probmass::VU
    _closed_left::Bool
    _mean::T
    _mode::T
    _var::T
end

export UvBinnedDist

function UvBinnedDist{T}(h::Histogram{<:Real,1}) where {T<:Real}
    nh = normalize(h)

    edge = T.(first(nh.edges))
    closed_left = nh.closed == :left

    bin_pdf = nh.weights
    bin_probmass = bin_pdf .* _bin_widths(edge)
    bin_probmass .*= inv(sum(bin_probmass))

    edge_cdf = pushfirst!(cumsum(bin_probmass), 0)
    @assert last(edge_cdf) ≈ 1
    edge_cdf[end] = 1

    mean_est_tpl = _mean(nh)
    mean_est = T(first(mean_est_tpl))
    mode_est = T(first(_mode(nh)))
    var_est = T(first(_var(nh, mean_est_tpl)))

    return UvBinnedDist(
        edge, edge_cdf, bin_pdf, bin_probmass, closed_left,
        mean_est, mode_est, var_est
    )
end

UvBinnedDist(h::Histogram{<:Real,1}) = UvBinnedDist{float(eltype(first(h.edges)))}(h)

function Adapt.adapt_structure(to, d::UvBinnedDist)
    UvBinnedDist(
        adapt(to, d._edge), adapt(to, d._edge_cdf),
        adapt(to, d._bin_pdf), adapt(to, d._bin_probmass),
        d._closed_left, d._mean, d._mode, d._var
    )
end

Histogram(d::UvBinnedDist) = Histogram((Array(d._edge),), Array(d._bin_pdf), (d._closed_left ? :left : :right), true)
Base.convert(::Type{Histogram}, d::UvBinnedDist) = Histogram(d)

Base.length(d::UvBinnedDist) = 1
Base.size(d::UvBinnedDist) = ()
Base.eltype(d::UvBinnedDist{T}) where {T} = T

Statistics.mean(d::UvBinnedDist) = d._mean
StatsBase.mode(d::UvBinnedDist) = d._mode
Statistics.var(d::UvBinnedDist) = d._var

Distributions.minimum(d::UvBinnedDist) = first(d._edge)
Distributions.maximum(d::UvBinnedDist) = last(d._edge)

function Distributions.pdf(d::UvBinnedDist{T,U}, x::Real) where {T,U}
    if insupport(d, x)
        i = _find_bin(d._edge, d._closed_left, x)
        i <= lastindex(d._bin_pdf) ? convert(U, d._bin_pdf[i]) : zero(U)
    else
        zero(U)
    end
end

Distributions.logpdf(d::UvBinnedDist, x::Real) = log(pdf(d, x))

function Distributions.cdf(d::UvBinnedDist{T}, x::Real) where T
    if x <= minimum(d)
        zero(T)
    elseif x >= maximum(d)
        one(T)
    else
        _linear_interpol(d._edge, d._edge_cdf, x)
    end
end

function Distributions.quantile(d::UvBinnedDist, x::Real)
    0 <= x <= 1 || throw(DomainError(x, "quantile requires value between 0 and 1"))
    _linear_interpol(d._edge_cdf, d._edge, x)
end

function Random.rand(rng::AbstractRNG, d::UvBinnedDist{T}) where T
    u = rand(T)
    @assert axes(d._edge) == axes(d._edge_cdf)
    i_lo, i_hi = _find_idxs_lohi(d._edge_cdf, u)
    x_lo, x_hi = d._edge[i_lo], d._edge[i_hi]
    _rand_uniform(rng, T, x_lo, x_hi)
end


text:
module CUDAnative

using CUDAdrv

using LLVM
using LLVM.Interop

using Adapt

using Pkg
using Libdl

const ext = joinpath(@__DIR__, "..", "deps", "ext.jl")
isfile(ext) || error("CUDAnative.jl has not been built, please run Pkg.build(\"CUDAnative\").")
include(ext)
if !configured
    # default (non-functional) values for critical variables,
    # making it possible to _load_ the package at all times.
    const target_support = [v"2.0"]
    const cuda_driver_version = v"5.5"
end

include("utils.jl")

# needs to be loaded _before_ the compiler infrastructure, because of generated functions
include(joinpath("device", "tools.jl"))
include(joinpath("device", "pointer.jl"))
include(joinpath("device", "array.jl"))
include(joinpath("device", "intrinsics.jl"))
include(joinpath("device", "libdevice.jl"))
include(joinpath("device", "runtime.jl"))

include("compiler.jl")
include("execution.jl")
include("reflection.jl")

include("deprecated.jl")

include("init.jl")

end


text:
module SkyRoom

export main

using Dates, WGLMakie, JSServe, ImageCore, FilePathsBase, CSV, Pkg.TOML, Tar, FileIO, ImageMagick, Observables, PyCall, LibSerialPort, Missings#, MemoryHunter
using FilePathsBase: /
using JSServe.DOM
using JSServe: @js_str

# function print_sizes()
#     mem = Pair{String, Int}[]
#     for m in values(Base.loaded_modules), vs in names(m, all = true)
#         if isdefined(m, vs)
#             v = getfield(m, vs)
#             x = Base.summarysize(v)
#             if x > 10^6
#                 push!(mem, string(m, "; ", vs, ": ") => x)
#             end
#         end
#     end
#     sort!(mem, by = last, rev = true)
#     n = length(mem)
#     if n > 30
#         deleteat!(mem, 31:n)
#     end
#     println("Memory usage:")
#     for (txt, i) in mem
#         println(txt, Base.format_bytes(i))
#     end
#     println("")
# end

const t₀ = Ref{DateTime}()
const datadir = Ref{PosixPath}()
const tempdatadir = Ref{PosixPath}()
const nicolas = Ref{Bool}()
const baudrate = Ref{Int}()
const top_rpm = Ref{Int}()

function __init__()
    t₀[] = now()
    datadir[] = p"/home/pi/mnt/data"
    isdir(datadir[]) || mkpath(datadir[])
    tempdatadir[] = p"/home/pi/mnt/tmp"
    isdir(tempdatadir[]) || mkpath(tempdatadir[])
    nicolas[] = Base.Libc.gethostname() == "nicolas"
    baudrate[] = 9600
    top_rpm[] = 12650
    py"""
    import picamera
    import io

    def PC():
        return picamera.PiCamera()

    def PIO():
        return io.BytesIO()
    """

end

include("cobs.jl")
include("abstractarduinos.jl")
include("leds.jl")
# if nicolas[]
include("winds.jl")
# end

get_rpms(::Nothing) = nothing
record(::Nothing, _) = nothing
Base.isopen(::Nothing) = true
Base.close(::Nothing) = nothing

include("camera.jl")

function initialize()
    fan_ports = ["/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_957353530323510141D0-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_95635333930351917172-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_95735353032351010260-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_55838323435351213041-if00", "/dev/serial/by-id/usb-Arduino__www.arduino.cc__0043_957353530323514121D0-if00"]
    allwind = nicolas[] ? AllWind([FanArduino(id, port) for (id, port) in enumerate(fan_ports) if isconnected(port)], 1)  : nothing
    led_arduino = LEDArduino()
    framerate = 30
    camera = PiCamera(framerate, 67, 67, 4)

    data = Observable((frame = snap(camera), trpms = get_rpms(allwind)))
    task = @async while isopen(allwind) && isopen(camera)
        # try
            data[] = (; frame = snap(camera), trpms = get_rpms(allwind))
            sleep(0.01)
        #=catch e
            @show now()
            MemoryHunter.print_sizes()
            @show e
            throw(e)
        end=#
    end
    return allwind, led_arduino, camera, data, task
end

function plotrpm(trpms)
    if isnothing(trpms[])
        return nothing
    end
    rpms = map(trpms) do (_, x)
        collect(Missings.replace(Iterators.flatten(x), NaN))
    end
    nports = length(last(trpms[]))
    rpmplt_cont = (colors = repeat(1:nports, inner = [3]), x = vcat(((i - 1)*(nports - 1) + 1 : (i - 1)*(nports - 1) + 3 for i in 1:nports)...), y = top_rpm[]*ones(3nports), resolution = (540, round(Int, 3nports + 3*540/(3nports+4))))
    rpmplot = Scene(show_axis = false, resolution = rpmplt_cont.resolution)
    barplot!(rpmplot, rpmplt_cont.x, rpmplt_cont.y, color = :white, strokecolor = :black, strokewidth = 1)
    barplot!(rpmplot, rpmplt_cont.x, rpms, color = rpmplt_cont.colors, strokecolor = :transparent, strokewidth = 0)
    disconnect!(AbstractPlotting.camera(rpmplot))
    return rpmplot
end

parse2wind(windrow) = [Wind(id, v) for (id, v) in enumerate(windrow)]
parse2stars(starsrow) = [Star(v...) for v in Iterators.partition(starsrow, 4) if !any(ismissing, v)]

function parse2one(setup_file)
    rs = NamedTuple{(:label, :stars),Tuple{String,Array{Star,1}}}[]
    for nt in CSV.Rows(setup_file, header = 1:2)
        if !ismissing(nt.setup_label)
            t = Tuple(nt)
            stars = parse2stars(t[2:end])
            push!(rs, (label = t[1], stars = stars))
        end
    end
    rs
end

function parse2both(setup_file)
    rs = NamedTuple{(:label, :fans, :stars),Tuple{String,Array{Wind,1},Array{Star,1}}}[]
    for nt in CSV.Rows(setup_file, header = 1:2)
        if !ismissing(nt.setup_label)
            t = Tuple(nt)
            fans = parse2wind(t[2:6])
            stars = parse2stars(t[7:end])
            push!(rs, (label = t[1], fans = fans, stars = stars))
        end
    end
    rs
end

function get_setups()
    setupsurl = nicolas[] ? "https://docs.google.com/spreadsheets/d/e/2PACX-1vQNLWhLfp_iuW68j7SM6Px8ysTmbrfmrP_7ipXK9BkfzBgfqn3Mj7ra177mZyHlY5NLA3SDtfYNTROv/pub?gid=0&single=true&output=csv" : "https://docs.google.com/spreadsheets/d/e/2PACX-1vSfv92ymTJjwdU-ft9dgglOOnxPVWwtk6gFIVSocHM3jSfHkjYk-mtEXl3g96-735Atbk1LBRt-8lAY/pub?gid=0&single=true&output=csv"
    setup_file = download(setupsurl)
    nicolas[] ? parse2both(setup_file) : parse2one(setup_file)
end

_fun(::Nothing, _) = nothing
_fun(allwind, setup) = for a in allwind.arduinos
    a.pwm[] = setup.fans[a.id].pwm
end

function button(allwind, led_arduino, setup, setuplog)
    button_class = "bg-blue-500 hover:bg-blue-700 text-white font-bold py-2 px-4 rounded"
    b = JSServe.Button(string(setup.label), class = button_class)
    on(b) do _
        _fun(allwind, setup)
        led_arduino.pwm[] = parse2arduino(setup.stars)
        push!(setuplog, now() => todict(setup))
    end
    return b
end

typedict(x::T) where {T} = Dict(fn=>getfield(x, fn) for fn ∈ fieldnames(T))

function todict(setup)
    x = Dict{Symbol, Any}(pairs(setup))
    if haskey(x, :fans)
        x[:fans] = typedict.(x[:fans])
    end
    x[:stars] = typedict.(x[:stars])
    return x
end

stop_record(::Nothing) = nothing
stop_record(allwind) = close(allwind.io)

function record(tf, allwind, camera, timestamp, setuplog)
    if tf
        timestamp[] = string(now())
        folder = datadir[] / timestamp[]
        !isdir(folder) && mkdir(folder)
        camera.cam.recording && camera.cam.stop_recording()
        camera.cam.start_recording(string(folder / "video.h264"))
        record(allwind, folder)
        deleteat!(setuplog, 1:length(setuplog) - 1)
    else
        camera.cam.stop_recording()
        stop_record(allwind)
    end
end

function save(timestamp, beetleid, comment, setuplog, left2backup, msg)
    folder = datadir[] / timestamp[]
    if !isdir(folder)
        msg[] = "You haven't recorded a video yet, there is nothing to save"
        return nothing
    end
    msg[] = "Saving"
    md = Dict()
    md["timestamp"] = timestamp[]
    md["beetleid"] = beetleid[]
    md["comment"] = comment[]
    md["setuplog"] = Dict(string(t) => v for (t,v) in setuplog)
    open(folder / "metadata.toml", "w") do io
        TOML.print(io, md)
    end
    left2backup[] += 1
    comment[] = ""
    beetleid[] = ""
    timestamp[] = "garbage"
    msg[] = "Saved"
    return nothing
end

function backup(left2backup, msg)
    msg[] = "Started backing up..."
    for folder in readpath(datadir[])
        name = basename(folder)
        tmp = tempdatadir[] / name * ".tar"
        Tar.create(string(folder), string(tmp))
        rm(folder, recursive = true)
        bucket = nicolas[] ? "nicolas-cage-skyroom" : "top-floor-skyroom2"
        run(`aws s3 mv $tmp s3://$bucket/$name.tar --quiet`)
        left2backup[] = length(readdir(datadir[]))
    end
    msg[] = "Finished backing up!"
end

function handler(allwind, led_arduino, camera, data, session, request)

    data2 = Observable(data[])
    on(session, data) do x
        data2[] = x
    end

    msg = Observable("")
    on(msg) do x
        if !isempty(x)
            @async begin
                sleep(3)
                msg[] = ""
            end
        end
    end

    frame = map(data2) do x
        x.frame
    end
    trpms = map(data2) do x
        x.trpms
    end

    rpmplot = plotrpm(trpms)
    frameplot = image(frame, scale_plot = false, show_axis = false)
    disconnect!(WGLMakie.camera(frameplot))

    setuplog = []
    timestamp = Ref("")
    recording = JSServe.Checkbox(false)
    on(x -> record(x, allwind, camera, timestamp, setuplog), recording)

    text_class = "border py-2 px-3 text-grey-darkest"
    comment = JSServe.TextField("", class = text_class)
    beetleid = JSServe.TextField("", class = text_class)

    button_class = "bg-blue-500 hover:bg-blue-700 text-white font-bold py-2 px-4 rounded"
    saving = JSServe.Button("Save", class = button_class)
    on(saving) do _
        if recording[]
            recording[] = false
            sleep(0.1)
        end
        save(timestamp, beetleid, comment, setuplog, left2backup, msg)
    end

    backingup = JSServe.Button("Backup", class = button_class)
    left2backup = Observable(length(readdir(datadir[])))
    on(_ -> backup(left2backup, msg), backingup)

    setups = get_setups()
    buttons = [button(allwind, led_arduino, setup, setuplog) for setup in setups]

    # print_sizes()

    grid_class = "grid auto-cols-max grid-cols-2 md:grid-cols-3 lg:grid-cols-5 gap-4"
    return DOM.div(JSServe.TailwindCSS,
                   DOM.div(rpmplot),
                   DOM.div(frameplot),
                   DOM.div(buttons..., class = grid_class),
                   DOM.div("Record ", recording),
                   DOM.div("Beetle ID ", beetleid),
    DOM.div("Comment ", comment),
    DOM.div(saving),
    DOM.div(backingup, left2backup, " runs left to backup"),
    DOM.div(msg),
    class = "grid grid-cols-1 gap-4"
   )
end

function main()
    allwind, led_arduino, camera, data, task = initialize()
    JSServe.Server((session, request) -> handler(allwind, led_arduino, camera, data, session, request), "0.0.0.0", 8082)
end

end # module

# PyCall, Dates, WGLMakie, JSServe, ImageCore, FilePathsBase, CSV, DataFrames, HTTP, Pkg.TOML, Tar, FileIO, ImageMagick, LibSerialPort, Observables


text:
using Polynomials

function _L(i, ord)
    ord % 2 == 1 || throw(DomainError("ord=$ord must be odd"))
    d = div(ord,2)
    result = Polynomials.Polynomial([big(1 // 1)])
    for j = -d:d+1
        if j != i
            result *= Polynomials.Polynomial([-j, 1 // 1] .// (i - j))
        end
    end
    return result
end

function _Lprim(i, ord)
    result = big(0//1)
    d = div(ord,2)
    for j = -d:d+1
        if i != j
            result += 1//(i-j)
        end
    end
    return result
end
_K(i,ord)=_L(i,ord)^2*Polynomial([-i,1//1])
_H(i,ord)=_L(i,ord)^2*(1-2*_Lprim(i,ord)*Polynomial([-i,1//1]))
# function bplusbis(i,ord)
#     i != 0 || throw(DomainError("i=$i must be different of zero"))
#     result=big(1//1)
#     for j=-d:d+1
#         if i != j
#             if j != 0
#                 result *= (-j)
#             end
#             result /= (i-j)
#         end
#     end
#     return result
# end
# ref function _bplus(i,ord)
#     d = div(ord,2)
#     res= prod([big(-j//1) for j=-d:d+1 if j!=0 && j!=i])
#     return res * prod([big(1//(i-j)) for j=-d:d+1 if j!=i])
# end
function _bplus(i, rplus, splus)
    res= prod([big(-j//1) for j=rplus:splus if j!=0 && j!=i])
    return res * prod([big(1//(i-j)) for j=rplus:splus if j!=i])
end
function _bminus(i, rminus, sminus)
    return -_bplus(-i, -sminus, -rminus)
end
struct PrecalHermite{ord,d}
    L::Vector{Polynomial{Rational{BigInt}}}
    Lprim::Vector{Rational{BigInt}}
    K::Vector{Polynomial{Rational{BigInt}}}
    H::Vector{Polynomial{Rational{BigInt}}}
    bplus::Vector{Rational{BigInt}}
    bminus::Vector{Rational{BigInt}}
    rplus::Int
    splus::Int
    rminus::Int
    sminus::Int
    function PrecalHermite(ord ; flbis=false)
        ord % 2 == 1 || throw(DomainError("ord=$ord must be odd"))
        d = div(ord,2)
        L = map(i-> _L(i,ord) , -d:d+1)
        Lprim = map(i-> _Lprim(i,ord) , -d:d+1)
        K = map(i-> _K(i,ord) , -d:d+1)
        H = map(i-> _H(i,ord) , -d:d+1)
        if flbis
            rplus = -d-1
            splus = d
        else
            rplus = -d
            splus = d+1
        end
        rminus = -splus
        sminus = -rplus

# ref        bplus = map(i-> i != 0 ? _bplus(i,ord) : 0 , -d:d+1)
        bplus = map(i-> i != 0 ? _bplus(i, rplus, splus) : 0 , rplus:splus)
# ref        bminus = map(i-> i != 0 ? _bminus(i,ord) : 0 , -d-1:d)
        bminus = map(i-> i != 0 ? _bminus(i, rminus, sminus) : 0 , rminus:sminus)
        bplus[1-rplus] = -sum(bplus)
        bminus[1-rminus] = -sum(bminus)
        return new{ord,d}(L,Lprim,K,H,bplus,bminus, rplus, splus, rminus, sminus)
    end
end
L(ph::PrecalHermite{ord,d}, i) where{ord,d} =ph.L[d+1+i]
Lprim(ph::PrecalHermite{ord,d}, i) where{ord,d} =ph.Lprim[d+1+i]
K(ph::PrecalHermite{ord,d}, i) where{ord,d} =ph.K[d+1+i]
H(ph::PrecalHermite{ord,d}, i) where{ord,d} =ph.H[d+1+i]
bplus(ph::PrecalHermite{ord,d}, i) where{ord,d} =ph.bplus[1-ph.rplus+i]
bminus(ph::PrecalHermite{ord,d}, i) where{ord,d} =ph.bminus[1-ph.rminus+i]

struct Hermite{T, edge, order} <: AbstractInterpolation{T,edge,order}
    tabfct::Vector{Polynomial{T}}
    function Hermite(order::Int, T::DataType = Float64; edge::EdgeType = CircEdge, flbis=false)
        if flbis
            order % 4 == 3 || throw(DomainError("order=$order modulo 4 must equal to 3"))
            ord = div(order, 2)
        else
            order % 4 == 1 || throw(DomainError("order=$order modulo 4 must equal to 1"))
            ord = div(order, 2)+1
        end
        decal = div(order+1,2)
        d = div(ord,2)
        ph = PrecalHermite(ord, flbis=flbis)
        tabfct = map( x->Polynomial(Rational{BigInt}[0]), 1:(order+1))
        for i=-d:d+1
            tabfct[decal+i] += H(ph,i)
        end
        for i=-d:0
            for k=ph.rplus:ph.splus
                tabfct[decal+i+k] += bplus(ph,k)*K(ph,i)
            end
        end
        for i=1:d+1
            for k=ph.rminus:ph.sminus
                tabfct[decal+i+k] += bminus(ph,k)*K(ph,i)
            end
        end
        return new{T,edge,order}(tabfct)
    end
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Many auxiliary functions for ground repeating, Sun-synchronous (GRSS) orbit
#   computations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export adjacent_track_distance_grss, adjacent_track_angle_grss
export compute_ss_orbit_by_num_rev_per_day
export list_ss_orbits_by_rep_period, sort_list_ss_orbits_by_height

"""
    adjacent_track_distance_grss(T::Number, i::Number, To::Int, lat::Number)

Compute the distance between adjacent ground tracks [m] at a given latitude
`lat` [rad] for a ground repeating, Sun-synchronous orbit with period `T` [s],
inclination `i` [rad], and orbit cycle `To` [days].

# Remarks

The functions **does not** check if the orbit is a GRSS orbit.

"""
function adjacent_track_distance_grss(T::Number, i::Number, To::Int, lat::Number)
    # Angle between two adjacent traces.
    theta = (T/To)*pi/43200.0

    # Angle of swath.
    beta = asin(sin(theta)*sin(i))

    # Minimum swath.
    beta*R0*cos(lat)
end

"""
    adjacent_track_distance_grss(a::Number, e::Number, i::Number, To::Int, lat::Number)

Compute the distance between adjacent ground tracks [m] at a given latitude
`lat` [rad] for a ground repeating, Sun-synchronous orbit with semi-major axis
`a` [m], eccentricity `e`, inclination `i` [rad], and orbit cycle `To` [days].

# Remarks

The functions *does not* check if the orbit is a GRSS orbit.

"""
function adjacent_track_distance_grss(a::Number, e::Number, i::Number, To::Int,
                                      lat::Number)
    # Orbit period.
    T = period(a, e, i, :J2)

    # Compute the distance between adjacent ground tracks.
    adjacent_track_distance_grss(T, i, To, lat)
end

"""
    adjacent_track_angle_grss(h::Number, T::Number, i::Number, To::Int, lat::Number)

Compute the angle between two adjacent ground tracks [rad] in a given latitude
`lat` [rad] measured from the satellite position for a ground repeating,
Sun-synchronous orbit with altitude in the Equator `h` [m], period `T` [s],
inclination `i` [rad], and orbit cycle `To` [days].

# Remarks

The functions **does not** check if the orbit is a GRSS orbit.

"""
function adjacent_track_angle_grss(h::Number, T::Number, i::Number, To::Int,
                                   lat::Number)
    # Angle between two adjacent traces.
    theta = (T/To)*pi/43200.0

    # Angle of swath.
    beta = asin(sin(theta)*sin(i))

    # Get the Earth radius in the plane perpendicular to the Equatorial plane.
    Rp = R0*cos(lat)

    # Compute the angle between the two ground tracks.
    b = sqrt(Rp^2+(Rp+h)^2-2*Rp*(Rp+h)*cos(beta/2.0))
    asin(Rp/b*sin(beta/2.0))
end

"""
    adjacent_track_angle_grss(h::Number, a::Number, e::Number, i::Number, To::Int, lat::Number)

Compute the angle between two adjacent ground tracks [rad] in a given latitude
`lat` [rad] measured from the satellite position for a ground repeating,
Sun-synchronous orbit with altitude in the Equator `h` [m], semi-major axis `a`
[m], eccentricity `e`, inclination `i` [rad], and orbit cycle `To` [days].

# Remarks

The functions *does not* check if the orbit is a GRSS orbit.

"""
function adjacent_track_angle_grss(h::Number, a::Number, e::Number, i::Number,
                                   To::Int, lat::Number)
    # Period (J2).
    T = period(a, e, i, :J2)

    # Compute the minimum half FOV.
    adjacent_track_angle_grss(h, T, i, To, lat)
end

"""
    compute_ss_orbit_by_num_rev_per_day(numRevPD::Number, e::Number)

Compute the Sun-synchronous orbit given the number of revolutions per day
`numRevPD` and the eccentricity `e`.

# Returns

* The semi-major axis [m].
* The inclination [rad].
* The residues of the two functions.
* A boolean variable that indicates if the numerical algorithm converged.

"""
function compute_ss_orbit_by_num_rev_per_day(numRevPD::Number, e::Number)
    compute_ss_orbit_by_ang_vel(numRevPD*2*pi/86400.0, e)
end

"""
    list_ss_orbits_by_rep_period(minRep::Int, maxRep::Int, minAlt::Number=-1.0, maxAlt::Number=-1.0, e::Number=0.0)

Compute a list of repeating Sun-synchronous orbits.

# Args

* `minRep`: Minimum repetition time of the orbit [days].
* `maxRep`: Maximum repetition time of the orbit [days].
* `minAlt`: Minimum altitude of the orbits on the list [m].
* `maxAlt`: Minimum altitude of the orbits on the list [m].
* `e`: Eccentricity.

# Returns

A matrix containing the orbits found with the following format:

    Semi-major axis [m] | Altitude [m] | Inclination [rad] | Period [s] | Int | Num | Den
    --------------------|--------------|-------------------|------------|-----|-----|-----

in which the period is Int + Num/Den.

# Remarks

If `minAlt` or `maxAlt` is < 0.0, then the altitude will not be checked when a
orbit is added to the list.

"""
function list_ss_orbits_by_rep_period(minRep::Int,         maxRep::Int,
                                      minAlt::Number=-1.0, maxAlt::Number=-1.0,
                                      e::Number=0.0)
    # Check if the arguments are valid.
    minRep <= 0       && throw(ArgumentError("The minimum repetition time must be greater than 0."))
    maxRep <= 0       && throw(ArgumentError("The maximum repetition time must be greater than 0."))
    minRep > maxRep   && throw(ArgumentError("The minimum repetition time must be smaller or equal to the maximum repetition time."))
    !( 0. <= e < 1. ) && throw(ArgumentError("The eccentricity must be within the interval 0 <= e < 1."))

    # Matrix to store the available orbits.
    ss_orbits = Matrix{Float64}(undef, 0, 7)

    # Integer part of the number of orbits per day.
    intNumOrb = (13, 14, 15, 16, 17)

    # Check if the altitude interval must be verified.
    check_altitude = (minAlt > 0) && (maxAlt > 0)
    check_altitude && (minAlt > maxAlt) && throw(ArgumentError("The minimum altitude must be lower than the maximum altitude."))

    # Loop for the possible repetition times.
    for den = minRep:maxRep
        for num = 0:den-1

            # Check if the fraction num/den is irreducible.
            if gcd(num, den) == 1

                # Loop through the integer parts.
                for ino in intNumOrb
                    addOrbit = false

                    # Number of revolutions per day.
                    numRevPD = ino + num/den

                    (a, i, f1r, f2r, converged) =
                        compute_ss_orbit_by_num_rev_per_day(numRevPD, e)

                    # Check if the orbit is valid.
                    !(@check_orbit(a,e)) && continue

                    # Check if the altitude interval must be verified.
                    if check_altitude
                        if (minAlt < a-R0) && (a-R0 < maxAlt) && (converged)
                            addOrbit = true
                        end
                    else
                        addOrbit = converged
                    end

                    # Check if the orbit must be added to the list.
                    if addOrbit
                        # Compute the period of the orbit considering
                        # perturbations (J2).
                        T = period(a, e, i, :J2)

                        ss_orbits = vcat(ss_orbits, [a a-R0 i T ino num den])
                    end
                end
            end
        end
    end

    # Return the list of orbits.
    ss_orbits
end

"""
    sort_list_ss_orbits_by_height(ss_orbits::Matrix)

Sort the list of Sun-synchronous orbits `ss_orbits` (see
`list_ss_orbits_by_rep_period`) by height and return a new matrix.

"""
sort_list_ss_orbits_by_height(ss_orbits::Matrix) =
    sortslices(ss_orbits, dims=1, by=x->x[1])


text:
function Initi_model_struct(H::T,Hc::T,N::T,Wf::T,B::T,αf::T,β₁::T,β₂::T,
    Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    cons = Constants()
    env = EnvironmentStruct()
    kinetic = PhotoKineticRates()
    photo = PhotoPar(kinetic,env.Tₐ)
    hydPar = HydraulicsPar(;Kₓₗ₀=Kₓₗ₀,i=i)
    treepar = TreePar(;αf=αf,β₁=β₁,β₂=β₂,rₘ=rₘ,Nₛ=Nₛ,a_Jmax=a_Jmax,b_Jmax=b_Jmax)

    Hs = H-Hc
    As = Wf/treepar.αf
    Ww = treepar.ρw*As*(treepar.β₁*H+treepar.β₂*Hs)
    treesize = TreeSize(Wf,Ww,H,Hs,As,B,N)

    model = CCPHStruct(cons,env,treepar,treesize,photo,hydPar)

    return model,kinetic
end

function Initi_model_struct(j::Integer,data::RoData,
    αf::T,β₁::T,β₂::T,Nₛ::T,rₘ::T,a_Jmax::T,b_Jmax::T,i::T,Kₓₗ₀::T) where {T<:Float64}
    data_ind = Find_data_ind(j)
    H = data.H[data_ind]
    N = data.N[data_ind]
    Wf = data.Wf[data_ind]
    B = data.B[data_ind]
    Hc = data.Hc
    model,kinetic = Initi_model_struct(H,Hc,N,Wf,B,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

    return model,kinetic
end

function Run_RO_C_F_CCPH(par::Array{Float64,1},weatherts_F::WeatherTS,
    weatherts_C::WeatherTS,data_F::RoData,data_C::RoData;sim_steps::Integer=83)

    αf,β₁,β₂,b_Jmax = 460.0,1.27,-0.27,0.0
    Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i = par[1:6]

    #Fertilized stand
    GPP_model_F = zeros(sim_steps+1)
    EC_model_F = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_F,αf,β₁,β₂,Nₛ,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_F[j],EC_model_F[j] = CalcModelOutput!(j,model,weatherts_F,kinetic)
    end

    #Control stand
    GPP_model_C = zeros(sim_steps+1)
    EC_model_C = zeros(sim_steps+1)
    for j = 1:sim_steps+1

        model,kinetic = Initi_model_struct(j,data_C,αf,β₁,β₂,Nₛ_C,rₘ,a_Jmax,b_Jmax,i,Kₓₗ₀)

        GPP_model_C[j],EC_model_C[j] = CalcModelOutput!(j,model,weatherts_C,kinetic)
    end

    return GPP_model_F,GPP_model_C,EC_model_F,EC_model_C
end

function calc_logP_term(ydata::T,ymodel::T,a::T,b::T) where {T<:Float64}
    e = abs(ydata-ymodel)
    σ = a+b*ymodel
    σ>0.0||error("Negative variance")
    return e/σ+log(2*σ)
end

function Calc_logP(
    GPP_model_F::Array{T,1},Ec_model_F::Array{T,1},GPP_data_F::Array{T,1},Ec_data_F::Array{T,1},
    GPP_model_C::Array{T,1},Ec_model_C::Array{T,1},GPP_data_C::Array{T,1},Ec_data_C::Array{T,1},
    a_GPP::T,b_GPP::T,a_Ec::T,b_Ec::T) where {T<:Float64}

    logP = 0.0
    for i = 1:length(GPP_data_F)
        logP += calc_logP_term(GPP_data_F[i],GPP_model_F[i],a_GPP,b_GPP)
        logP += calc_logP_term(GPP_data_C[i],GPP_model_C[i],a_GPP,b_GPP)

        logP += calc_logP_term(Ec_data_F[i],Ec_model_F[i],a_Ec,b_Ec)
        logP += calc_logP_term(Ec_data_C[i],Ec_model_C[i],a_Ec,b_Ec)
    end

    return  -logP
end

function Post_distri_RO_C_F_CCPH(par::Array{Float64,1},RO_data::RO_raw_data,
    GPP_data_F_in::Array{Float64,1},Ec_data_F_in::Array{Float64,1},
    GPP_data_C_in::Array{Float64,1},Ec_data_C_in::Array{Float64,1})

        try
            X0,τ,Smax,τ_C = par[7:10]
            weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
            weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

            data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

            GPP_model_F,GPP_model_C,EC_model_F,EC_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

            a_GPP,b_GPP,a_Ec,b_Ec = par[11:14]

            logP = Calc_logP(
            GPP_model_F,EC_model_F,GPP_data_F_in,Ec_data_F_in,
            GPP_model_C,EC_model_C,GPP_data_C_in,Ec_data_C_in,
            a_GPP,b_GPP,a_Ec,b_Ec)

            return post = logP

        catch err
            println("Parameters Error: ", err)
            return post = -Inf
        end
end

function getmodeldata(par::Array{Float64,1},RO_data::RO_raw_data)
    X0,τ,Smax,τ_C = par[7:10]
    weatherts_F,GPP_data_F = create_weather_struct_RO(RO_data;stand_type="Fertilized",X0=X0,τ=τ,Smax=Smax)
    weatherts_C,GPP_data_C = create_weather_struct_RO(RO_data;stand_type="Control",X0=X0,τ=τ_C,Smax=Smax)

    data_F,data_C = Create_RoData_C_F(GPP_data_F,GPP_data_C,weatherts_F,weatherts_C)

    GPP_model_F,GPP_model_C,EC_model_F,EC_model_C = Run_RO_C_F_CCPH(par,weatherts_F,weatherts_C,data_F,data_C)

    return GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,weatherts_F,weatherts_C
end

function generatedatafromparameter(par::Array{Float64,1},RO_data::RO_raw_data)

    GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,weatherts_F,weatherts_C = getmodeldata(par,RO_data)

    a_GPP,b_GPP,a_Ec,b_Ec = par[11:14]

    error_GPP_F = []
    error_GPP_C = []
    error_Ec_F = []
    error_Ec_C = []

    for i = 1:length(GPP_model_F)
        push!(error_GPP_F,rand(Laplace(0.0, a_GPP+b_GPP*GPP_model_F[i])))
        push!(error_GPP_C,rand(Laplace(0.0, a_GPP+b_GPP*GPP_model_C[i])))
        push!(error_Ec_F,rand(Laplace(0.0, a_Ec+b_Ec*EC_model_F[i])))
        push!(error_Ec_C,rand(Laplace(0.0, a_Ec+b_Ec*EC_model_C[i])))
    end

    return GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,error_GPP_F,error_GPP_C,error_Ec_F,error_Ec_C,weatherts_F,weatherts_C
end

function run_test()
    RO_data = Load_RO_data()

    par = [0.011976551772617215, 25.42466609569836, 0.5009242298864488, 0.023503724920524185,
    0.010452800544803813,1.0,-1.312709470268466, 3.3784814235281933, 17.510523050796454, 9.566445194647766,
    0.237,0.0845,0.0688,0.146]

    GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,
    error_GPP_F,error_GPP_C,error_Ec_F,error_Ec_C,weatherts_F,weatherts_C = generatedatafromparameter(par,RO_data)

    GPP_data_F = GPP_model_F.+error_GPP_F
    GPP_data_C = GPP_model_C.+error_GPP_C
    EC_data_F = EC_model_F.+error_Ec_F
    EC_data_C = EC_model_C.+error_Ec_C

    plot(weatherts_F.date,GPP_model_F)
    pl1 = plot!(weatherts_F.date,GPP_model_F.+error_GPP_F)
    plot(weatherts_F.date,EC_model_F)
    pl2 = plot!(weatherts_F.date,EC_model_F.+error_Ec_F)

    plot(weatherts_C.date,GPP_model_C)
    pl3 = plot!(weatherts_C.date,GPP_model_C.+error_GPP_C)
    plot(weatherts_C.date,EC_model_C)
    pl4 = plot!(weatherts_C.date,EC_model_C.+error_Ec_C)

    @show GPP_R²_F =  calcR²(GPP_model_F.+error_GPP_F,GPP_model_F)
    @show GPP_cor_F = cor(GPP_model_F.+error_GPP_F,GPP_model_F)

    @show Ec_R²_F =  calcR²(EC_model_F.+error_Ec_F,EC_model_F)
    @show Ec_cor_F = cor(EC_model_F.+error_Ec_F,EC_model_F)

    @show GPP_R²_C =  calcR²(GPP_model_C.+error_GPP_C,GPP_model_C)
    @show GPP_cor_C = cor(GPP_model_C.+error_GPP_C,GPP_model_C)

    @show Ec_R²_C =  calcR²(EC_model_C.+error_Ec_C,EC_model_C)
    @show Ec_cor_C = cor(EC_model_C.+error_Ec_C,EC_model_C)

    plot(pl1,pl2,pl3,pl4,layout=(2,2),legends=false)

    #Nₛ,rₘ,a_Jmax,Kₓₗ₀,Nₛ_C,i,X0,τ,Smax,τ_C,a_GPP,b_GPP,a_Ec,b_Ec
    ranges = [(0.001,0.1),(10.0,40.0),(0.01,1.0),(0.001,0.1),(0.001,0.1),
    (0.1,6.0),(-18.0,-0.9),(1.0,15.0),(10.0,25.0),(1.0,15.0),
    (0.0001,5.0),(0.0001,3.0),(0.0001,5.0),(0.0001,3.0)]

    res = BlackBoxOptim.bboptimize(x::Array{Float64,1}->-Post_distri_RO_C_F_CCPH(x::Array{Float64,1},RO_data,
    GPP_data_F,EC_data_F,GPP_data_C,EC_data_C)
    ; SearchRange = ranges,PopulationSize = 50,MaxSteps=2000)

    @show xopt = BlackBoxOptim.best_candidate(res)

    @show -Post_distri_RO_C_F_CCPH(par,RO_data,GPP_data_F,EC_data_F,GPP_data_C,EC_data_C)

    GPP_model_F,GPP_model_C,EC_model_F,EC_model_C,weatherts_F,weatherts_C = getmodeldata(xopt,RO_data)

    @show GPP_R²_F =  calcR²(GPP_data_F,GPP_model_F)
    @show GPP_cor_F = cor(GPP_data_F,GPP_model_F)
    @show gₛ_R²_F =  calcR²(EC_data_F,EC_model_F)
    @show gₛ_cor_F = cor(EC_data_F,EC_model_F)

    @show GPP_R²_C =  calcR²(GPP_data_C,GPP_model_C)
    @show GPP_cor_C = cor(GPP_data_C,GPP_model_C)
    @show gₛ_R²_C =  calcR²(EC_data_C,EC_model_C)
    @show gₛ_cor_C = cor(EC_data_C,EC_model_C)

    a_GPP,b_GPP,a_Ec,b_Ec = xopt[11:14]

    σ_GPP_F = a_GPP.+GPP_model_F*b_GPP
    σ_GPP_C = a_GPP.+GPP_model_C*b_GPP

    σ_EC_F = a_Ec.+b_Ec*EC_model_F
    σ_EC_C = a_Ec.+b_Ec*EC_model_C

    c = -log(0.025*2) #Use for calculating 95% credible intervals

    plot(weatherts_F.date,GPP_data_F,label="Data",ylabel="GPP")
    plot!(weatherts_F.date,GPP_model_F+c*σ_GPP_F)
    plot!(weatherts_F.date,GPP_model_F-c*σ_GPP_F)
    pl1 = plot!(weatherts_F.date,GPP_model_F,label="Model")
    plot(weatherts_C.date,GPP_data_C,label="Data",ylabel="GPP")
    plot!(weatherts_C.date,GPP_model_C+c*σ_GPP_C)
    plot!(weatherts_C.date,GPP_model_C-c*σ_GPP_C)
    pl2 = plot!(weatherts_C.date,GPP_model_C,label="Model")

    plot(weatherts_F.date,EC_data_F,label="Data",ylabel="E_C")
    plot!(weatherts_F.date,EC_model_F+c*σ_EC_F)
    plot!(weatherts_F.date,EC_model_F-c*σ_EC_F)
    pl3 = plot!(weatherts_F.date,EC_model_F,label="Model")
    plot(weatherts_C.date,EC_data_C,label="Data",ylabel="E_C")
    plot!(weatherts_C.date,EC_model_C+c*σ_EC_C)
    plot!(weatherts_C.date,EC_model_C-c*σ_EC_C)
    pl4 = plot!(weatherts_C.date,EC_model_C,label="Model")

    plot(pl1,pl3,pl2,pl4,layout=(2,2),legends=false)
end

run_test()


text:
using SumBassSolver
using Random

f = [x-> 1/(1 + exp(-x)), x->2/(1 + exp(-x/2))]
df = [x->exp(-x)/(1 + exp(-x))^2, x-> exp(-x/2)/(1 + exp(-x/2))^2]
z = [0,0]
F = [-1,-1]
T = [2,2]
Pi = 2

sumSigmoidSolverSB(f,df,z,F,T,Pi)

Random.seed!(121)
# Generate random data. For implementation input proper data
n=50
p = 0.003*rand(n)
q = 0.3*rand(n)
m = 1000 .+ 500*rand(n)
c0 = m .* rand(n) ./ 2
alpha = rand(n)
F = zeros(n)
T = fill(100.0,n)
Pi = 2500

f = Vector{Function}(undef,n)
df = Vector{Function}(undef,n)
z = Vector{Float64}(undef,n)
df_inv = Vector{Function}(undef,n)
for i in 1:n
    f[i] = t -> alpha[i] * Bass(t,p[i],q[i],m[i],c0[i])
    df[i] = t -> alpha[i] * Bass_derivative(t,p[i],q[i],m[i],c0[i])
    z[i] = Bass_inflection(p[i],q[i],m[i],c0[i])
    df_inv[i] = t -> Bass_derivative_inverse(t/alpha[i],p[i],q[i],m[i],c0[i])[1]
end

x1,v1 = sumSigmoidSolverSB(f,df,z,F,T,Pi)
x2,v2 = sumSigmoidSolverSB(f,df,z,F,T,Pi,df_inv)
x3,v3 = @time sumBassSolver(p,q,m,c0,alpha,F,T,Pi,alg="SB",localopt=true)

max_stars_repo_name: andres-alban/res-alloc-bass-demands

max_issues_repo_name: andres-alban/res-alloc-bass-demands

max_forks_repo_name: andres-alban/res-alloc-bass-demands


text:
module BLASBenchmarksCPU

# BLAS libs (& Libdl)
@static if Sys.ARCH === :x86_64
using MKL_jll
end
using OpenBLAS_jll, blis_jll#, Libdl
# Julia BLAS
using Tullio, Octavian, Gaius

# utils: LoopVectorization for Tullio.jl, VectorizationBase for info
using LoopVectorization, VectorizationBase
using VectorizationBase: num_cores, align
using RecursiveFactorization

using Random

# Adjoint
using LinearAlgebra

# Utils
using BenchmarkTools, ProgressMeter

# Plotting & presenting results
using Cairo
using Fontconfig, Gadfly, Colors, DataFrames

export benchmark_result_type
export benchmark_result_df
export benchmark_result_threaded
export logspace
export plot
export runbench

# BLIS
export gemmblis!
export blis_set_num_threads

# Octavian.jl
export matmul!

# OpenBLAS
export gemmopenblas!
export openblas_set_num_threads

# MKL
export gemmmkl!, gemmmkl_direct!
export mkl_set_num_threads

# set threads
@static if Sys.ARCH === :x86_64
const libMKL = MKL_jll.libmkl_rt # more convenient name
function mkl_set_num_threads(N::Integer)
    ccall((:MKL_Set_Num_Threads,libMKL), Cvoid, (Int32,), N % Int32)
end
end
const libOPENBLAS = OpenBLAS_jll.libopenblas # more convenient name
function openblas_set_num_threads(N::Integer)
    ccall((:openblas_set_num_threads64_,libOPENBLAS), Cvoid, (Int64,), N)
end

const libBLIS = blis_jll.blis # more convenient name
function blis_set_num_threads(N::Integer)
    ccall((:bli_thread_set_num_threads,libBLIS), Cvoid, (Int32,), N)
end
function blis_get_num_threads(N::Integer)
    ccall((:bli_thread_get_num_threads,libBLIS), Int32, ())
end

include("ccallblas.jl")
include("benchconfig.jl")
include("runbenchmark.jl")
include("plotting.jl")

function __init__()
    Sys.ARCH === :x86_64 && mkl_set_num_threads(num_cores())
    openblas_set_num_threads(num_cores())
    blis_set_num_threads(num_cores())
end

end


text:
#=
Examples of strings.
=#

# String as a double quote
s1 = "This is a string"
println(s1)

# String as a triple double quote
s2 = """This is "another" string"""
println(s2)

# Indexing a string
println("s2[begin] ", s2[begin])
println("s2[end] ", s2[end])
println("s2[1] ", s2[1])
println("s2[2:5] ", s2[2:6])

# Combining strings with $
s3 = "Hello"
s4 = "Julia"
println("$s3 $s4 programming!")

# Combining strings with string()
s5 = string(s3, " ", s4, " programming!!")
println(s5)

# Combining strings with *
s6 = s3 * " " * s4 * " programming!!!"
println(s6)

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
struct RandomQs <: QueryStrategy end

function qs_score(qs::RandomQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    return rand(size(x, 2))
end


text:
using FastLapackInterface
using LinearAlgebra
using Test

include("LinSolveAlgo_test.jl")
include("SchurAlgo_test.jl")
include("QrAlgo_test.jl")


text:
#=
Running order:
1. refresh_norsk_dt_210423()
2. cv_210510()
=#
using DataFrames, JLD2, CSV, GZip, LinearAlgebra, SparseArrays, Statistics, CodecZlib
import RD2:dat_dir

function ccwcv(yt, wt, yv, Gi, X)
    nt = length(yt)
    nv = length(yv)
    nx = size(X)[2] + 1

    #Ri = diagm(wt)
    XR = X .* wt
    Z  = sparse([I zeros(nt, nv)])
    ZR = Z .* wt
    lhs = [ XR'X XR'Z
            ZR'X ZR'Z + Gi ] # as va = ve = 1, h2=0.5.
    rhs = [ XR'yt
            ZR'yt ]
    bv = lhs \ rhs
    cc = cor(yv, bv[nt+nx:end])
    return cc, bv[nt+nx:end]
end

function cccv(yt, yv, Gi, X)
    nt = length(yt)
    nv = length(yv)
    nx = size(X)[2] + 1

    Z = sparse([I zeros(nt, nv)])
    lhs = [ X'X X'Z
            Z'X Z'Z + Gi ]      # as va = ve = 1, h2=0.5.
    rhs = [ X'yt
            Z'yt ]
    bv = lhs \ rhs
    cc = cor(yv, bv[nt+nx:end])
    return cc, bv[nt+nx:end]
end

"""
    function mtwcv(y1, y2, yv, wt, ct, Gi, X)
---
## `M`ultiple `t`rait `w`eighted `c`ross `v`alidation.

- `y1`: the 1st training set
- `y2`: the 2nd training set
- `yv`: the validation set
- `wt`: weight for training values
- `ct`: correlation between the traits
- `Gi`: inverse of GRM
- `X`: fixed effects matrix
"""
function mtwcv(y1, y2, yv, wt, ct, Gi, X, grp)
    XR = X .* wt
    yt = [y1; y2]

    n1 = length(y1)
    n2 = length(y2)
    nt = n1 + n2
    nv = length(yv)
    nn = nt + nv
    nx = size(X)[2]

    Z = sparse(zeros(nt, 2nn))
    Z[1:n1, 1:n1] = I(n1)
    Z[n1+1:n1+n2, nn+n1+1:nn+n1+n2] = I(n2)
    ZR = Z .* wt

    T = inv([1 ct; ct 1])
    lhs = [ XR'X XR'Z
            ZR'X ZR'Z + kron(T, Gi) ]
    rhs = [ XR'yt
            ZR'yt ]

    ebv = lhs \ rhs             # all ebv
    fra = (grp == 1) ?
        nx + n1 + n2 :
        nx + n1 + n2 + nn

    vbv = ebv[fra+1:fra+nv]     # ebv for the validation set
    cc = cor(yv, vbv)
    return cc, vbv
end

function mtcv(y1, y2, yv, ct, Gi, X, grp)
    yt = [y1; y2]

    n1 = length(y1)
    n2 = length(y2)
    nt = n1 + n2
    nv = length(yv)
    nn = nt + nv
    nx = size(X)[2]

    Z = sparse(zeros(nt, 2nn))
    Z[1:n1, 1:n1] = I(n1)
    Z[n1+1:n1+n2, nn+n1+1:nn+n1+n2] = I(n2)

    T = inv([1 ct; ct 1])
    lhs = [ X'X X'Z
            Z'X Z'Z + kron(T, Gi) ]
    rhs = [ X'yt
            Z'yt ]

    ebv = lhs \ rhs             # all ebv
    fra = (grp == 1) ?
        nx + n1 + n2 :
        nx + n1 + n2 + nn
    vbv = ebv[fra+1:fra+nv]     # ebv for the validation set
    cc = cor(yv, vbv)
    return cc, vbv
end

function cv_G(d1, d2, ct, vs, G, trait, erc, X; grp = 1)
    y1 = d1[:, trait]
    y2 = d2[:, trait]

    w1 = d1[:, erc]
    w2 = d2[:, erc]
    wt = [w1; w2]

    yv = vs[:, trait]

    i1 = d1[:, :ix]
    i2 = d2[:, :ix]
    i3 = vs[:, :ix]
    id = [i1; i2; i3]

    gi = inv(G[id, id])

    c1, _ = mtwcv(y1, y2, yv, wt, ct, gi, X, grp)
    c2, _ = mtcv(y1, y2, yv, ct, gi, X, grp)
    return c1, c2
end

function cv_G(ts, vs, G, trait, erc, X)
    yt = ts[:, trait]
    wt = ts[:, erc]

    yv = vs[:, trait]

    it = ts[:, :ix]
    iv = vs[:, :ix]
    id = [it; iv]
    gi = inv(G[id, id])

    c1, _ = ccwcv(yt, wt, yv, gi, X)
    c2, _ = cccv(yt, yv, gi, X)
    return c1, c2
end

function make_X(len...)
    nc = length(len)
    nr = sum(len)
    X = zeros(nr, nc)
    ac = 0
    for i in 1:nc
        X[ac+1:ac+len[i], i] .=1
        ac += len[i]
    end
    X
end

function mk_cv_tbl(dts, dvs, gts, gvs, nts, nvs, G, trait, erc)
    # target is each country in turn
    rst = zeros(6, 7)           # 6 = w/nw by (DGN); 7 = D G N DG DN GN DGN
    ct = 0.9
    println("correlation = $ct")

    @info "Prediction of one country with one of 3 coutries"
    # D → D
    X = make_X(nrow(dts))
    rst[1, 1], rst[4, 1] = cv_G(dts, dvs, G, trait, erc, X)
    # G → D
    X = make_X(nrow(gts))
    rst[1, 2], rst[4, 2] = cv_G(gts, dvs, G, trait, erc, X)
    # N → D
    X = make_X(nrow(nts))
    rst[1, 3], rst[4, 3] = cv_G(nts, dvs, G, trait, erc, X)
    # D → G
    X = make_X(nrow(dts))
    rst[2, 1], rst[5, 1] = cv_G(dts, gvs, G, trait, erc, X)
    # G → G
    X = make_X(nrow(gts))
    rst[2, 2], rst[5, 2] = cv_G(gts, gvs, G, trait, erc, X)
    # N → G
    X = make_X(nrow(nts))
    rst[2, 3], rst[5, 3] = cv_G(nts, gvs, G, trait, erc, X)
    # D → G
    X = make_X(nrow(dts))
    rst[3, 1], rst[6, 1] = cv_G(dts, nvs, G, trait, erc, X)
    # G → G
    X = make_X(nrow(gts))
    rst[3, 2], rst[6, 2] = cv_G(gts, nvs, G, trait, erc, X)
    # N → G
    X = make_X(nrow(nts))
    rst[3, 3], rst[6, 3] = cv_G(nts, nvs, G, trait, erc, X)

    @info "Prediction of one country with two of 3 coutries"
    dgts = vcat(dts, gts)
    dnts = vcat(dts, nts)
    gnts = vcat(gts, nts)
    # D + G → D
    X = make_X(nrow(dts), nrow(gts))
    rst[1, 4], rst[4, 4] = cv_G(dgts, dvs, G, trait, erc, X)
    # D + N → D
    X = make_X(nrow(dts), nrow(nts))
    rst[1, 5], rst[4, 5] = cv_G(dts, nts, ct, dvs, G, trait, erc, X)
    # G + N → D
    X = make_X(nrow(gts), nrow(nts))
    rst[1, 6], rst[4, 6] = cv_G(gts, nts, ct, dvs, G, trait, erc, X)
    # D + G → G
    X = make_X(nrow(dts), nrow(gts))
    rst[2, 4], rst[5, 4] = cv_G(dgts, gvs, G, trait, erc, X)
    # D + N → G
    X = make_X(nrow(dts), nrow(nts))
    rst[2, 5], rst[5, 5] = cv_G(dts, nts, ct, gvs, G, trait, erc, X)
    # G + N → G
    X = make_X(nrow(gts), nrow(nts))
    rst[2, 6], rst[5, 6] = cv_G(gts, nts, ct, gvs, G, trait, erc, X)
    # D + G → N
    X = make_X(nrow(dts), nrow(gts))
    rst[3, 4], rst[6, 4] = cv_G(dts, gts, ct, nvs, G, trait, erc, X, grp = 2)
    # N + D → N
    X = make_X(nrow(nts), nrow(dts))
    rst[3, 5], rst[6, 5] = cv_G(nts, dts, ct, nvs, G, trait, erc, X)
    # N + G → N
    X = make_X(nrow(nts), nrow(gts))
    rst[3, 6], rst[6, 6] = cv_G(nts, gts, ct, nvs, G, trait, erc, X)

    @info "Prediction of one country with all 3 coutries"
    #dgnts = vcat(dts, gts, nts)
    # D + G + N → D
    X = make_X(nrow(dts), nrow(gts), nrow(nts))
    rst[1, 7], rst[4, 7] = cv_G(dgts, nts, ct, dvs, G, trait, erc, X)
    # D + G + N → G
    X = make_X(nrow(dts), nrow(gts), nrow(nts))
    rst[2, 7], rst[5, 7] = cv_G(dgts, nts, ct, gvs, G, trait, erc, X)
    # D + G + N → N
    X = make_X(nrow(dts), nrow(gts), nrow(nts))
    rst[3, 7], rst[6, 7] = cv_G(dgts, nts, ct, nvs, G, trait, erc, X, grp = 2)

    round.(rst, sigdigits = 2)
end

function make_tv(ts, cs, bv, drp, erc)
    t1 = filter(row ->row.st == 't', cs)
    df1 = select(innerjoin(t1, ts, on = :ID), :ix, drp, erc)
    dft = filter(row -> row[drp] < 999, df1)
    t2  = filter(row -> row.st == 'v', cs)
    dfv = select(innerjoin(t2, bv, on = :ID), :ix, drp)
    return dft, dfv
end

function cv_210621()
    @load "$dat_dir/jld/cv-setup.jld" dcs gcs ncs
    @load "$dat_dir/jld/drp-training-2021-04-23.jld" dts gts nts
    @load "$dat_dir/jld/ebv_202105.jld"  dbv gbv nbv
    @load "$dat_dir/jld/GAID.jld"     G A ID

    begin                       # The milk trait
        dt, dv = make_tv(dts, dcs, dbv, :Milk, :Milk_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :Milk, :Milk_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :Milk, :Milk_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Milk, :Milk_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
    begin                       # The protein trait
        dt, dv = make_tv(dts, dcs, dbv, :Protein, :Protein_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :Protein, :Protein_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :Protein, :Protein_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Protein, :Protein_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
    begin                       # The fat trait
        dt, dv = make_tv(dts, dcs, dbv, :Fat, :Fat_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :Fat, :Fat_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :Fat, :Fat_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :Fat, :Fat_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
    begin                       # The SCS trait
        dt, dv = make_tv(dts, dcs, dbv, :SCS, :SCS_ERC)
        gt, gv = make_tv(gts, gcs, gbv, :SCS, :SCS_ERC)
        nt, nv = make_tv(nts, ncs, nbv, :SCS, :SCS_ERC)
        rst = mk_cv_tbl(dt, dv, gt, gv, nt, nv, G, :SCS, :SCS_ERC)
        for row in eachrow(rst)
            println("| ", join(row, " | "), " |")
        end
    end
end

#=
proto type from Theo.

function mt_am(y, RI, ifix, itr, id, GRMI, GI)
    # y = n*1 vector of records across traits
    # RI = n*1 vector of weights of records (inverse of residual variances)
    # ifix = n*1 indicates fixed effect class variables across traits (i.e. if ifix=1 records belong to the same fixed effect class irrespective trait)
    # itr = n*1 indicates trait of record
    # id = n*1 indicates ID of animal (agrees with row/column in GRM)
    # GRMI = q*q inverse of matrix of genomic relationships
    # GI = ntrait*ntrait inverse of matrix of genetic (co)variances of the traits
    ntrait = size(GI, 1)
    nanim = size(GRMI, 1)
    n = size(y, 1)
    X = zeros(n, maximum(ifix))
    XR = zero(n, maximum(ffix))
    for i = 1:n
        X[i,ifix[i]] = 1
        XR[i,ifix[i]] = RI[i]
    end
    Z = zeros(n, nanim * ntrait)
    ZR = zeros(n, nanim*ntrait)
    for i = 1:n
        iid = (itr[i] - 1) * nanim + id[i]
        Z[i, iid] = 1
        ZR[i, iid] = RI[i]
    end
    LHS = [XR'X XR'Z
           ZR'X ZR'Z + kron(GI, GRMI)]
    RHS = [XR'y
           ZR'y]
    EBV = LHS\RHS
    return EBV #size(EBV) = nanim*ntrait; EBV' = [ebv_tr1 ebv_tr2 ebv_tr3 ....]
end
=#


text:
using Javis
using Documenter

makedocs(;
    modules = [Javis],
    authors = "Ole Kröger <o.kroeger@opensourc.es> and contributors",
    repo = "https://github.com/Wikunia/Javis.jl/blob/{commit}{path}#L{line}",
    sitename = "Javis.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://Wikunia.github.io/Javis.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Tutorials" => [
            "tutorials.md",
            "tutorials/tutorial_1.md",
            "tutorials/tutorial_2.md",
            "tutorials/tutorial_3.md",
            "tutorials/tutorial_4.md",
        ],
        "Mission" => "mission.md",
        "References" => "references.md",
        "Contributing" => "contributing.md",
    ],
)

deploydocs(; repo = "github.com/Wikunia/Javis.jl")


text:
#!/usr/bin/env julia
# Tim Sterne-Weiler 2015

const dir = abspath( splitdir(@__FILE__)[1] )
const ver = chomp(readline(open(dir * "/../VERSION")))

start = time_ns()
println( stderr, "Whippet $ver loading and compiling... " )

using Libz
using ArgParse
using StatsBase
using BioSequences
using Combinatorics

push!( LOAD_PATH, dir * "/../../src" )
using Whippet

function parse_cmd()
  s = ArgParseSettings()
  # TODO finish options...
  @add_arg_table s begin
    "--index", "-x"
      help = "Prefix or full-name of index 'dir/prefix' (default Whippet/index/graph)"
      arg_type = String
      default  = fixpath( "$(dir)/../../index/graph" )
    "--out", "-o"
      help = "Where should the gzipped output go 'dir/prefix'?"
      arg_type = String
      default  = fixpath( "$(dir)/../../simul" )
    "--num-nodes", "-n"
      help = "Maximum number of consecutive nodes to output sliding combinatorial paths through"
      arg_type = Int64
      default  = 8
    "--num-genes", "-g"
      help = "Choose the first -n genes from the index to simulate (starting at offset -i)."
      arg_type = Int64
      default  = 100000
    "--offset", "-i"
      help = "Offset of gene number to start simulating from"
      arg_type = Int64
      default  = 1
    "--theoretical"
      action = :store_true
    "--verbose"
      help = "Print progress messages.."
      action = :store_true
  end
  return parse_args(s)
end

function main()

   args = parse_cmd()

   println(stderr, " $( round( (time_ns()-start)/1e9, digits=6 ) ) seconds" )

   println(stderr, "Loading splice graph index... $( args["index"] ).jls")
   @timer const lib = open(deserialize, "$( args["index"] ).jls")

   println(stderr, "Simulating combinatorial transcripts..")
   @timer simulate_genes( lib, output=args["out"],
                               gene_num=args["num-genes"],
                               offset=args["offset"],
                               node_num=args["num-nodes"],
                               verbose=args["verbose"],
                               theoretical=args["theoretical"] )
end

mutable struct SimulTranscript
   seq::Whippet.SGSequence
   nodes::Vector{UInt}
end

struct SimulGene
   trans::Vector{SimulTranscript}
   gene::String
end

istxstart( edge::EdgeType ) = edge == EDGETYPE_SL || edge == EDGETYPE_LS ? true : false
istxstop(  edge::EdgeType ) = edge == EDGETYPE_SR || edge == EDGETYPE_RS ? true : false

function collect_txstarts( sg::SpliceGraph )
   starts = Int[]
   for i in 1:length(sg.nodelen)
      if istxstart( sg.edgetype[i] )
         push!( starts, i )
      end
   end
   starts
end

function collect_txstops( sg::SpliceGraph )
   stops = Int[]
   for i in 2:length(sg.edgetype)
      if istxstop( sg.edgetype[i] )
         push!( stops, i-1 )
      end
   end
   stops
end

function combinatorial_paths( sg::SpliceGraph; max_nodes=10, verbose=false, theoretical=false )
   verbose && println("collecting starts & stops")
   starts = collect_txstarts( sg )
   stops  = collect_txstops( sg )
   paths  = Vector{Vector{Int16}}()
   verbose && println("going through $(length(starts)) starts and $(length(stops)) stops")
   isocnt = 0
   junc = Dict{Tuple{Int,Int},Int}()

   function add_juncs!( dict, path )
      for i in 1:length(path)-1
         dict[(path[i],path[i+1])] = 1
      end
   end

   function add_combinations!( paths, i, j, sub_range, val )
       for m in combinations( sub_range, val )
          for iv in reverse(i)
             pushfirst!( m, iv )
          end
          for jv in j
             push!( m, jv )
          end
          if isvalid_path( m, sg )
             isocnt += 1
             if theoretical
                add_juncs!( junc, m )
             else
                push!( paths, map(x->convert(Int16,x), m) )
             end
          end
       end
   end

   for i in starts, j in stops
      if i == j
         push!( paths, Int16[i] )
      elseif i == j - 1
         push!( paths, Int16[i,j] )
      elseif i < j
         range = i+1:j-1
         for k in 0:length(range)
            verbose && println(" going through combinations of $range, $k")
            if length(range) >= max_nodes
               sub_range = first(range):(first(range)+(length(range)-max_nodes))
               for r in sub_range
                  verbose && println(" adding sub_combinations of $(r:r+max_nodes-1) with $(i:r-1) and $(r+max_nodes:j) flanking")
                  add_combinations!( paths, i:r-1, r+max_nodes:j, r:r+max_nodes-1, k )
               end
            else
               add_combinations!( paths, [i], [j], range, k )
            end
         end
      end
   end
   paths, isocnt, length(values(junc))
end

function isvalid_path( path::Vector{Int}, sg::SpliceGraph )
   for i in 1:length(path)-1
      if path[i]+1 == path[i+1] # adjacent nodes
         if sg.edgetype[path[i+1]] in (EDGETYPE_LS, EDGETYPE_SR)
            return false
         end
      else
         if !(sg.edgetype[path[i]+1] in (EDGETYPE_LS, EDGETYPE_LR, EDGETYPE_LL) &&
              sg.edgetype[path[i+1]] in (EDGETYPE_SR, EDGETYPE_LR, EDGETYPE_RR))
            return false
         end
      end
   end
   true
end

function simulate_genes( lib; output="simul_genes", gene_num=length(lib.graphs),
                              offset=1, node_num=10, verbose=false, theoretical=false )
   fastaout = open( output * ".fa.gz", "w" )
   gtfout   = open( output * ".gtf.gz", "w" )
   fastastr = ZlibDeflateOutputStream( fastaout )
   gtfstr = ZlibDeflateOutputStream( gtfout )

   isocnt  = 0
   junccnt = 0

   for g in offset:min(offset+gene_num-1, length(lib.graphs)) #sample( 1:length(lib.graphs), min( length(lib.graphs), gene_num ), replace=false )
      sgene = SimulGene( Vector{SimulTranscript}(), lib.names[g] )
      iso,junc = simulate_transcripts( fastastr, gtfstr, sgene, lib.graphs[g], lib.info[g],
                                             node_num=node_num, verbose=verbose, theoretical=theoretical )
      #output_nodes( nodesstr, lib.names[g], lib.info[g], lib.graphs[g] )
      isocnt  += iso
      junccnt += junc
   end
   println(stderr, "$isocnt Unique Isoforms Simulated...")
   println(stderr, "$junccnt Unique Exon-exon Junctions..")
   close( fastastr )
   close( fastaout )
   close( gtfstr )
   close( gtfout )
end

function simulate_transcripts( fastream, gtfstream, simul::SimulGene, sg::SpliceGraph, info::GeneInfo;
                               node_num=10, verbose=false, theoretical=false )
   verbose && println("simulating paths for length $(length(sg.nodelen))")
   paths,isocnt,junccnt = combinatorial_paths( sg, max_nodes=node_num, verbose=verbose, theoretical=theoretical )
   verbose && println("length $(length(paths))")
   theoretical && return isocnt,junccnt
   for s in paths
      txid = simul.gene * "_" * join( s, "-" )
      output_transcript( fastream, s, sg, header=txid )
      output_gtf( gtfstream, s, sg, info, gene_id=simul.gene, transcript_id=txid )
   end
   isocnt, junccnt
end

function output_transcript( stream, path, sg::SpliceGraph; header::String="PREFIX" )
   write( stream, '>' )
   write( stream, header )
   write( stream, "\n" )
   for n in path
      noderange = sg.nodeoffset[n]:(sg.nodeoffset[n]+sg.nodelen[n]-1)
      for i in sg.seq[ noderange ]
         write( stream, convert(Char, i) )
      end
   end
   write( stream, "\n" )
end

function output_gtf( stream, path, sg::SpliceGraph, info::GeneInfo; gene_id::String="GENE", transcript_id::String="TXID" )
   start,stop = 0,0
   i = 1
   while i <= length(path)
      if i < length(path) && (sg.nodecoord[path[i+1]] == sg.nodecoord[path[i]]+sg.nodelen[path[i]] ||
                              sg.nodecoord[path[i+1]]+sg.nodelen[path[i+1]] == sg.nodecoord[path[i]])
         if info.strand
            start = start > 0 ? start : sg.nodecoord[path[i]]
            stop  = sg.nodecoord[path[i+1]]+sg.nodelen[path[i+1]]-1
         else
            start = sg.nodecoord[path[i+1]]
            stop  = stop > 0 ? stop : sg.nodecoord[path[i]]+sg.nodelen[path[i]]-1
         end
      else
         start = start > 0 ? start : Int(sg.nodecoord[path[i]])
         stop  = stop  > 0 ? stop  : Int(sg.nodecoord[path[i]]+sg.nodelen[path[i]]-1)

         write( stream, info.name * "\tWHIPPET_COMBINATORIAL\texon\t" )
         write( stream, string(start) * "\t" )
         write( stream, string(stop) * "\t\.\t" )
         write( stream, info.strand ? '+' : '-' )
         write( stream, "\t\.\tgene_id \"" * gene_id * "\"\; " )
         write( stream, "transcript_id \"" * transcript_id * "\"\;\n" )

         start,stop = 0,0
      end
      i += 1
   end
end

function output_nodes( stream, gname::String, info, sg::SpliceGraph)
   for n in 1:length(sg.nodelen)
      coord = "$(info.name):$(sg.nodecoord[n])-$(sg.nodecoord[n] + sg.nodelen[n] - 1)"
      write( stream, "$gname\t$n\t$coord\t$(info.gene)\n" )
   end
end

@timer main()