text:
##################################################
# Go into (Jr,L) space
##################################################

"""
    radius_s_bounds_action(Jr,L,[nbAvr=nbAvr_default])

Computes the bounds of the orbit with action parameters `(L, Jr)`.
Bounds are in the variable `s`, defined as `s^2 = 1 + x^2`, where `x = r/b`.

# Remarks
- Always takes in allowed orbital parameters.
- Uses `PolynomialRoots`.

# Arguments
- `Jr::Float64`: Radial action parameter.
- `L::Float64`: Angular momentum parameter.
- `nbAvr::Int64`: Orbit average sampling number.
"""
function radius_s_bounds_action(Jr::Float64, L::Float64, nbAvr::Int64=nbAvr_default)

    E = _E_from_Jr(Jr,L,nbAvr)

    tE = _tE(E)
    tL = _tL(L)
    if (L != 0.0)
        rts = sort(real(roots([1,-(tE-tL^2/2),-1,tE])))
        return Float64(rts[2]), Float64(rts[3])
    else
        return 1.0, 1/tE
    end

end

"""
    halfperiod_action(Jr,L,[nbAvr=nbAvr_default])

Computes half the radial period of the orbit with action parameters `(L, Jr)`.

# Remarks
- Well-defined everywhere by using an effective anomaly.

# Arguments
- `Jr::Float64`: Radial action parameter.
- `L::Float64`: Angular momentum parameter.
- `nbAvr::Int64`: Orbit average sampling number.
"""
function halfperiod_action(Jr::Float64, L::Float64, nbAvr::Int64=nbAvr_default)

    sp, sa = radius_s_bounds_action(Jr,L)
    halfperiod = 0.0

    for iu=1:nbAvr
        uloc = -1+2*(iu-0.5)/nbAvr
        jac_loc = dth1du(uloc,sp,sa)
        halfperiod += jac_loc
    end

    halfperiod *= (1/_Omega0)*(2/nbAvr)

    return halfperiod
end

"""
    PHIeff(s, L)

Compute the effective potential psi_eff in the variable `s`.
Variable `s` is defined as `s^2 = 1 + x^2`, where `x = r/b`.

# Arguments
- `s::Float64`: Radial parameter.
- `L::Float64`: Angular momentum parameter.
"""
function PHIeff(s::Float64, L::Float64)
    return _E0/s + L^2/(2*_b^2*(s^2-1))
end

"""
    dPHIeff(s, L)

Compute the derivatives of `PHIeff`, in the variables `s` and `L`.
Return a tuple `(dPHIeff/ds, dPHIeff/dL)`.

# Arguments
- `s::Float64`: Radial parameter.
- `L::Float64`: Angular momentum parameter.
"""
function dPHIeff(s::Float64, L::Float64)
    ds = -_E0/s^2 - s*L^2/(_b^2*(s^2-1)^2)
    dL = L/(_b^2*(s^2-1))
    return ds, dL
end

"""
    _fAnomaly(u)

Function f(u) used to define the effective anomaly of the Plummer cluster.

# Arguments
- `u::Float64`: Effective anomaly.
"""
function _fAnomaly(u::Float64)
    return u*(1.5-u^2/2)
end

"""
    _djacdsp(u,sp,sa)

`sp`-derivative of the non-dimensional Jacobian of the transformation from angles `θ` to effective anomaly `u`.

# Arguments
- `u::Float64`: Effective anomaly.
- `sp::Float64`: Pericenter in the variable `s`.
- `sa::Float64`: Apocenter in the variable `s`.
"""
function _djacdsp(u::Float64, sp::Float64, sa::Float64)

    num = (3 *(4 *sp^3 *(-1 + u)^2 *(12 - 3 *u^2 + 2 *u^3 + u^4) +
           sa *sp^2 *(108 - 120 *u - 27 *u^2 + 40 *u^3 + 18 *u^4 - 3 *u^6 +
           3 *sp^2 *(-1 + u)^4 *(2 + u)^2) -
           2 *sa^2 *sp *(-6 - 3 *u + u^3) *(6 - 3 *u + u^3 +
           sp^2 *(2 - 3 *u + u^3)) -
           sa^3 *(-2 + u) *(1 + u)^2 *(6 + 3* u - u^3 + sp^2 *(6 - 3 *u + u^3))))

    den = (8 *sqrt(2)* sp *(sa +
           sp) *sqrt(-(((-4 + u^2) *(sa^2 *sp *(-2 + u) *(1 + u)^2 -
           sp *(6 - 3 *u + u^3) +
           sa *(-6 - 3 *u + u^3 - sp^2 *(2 - 3* u + u^3))))/(
           sa *sp *(sa + sp) *(sa *(-2 + u) *(1 + u)^2 -
           sp *(2 - 3 *u + u^3)))))* (sa^2* sp *(2 + 3 *u - u^3) +
           sp *(6 - 3* u + u^3) + sa *(6 + 3*u - u^3 + sp^2 *(2 - 3 *u + u^3))))

    return num/den
end

"""
    _djacdsa(u,sp,sa)

`sa`-derivative of the non-dimensional Jacobian of the transformation from angles `θ` to effective anomaly `u`.

# Arguments
- `u::Float64`: Effective anomaly.
- `sp::Float64`: Pericenter in the variable `s`.
- `sa::Float64`: Apocenter in the variable `s`.
"""
function _djacdsa(u::Float64, sp::Float64, sa::Float64)

    num = (3 *(-4 + u^2) *(3 *sa^4 *sp *(-2 + u)^2 *(1 + u)^4 +
           sp^3 *(-1 + u)^2 *(12 - 3 *u^2 + 2 *u^3 + u^4) -
           2 *sa *sp^2 *(-36 + 9 *u^2 - 6 *u^4 + u^6) -
           2 *sa^3 *(-2 + u) *(1 + u)^2 *(12 + 6 *u - 2 *u^3 +
           sp^2 *(6 - 3 *u + u^3)) +
           sa^2 *sp *(108 + 120 *u - 27 *u^2 - 40*u^3 + 18 *u^4 - 3 *u^6 -
           sp^2 *(-12 + 12 *u + 9 *u^2 - 4 *u^3 - 6 *u^4 + u^6))))

    den = (8 *sqrt(2)*
           sa^2 *sp *(sa + sp)^2 *(sa *(-2 + u) *(1 + u)^2 -
           sp *(2 - 3 *u + u^3)) *(-(((-4 + u^2) *(sa^2 *sp *(-2 + u) *(1 + u)^2 -
           sp *(6 - 3* u + u^3) +
           sa *(-6 - 3* u + u^3 - sp^2 *(2 - 3 *u + u^3))))/(
           sa *sp* (sa + sp) *(sa *(-2 + u) *(1 + u)^2 - sp *(2 - 3 *u + u^3)))))^(3/2))

    return num/den
end

"""
    djac_and_ds(u,sp,sa,L)

Gradient the non-dimensional Jacobian of the transformation from angles `θ` to effective anomaly `u`.
Partial derivatives are in `E` and `L` variables.
Also compute the `L`-derivative of the `s`-bounds of motions.
Return a tuple `(djac/dE, djac/dL, ds/dL)`.

# Arguments
- `u::Float64`: Effective anomaly.
- `sp::Float64`: Pericenter in the variable `s`.
- `sa::Float64`: Apocenter in the variable `s`.
- `L::Float64`: Angular momentum parameter.
"""
function djac_and_ds(u::Float64, sp::Float64, sa::Float64, L::Float64)
    dPHIdsp, dPHIdLp = dPHIeff(sp,L)
    dPHIdsa, dPHIdLa = dPHIeff(sa,L)
    dspdE = 1/dPHIdsp
    dsadE = 1/dPHIdsa
    dspdL = -dPHIdLp/dPHIdsp
    dsadL = -dPHIdLa/dPHIdsa

    djacdsp = _djacdsp(u,sp,sa)
    djacdsa = _djacdsa(u,sp,sa)

    dsdsp = (1-_fAnomaly(u))/2
    dsdsa = (1+_fAnomaly(u))/2

    djacdE = dsadE*djacdsa + dspdE*djacdsp
    djacdL = dsadL*djacdsa + dspdL*djacdsp

    dsdL = dsadL*dsdsa+dspdL*dsdsp

    return djacdE, djacdL, dsdL
end

"""
    avr_action_coefficients(Jr,L,m_field,[nbK=nbK_default,nbAvr=nbAvr_default,m_test=m_field])

Orbit-averaged diffusion coefficients in action space `(L,Jr)`.
Returns a 5-uple `(avrDJr, avrDL, avrDJrJr, avrDJrL, avrDLL)`.

# Arguments
- `Jr::Float64`: Radial action parameter.
- `L::Float64`: Angular momentum parameter.
- `m_field::Float64`: Mass of field star.
- `nbK::Int64`: Sampling number of Rosenbluth integrals.
- `nbAvr::Int64`: Sampling number of orbit average.
- `m_test::Float64`: Mass of test star.
"""
function avr_action_coefficients(Jr::Float64, L::Float64, m_field::Float64,
                                nbK::Int64=nbK_default, nbAvr::Int64=nbAvr_default,
                                m_test::Float64=m_field)

    E = _E_from_Jr(Jr,L,nbAvr)
    sp, sa = radius_s_bounds_action(Jr,L,nbAvr)

    a = sma_eff(sp,sa)
    ecc = ecc_eff(sp,sa)

    avrDE = 0.0
    avrDL = 0.0
    avrDEE = 0.0
    avrDEL = 0.0
    avrDLL = 0.0
    halfperiod = 0.0

    dJrdE = 0.0
    dJrdL = 0.0
    d2JrdE2 = 0.0
    d2JrdEL = 0.0
    d2JrdL2 = 0.0

    for iu=1:nbAvr
        uloc = -1+2*(iu-0.5)/nbAvr
        rloc = r_eff_anomaly(uloc,a,ecc)
        trloc = _tr(rloc)
        sloc = sqrt(trloc^2+1)

        jac_loc = dth1du(uloc,sp,sa)
        djacdE_loc, djacdL_loc, dsdL_loc = djac_and_ds(uloc,sp,sa,L)

        dEloc, dE2loc, dLloc, dEdLloc, dL2loc = localOrbitChange(rloc,E,L,m_field,nbK,m_test)

        avrDE += jac_loc*dEloc
        avrDL += jac_loc*dLloc
        avrDEE += jac_loc*dE2loc
        avrDEL += jac_loc*dEdLloc
        avrDLL += jac_loc*dL2loc
        halfperiod += jac_loc

        dJrdE += jac_loc
        dJrdL += jac_loc/rloc^2
        d2JrdE2 += djacdE_loc
        d2JrdEL += djacdL_loc
        d2JrdL2 += (jac_loc/trloc^2+L/trloc^4*
                (djacdL_loc*trloc^2-2*jac_loc*sloc*dsdL_loc))

    end

    avrDE /= halfperiod
    avrDL /= halfperiod
    avrDEE /= halfperiod
    avrDEL /= halfperiod
    avrDLL /= halfperiod

    dJrdE *= (2/nbAvr)*(1/_Omega0)/PI
    dJrdL *= (2/nbAvr)*(-L/PI)*(1/_Omega0)
    d2JrdE2 *= (2/nbAvr)*(1/_Omega0)/PI
    d2JrdEL *= (2/nbAvr)*(1/_Omega0)/PI
    d2JrdL2 *= (2/nbAvr)*(-1/PI)/_L0

    avrDJr = dJrdE*avrDE+dJrdL*avrDL+(1/2)*d2JrdE2*avrDEE+(1/2)*d2JrdL2*avrDLL+d2JrdEL*avrDEL
    avrDJrJr = dJrdE^2*avrDEE+dJrdL^2*avrDLL+2*dJrdE*dJrdL*avrDEL
    avrDJrL = dJrdE*avrDEL+dJrdL*avrDLL

    return avrDJr, avrDL, avrDJrJr, avrDJrL, avrDLL
end

"""
    _Jr(E,L,[nbAvr=nbAvr_default])

Computes the radial action of an orbit with energy `E` and angular momentum `L`.

# Remarks
- Can be seen as an orbit average of the square radial velocity.

# Arguments
- `E::Float64`: Energy parameter.
- `L::Float64`: Angular momentum parameter.
- `nbAvr::Int64`: Sampling number of orbit average.
"""
function _Jr(E::Float64, L::Float64, nbAvr::Int64=nbAvr_default)
    Jr = 0.0

    sp, sa = radius_s_bounds(E,L)
    a = sma_eff(sp,sa)
    ecc = ecc_eff(sp,sa)

    for iu=1:nbAvr
        uloc = -1+2*(iu-0.5)/nbAvr
        rloc = r_eff_anomaly(uloc,a,ecc)
        jac_loc = dth1du(uloc,sp,sa)
        vrSqloc = radialVelocitySq(rloc,E,L)

        Jr += jac_loc*vrSqloc
    end

    Jr *= (2/nbAvr)*(1/_Omega0)/PI
    return Jr
end

"""
    _E_from_Jr(Jr,L,[nbAvr=nbAvr_default,eps=4.0*10^(-15)])

Computes the energy of an orbit with radial action `Jr` and angular momentum `L`.

# Remarks
- Inverse bijection of the function `E -> _Jr(E,L)` at fixed `L`.
- Computed by a modified dichotomy.

# Arguments
- `Jr::Float64`: Radial action parameter.
- `L::Float64`: Angular momentum parameter.
- `nbAvr::Tnt64`: Sampling number of orbit average.
- `eps::Tnt64`: Required precision for the dichotomy.
"""
function _E_from_Jr(Jr::Float64, L::Float64, nbAvr::Int64=nbAvr_default, eps::Float64=4.0*10^(-15))
    En = -0.5 # E0
    delta = 0.5 # initial precision
    JrE = 0.0
    if (Jr==0.0 && L==0.0)
        return -1.0
    else
        while (delta>eps) # while E is not a solution precision enough
            delta /= 2
            if (L > Lc(_tE(En))) # if (En,L) is not an orbit, we move towards the accepted (E,L) region
                En += delta
            else # if within the accepted (E,L) region, use the fact that Jr(E,L) increases with E
                JrE = _Jr(En,L,nbAvr)
                if (JrE < Jr)
                    En += delta
                else
                    En -= delta
                end
            end
        end
        return En
    end
end

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
module spt

using ArgCheck
using Clustering
using CSV
using DataFrames
using Distributions
using LinearAlgebra
using LsqFit
using PyPlot
using Statistics
using StatsBase
using Glob
using HypergeometricFunctions
# using JuMP
using Random
using QHull
using DataFramesMeta
using Base.Threads

import Random: AbstractRNG, GLOBAL_RNG
import StatsBase: sem
import Distributions:
    quantile, minimum, maximum, pdf, cdf, rand, logpdf, @check_args, @distr_support
import MLJBase: int

# exportしたい関数一覧
export
    # angle.jl
    moving_angle,
    spiff,
    anisotropy_delta_t,
    anisotropy_mean_displacement,
    # displacement.jl
    Diffusion,
    suffstats,
    fit_mle,
    # msd.jl
    # mean_square_disaplcement,
    # plot_msd,
    # messages.jl
    forward!,
    backward!,
    posterior!,
    update_ξ!,
    # mle.jl
    fit_baumwelch,
    # example_mle,
    # likelihoods.jl
    likelihood!,
    loglikelihood!,
    # utilities.jl
    create_prior,
    preproccsing!,
    add_noise!,
    # data2matrix,
    nomapround,
    label_mean_displacement!,
    group_files,
    # viterbi.jl
    viterbi,
    viterbi!,
    # velocity_autocorrelation.jl
    vacf,
    # ctrw.jl
    ctrw,
    # dynamical_functional.jl
    dynamical_functional,
    # convex_hull.jl
    convex_hull,
    # gyration.jl
    gyration,
    # mme.jl
    # mme,
    # moment,
    # displacement,
    # squared_displacement,
    # ergodicity_test.jl
    dynamical_functional,
    ergodicity_estimator,
    # features.jl
    asymmetry

# 読み込みたいjlファイル
# include("angle.jl")
include("displacement.jl")
include("angle.jl")
include("mle.jl")
include("mle_init.jl")
# include("msd.jl")
include("ctrw.jl")
include("messages.jl")
# include("mme.jl")
include("viterbi.jl")
include("likelihoods.jl")
include("utilities.jl")
include("velocity_autocorrelation.jl")
include("preprocessing.jl")
include("ergodicity_test.jl")
include("convex_hull.jl")
include("gyration.jl")
include("features.jl")

# To be removed in a future version
# ---------------------------------
# export
# n_parameters,

# @deprecate n_parameters(hmm) nparams(hmm)

end # module

ext: jl

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
module TestTransducers

import LoadAllPackages
import PerformanceTestTools
using Distributed: addprocs, nworkers
using Test

if get(ENV, "CI", "false") == "true"
    addprocs(3)

    # Tests in `PerformanceTestTools.@include_foreach` might cause
    # pre-compilation errors as two processes try to compile packages
    # at the same time.  This can happen when the tests are run via
    # `Pkg.test`.  Doing this after `addprocs` to workaround a quirk
    # in Distributed.jl.
    LoadAllPackages.loadall()
end
@info "Testing with:" nworkers()

@testset "$file" for file in sort([file for file in readdir(@__DIR__) if
                                   match(r"^test_.*\.jl$", file) !== nothing])

    if file == "test_doctest.jl"
        if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
            @info "Skipping doctests on PkgEval."
            continue
        elseif !(v"1.6" <= VERSION < v"1.7-")
            @info "Skipping doctests on Julia $VERSION."
            continue
        end
    end

    VERSION < v"1.3" && file == "test_nondeterministic_threading.jl" && continue

    # Do not support `Broadcasting` in Julia 1.0
    VERSION < v"1.1" && file == "test_broadcasting.jl" && continue

    # BlockArrays.jl 0.11 only support Julia 1.1 and above.
    VERSION < v"1.1" && file == "test_blockarrays.jl" && continue

    include(file)
end

PerformanceTestTools.@include_foreach(
    "threads/runtests.jl",
    [nothing, ["JULIA_NUM_THREADS" => Threads.nthreads() > 1 ? "1" : "2"]],
    parallel = true,
)

end  # module


text:
@testset "exceptional $T" for T in (Float32, Float64)

@testset "exceptional $xatan" for xatan in (SLEEF.atan_fast, SLEEF.atan)

    @test xatan(T(0.0),  T(-0.0)) === T(pi)
    @test xatan(T(-0.0), T(-0.0)) === -T(pi)
    @test ispzero(xatan(T(0.0), T(0.0)))
    @test isnzero(xatan(T(-0.0), T(0.0)))
    @test xatan( T(Inf), -T(Inf)) === T(3*pi/4)
    @test xatan(-T(Inf), -T(Inf)) === T(-3*pi/4)
    @test xatan( T(Inf),  T(Inf))  === T(pi/4)
    @test xatan(-T(Inf),  T(Inf))  === T(-pi/4)

    y = T(0.0)
    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5]
    for x in xa
        @test xatan(y,x) === T(pi)
    end

    y = T(-0.0)
    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5]
    for x in xa
        @test xatan(y,x) === T(-pi)
    end

    ya = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5]
    xa = T[T(0.0), T(-0.0)]
    for x in xa, y in ya
        @test xatan(y,x) === T(-pi/2)
    end

    ya = T[100000.5, 100000, 3, 2.5, 2, 1.5, 1.0, 0.5]
    xa = T[T(0.0), T(-0.0)]
    for x in xa, y in ya
        @test xatan(y,x) === T(pi/2)
    end

    y = T(Inf)
    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5, -0.0, +0.0, 0.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    for x in xa
        @test xatan(y,x) === T(pi/2)
    end

    y = T(-Inf)
    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5, -0.0, +0.0, 0.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    for x in xa
        @test xatan(y,x) === T(-pi/2)
    end

    ya = T[0.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    x = T(Inf)
    for y in ya
        @test ispzero(xatan(y,x))
    end

    ya = T[-0.5, -1.5, -2.0, -2.5, -3.0, -100000, -100000.5]
    x = T(Inf)
    for y in ya
        @test isnzero(xatan(y,x))
    end

    ya = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5, -0.0, +0.0, 0.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5, NaN]
    x = T(NaN)
    for y in ya
        @test isnan(xatan(y,x))
    end

    y = T(NaN)
    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5, -0.0, +0.0, 0.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5, NaN]
    for x in xa
        @test isnan(xatan(y,x))
    end

end # denormal/nonumber atan

@testset "exceptional xpow" begin

    @test SLEEF.pow(T(1), T(NaN))   === T(1)
    @test SLEEF.pow( T(NaN), T(0))  === T(1)
    @test SLEEF.pow(-T(1), T(Inf))  === T(1)
    @test SLEEF.pow(-T(1), T(-Inf)) === T(1)

    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5]
    ya = T[-100000.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 100000.5]
    for x in xa, y in ya
        @test isnan(SLEEF.pow(x,y))
    end

    x = T(NaN)
    ya = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    for y in ya
        @test isnan(SLEEF.pow(x,y))
    end

    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5, -0.0, +0.0, 0.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    y = T(NaN)
    for x in xa
        @test isnan(SLEEF.pow(x,y))
    end

    x = T(0.0)
    ya = T[1, 3, 5, 7, 100001]
    for y in ya
        @test ispzero(SLEEF.pow(x,y))
    end

    x = T(-0.0)
    ya = T[1, 3, 5, 7, 100001]
    for y in ya
        @test isnzero(SLEEF.pow(x,y))
    end

    xa = T[0.0, -0.0]
    ya = T[0.5, 1.5, 2.0, 2.5, 4.0, 100000, 100000.5]
    for x in xa, y in ya
        @test ispzero(SLEEF.pow(x,y))
    end

    xa = T[-0.999, -0.5, -0.0, +0.0, +0.5, +0.999]
    y = T(-Inf)
    for x in xa
        @test SLEEF.pow(x,y) === T(Inf)
    end

    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    y = T(-Inf)
    for x in xa
        @test ispzero(SLEEF.pow(x,y))
    end

    xa = T[-0.999, -0.5, -0.0, +0.0, +0.5, +0.999]
    y = T(Inf)
    for x in xa
        @test ispzero(SLEEF.pow(x,y))
    end

    xa = T[-100000.5, -100000, -3, -2.5, -2, -1.5, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    y = T(Inf)
    for x in xa
        @test SLEEF.pow(x,y) === T(Inf)
    end

    x = T(-Inf)
    ya = T[-100001, -5, -3, -1]
    for y in ya
        @test isnzero(SLEEF.pow(x,y))
    end

    x = T(-Inf)
    ya = T[-100000.5, -100000, -4, -2.5, -2, -1.5, -0.5]
    for y in ya
        @test ispzero(SLEEF.pow(x,y))
    end

    x = T(-Inf)
    ya = T[1, 3, 5, 7, 100001]
    for y in ya
        @test SLEEF.pow(x,y) === T(-Inf)
    end

    x = T(-Inf)
    ya = T[0.5, 1.5, 2, 2.5, 3.5, 4, 100000, 100000.5]
    for y in ya
        @test SLEEF.pow(x,y) === T(Inf)
    end

    x = T(Inf)
    ya = T[-100000.5, -100000, -3, -2.5, -2, -1.5, -1.0, -0.5]
    for y in ya
        @test ispzero(SLEEF.pow(x,y))
    end

    x = T(Inf)
    ya = T[0.5, 1, 1.5, 2.0, 2.5, 3.0, 100000, 100000.5]
    for y in ya
        @test SLEEF.pow(x,y) === T(Inf)
    end

    x = T(0.0)
    ya = T[-100001, -5, -3, -1]
    for y in ya
        @test SLEEF.pow(x,y) ===  T(Inf)
    end

    x = T(-0.0)
    ya = T[-100001, -5, -3, -1]
    for y in ya
        @test SLEEF.pow(x,y) ===  T(-Inf)
    end

    xa = T[0.0, -0.0]
    ya = T[-100000.5, -100000, -4, -2.5, -2, -1.5, -0.5]
    for x in xa, y in ya
        @test SLEEF.pow(x,y) === T(Inf)
    end

    xa = T[1000, -1000]
    ya = T[1000, 1000.5, 1001]
    for x in xa, y in ya
        @test cmpdenorm(SLEEF.pow(x,y), Base.:^(BigFloat(x),BigFloat(y)))
    end

end # denormal/nonumber pow

fun_table = Dict(SLEEF.sin_fast => Base.sin, SLEEF.sin => Base.sin)
@testset "exceptional $xtrig" for (xtrig, trig) in fun_table
    xa = T[NaN, -0.0, 0.0, Inf, -Inf]
    for x in xa
        @test cmpdenorm(xtrig(x), trig(BigFloat(x)))
    end
end

fun_table = Dict(SLEEF.cos_fast => Base.cos, SLEEF.cos => Base.cos)
@testset "exceptional $xtrig" for (xtrig, trig) in fun_table
    xa = T[NaN, -0.0, 0.0, Inf, -Inf]
    for x in xa
        @test cmpdenorm(xtrig(x), trig(BigFloat(x)))
    end
end

@testset "exceptional sin in $xsincos"for xsincos in (SLEEF.sincos_fast, SLEEF.sincos)
    xa = T[NaN, -0.0, 0.0, Inf, -Inf]
    for x in xa
        q = xsincos(x)[1]
        @test cmpdenorm(q, Base.sin(BigFloat(x)))
    end
end

@testset "exceptional cos in $xsincos"for xsincos in (SLEEF.sincos_fast, SLEEF.sincos)
    xa = T[NaN, -0.0, 0.0, Inf, -Inf]
    for x in xa
        q = xsincos(x)[2]
        @test cmpdenorm(q, Base.cos(BigFloat(x)))
    end
end

@testset "exceptional $xtan" for xtan in (SLEEF.tan_fast, SLEEF.tan)
    xa = T[NaN, Inf, -Inf, -0.0, 0.0, pi/2, -pi/2]
    for x in xa
        @test cmpdenorm(xtan(x), Base.tan(BigFloat(x)))
    end
end

fun_table = Dict(SLEEF.asin => Base.asin, SLEEF.asin_fast => Base.asin, SLEEF.acos => Base.acos, SLEEF.acos_fast => Base.acos)
@testset "exceptional $xatrig" for (xatrig, atrig) in fun_table
    xa = T[NaN, Inf, -Inf, 2, -2, 1, -1, -0.0, 0.0]
    for x in xa
        @test cmpdenorm(xatrig(x), atrig(BigFloat(x)))
    end
end

@testset "exceptional $xatan" for xatan in (SLEEF.atan, SLEEF.atan_fast)
    xa = T[NaN, Inf, -Inf, -0.0, 0.0]
    for x in xa
        @test cmpdenorm(xatan(x), Base.atan(BigFloat(x)))
    end
end

@testset "exceptional exp" begin
    xa = T[NaN, Inf, -Inf, 10000, -10000]
    for x in xa
        @test cmpdenorm(SLEEF.exp(x), Base.exp(BigFloat(x)))
    end
end

@testset "exceptional sinh" begin
    xa = T[NaN, 0.0, -0.0, Inf, -Inf, 10000, -10000]
    for x in xa
        @test cmpdenorm(SLEEF.sinh(x), Base.sinh(BigFloat(x)))
    end
end

@testset "exceptional cosh" begin
    xa = T[NaN, 0.0, -0.0, Inf, -Inf, 10000, -10000]
    for x in xa
        @test cmpdenorm(SLEEF.cosh(x),  Base.cosh(BigFloat(x)))
    end
end

@testset "exceptional tanh" begin
    xa = T[NaN, 0.0, -0.0, Inf, -Inf, 10000, -10000]
    for x in xa
        @test cmpdenorm(SLEEF.tanh(x),  Base.tanh(BigFloat(x)))
    end
end

@testset "exceptional asinh" begin
    xa = T[NaN, 0.0, -0.0, Inf, -Inf, 10000, -10000]
    for x in xa
        @test cmpdenorm(SLEEF.asinh(x), Base.asinh(BigFloat(x)))
    end
end

@testset "exceptional acosh" begin
    xa = T[NaN, 0.0, -0.0, 1.0, Inf, -Inf, 10000, -10000]
    for x in xa
        @test cmpdenorm(SLEEF.acosh(x), Base.acosh(BigFloat(x)))
    end
end

@testset "exceptional atanh" begin
    xa = T[NaN, 0.0, -0.0, 1.0, -1.0, Inf, -Inf, 10000, -10000]
    for x in xa
        @test cmpdenorm(SLEEF.atanh(x), Base.atanh(BigFloat(x)))
    end
end

@testset "exceptional $xcbrt" for xcbrt = (SLEEF.cbrt, SLEEF.cbrt_fast)
    xa = T[NaN, Inf, -Inf, 0.0, -0.0]
    for x in xa
        @test cmpdenorm(SLEEF.cbrt(x), Base.cbrt(BigFloat(x)))
    end
end

@testset "exceptional exp2" begin
    xa = T[NaN, Inf, -Inf]
    for x in xa
        @test cmpdenorm(SLEEF.exp2(x), Base.exp2(BigFloat(x)))
    end
end

@testset "exceptional exp10" begin
    xa = T[NaN, Inf, -Inf]
    for x in xa
        @test cmpdenorm(SLEEF.exp10(x), Base.exp10(BigFloat(x)))
    end
end

@testset "exceptional expm1" begin
    xa = T[NaN, Inf, -Inf, 0.0, -0.0]
    for x in xa
        @test cmpdenorm(SLEEF.expm1(x), Base.expm1(BigFloat(x)))
    end
end

@testset "exceptional $xlog" for xlog in (SLEEF.log, SLEEF.log_fast)
    xa = T[NaN, Inf, -Inf, 0, -1]
    for x in xa
        @test cmpdenorm(xlog(x), Base.log(BigFloat(x)))
    end
end

@testset "exceptional log10" begin
    xa = T[NaN, Inf, -Inf, 0, -1]
    for x in xa
        @test cmpdenorm(SLEEF.log10(x), Base.log10(BigFloat(x)))
    end
end

@testset "exceptional log2" begin
    xa = T[NaN, Inf, -Inf, 0, -1]
    for x in xa
        @test cmpdenorm(SLEEF.log2(x), Base.log2(BigFloat(x)))
    end
end

@testset "exceptional log1p" begin
    xa = T[NaN, Inf, -Inf, 0.0, -0.0, -1.0, -2.0]
    for x in xa
        @test cmpdenorm(SLEEF.log1p(x), Base.log1p(BigFloat(x)))
    end
end

@testset "exceptional ldexp" begin
    for i = -10000:10000
        a = SLEEF.ldexp(T(1.0), i)
        b = Base.ldexp(BigFloat(1.0), i)
        @test (isfinite(b) && a == b || cmpdenorm(a,b))
    end
end

@testset "exceptional ilogb" begin
    @test SLEEF.ilogb(+T(Inf)) == SLEEF.INT_MAX
    @test SLEEF.ilogb(-T(Inf)) == SLEEF.INT_MAX
    @test SLEEF.ilogb(+T(0.0)) == SLEEF.FP_ILOGB0
    @test SLEEF.ilogb(-T(0.0)) == SLEEF.FP_ILOGB0
    @test SLEEF.ilogb( T(NaN)) == SLEEF.FP_ILOGBNAN
end

end #exceptional


text:
using StaticArrays

@test begin
    x0 = SVector(0., 0.)
    x1 = SVector(1., 1.)

    curve = Curve([x0, x1])
    is_consistent(curve)
end

@test begin
    x0 = SVector(0., 0.)
    x1 = SVector(1., 1.)
    x2 = SVector(1., 1.)

    curve = Curve([x0, x1, x2])
    ! is_consistent(curve)
end

@test begin
    x0 = SVector(0., 0.)
    x1 = SVector(0., 0.)

    curve = Curve([x0, x1])
    !is_consistent(curve)
end

@test begin
    x0 = SVector(0., 0.)
    x1 = SVector(1., 1.)
    x2 = SVector(2., 2.)

    curve = Curve([x0, x1, x2], [(1, 2), (2, 3), (3, 1)])
    is_consistent(curve)
end

@test begin
    x0 = SVector(0., 0.)
    x1 = SVector(1., 1.)
    x2 = SVector(2., 2.)

    curve = Curve([x0, x1, x2], [(1, 1), (2, 3), (3, 1)])
    !is_consistent(curve)
end

@test begin
    x0 = SVector(0., 0.)
    x1 = SVector(1., 1.)
    x2 = SVector(2., 2.)

    curve = Curve([x0, x1, x2], [(1, 9), (2, 3), (3, 1)])
    !is_consistent(curve)
end

# IO
@test begin
    ext = ".vtu"
    all(length(Curve(vtu_path*ext).points) > 0 for vtu_path in ("mesh2d", "mesh2din3d", "tet_mesh"))
end

@test begin
    ext = ".xml"
    all(length(Curve(vtu_path*ext).points) > 0 for vtu_path in ("mesh2d", "mesh2din3d", "tet_mesh"))
end

@test begin
    ext = ".msh"
    all(length(Curve(vtu_path*ext).points) > 0 for vtu_path in ("mesh2d", "mesh2din3d", "tet_mesh"))
end


text:
"""
    mincut(flow_graph::Graphs.IsDirected, source::Integer, target::Integer, capacity_matrix::AbstractMatrix, algorithm::AbstractFlowAlgorithm)

Compute the min-cut between `source` and `target` for the given graph.
First computes the maxflow using `algorithm` and then builds the partition of the residual graph
Returns a triplet `(part1, part2, flow)` with the partition containing the source, the partition containing the target (the rest) and the min-cut(max-flow) value
"""
function mincut(
        flow_graph::Graphs.DiGraph,             # the input graph
        source::Integer,                       # the source vertex
        target::Integer,                       # the target vertex
        capacity_matrix::AbstractMatrix,       # edge flow capacities
        algorithm::AbstractFlowAlgorithm       # keyword argument for algorithm
    )
    flow, flow_matrix = maximum_flow(flow_graph, source, target, capacity_matrix, algorithm)
    residual_matrix = spzeros(Graphs.nv(flow_graph),Graphs.nv(flow_graph))
    for edge in Graphs.edges(flow_graph)
        residual_matrix[edge.src,edge.dst] = max(0.0, capacity_matrix[edge.src,edge.dst] - flow_matrix[edge.src,edge.dst])
        residual_matrix[edge.dst,edge.src] = max(0.0, capacity_matrix[edge.dst,edge.src] - flow_matrix[edge.dst,edge.src])
    end
    part1 = typeof(source)[]
    queue = [source]
    while !isempty(queue)
        node = pop!(queue)
        push!(part1, node)
        dests = [dst for dst in 1:Graphs.nv(flow_graph) if residual_matrix[node,dst]>0.0 && dst ∉ part1 && dst ∉ queue]
        append!(queue, dests)
    end
    part2 = [node for node in 1:Graphs.nv(flow_graph) if node ∉ part1]
    return (part1, part2, flow)
end


text:
module PortAudio

using libportaudio_jll, SampledSignals

import Base: eltype, show
import Base: close, isopen
import Base: read, read!, write, flush

import LinearAlgebra
import LinearAlgebra: transpose!

export PortAudioStream

include("libportaudio.jl")

# This size is in frames

# data is passed to and from portaudio in chunks with this many frames, because
# we need to interleave the samples
const CHUNKSIZE=128

function versioninfo(io::IO=stdout)
    println(io, Pa_GetVersionText())
    println(io, "Version: ", Pa_GetVersion())
end

mutable struct PortAudioDevice
    name::String
    hostapi::String
    maxinchans::Int
    maxoutchans::Int
    defaultsamplerate::Float64
    idx::PaDeviceIndex
    lowinputlatency::Float64
    lowoutputlatency::Float64
    highinputlatency::Float64
    highoutputlatency::Float64
end

PortAudioDevice(info::PaDeviceInfo, idx) = PortAudioDevice(
        unsafe_string(info.name),
        unsafe_string(Pa_GetHostApiInfo(info.host_api).name),
        info.max_input_channels,
        info.max_output_channels,
        info.default_sample_rate,
        idx,
        info.default_low_input_latency,
        info.default_low_output_latency,
        info.default_high_input_latency,
        info.default_high_output_latency)

function devices()
    ndevices = Pa_GetDeviceCount()
    infos = PaDeviceInfo[Pa_GetDeviceInfo(i) for i in 0:(ndevices - 1)]
    PortAudioDevice[PortAudioDevice(info, idx-1) for (idx, info) in enumerate(infos)]
end

# not for external use, used in error message printing
devnames() = join(["\"$(dev.name)\"" for dev in devices()], "\n")

##################
# PortAudioStream
##################

mutable struct PortAudioStream{T}
    samplerate::Float64
    latency::Float64
    stream::PaStream
    warn_xruns::Bool
    sink # untyped because of circular type definition
    source # untyped because of circular type definition

    # this inner constructor is generally called via the top-level outer
    # constructor below

    # TODO: pre-fill outbut buffer on init
    # TODO: recover from xruns - currently with low latencies (e.g. 0.01) it
    # will run fine for a while and then fail with the first xrun.
    # TODO: figure out whether we can get deterministic latency...
    function PortAudioStream{T}(indev::PortAudioDevice, outdev::PortAudioDevice,
                                inchans, outchans, sr,
                                latency, warn_xruns) where {T}
        inchans = inchans == -1 ? indev.maxinchans : inchans
        outchans = outchans == -1 ? outdev.maxoutchans : outchans
        inparams = (inchans == 0) ?
            Ptr{Pa_StreamParameters}(0) :
            Ref(Pa_StreamParameters(indev.idx, inchans, type_to_fmt[T], latency, C_NULL))
        outparams = (outchans == 0) ?
            Ptr{Pa_StreamParameters}(0) :
            Ref(Pa_StreamParameters(outdev.idx, outchans, type_to_fmt[T], latency, C_NULL))
        this = new(sr, latency, C_NULL, warn_xruns)
        # finalizer(close, this)
        this.sink = PortAudioSink{T}(outdev.name, this, outchans)
        this.source = PortAudioSource{T}(indev.name, this, inchans)
        this.stream = suppress_err() do
            Pa_OpenStream(inparams, outparams, sr, 0, paNoFlag,
                          nothing, nothing)
        end

        Pa_StartStream(this.stream)

        this
    end
end

defaultlatency(devices...) = maximum(d -> max(d.highoutputlatency, d.highinputlatency), devices)

# this is the top-level outer constructor that all the other outer constructors end up calling
"""
    PortAudioStream(inchannels=2, outchannels=2; options...)
    PortAudioStream(duplexdevice, inchannels=2, outchannels=2; options...)
    PortAudioStream(indevice, outdevice, inchannels=2, outchannels=2; options...)

Audio devices can either be `PortAudioDevice` instances as returned
by `PortAudio.devices()`, or strings with the device name as reported by the
operating system. If a single `duplexdevice` is given it will be used for both
input and output. If no devices are given the system default devices will be
used.

Options:

* `eltype`:       Sample type of the audio stream (defaults to Float32)
* `samplerate`:   Sample rate (defaults to device sample rate)
* `latency`:      Requested latency. Stream could underrun when too low,
                  consider using provided device defaults
* `warn_xruns`:   Display a warning if there is a stream overrun or underrun,
                  which often happens when Julia is compiling, or with a
                  particularly large GC run. This can be quite verbose so is
                  false by default.
"""
function PortAudioStream(indev::PortAudioDevice, outdev::PortAudioDevice,
        inchans=2, outchans=2; eltype=Float32, samplerate=-1,
        latency=defaultlatency(indev, outdev), warn_xruns=false)
    if samplerate == -1
        sampleratein = indev.defaultsamplerate
        samplerateout = outdev.defaultsamplerate
        if inchans > 0 && outchans > 0 && sampleratein != samplerateout
            error("""
            Can't open duplex stream with mismatched samplerates (in: $sampleratein, out: $samplerateout).
                   Try changing your sample rate in your driver settings or open separate input and output
                   streams""")
        elseif inchans > 0
            samplerate = sampleratein
        else
            samplerate = samplerateout
        end
    end
    PortAudioStream{eltype}(indev, outdev, inchans, outchans, samplerate, latency, warn_xruns)
end

# handle device names given as streams
function PortAudioStream(indevname::AbstractString, outdevname::AbstractString, args...; kwargs...)
    indev = nothing
    outdev = nothing
    for d in devices()
        if d.name == indevname
            indev = d
        end
        if d.name == outdevname
            outdev = d
        end
    end
    if indev == nothing
        error("No device matching \"$indevname\" found.\nAvailable Devices:\n$(devnames())")
    end
    if outdev == nothing
        error("No device matching \"$outdevname\" found.\nAvailable Devices:\n$(devnames())")
    end

    PortAudioStream(indev, outdev, args...; kwargs...)
end

# if one device is given, use it for input and output, but set inchans=0 so we
# end up with an output-only stream
function PortAudioStream(device::PortAudioDevice, inchans=2, outchans=2; kwargs...)
    PortAudioStream(device, device, inchans, outchans; kwargs...)
end
function PortAudioStream(device::AbstractString, inchans=2, outchans=2; kwargs...)
    PortAudioStream(device, device, inchans, outchans; kwargs...)
end

# use the default input and output devices
function PortAudioStream(inchans=2, outchans=2; kwargs...)
    inidx = Pa_GetDefaultInputDevice()
    indevice = PortAudioDevice(Pa_GetDeviceInfo(inidx), inidx)
    outidx = Pa_GetDefaultOutputDevice()
    outdevice = PortAudioDevice(Pa_GetDeviceInfo(outidx), outidx)
    PortAudioStream(indevice, outdevice, inchans, outchans; kwargs...)
end

# handle do-syntax
function PortAudioStream(fn::Function, args...; kwargs...)
    str = PortAudioStream(args...; kwargs...)
    try
        fn(str)
    finally
        close(str)
    end
end

function close(stream::PortAudioStream)
    if stream.stream != C_NULL
        Pa_StopStream(stream.stream)
        Pa_CloseStream(stream.stream)
        stream.stream = C_NULL
    end

    nothing
end

isopen(stream::PortAudioStream) = stream.stream != C_NULL

SampledSignals.samplerate(stream::PortAudioStream) = stream.samplerate
eltype(stream::PortAudioStream{T}) where T = T

read(stream::PortAudioStream, args...) = read(stream.source, args...)
read!(stream::PortAudioStream, args...) = read!(stream.source, args...)
write(stream::PortAudioStream, args...) = write(stream.sink, args...)
write(sink::PortAudioStream, source::PortAudioStream, args...) = write(sink.sink, source.source, args...)
flush(stream::PortAudioStream) = flush(stream.sink)

function show(io::IO, stream::PortAudioStream)
    println(io, typeof(stream))
    println(io, "  Samplerate: ", samplerate(stream), "Hz")
    if nchannels(stream.sink) > 0
        print(io, "\n  ", nchannels(stream.sink), " channel sink: \"", name(stream.sink), "\"")
    end
    if nchannels(stream.source) > 0
        print(io, "\n  ", nchannels(stream.source), " channel source: \"", name(stream.source), "\"")
    end
end

##################################
# PortAudioSink & PortAudioSource
##################################

# Define our source and sink types
for (TypeName, Super) in ((:PortAudioSink, :SampleSink),
                          (:PortAudioSource, :SampleSource))
    @eval mutable struct $TypeName{T} <: $Super
        name::String
        stream::PortAudioStream{T}
        chunkbuf::Array{T, 2}
        nchannels::Int

        function $TypeName{T}(name, stream, channels) where {T}
            # portaudio data comes in interleaved, so we'll end up transposing
            # it back and forth to julia column-major
            chunkbuf = zeros(T, channels, CHUNKSIZE)
            new(name, stream, chunkbuf, channels)
        end
    end
end

SampledSignals.nchannels(s::Union{PortAudioSink, PortAudioSource}) = s.nchannels
SampledSignals.samplerate(s::Union{PortAudioSink, PortAudioSource}) = samplerate(s.stream)
eltype(::Union{PortAudioSink{T}, PortAudioSource{T}}) where {T} = T
function close(s::Union{PortAudioSink, PortAudioSource})
    throw(ErrorException("Attempted to close PortAudioSink or PortAudioSource.
                          Close the containing PortAudioStream instead"))
end
isopen(s::Union{PortAudioSink, PortAudioSource}) = isopen(s.stream)
name(s::Union{PortAudioSink, PortAudioSource}) = s.name

function show(io::IO, ::Type{PortAudioSink{T}}) where T
    print(io, "PortAudioSink{$T}")
end

function show(io::IO, ::Type{PortAudioSource{T}}) where T
    print(io, "PortAudioSource{$T}")
end

function show(io::IO, stream::T) where {T <: Union{PortAudioSink, PortAudioSource}}
    print(io, nchannels(stream), "-channel ", T, "(\"", stream.name, "\")")
end

function SampledSignals.unsafe_write(sink::PortAudioSink, buf::Array, frameoffset, framecount)
    nwritten = 0
    while nwritten < framecount
        n = min(framecount-nwritten, CHUNKSIZE)
        # make a buffer of interleaved samples
        transpose!(view(sink.chunkbuf, :, 1:n),
                   view(buf, (1:n) .+ nwritten .+ frameoffset, :))
        # TODO: if the stream is closed we just want to return a
        # shorter-than-requested frame count instead of throwing an error
        Pa_WriteStream(sink.stream.stream, sink.chunkbuf, n, sink.stream.warn_xruns)
        nwritten += n
    end

    nwritten
end

function SampledSignals.unsafe_read!(source::PortAudioSource, buf::Array, frameoffset, framecount)
    nread = 0
    while nread < framecount
        n = min(framecount-nread, CHUNKSIZE)
        # TODO: if the stream is closed we just want to return a
        # shorter-than-requested frame count instead of throwing an error
        Pa_ReadStream(source.stream.stream, source.chunkbuf, n, source.stream.warn_xruns)
        # de-interleave the samples
        transpose!(view(buf, (1:n) .+ nread .+ frameoffset, :),
                   view(source.chunkbuf, :, 1:n))

        nread += n
    end

    nread
end

function suppress_err(dofunc::Function)
    nullfile = @static Sys.iswindows() ? "nul" : "/dev/null"
    open(nullfile, "w") do io
        redirect_stderr(dofunc, io)
    end
end

function __init__()
    if Sys.islinux()
        envkey = "ALSA_CONFIG_DIR"
        if envkey ∉ keys(ENV)
            searchdirs = ["/usr/share/alsa",
                          "/usr/local/share/alsa",
                          "/etc/alsa"]
            confdir_idx = findfirst(searchdirs) do d
                isfile(joinpath(d, "alsa.conf"))
            end
            if confdir_idx === nothing
                throw(ErrorException(
                    """
                    Could not find ALSA config directory. Searched:
                    $(join(searchdirs, "\n"))

                    if ALSA is installed, set the "ALSA_CONFIG_DIR" environment
                    variable. The given directory should have a file "alsa.conf".

                    If it would be useful to others, please file an issue at
                    https://github.com/JuliaAudio/PortAudio.jl/issues
                    with your alsa config directory so we can add it to the search
                    paths.
                    """))
            end
            confdir = searchdirs[confdir_idx]
            ENV[envkey] = confdir
        end
    end
    # initialize PortAudio on module load. libportaudio prints a bunch of
    # junk to STDOUT on initialization, so we swallow it.
    # TODO: actually check the junk to make sure there's nothing in there we
    # don't expect
    suppress_err() do
        Pa_Initialize()
    end

    atexit() do
        Pa_Terminate()
    end
end

end # module PortAudio


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Roots
import DataFrames
import OpticSim.Repeat

#Really should extend AbstractLattice to support cosets, then wouldn't need this ad hoc stuff. TODO
colorbasis(::Repeat.HexBasis1) = SMatrix{2,2}(2, -1, 1, 1)
colorbasis(::Repeat.HexBasis3) = SMatrix{2,2}(2, -1, 1, 1)
colororigins(::Repeat.HexBasis1) = ((0, 0), (-1, 0), (-1, 1))
colororigins(::Repeat.HexBasis3) = ((0, 0), (0, -1), (1, -1))

""" For lenslets arranged in a hexagonal pattern cross talk between lenslets can be reduced by arranging the color of the lenslets so that each lenslet is surrounded only by lenslets of a different color. This function computes the color to assign to any lattice point to ensure this property holds."""
function pointcolor(point, cluster::Repeat.AbstractLatticeCluster)
    latticematrix = colorbasis(Repeat.elementbasis(cluster))
    origins = colororigins(Repeat.elementbasis(cluster))
    colors = zip(origins, ("red", "green", "blue"))
    for (origin, color) in colors
        if nothing !== Repeat.latticepoint(latticematrix, origin, point)
            return color
        end
    end
end
export pointcolor

defaultclusterproperties() = (mtf = .2, minfnumber = 2.0, cyclesperdegree = 11, λ = 530nm, pixelpitch = .9μm)
export defaultclusterproperties

"""maximum allowable display for visibility reasons"""
maxlensletdisplaysize() = (250μm, 250μm)

const ρ_quartervalue = 2.21509 # value of ρ at which the airy disk function has magnitude .25
export ρ_quartervalue
const ρ_zerovalue = 3.832 # value of ρ at which the airy disk function has magnitude 0

"""given pixelpitch and angular subtense (in degrees) of pixel returns focal length"""
function compute_focal_length(fov,lensletdiameter, minfnumber, maxdisplaysize)
    maxangle = deg2rad(max(fov...)) #leaving angles as degrees caused trouble for not obvious reason

    tempfl = uconvert(mm, maxdisplaysize/(2*tan(maxangle/2.0)))
    fnum = tempfl/lensletdiameter
    return tempfl*minfnumber/fnum
end
export compute_focal_length

pixelsperdegree(focal_length,pixelpitch) = 1/(2.0*atand(uconvert(Unitful.NoUnits,pixelpitch/(2.0*focal_length))))
export pixelsperdegree

function lensletdisplaysize(fov,focal_length)
    uconvert.(μm,2*focal_length .* tand.(fov./2))
end
export lenseletdisplaysize

"""diffraction limited response for a circular aperture, normalized by maximum cutoff frequency"""
function mtfcircular(freq, freqcutoff)
    s =  freq / freqcutoff
    return 2 / π * (acos(s) - s * ((1 - s^2)^.5))
end
export mtfcircular

"""
# Returns the diffraction limit frequency in cycles/deg. At this frequeny the response of the system is zero.

## Derivation from the more common formula relating f number and diffraction limit.

focal length = 𝒇𝒍
diffraction cutoff frequency,fc, in cycles/mm = 1/λF# = diameter/λ*𝒇𝒍
cutoff wavelength, Wc, = 1/cutoff frequency = λ*𝒇𝒍/diameter

angular wavelength, θc, radians/cycle, corresponding to cutoff wavelength:

θ ≈ tanθ for small θ
θc corresponding to Wc:  θc ≈ tanθ = Wc/𝒇𝒍
Wc = θc*𝒇𝒍

from equation for Wc:

λ*𝒇𝒍/diameter = Wc = θc*𝒇𝒍
θc = λ/diameter
cycles/rad = 1/θc = diameter/λ
"""
diffractionlimit(λ::Unitful.Length,diameter::Unitful.Length) = uconvert(Unitful.NoUnits, diameter / λ) / rad2deg(1)
export diffractionlimit

"""computes the diameter of a diffraction limited lens that has response mtf at frequency cyclesperdeg"""
function diameter_for_cycles_deg(mtf, cyclesperdeg, λ::Unitful.Length)
    cyclesperrad = rad2deg(1) * cyclesperdeg
    f(x) = mtfcircular(x, 1.0) - mtf
    normalizedfrequency = find_zero(f, (0.0, 1.0))
    fc = cyclesperrad / normalizedfrequency
    return uconvert(mm, λ * fc)
end
export diameter_for_cycles_deg

airydisk(ρ) = (2 * SpecialFunctions.besselj1(ρ) / ρ)^2

"""The f# required for the first zero of the airy diffraction disk to be at the next sample point"""
diffractionfnumber(λ::Unitful.Length,pixelpitch::Unitful.Length,indexofrefraction) = uconvert(Unitful.NoUnits, pixelpitch / (2.44λ / indexofrefraction))
export diffractionfnumber

"""returns ρ, the normalized distance, at which the airy disk will have the value airyvalue"""
function ρatairyvalue(airyvalue)
    ρ = 1e-5 # start at small angle
    while airydisk(ρ) > airyvalue
        ρ += 1e-5 # numerical precision not an issue here
    end
    return ρ
end
export ρatairyvalue

""" Spacing between lenslets which guarantess that for any pixel visible in all lenslets every point on the eyebox plane is covered. This is closest packing of circles."""
closestpackingdistance(pupildiameter::Unitful.Length) = pupildiameter * cosd(30)
export closestpackingdistance

"""Tries clusters of various sizes to choose the largest one which fits within the eye pupil. Larger clusters allow for greater reduction of the fov each lenslet must cover so it returns the largest feasible cluster"""
function choosecluster(pupildiameter::Unitful.Length, lensletdiameter::Unitful.Length)
    clusters = (hex3RGB, hex4RGB, hex7RGB , hex9RGB, hex12RGB,hex19RGB) #, hex37RGB) # hex37RGB leave out for now. Leads to designs with thousands of small lenslets. May not be practical.
    pupildiameter
    ratio = 0.0
    clusterindex = 0
    scale = 0

    for clusterfunc in clusters
        cluster = clusterfunc() #create an instance of the cluster type with default unit scale
        scale = ustrip(mm,lensletdiameter)/euclideandiameter(elementbasis(cluster)) #scale the element basis of the cluster to match lenslet diameter
        scaledcluster = clusterfunc(scale) #make a new instance of the cluster type scaled so the element basis has diameter equal to lenslet diameter

        temp = ustrip(mm,pupildiameter) / (euclideandiameter(scaledcluster))

        if temp >= 1.0
            ratio = temp
            clusterindex += 1
        else
            break
        end
    end

    maxcluster = clusters[clusterindex](ratio*scale)
    @assert isapprox(Repeat.euclideandiameter(maxcluster), ustrip(mm,pupildiameter))

    scaledlensletdiameter = lensletdiameter*ratio
    return (cluster = maxcluster, lensletdiameter = scaledlensletdiameter, packingdistance = pupildiameter * ustrip(mm, lensletdiameter))
end
export choosecluster

"""Computes the largest feasible cluster size meeting constraints"""
function choosecluster(pupildiameter::Unitful.Length, λ::Unitful.Length, mtf, cyclesperdeg::T) where {T <: Real}
    diam = diameter_for_cycles_deg(mtf, cyclesperdeg, λ)
    return choosecluster(pupildiameter, diam) # use minimum diameter for now.
end

"""This computes the approximate size of the entire display, not the individual lenslet displays."""
sizeofdisplay(fov,eyerelief) = @. 2 * tand(fov / 2) * eyerelief
export sizeofdisplay

"""Computes the number of lenslets required to create a specified field of view at the given eyerelief"""
function numberoflenslets(fov, eyerelief::Unitful.Length, lensletdiameter::Unitful.Length)
    lensletarea = π * (lensletdiameter / 2)^2
    dispsize = sizeofdisplay(fov, eyerelief)
     return dispsize[1] * dispsize[2] / lensletarea
end
export numberoflenslets

"""angular size of the eyebox when viewed from distance eyerelief"""
eyeboxangles(eyebox,eyerelief) = @. atand(uconvert(Unitful.NoUnits, eyebox / eyerelief))
export eyeboxangles

"""computes how the fov can be subdivided among lenslets based on cluster size. Assumes the horizontal size of the eyebox is larger than the vertical so the larger number of subdivisions will always be the first number in the returns Tuple."""
function anglesubdivisions(pupildiameter::Unitful.Length, λ::Unitful.Length, mtf, cyclesperdegree;RGB=true)
    cluster, _ = choosecluster(pupildiameter, λ, mtf, cyclesperdegree)
    numelements = Repeat.clustersize(cluster)
    if numelements == 37
        return RGB ? (4,3) : (6,6)
    elseif numelements == 19
        return RGB ? (3, 2) : (5, 3)
    elseif numelements == 12
        return RGB ? (2, 2) : (4, 3)
   elseif numelements == 9
        return RGB ? (3, 1) : (3, 3)
    elseif numelements == 7
        return RGB ? (3, 1) : (3, 2)
    elseif numelements == 4 || numelements == 3
        return RGB ? (1, 1) : (3, 1)
    else throw(ErrorException("this should never happen"))
    end
end
export anglesubdivisions

"""Computes the approximate fov required of each lenslet for the given constraints. This is strictly correct only for a lenslet centered in front of the eyebox, but the approximation is good enough for high level analysis"""
function lensletangles(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length; clusterproperties=defaultclusterproperties(), RGB=true)
    return eyeboxangles(eyebox, eyerelief) ./ anglesubdivisions(pupildiameter, clusterproperties.λ, clusterproperties.mtf, clusterproperties.cyclesperdegree, RGB=RGB)
end
export lensletangles

"""Multilens displays tradeoff pixel redundancy for a reduction in total track of the display, by using many short focal length lenses to cover the eyebox. This function computes the ratio of pixels in the multilens display vs. a conventional display of the same nominal resolution"""
function pixelredundancy(fov, eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length, ppd; RGB=true)
    lensprops = defaultclusterproperties()
    clusterdata = choosecluster(pupildiameter, lensprops.λ, lensprops.mtf, lensprops.cyclesperdegree)
    nominalresolution = ustrip.(°, fov) .* ppd #remove degree units so pixel redundancy doesn't have units of °^-2
    angles = lensletangles(eyerelief, eyebox, pupildiameter, RGB=RGB)
    pixelsperlenslet = angles .* ppd
    numlenses = numberoflenslets(fov, eyerelief, clusterdata.lensletdiameter)
    return (numlenses * pixelsperlenslet[1] * pixelsperlenslet[2]) / ( nominalresolution[1] * nominalresolution[2])
end
export pixelredundancy

label(color) = color ? "RGB" : "Monochrome"

"""generates a contour plot showing pixel redundancy as a function of ppd and pupil diameter"""
function redundancy_ppdvspupildiameter()
    x = 20:2:45
    y = 3.0:.05:4

    RGB = true
    Plots.plot(Plots.contour(x, y, (x, y) -> pixelredundancy((50, 35), 18mm, (10mm, 6mm), y * mm, x, RGB=RGB), fill=true, xlabel="pixels per degree", ylabel="pupil diameter", legendtitle="pixel redundancy", title="$(label(RGB)) lenslets"))
end
export redundancy_ppdvspupildiameter

"""generates a contour plot of lenslet display size as a function of ppd and pupil diameter"""
function displaysize_ppdvspupildiameter()
    x = 20:2:45
    y = 3.0:.05:4
    RGB = true

    Plots.plot(Plots.contour(x, y, (x, y) -> maximum(ustrip.(μm, lensletdisplaysize((50, 35), 18mm, (10mm, 6mm), y * mm, x, RGB=RGB))), fill=true, xlabel="pixels per degree", ylabel="pupil diameter", legendtitle="display size μm", title="$(label(RGB)) lenslets"))
end
export displaysize_ppdvspupildiameter

""" Compute high level properties of a lenslet display system. Uses basic geometric and optical constraints to compute approximate values which should be within 10% or so of a real physical system.

`eyerelief` is in mm: `18mm` not `18`.

`eyebox` is a 2 tuple of lengths, also in mm, that specifies the size of the rectangular eyebox, assumed to lie on vertex of the cornea when the eye is looking directly forward.

`fov` is the field of view of the display as seen by the eye, in degrees (no unit type necessary here).

`pupildiameter` is the diameter of the eye pupil in mm.

`mtf` is the MTF response of the system at the angular frequency `cyclesperdegree`.

`pixelsperdegree` is the number of pixels per degree on the display, as seen by the eye through the lenslets. This is not a specification of the optical resolution of the system (use `mtf` and `cyclesperdegree` for that). It is a physical characteristic of the display.

Example:
```
julia> system_properties(18mm,(10mm,9mm),(55°,45°),4.0mm,.2,11.0)
Dict{Symbol, Any} with 14 entries:
  :lenslet_diameter     => 0.742781 mm
  :lenslet_display_size => (251.708 μm, 345.939 μm)
  :fnumber              => 2.0
  :number_lenslets      => 644.903
  :total_silicon_area   => 56.1555 mm^2
  :subdivisions         => (3, 2)
  :pixel_redundancy     => 33.5192
  :diffraction_limit    => 24.4603
  :eyebox_angles        => (29.0546, 26.5651)
  :pixels_per_degree    => 28.8088
  :lenslet_fov          => (9.68487, 13.2825)
  :display_size         => (18.7404 mm, 14.9117 mm)
  :cluster_data         => (cluster = ClusterWithProperties{19, 2, Float64}(Lat…
  :focal_length         => 1.48556 mm

```
"""
function system_properties(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, fov, pupildiameter::Unitful.Length, mtf, cyclesperdegree; minfnumber=2.0,RGB=true,λ=530nm,pixelpitch=.9μm, maxdisplaysize = 250μm)::Dict{Symbol,Any}
    diameter = diameter_for_cycles_deg(mtf, cyclesperdegree, λ)
    clusterdata = choosecluster(pupildiameter, diameter)
    difflimit = diffractionlimit(λ, clusterdata.lensletdiameter)
    numlenses = numberoflenslets(fov, eyerelief, clusterdata.lensletdiameter)
    redundancy = pixelredundancy(fov, eyerelief, eyebox, pupildiameter, difflimit, RGB=RGB)
    subdivisions = anglesubdivisions(pupildiameter, λ, mtf, cyclesperdegree, RGB=RGB)
    eyebox_angles = eyeboxangles(eyebox,eyerelief)
    angles = lensletangles(eyerelief, eyebox, pupildiameter, clusterproperties=(mtf = mtf, minfnumber = minfnumber, cyclesperdegree = cyclesperdegree, λ = λ, pixelpitch = pixelpitch))

    focal_length = compute_focal_length(angles,clusterdata.lensletdiameter,minfnumber,maxdisplaysize)
    dispsize = lensletdisplaysize(angles,focal_length)
    fnumber = focal_length/clusterdata.lensletdiameter
    siliconarea = uconvert(mm^2,numlenses  * dispsize[1]*dispsize[2])
    fulldisplaysize = sizeofdisplay(fov,eyerelief)

    pixels_per_degree = pixelsperdegree(focal_length,pixelpitch)
    return Dict(:cluster_data => clusterdata, :lenslet_diameter => clusterdata.lensletdiameter,:pixels_per_degree => pixels_per_degree, :diffraction_limit => difflimit, :fnumber => fnumber, :focal_length => focal_length, :display_size => fulldisplaysize, :lenslet_display_size => dispsize, :total_silicon_area => siliconarea, :number_lenslets => numlenses, :pixel_redundancy => redundancy, :eyebox_angles => eyebox_angles, :lenslet_fov => angles, :subdivisions => subdivisions)
end
export system_properties

"""prints system properties nicely"""
function printsystemproperties(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, fov, pupildiameter::Unitful.Length, mtf, cyclesperdegree; minfnumber=2.0,RGB=true,λ=530nm,pixelpitch=.9μm,maxdisplaysize = 350μm)
    println("Input parameters")
    println()
    println("pixel pitch $pixelpitch")
    println("eye relief = $eyerelief")
    println("eye box = $eyebox")
    println("fov = $(fov)")
    println("pupil diameter = $pupildiameter")
    println("mtf = $mtf @ $cyclesperdegree cycles/°")
    println()
    printsystemproperties(system_properties(eyerelief, eyebox, fov, pupildiameter, mtf, cyclesperdegree, minfnumber = minfnumber,RGB=RGB,λ=λ,pixelpitch=pixelpitch,maxdisplaysize = maxdisplaysize))
end
export printsystemproperties

function printsystemproperties(props)
    println("Output values")
    println()

    for (key,value) in pairs(props)
        if key == :cluster_data
            println("cluster = $(typeof(props[:cluster_data][:cluster]))")
        elseif key == :diffraction_limit
                println("$key = $value cycles/°")
        else
                println("$key = $value")
        end
    end

    println("occlusion ratio  $(π*uconvert(Unitful.NoUnits,*(props[:lenslet_display_size]...)/(props[:lenslet_diameter])^2))")
    cluster = props[:cluster_data][:cluster]
    clusterdiameter = Repeat.euclideandiameter(cluster)
    println("cluster diameter (approx): $(clusterdiameter)")
end


text:
using LinearAlgebra
using Plots

include("cheb.jl")
include("parameters.jl")
include("geometry.jl")
include("profile.jl")
include("build_A.jl")
include("compute_spectrum.jl")
include("mesh.jl")
include("plot_fields.jl")

params=parameters(
     H = 1,
    Lϕ = π/100,
  ϕmin = π/4,
   TwoΩ = 1,
      g = 10,
      a = 1,
    Lj = π/2000,
    Uj = 1.0,
    Nϕ = 100,
    dk = 1e-2,
    kₘ = 3
    );

params.Fr = (params.TwoΩ * params.Lj)^2/(params.g * params.H);
params.Ro =  params.Uj/(params.TwoΩ * params.Lj)

print("\n")
print("Linear-Stability-Caluculator\n")
print("============================\n\n")
print("Repo:      Linear-Stability-Calculators\n")
print("Code:      ShallowWater/Julia/Spherical/linear_stability_shallow_water.jl\n")
print("Model:     Rotating Shallow Water\n")
print("Geometry:  Spherical\n")
print("Structure: Bickley Jet\n")
print("\n")
print("Nondimensional Parameters\n")
print("=========================\n")
print("Fr = ", params.Fr, "\n")
print("Ro = ", params.Ro, "\n")

# Wavenumbers
ks = collect(params.dk:params.dk:params.kₘ) / params.Lj;
Nk = length(ks);

# Grid and basic state
Dϕ, Dϕ2, ϕ = geometry(params);
      U,  E = profile(  ϕ, params);

plot_basic_state(ϕ, U, E, "basic_state.png")

# initialize fields to store
Nmodes = 2
     σ = zeros(Nmodes, Nk);
     ω = zeros(Nmodes, Nk);
σmodes = zeros(ComplexF64, 3*params.Nϕ+1, Nmodes, Nk);

# Compute growth rates
for cnt in 1:Nk
    local k = ks[cnt]
    local A = build_A(k, U, E, Dϕ, Dϕ2, ϕ, params);
    local σ[:,cnt], ω[:,cnt], σmodes[:,:,cnt] = compute_spectrum(A, k, params, Nmodes);
end

plot_growth_rates(ks, σ, Nmodes, "growth_rates_spherical_Bickley_jet.png")

#print("σ = ", σ[1,:],"\n")

mode_number = 1;
      σ_max = maximum(σ[mode_number,:]);
    k_index = sortperm(σ[mode_number,:],rev=true)[1];
          k = ks[k_index]; # pick wavenumber

plot_1D_streamfunction(k_index, k, ϕ, σmodes, mode_number, "modes_1D_streamfunction.png")
plot_2D_streamfunction(k_index, k, ϕ, σmodes, mode_number, "modes_2D_streamfunction.png")
plot_2D_vorticity(Dϕ,  k_index, k, ϕ, σmodes, mode_number, "modes_2D_vorticity.png")


text:
@testset "Test flux limiters" begin
    include(joinpath("VlasovSolver", "test_flux_limiters.jl"))
end

@testset "Test Vlasov solvers" begin
    include(joinpath("VlasovSolver", "test_1d_advection.jl"))
end


text:
# Autogenerated wrapper script for Fontconfig_jll for i686-linux-musl
export fc_cache, fc_cat, fc_conflist, fc_list, fc_match, fc_pattern, fc_query, fc_scan, fc_validate, fonts_conf, libfontconfig

using FreeType2_jll
using Bzip2_jll
using Zlib_jll
using Libuuid_jll
using Expat_jll
JLLWrappers.@generate_wrapper_header("Fontconfig")
JLLWrappers.@declare_executable_product(fc_cache)
JLLWrappers.@declare_executable_product(fc_cat)
JLLWrappers.@declare_executable_product(fc_conflist)
JLLWrappers.@declare_executable_product(fc_list)
JLLWrappers.@declare_executable_product(fc_match)
JLLWrappers.@declare_executable_product(fc_pattern)
JLLWrappers.@declare_executable_product(fc_query)
JLLWrappers.@declare_executable_product(fc_scan)
JLLWrappers.@declare_executable_product(fc_validate)
JLLWrappers.@declare_file_product(fonts_conf)
JLLWrappers.@declare_library_product(libfontconfig, "libfontconfig.so.1")
function __init__()
    JLLWrappers.@generate_init_header(FreeType2_jll, Bzip2_jll, Zlib_jll, Libuuid_jll, Expat_jll)
    JLLWrappers.@init_executable_product(
        fc_cache,
        "bin/fc-cache",
    )

    JLLWrappers.@init_executable_product(
        fc_cat,
        "bin/fc-cat",
    )

    JLLWrappers.@init_executable_product(
        fc_conflist,
        "bin/fc-conflist",
    )

    JLLWrappers.@init_executable_product(
        fc_list,
        "bin/fc-list",
    )

    JLLWrappers.@init_executable_product(
        fc_match,
        "bin/fc-match",
    )

    JLLWrappers.@init_executable_product(
        fc_pattern,
        "bin/fc-pattern",
    )

    JLLWrappers.@init_executable_product(
        fc_query,
        "bin/fc-query",
    )

    JLLWrappers.@init_executable_product(
        fc_scan,
        "bin/fc-scan",
    )

    JLLWrappers.@init_executable_product(
        fc_validate,
        "bin/fc-validate",
    )

    JLLWrappers.@init_file_product(
        fonts_conf,
        "etc/fonts/fonts.conf",
    )

    JLLWrappers.@init_library_product(
        libfontconfig,
        "lib/libfontconfig.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
    ENV["FONTCONFIG_FILE"] = get(ENV, "FONTCONFIG_FILE", fonts_conf)
    ENV["FONTCONFIG_PATH"] = get(ENV, "FONTCONFIG_PATH", dirname(ENV["FONTCONFIG_FILE"]))

end  # __init__()


text:
export generate_ref_dataset, train_iceflow_UDE
export predict_A̅, A_fake

"""
    generate_ref_dataset(temp_series, H₀)

Generate reference dataset based on the iceflow PDE
"""
function generate_ref_dataset(gdirs_climate, tspan, solver = Ralston())

    # Perform reference simulation with forward model
    println("Running forward PDE ice flow model...\n")
    # Run batches in parallel
    gdirs = gdirs_climate[2]
    longterm_temps = gdirs_climate[3]
    A_noises = randn(rng_seed(), length(gdirs)) .* 10e-18
    refs = @showprogress pmap((gdir, longterm_temp, A_noise) -> batch_iceflow_PDE(gdir, longterm_temp, A_noise, tspan, solver), gdirs, longterm_temps, A_noises)

    # Split into different vectors
    H_refs, V̄x_refs, V̄y_refs = [],[],[]
    for ref in refs
        push!(H_refs, ref["H"])
        push!(V̄x_refs, ref["Vx"])
        push!(V̄y_refs, ref["Vy"])
    end

    return H_refs, V̄x_refs, V̄y_refs
end

"""
    batch_iceflow_PDE(climate, gdir, context)

Solve the Shallow Ice Approximation iceflow PDE for a given temperature series batch
"""
function batch_iceflow_PDE(gdir, longterm_temp, A_noise, tspan, solver)
    println("Processing glacier: ", gdir.rgi_id)

    context, H = build_PDE_context(gdir, longterm_temp, A_noise, tspan)
    refs = simulate_iceflow_PDE(H, context, solver)

    return refs
end

"""
    simulate_iceflow_PDE(H, context, solver)

Make forward simulation of the SIA PDE.
"""
function simulate_iceflow_PDE(H, context, solver)
    tspan = context.x[22]
    iceflow_prob = ODEProblem(iceflow!,H,tspan,context)
    iceflow_sol = solve(iceflow_prob, solver,
                    reltol=1e-6, save_everystep=false,
                    progress=true, progress_steps = 10)
    # Compute average ice surface velocities for the simulated period
    H_ref = iceflow_sol.u[end]
    temps = context.x[7]
    V̄x_ref, V̄y_ref = avg_surface_V(context, H_ref, mean(temps), "PDE") # Average velocity with average temperature
    refs = Dict("Vx"=>V̄x_ref, "Vy"=>V̄y_ref, "H"=>H_ref)
    return refs
end

"""
    train_iceflow_UDE(H₀, UA, θ, train_settings, PDE_refs, temp_series)

Train the Shallow Ice Approximation iceflow UDE
"""
function train_iceflow_UDE(gdirs_climate, tspan, train_settings, PDE_refs, θ_trained=[], solver = ROCK4(), loss_history=[])
    if length(θ_trained) == 0
        global current_epoch = 1 # reset epoch count
    end
    optimizer = train_settings[1]
    epochs = train_settings[2]
    UA, θ = get_NN(θ_trained)
    gdirs = gdirs_climate[2]
    H_refs = PDE_refs["H_refs"]
    Vx_refs = PDE_refs["V̄x_refs"]
    Vy_refs = PDE_refs["V̄y_refs"]
    # Build context for all the batches before training
    println("Building context...")
    context_batches = pmap((gdir, H_ref, Vx_ref, Vy_ref) -> build_UDE_context(gdir, H_ref, Vx_ref, Vy_ref, tspan), gdirs, H_refs, Vx_refs, Vy_refs)
    loss(θ) = loss_iceflow(θ, UA, gdirs_climate, context_batches, PDE_refs, solver) # closure

    println("Training iceflow UDE...")
    temps = gdirs_climate[3]
    A_noise = randn(rng_seed(), length(gdirs)).*6e-18
    cb(θ, l, UA) = callback(θ, l, UA, temps, A_noise)
    iceflow_trained = DiffEqFlux.sciml_train(loss, θ, optimizer, cb=cb, maxiters = epochs)
    # iceflow_trained = DiffEqFlux.sciml_train(loss, θ, optimizer, cb=callback, maxiters = epochs)

    return iceflow_trained, UA
end

callback = function (θ, l, UA, temps, A_noise) # callback function to observe training
    println("Epoch #$current_epoch - Loss $loss_type: ", l)

    avg_temps = [mean(temps[i]) for i in 1:length(temps)]
    p = sortperm(avg_temps)
    avg_temps = avg_temps[p]
    pred_A = predict_A̅(UA, θ, collect(-23:1:0)')
    pred_A = [pred_A...] # flatten
    true_A = A_fake(avg_temps, A_noise[p], noise)

    Plots.scatter(avg_temps, true_A, label="True A")
    plot_epoch = Plots.plot!(-23:1:0, pred_A, label="Predicted A",
                        xlabel="Long-term air temperature (°C)",
                        ylabel="A", ylims=(minA,maxA),
                        legend=:topleft)
    Plots.savefig(plot_epoch,joinpath(root_plots,"training","epoch$current_epoch.png"))
    global current_epoch += 1
    push!(loss_history, l)

    false
end

"""
    loss_iceflow(θ, context, UA, PDE_refs::Dict{String, Any}, temp_series)

Loss function based on glacier ice velocities and/or ice thickness
"""
function loss_iceflow(θ, UA, gdirs_climate, context_batches, PDE_refs::Dict{String, Any}, solver)
    H_V_preds = predict_iceflow(θ, UA, gdirs_climate, context_batches, solver)

    # Compute loss function for the full batch
    l_Vx, l_Vy, l_H = 0.0, 0.0, 0.0
    for i in 1:length(H_V_preds)

        # Get ice thickness from the reference dataset
        H_ref = PDE_refs["H_refs"][i]
        # Get ice velocities for the reference dataset
        Vx_ref = PDE_refs["V̄x_refs"][i]
        Vy_ref = PDE_refs["V̄y_refs"][i]
        # Get ice thickness from the UDE predictions
        H = H_V_preds[i][1]
        # Get ice velocities prediction from the UDE
        V̄x_pred = H_V_preds[i][2]
        V̄y_pred = H_V_preds[i][3]

        if random_sampling_loss
            # sample random indices for which V_ref is non-zero
            n_sample, n_counts = 50, 0
            nxy = length(H_ref[H_ref .!= 0.0])

            while n_counts < n_sample
                j = rand(1:nxy-2)
                H_ref_f = H_ref[H_ref .!= 0.0]
                Vx_ref_f = Vx_ref[Vx_ref .!= 0.0]
                Vy_ref_f = Vy_ref[Vy_ref .!= 0.0]
                H_f = H[H_ref .!= 0.0]
                V̄x_pred_f = V̄x_pred[Vx_ref .!= 0.0]
                V̄y_pred_f = V̄y_pred[Vy_ref .!= 0.0]

                if norm_loss
                    normH = H_ref_f[j]  .+ ϵ
                    normVx = Vx_ref_f[j] .+ ϵ
                    normVy = Vy_ref_f[j] .+ ϵ
                else
                    normH, normH, normVy = 1, 1, 1
                end
                l_H += Flux.Losses.mse(H_f[j]  ./normH, H_ref_f[j]  ./normH; agg=mean)
                l_Vx += Flux.Losses.mse(V̄x_pred_f[j]  ./normVx, Vx_ref_f[j]  ./normVx; agg=mean)
                l_Vy += Flux.Losses.mse(V̄y_pred_f[j]  ./normVy, Vy_ref_f[j]  ./normVy; agg=mean)

                n_counts += 1
            end
            l_H = l_H / n_sample
            l_Vx = l_Vx / n_sample
            l_Vy = l_Vy / n_sample
        else
            # Classic loss function with the full matrix
            if scale_loss
                normH = H_ref[H_ref .!= 0.0] .+ ϵ
                normVx = Vx_ref[Vx_ref .!= 0.0] .+ ϵ
                normVy = Vy_ref[Vy_ref .!= 0.0] .+ ϵ
                l_H += Flux.Losses.mse(H[H_ref .!= 0.0] ./normH, H_ref[H_ref.!= 0.0] ./normH; agg=mean)
                l_Vx += Flux.Losses.mse(V̄x_pred[Vx_ref .!= 0.0] ./normVx, Vx_ref[Vx_ref.!= 0.0] ./normVx; agg=mean)
                l_Vy += Flux.Losses.mse(V̄y_pred[Vy_ref .!= 0.0] ./normVy, Vy_ref[Vy_ref.!= 0.0] ./normVy; agg=mean)
            else
                l_H += Flux.Losses.mse(H[H_ref .!= 0.0], H_ref[H_ref.!= 0.0]; agg=mean)
                l_Vx += Flux.Losses.mse(V̄x_pred[Vx_ref .!= 0.0], Vx_ref[Vx_ref.!= 0.0]; agg=mean)
                l_Vy += Flux.Losses.mse(V̄y_pred[Vy_ref .!= 0.0], Vy_ref[Vy_ref.!= 0.0]; agg=mean)
            end
        end
    end

    @assert (loss_type == "H" || loss_type == "V" || loss_type == "HV") "Invalid loss_type. Needs to be 'H', 'V' or 'HV'"
    if loss_type == "H"
        l_avg = l_H/length(PDE_refs["H_refs"])
    elseif loss_type == "V"
        l_avg = (l_Vx/length(PDE_refs["V̄x_refs"]) + l_Vy/length(PDE_refs["V̄y_refs"]))/2
    elseif loss_type == "HV"
        l_avg = (l_Vx/length(PDE_refs["V̄x_refs"]) + l_Vy/length(PDE_refs["V̄y_refs"]) + l_H/length(PDE_refs["H_refs"]))/3
    end
    return l_avg, UA
end

"""
    predict_iceflow(θ, UA, context, temp_series)

Makes a prediction of glacier evolution with the UDE for a given temperature series
"""
function predict_iceflow(θ, UA, gdirs_climate, context_batches, solver)

    # Train UDE in parallel
    # gdirs_climate = (dates, gdirs, longterm_temps, annual_temps)
    longterm_temps = gdirs_climate[3]
    H_V_pred = pmap((context, longterm_temps_batch) -> batch_iceflow_UDE(θ, UA, context, longterm_temps_batch, solver), context_batches, longterm_temps)
    return H_V_pred
end

"""
    batch_iceflow_UDE(θ, H, climate, context)

Solve the Shallow Ice Approximation iceflow UDE for a given temperature series batch
"""
function batch_iceflow_UDE(θ, UA, context, longterm_temps_batch, solver)
    # Retrieve long-term temperature series
    H = context[3]
    tspan = context[6]
    iceflow_UDE_batch(H, θ, t) = iceflow_NN(H, θ, t, UA, context, longterm_temps_batch) # closure
    iceflow_prob = ODEProblem(iceflow_UDE_batch,H,tspan,θ)
    iceflow_sol = solve(iceflow_prob, solver, u0=H, p=θ,
                    reltol=1e-6, save_everystep=false,
                    progress=true, progress_steps = 10)
    # Get ice velocities from the UDE predictions
    H_pred = iceflow_sol.u[end]
    V̄x_pred, V̄y_pred = avg_surface_V(context, H_pred, mean(longterm_temps_batch), "UDE", θ, UA) # Average velocity with average temperature
    H_V_pred = (H_pred, V̄x_pred, V̄y_pred)
    return H_V_pred
end

"""
    iceflow!(dH, H, context,t)

Runs a single time step of the iceflow PDE model in-place
"""
function iceflow!(dH, H, context,t)
    # Unpack parameters
    #A, B, S, dSdx, dSdy, D, temps, dSdx_edges, dSdy_edges, ∇S, Fx, Fy, Vx, Vy, V, C, α, current_year
    current_year = Ref(context.x[18])
    A = Ref(context.x[1])
    t₁ = context.x[22][end]
    A_noise = context.x[23]

    # Get current year for MB and ELA
    year = floor(Int, t) + 1
    if year != current_year[] && year <= t₁
        temp = Ref{Float64}(context.x[7][year])
        A[] .= A_fake(temp[], A_noise, noise)
        current_year[] .= year
    end

    # Compute the Shallow Ice Approximation in a staggered grid
    SIA!(dH, H, context)
end

"""
    iceflow_NN(H, θ, t, context, temps, UA)

Runs a single time step of the iceflow UDE model
"""
function iceflow_NN(H, θ, t, UA, context, temps)

    year = floor(Int, t) + 1
    t₁ = context[6][end]
    if year <= t₁
        temp = temps[year]
    else
        temp = temps[year-1]
    end
    A = predict_A̅(UA, θ, [temp]) # FastChain prediction requires explicit parameters

    # Compute the Shallow Ice Approximation in a staggered grid
    return SIA(H, A, context)
end

"""
    SIA!(dH, H, context)

Compute an in-place step of the Shallow Ice Approximation PDE in a forward model
"""
function SIA!(dH, H, context)
    # Retrieve parameters
    #[A], B, S, dSdx, dSdy, D, copy(temp_series[1]), dSdx_edges, dSdy_edges, ∇S, Fx, Fy, Vx, Vy, V, C, α, [current_year], nxy, Δxy
    A = context.x[1]
    B = context.x[2]
    S = context.x[3]
    dSdx = context.x[4]
    dSdy = context.x[5]
    D = context.x[6]
    dSdx_edges = context.x[8]
    dSdy_edges = context.x[9]
    ∇S = context.x[10]
    Fx = context.x[11]
    Fy = context.x[12]
    Δx = context.x[20][1]
    Δy = context.x[20][2]

    # Update glacier surface altimetry
    S .= B .+ H

    # All grid variables computed in a staggered grid
    # Compute surface gradients on edges
    dSdx .= diff_x(S) ./ Δx
    dSdy .= diff_y(S) ./ Δy
    ∇S .= (avg_y(dSdx).^2 .+ avg_x(dSdy).^2).^((n - 1)/2)

    Γ = 2.0 * A * (ρ * g)^n / (n+2.0) # 1 / m^3 s
    D .= Γ .* avg(H).^(n + 2.0) .* ∇S

    # Compute flux components
    dSdx_edges .= diff_x(S[:,2:end - 1]) ./ Δx
    dSdy_edges .= diff_y(S[2:end - 1,:]) ./ Δy
    Fx .= .-avg_y(D) .* dSdx_edges
    Fy .= .-avg_x(D) .* dSdy_edges

    #  Flux divergence
    inn(dH) .= .-(diff_x(Fx) ./ Δx .+ diff_y(Fy) ./ Δy) # MB to be added here
end

"""
    SIA(H, A, context)

Compute a step of the Shallow Ice Approximation UDE in a forward model. Allocates memory.
"""
function SIA(H, A, context)
    # Retrieve parameters
    # context = (B, H₀, H, nxy, Δxy)
    B = context[1]
    Δx = context[5][1]
    Δy = context[5][2]

    # Update glacier surface altimetry
    S = B .+ H

    # All grid variables computed in a staggered grid
    # Compute surface gradients on edges
    dSdx = diff_x(S) ./ Δx
    dSdy = diff_y(S) ./ Δy
    ∇S = (avg_y(dSdx).^2.0 .+ avg_x(dSdy).^2.0).^((n - 1.0)/2.0)

    Γ = 2.0 * A * (ρ * g)^n / (n+2.0) # 1 / m^3 s
    D = Γ .* avg(H).^(n + 2.0) .* ∇S

    # Compute flux components
    dSdx_edges = diff_x(S[:,2:end - 1]) ./ Δx
    dSdy_edges = diff_y(S[2:end - 1,:]) ./ Δy
    Fx = .-avg_y(D) .* dSdx_edges
    Fy = .-avg_x(D) .* dSdy_edges

    #  Flux divergence
    @tullio dH[i,j] := -(diff_x(Fx)[pad(i-1,1,1),pad(j-1,1,1)] / Δx + diff_y(Fy)[pad(i-1,1,1),pad(j-1,1,1)] / Δy) # MB to be added here

    return dH
end

"""
    avg_surface_V(H, B, temp)

Computes the average ice velocity for a given input temperature
"""
function avg_surface_V(context, H, temp, sim, θ=[], UA=[])
    # context = (B, H₀, H, nxy, Δxy)
    B, H₀, Δx, Δy, A_noise = retrieve_context(context)

    # Update glacier surface altimetry
    S = B .+ (H₀ .+ H)./2.0 # Use average ice thickness for the simulated period

    # All grid variables computed in a staggered grid
    # Compute surface gradients on edges
    dSdx = diff_x(S) / Δx
    dSdy = diff_y(S) / Δy
    ∇S = (avg_y(dSdx).^2.0 .+ avg_x(dSdy).^2.0).^((n - 1.0)/2.0)

    @assert (sim == "UDE" || sim == "PDE") "Wrong type of simulation. Needs to be 'UDE' or 'PDE'."
    if sim == "UDE"
        A = predict_A̅(UA, θ, [temp]) # FastChain prediction requires explicit parameters
    elseif sim == "PDE"
        A = A_fake(temp, A_noise, noise)
    end
    Γꜛ = 2.0 * A * (ρ * g)^n / (n+1.0) # surface stress (not average)  # 1 / m^3 s
    D = Γꜛ .* avg(H).^(n + 1.0) .* ∇S

    # Compute averaged surface velocities
    Vx = - D .* avg_y(dSdx)
    Vy = - D .* avg_x(dSdy)

    return Vx, Vy

end

# Polynomial fit for Cuffey and Paterson data
A_f = fit(A_values[1,:], A_values[2,:]) # degree = length(xs) - 1

"""
    A_fake(temp, noise=false)

Fake law establishing a theoretical relationship between ice viscosity (A) and long-term air temperature.
"""
function A_fake(temp, A_noise=nothing, noise=false)
    # A = @. minA + (maxA - minA) * ((temp-minT)/(maxT-minT) )^2
    A = A_f.(temp) # polynomial fit
    if noise
        A = abs.(A .+ A_noise)
    end
    return A
end

"""
    predict_A̅(UA, θ, temp)

Predicts the value of A with a neural network based on the long-term air temperature.
"""
function predict_A̅(UA, θ, temp)
    return UA(temp, θ) .* 1e-17
end

# function fake_temp_series(t, means=Array{Float64}([0.0,-2.0,-3.0,-5.0,-10.0,-12.0,-14.0,-15.0,-20.0]))
#     temps, norm_temps, norm_temps_flat = [],[],[]
#     for mean in means
#         push!(temps, mean .+ rand(t).*1e-1) # static
#         append!(norm_temps_flat, mean .+ rand(t).*1e-1) # static
#     end

#     # Normalise temperature series
#     norm_temps_flat = Flux.normalise([norm_temps_flat...]) # requires splatting

#     # Re-create array of arrays
#     for i in 1:t₁:length(norm_temps_flat)
#         push!(norm_temps, norm_temps_flat[i:i+(t₁-1)])
#     end

#     return temps, norm_temps
# end

"""
    get_initial_geometry(glacier_gd)

Retrieves the initial glacier geometry (bedrock + ice thickness) for a glacier.
"""
function get_initial_geometry(gdir, smoothing=true)
    # Load glacier gridded data
    glacier_gd = xr.open_dataset(gdir.get_filepath("gridded_data"))
    H₀ = glacier_gd.consensus_ice_thickness.data # initial ice thickness conditions for forward model
    fillNaN!(H₀) # Fill NaNs with 0s to have real boundary conditions
    if smoothing
        smooth!(H₀)  # Smooth initial ice thickness to help the solver
    end
    H = deepcopy(H₀)
    B = glacier_gd.topo.data .- H₀ # bedrock

    nx = glacier_gd.y.size # glacier extent
    ny = glacier_gd.x.size # really weird, but this is inversed
    Δx = abs(gdir.grid.dx)
    Δy = abs(gdir.grid.dy)

    return H₀, H, B, (nx,ny), (Δx,Δy)
end

function build_PDE_context(gdir, longterm_temp, A_noise, tspan)
    # Determine initial geometry conditions
    H₀, H, B, nxy, Δxy = get_initial_geometry(gdir)
    # Initialize all matrices for the solver
    nx, ny = nxy
    S, dSdx, dSdy = zeros(Float64,nx,ny),zeros(Float64,nx-1,ny),zeros(Float64,nx,ny-1)
    dSdx_edges, dSdy_edges, ∇S = zeros(Float64,nx-1,ny-2),zeros(Float64,nx-2,ny-1),zeros(Float64,nx-1,ny-1)
    D, Fx, Fy = zeros(Float64,nx-1,ny-1),zeros(Float64,nx-1,ny-2),zeros(Float64,nx-2,ny-1)
    V, Vx, Vy = zeros(Float64,nx-1,ny-1),zeros(Float64,nx-1,ny-1),zeros(Float64,nx-1,ny-1)
    A = 2e-16
    α = 0                       # Weertman-type basal sliding (Weertman, 1964, 1972). 1 -> sliding / 0 -> no sliding
    C = 15e-14                  # Sliding factor, between (0 - 25) [m⁸ N⁻³ a⁻¹]

    # Gather simulation parameters
    current_year = 0
    context = ArrayPartition([A], B, S, dSdx, dSdy, D, longterm_temp, dSdx_edges, dSdy_edges, ∇S, Fx, Fy, Vx, Vy, V, C, α, [current_year], nxy, Δxy, H₀, tspan, A_noise)
    return context, H
end

function build_UDE_context(gdir, H_ref, Vx_ref, Vy_ref, tspan)
    H₀, H, B, nxy, Δxy = get_initial_geometry(gdir)

    # Tuple with all the temp series and H_refs
    context = (B, H₀, H, nxy, Δxy, tspan)

    return context
end

"""
    retrieve_context(context::Tuple)

Retrieves context variables for computing the surface velocities of the UDE.
"""
function retrieve_context(context::Tuple)
    B = context[1]
    H₀ = context[2]
    Δx = context[5][1]
    Δy = context[5][2]
    return B, H₀, Δx, Δy, nothing
end

"""
    retrieve_context(context::ArrayPartition)

Retrieves context variables for computing the surface velocities of the PDE.
"""
function retrieve_context(context::ArrayPartition)
    B = context.x[2]
    H₀ = context.x[21]
    Δx = context.x[20][1]
    Δy = context.x[20][2]
    A_noise = context.x[23]
    return B, H₀, Δx, Δy, A_noise
end

"""
    get_NN()

Generates a neural network.
"""
function get_NN(θ_trained)
    UA = FastChain(
        FastDense(1,3, x->softplus.(x)),
        FastDense(3,10, x->softplus.(x)),
        FastDense(10,3, x->softplus.(x)),
        FastDense(3,1, sigmoid_A)
    )
    # See if parameters need to be retrained or not
    if isempty(θ_trained)
        θ = initial_params(UA)
    else
        θ = θ_trained
    end
    return UA, θ
end

function sigmoid_A(x)
    minA_out = 8.5e-3 # /!\     # these depend on predict_A̅, so careful when changing them!
    maxA_out = 8.0
    return minA_out + (maxA_out - minA_out) / ( 1.0 + exp(-x) )
end

max_stars_repo_name: ODINN-SciML/ODINN

max_issues_repo_name: ODINN-SciML/ODINN

max_forks_repo_name: ODINN-SciML/ODINN


text:
module JuPOT

import JuMP
import DataFrames
# We extend length so it needs to be imported

import Base: length, show

export AbstractModel,
       optimize,
       AssetsCollection,
       SimpleMVO,
       RobustMVO,
       CVaRO,
       MinVarO,
       getAssetAndReturnsFromCSV,
       getReturns,
       setReturns,
       getCovariance,
       getCoVarForAssetPair,
       getVarForAsset,
       getReturnForAsset,
       setReturns,
       setCovariance,
       setCoVarForAssetPair,
       setVarForAsset,
       setReturnForAsset,
       getNames,
       setNames,
       exportModelResultsToCSV

include("core/utils.jl")
include("core/assetscollection.jl")
include("core/abstractmodel.jl")

include("models/simplemvo.jl")
include("models/robustmvo.jl")
include("models/minvaro.jl")
include("models/cvaro.jl")
end


text:
using OptimizationAlgorithms
const Optimization = OptimizationAlgorithms

# optimization of kernel hyperparameters of Gaussian Process
# could add optimization w.r.t. leave-one-out loss
# TODO: optimization of noise variance
# IDEA: pass inference method (Cholesky, Iterative, Stochastic, ...)
optimize(k, θ, x, y, σ²::Real) = optimize(k, θ, x, 𝒩(y, σ²))
function optimize(k, θ::AbstractVector, x::AbstractVector,
                  y::Union{AbstractVector, Multivariate})
    optimize!(k, copy(θ), x, y)
end
optimize!(k, θ, x, y, σ²::Real) = optimize!(k, θ, x, 𝒩(y, σ²))
function optimize!(k, θ::AbstractVector, x::AbstractVector,
                   y::Union{AbstractVector, Multivariate})
    val(θ) = nlml(k(θ), x, y) # make this "safe"?
    function valdir(θ)
        v, g = value_gradient(nlml, k, θ, x, y) # need to add tol if y is Multivariate
        g .*= -1
        v, g
    end
    # optimization based on L-BFGS with m stages
    m = 5
    D = Optimization.CustomDirection(val, valdir, θ)
    D = Optimization.LBFGS(D, θ, m, check = false)
    # scaling proposed Dection to unit length before doing line search,
    # avoids jumping to sub-optimal and very large parameter vectors
    D = Optimization.UnitDirection(D)
    # to accept a set, do line search for step size to ensure decrease
    D = Optimization.DecreasingStep(D, θ)
    # finally, find fixed point of the iterations
    _, t = Optimization.fixedpoint!(D, θ)
    return θ
end

# here, k(θ) is assumed to be a kernel function
# computes both value and gradient
function value_gradient(::typeof(nlml), k, θ::AbstractVector, x::AbstractVector,
                        y::Union{AbstractVector, Multivariate})
    val, push = pushforward(nlml, k(θ), x, y)

    ∇ = (x, y) -> reshape(ForwardDiff.gradient(z->k(z)(x, y), θ), :, 1)
    G = gramian(∇, x) # lazily evaluated matrix of kernel gradients
    G_matofvec = G.A
    G = matvec2vecmat(G_matofvec) # casts G to a vector of dense matrices. IDEA: make lazy!

    grad = similar(θ)
    @threads for i in eachindex(grad) # parallelizing here is efficient
        grad[i] = push(G[i])
    end
    return val, grad
end

function gradient(::typeof(nlml), k, θ::AbstractVector,
                    x::AbstractVector, y::Union{AbstractVector, Multivariate})
    value_gradient(nlml, k, θ, x, y)[2]
end
using CovarianceFunctions: Gramian
matvec2vecmat(A::Gramian{<:AbstractVecOrMat}) = matvec2vecmat(Matrix(A))
# convert matrix of vectors to vector of matrices
function matvec2vecmat(A::AbstractMatrix{<:AbstractVecOrMat})
    k = length(A[1])
    all(==(k), (length(Ai) for Ai in A)) || throw(DimensionMismatch("component vectors do not have the same size"))
    n, m = size(A)
    B = Vector{Matrix{eltype(A[1])}}(undef, k)
    for i in eachindex(B) # this is still inefficient, because it calculates the gradient multiple times
        B[i] = [a[i] for a in A]
    end
    return B
end

########################### optimizing noise variance ##########################
function value_gradient(::typeof(nlml), log_σ²::Real, K::Factorization, y::AbstractVector)
    n = length(y)
    Σ = exp(log_σ²)*I(n)
    W = Woodbury(Σ, K)
    val, push = pushforward(nlml, W, y)
    grad = push(Σ)
    return val, grad
end

function optimize(σ²::Real, k, x::AbstractVector, y::AbstractVector; tol::Real = tol)
    # pre-factorize kernel matrix
    n = length(x)
    K = gramian(k, x)
    K = factorize(K, tol = tol)
    function val(log_σ²::Real)
        Σ = exp(log_σ²)*I(n)
        nlml(Woodbury(Σ, K), y)
    end
    val(θ::AbstractVector) = val(θ[1])
    function valdir(log_σ²::Real)
        v, g = value_gradient(nlml, log_σ², K, y)
        g = -g
        v, g
    end
    function valdir(θ::AbstractVector)
        v, g = valdir(θ[1])
        v, [g]
    end
    θ = [log(σ²)]
    D = Optimization.CustomDirection(val, valdir, θ)
    D = Optimization.LBFGS(D, θ, 2, check = false) # does this reduce to secant?
    D = Optimization.UnitDirection(D)
    D = Optimization.DecreasingStep(D, θ)
    Optimization.fixedpoint!(D, θ)
    return σ² = exp(θ[1])
end


text:
module Collector

using Compat; import Compat.String

import ..Utilities:

    Utilities,
    isexpr,
    parsefile,
    Head,
    @H_str

import ..Runner:

    State,
    findmethods,
    findtuples,
    findvcats,
    exec,
    getvar,
    pushscope!,
    popscope!

import ..Formats

include("utilities.jl")
include("search.jl")
include("types.jl")
include("docstrings.jl")

end


text:
module ChainLadder

using DataFrames

export ClaimTriangle,
    latest_diagonal,
    IncrementalTriangle, CumulativeTriangle,
    link_ratio

abstract type TriangleType end

struct CumulativeTriangle <: TriangleType end
struct IncrementalTriangle <: TriangleType end

struct ClaimTriangle
    origin_indices
	development_indices
    values
    type::TriangleType

end

function ClaimTriangle(df,origin,development,values,type)
	sort!(df,[origin,order(development, rev=true)])
	vals = map(groupby(df,origin)) do subdf
		return reverse(subdf[:,values])
	end

	return ClaimTriangle(
			sort!(unique(df.origin)),
			sort!(unique(df.origin),rev=true),
			vals,
            type
	)
end

function latest_diagonal(t::ClaimTriangle)
	return [v[end] for v in t.values]
end

function link_ratio(::CumulativeTriangle,v::Vector)
	v[2:end] ./ v[1:end-1]
end

function link_ratio(::IncrementalTriangle,v::Vector)
	cumsum(v)[2:end] ./ v[1:end-1]
end

function link_ratio(t::ClaimTriangle)
    return [link_ratio(t.type,v) for v in t.values[1:end-1]]
end

function valuation_date(t::ClaimTriangle)
	return last(t.development_indices)
end
end


text:
# ---
# title: 813. Largest Sum of Averages
# id: problem813
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/largest-sum-of-averages/description/>
# hidden: true
# ---
#
# We partition a row of numbers `A` into at most `K` adjacent (non-empty)
# groups, then our score is the sum of the average of each group. What is the
# largest score we can achieve?
#
# Note that our partition must use every number in A, and that scores are not
# necessarily integers.
#
#
#
#     **Example:**
#     Input:
#     A = [9,1,2,3,9]
#     K = 3
#     Output: 20
#     Explanation:
#     The best choice is to partition A into [9], [1, 2, 3], [9]. The answer is 9 + (1 + 2 + 3) / 3 + 9 = 20.
#     We could have also partitioned A into [9, 1], [2], [3, 9], for example.
#     That partition would lead to a score of 5 + 2 + 6 = 13, which is worse.
#
#
#
#
# **Note:**
#
#   * `1 <= A.length <= 100`.
#   * `1 <= A[i] <= 10000`.
#   * `1 <= K <= A.length`.
#   * Answers within `10^-6` of the correct answer will be accepted as correct.
#
#
## @lc code=start
using LeetCode

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


text:
@testset "Discretization" begin
  @testset "FIST" begin
    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [2,4,5]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(2,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [2,4]

    𝒫 = Chain(P2[(0,0),(1,0),(1,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == [2,4]

    𝒫 = Chain(P2[(0,0),(1,1),(1,2),(0,0)])
    @test Meshes.ears(𝒫) == []

    𝒫 = Chain(P2[(0.443339268495331, 0.283757618605357),
                 (0.497822414616971, 0.398142813114205),
                 (0.770343126156527, 0.201815462842808),
                 (0.761236456732531, 0.330085709922366),
                 (0.985658085510286, 0.221530395507904),
                 (0.877899962498139, 0.325516131702896),
                 (0.561404274882782, 0.540334008885703),
                 (0.949459768187313, 0.396227653478068),
                 (0.594962560615951, 0.584927547374551),
                 (0.324208409133154, 0.607290684450708),
                 (0.424085089823892, 0.493532112641353),
                 (0.209843417261654, 0.590030658255966),
                 (0.27993878548962, 0.525162463476181),
                 (0.385557753911967, 0.322338556632868),
                 (0.443339268495331, 0.283757618605357)])
    @test Meshes.ears(𝒫) == [1,3,5,6,8,10,12,14]

    points = P2[(0,0),(1,0),(1,1),(2,1),(2,2),(1,2),(0,0)]
    connec = connect.([(4,5,6),(3,4,6),(3,6,1),(1,2,3)], Triangle)
    target = UnstructuredMesh(points[1:end-1], connec)
    poly = PolyArea(points)
    mesh = discretize(poly, FIST())
    @test mesh == target
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "poly1.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "poly2.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "poly3.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "poly4.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "poly5.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "smooth1.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    # TODO: hanging when T=Float32
    # poly = readpoly(T, joinpath(datadir, "smooth2.line"))
    # mesh = discretize(poly, FIST())
    # @test Set(vertices(poly)) == Set(vertices(mesh))
    # @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "smooth3.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "smooth4.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "smooth5.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "hole1.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "hole2.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "hole3.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "hole4.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2

    poly = readpoly(T, joinpath(datadir, "hole5.line"))
    mesh = discretize(poly, FIST())
    @test Set(vertices(poly)) == Set(vertices(mesh))
    @test nelements(mesh) == length(vertices(mesh)) - 2
  end
end


text:
using Documenter, SingularityAppBuilder

makedocs(;
    modules=[SingularityAppBuilder],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/SebastianM-C/SingularityAppBuilder.jl/blob/{commit}{path}#L{line}",
    sitename="SingularityAppBuilder.jl",
    authors="Sebastian Micluța-Câmpeanu <m.c.sebastian95@gmail.com>",
    assets=String[],
)


text:
function LegendreRuleFactory()
    return LegendreRuleFactory(())
end


text:
# Test basic STRIPS functionality in a gripper domain
path = joinpath(dirname(pathof(PDDL)), "..", "test", "strips")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == :gripper
@test domain.predicates[:room] == pddl"(room ?r)"

problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.name == Symbol("gripper-problem")
@test problem.objects == @pddl("rooma", "roomb", "ball1", "ball2", "left", "right")

# Test forward execution of plans
state = init_state(problem)
state = execute(pddl"(pick ball1 rooma left)", state, domain)
@test state[pddl"(carry ball1 left)"] == true
state = execute(pddl"(move rooma roomb)", state, domain)
@test state[pddl"(robbyat roomb)"] == true
state = execute(pddl"(drop ball1 roomb left)", state, domain)
@test state[pddl"(at ball1 roomb)"] == true

@test satisfy(problem.goal, state)[1] == true

state = init_state(problem)
plan = @pddl("(pick ball1 rooma left)",
             "(move rooma roomb)",
             "(drop ball1 roomb left)")
state = execute(plan, state, domain)
@test satisfy(problem.goal, state)[1] == true

# Test action availability
state = init_state(problem)
@test Set(available(state, domain)) == Set{Term}(@pddl(
    "(pick ball1 rooma right)", "(pick ball1 rooma left)",
    "(pick ball2 rooma right)", "(pick ball2 rooma left)",
    "(move rooma roomb)", "(move rooma rooma)"
))

# Test backward regression of plans
state = goal_state(problem)
state = regress(pddl"(drop ball1 roomb left)", state, domain)
@test state[pddl"(carry ball1 left)"] == true
state = regress(pddl"(move rooma roomb)", state, domain)
@test state[pddl"(robbyat rooma)"] == true
state = regress(pddl"(pick ball1 rooma left)", state, domain)
@test state[pddl"(at ball1 rooma)"] == true
@test issubset(state, init_state(problem))

# Test action relevance
state = goal_state(problem)
@test Set(relevant(state, domain)) == Set{Term}(@pddl(
    "(drop ball1 roomb left)", "(drop ball1 roomb right)",
    "(drop ball1 roomb ball1)", "(drop ball1 roomb ball2)",
    "(drop ball1 roomb rooma)", "(drop ball1 roomb roomb)"
))


text:
module Day1

export get_inputs, get_solution1, get_solution2

## Setup
using ..FileIO: read_simple

## Inputs (change these)
function get_inputs()
    test_input1 = test_input2 = [
        1721
        979
        366
        299
        675
        1456
    ]
    test_output1 = 514579
    test_output2 = 241861950

    data = read_simple(joinpath(@__DIR__, "input.txt")) #|> sort

    return (; test_input1, test_input2, test_output1, test_output2, data)
end

## Solution functions (change these)
# Part 1
function get_solution1(data)
    for i in eachindex(data)
        for j in i+1:length(data)
            these_data = (data[i], data[j])
            if sum(these_data)==2020
                return prod(these_data)
            end
        end
    end
end

# Part 2
function get_solution2(data)
    for i in eachindex(data)
        for j in i+1:length(data)
            for k in j+1:length(data)
                these_data = (data[i], data[j], data[k])
                if sum(these_data)==2020
                    return prod(these_data)
                end
            end
        end
    end
end

end