text:
"""
    jit(fun, static_argnums=tuple(); device, backend)

Sets up `fun` for just-in-time compilation with XLA.

Arguments:
  - `fun`: Function to be jitted. Should be a pure function, as side-effects may
    only be executed once. Its arguments and return value should be arrays,
    scalars, or (nested) standard Python containers (tuple/list/dict) thereof.
    Positional arguments indicated by `static_argnums` can be anything at
    all, provided they are hashable and have an equality operation defined.
    Static arguments are included as part of a compilation cache key, which is
    why hash and equality operators must be defined.
  - `static_argnums`: An int or tuple of ints specifying which positional
    arguments to treat as static (compile-time constant). Operations that only
    depend on static arguments will be constant-folded. Calling the jitted
    function with different values for these constants will trigger
    recompilation. If the jitted function is called with fewer positional
    arguments than indicated by `static_argnums` then an error is raised.
    Defaults to ().

  - `device`: This is an experimental feature and the API is likely to change.
    Optional, the Device the jitted function will run on. (Available devices
    can be retrieved via `jax.devices()`.) The default is inherited from
    XLA's DeviceAssignment logic and is usually to use ``jax.devices()[0]``.

  - `backend`: This is an experimental feature and the API is likely to change.
    Optional, a string representing the xla backend. 'cpu','gpu', or 'tpu'.

Returns:
  A wrapped version of ``fun``, set up for just-in-time compilation.
"""
jit(f; static_argnums=tuple(), device=nothing, backend=nothing) = begin
    jit(f, static_argnums...; device=device, backend=backend)
end

function jit(f, s, static_argnums...; device=nothing, backend=nothing)
    # Change from Julia's 1 based to Python's 0 based indexing
    s = tuple(s..., static_argnums...)
    _static_argnums = s .- 1
    return jax.jit(f, static_argnums=_static_argnums, device=device, backend=backend)
end

##

"""
    grad(f, argnums; has_aux=false, holomorphic=false)

Creates a function which evaluates the gradient of ``fun``.

Args:
  - fun: Function to be differentiated. Its arguments at positions specified by
    ``argnums`` should be arrays, scalars, or standard Python containers. It
    should return a scalar (which includes 0-dimensional arrays but not
    1-dimensional arrays of size 1.)
  - argnums: Optional, integer or sequence of integers. Specifies which
    positional argument(s) to differentiate with respect to (default 0).
  - has_aux: Optional, bool. Indicates whether ``fun`` returns a pair where the
    first element is considered the output of the mathematical function to be
    differentiated and the second element is auxiliary data. Default False.
  - holomorphic: Optional, bool. Indicates whether ``fun`` is promised to be
    holomorphic. Default False.

Returns:
  A function with the same arguments as ``fun``, that evaluates the gradient
  of ``fun``. If ``argnums`` is an integer then the gradient has the same
  shape and type as the positional argument indicated by that integer. If
  argnums is a tuple of integers, the gradient is a tuple of values with the
  same shapes and types as the corresponding arguments. If ``has_aux`` is True
  then a pair of (gradient, auxiliary_data) is returned.

For example:

>>> grad_tanh = jax.grad(jax.numpy.tanh)
>>> print(grad_tanh(0.2))
0.961043
"""
grad(f; argnums=tuple(), has_aux=false, holomorphic=false) = begin
    grad(f, static_argnums...; device=device, backend=backend)
end

function grad(f, s, argnums...; has_aux=false, holomorphic=false)
    # Change from Julia's 1 based to Python's 0 based indexing
    s = tuple(s..., argnums...)
    _argnums = s .- 1
    return jax.grad(f, argnums=_argnums, has_aux=has_aux, holomorphic=holomorphic)
end


text:
#!/usr/bin/env julia

ENV["JULIA_MATHKERNEL"] = "/usr/local/Wolfram/WolframEngine/13.0/Executables/WolframKernel"
ENV["JULIA_MATHLINK"]   = "/usr/local/Wolfram/WolframEngine/13.0/SystemFiles/Libraries/Linux-x86-64/libML$(Sys.WORD_SIZE)i4"

using Pkg
Pkg.instantiate()
Pkg.precompile()

using Klio

ext: jl


text:
__precompile__(true)

module Photon

using Printf
using Statistics
import Knet
import CUDA

include("core.jl")
include("layers/layers.jl")
include("metrics/metrics.jl")
include("losses/losses.jl")
include("optimisers/optimisers.jl")
include("train.jl")
include("callbacks/callbacks.jl")
include("data/data.jl")
include("zoo/zoo.jl")

@debug "Loaded Photon"

end


text:
using Compat
using FilePathsBase
using Compat.LinearAlgebra
using Compat.Test

@testset "FilePathsBase" begin

include("mode.jl")
include("path.jl")

end


text:
using Statistics
using Printf

function compute_simd(x::Array{Float64, 1}, y::Array{Float64, 1},
                 grads::Array{Float64, 3}, out::Array{Float64, 1})

    det = -(x[2] - x[3])*y[1] + (x[1] - x[3])*y[2] - (x[1] - x[2])*y[3]
    dv = 1.0 / (6.0 * det)
    ii = 1
    for i=1:6
        for j=1:i
            s = 0.0
            @simd for k=1:3
                s += grads[1, k, i] * grads[1, k, j] + grads[2, k, i] * grads[2, k, j]
            end
            out[ii] = dv * s
            ii += 1
        end
    end

    return out
end

function compute(x::Array{Float64, 1}, y::Array{Float64, 1},
                 grads::Array{Float64, 3}, out::Array{Float64, 1})

    det = -(x[2] - x[3])*y[1] + (x[1] - x[3])*y[2] - (x[1] - x[2])*y[3]
    dv = 1.0 / (6.0 * det)
    ii = 1
    for i=1:6
        for j=1:i
            s = 0.0
            for k=1:3
                s += grads[1, k, i] * grads[1, k, j] + grads[2, k, i] * grads[2, k, j]
            end
            out[ii] = dv * s
            ii += 1
        end
    end

    return out
end

const x = Float64[0., 1., 0.]
const y = Float64[0., 0., 1.]

a11 = -y[1] + y[3]
a12 =  y[1] - y[2]
a21 =  x[1] - x[3]
a22 = -x[1] + x[2]

gq = reshape(
    [-1.0,-1.0,1.0,1.0,-1.0,-1.0,1.0,0.0,1.0,0.0,-1.0,0.0,0.0,
     -1.0,0.0,1.0,0.0,1.0,0.0,-2.0,-2.0,-2.0,2.0,0.0,0.0,2.0,2.0,
     2.0,2.0,0.0,0.0,2.0,-2.0,-2.0,-2.0,0], 2, 3, 6)

grads = zeros(2, 3, 6)

for f=1:6
    for p=1:3
        grads[1, p, f] = a11 * gq[1, p, f] + a12 * gq[2, p, f]
        grads[2, p, f] = a21 * gq[1, p, f] + a22 * gq[2, p, f]
    end
end

result = zeros(21)

nb_runs = 1000

times = zeros(nb_runs)

for irun in 1:nb_runs
    times[irun] = @elapsed compute(x, y, grads, result)
end

time_without_simd = median(times)

@printf("Julia without simd %.3f µs\n", time_without_simd * 1e6)

for irun in 1:nb_runs
    times[irun] = @elapsed compute_simd(x, y, grads, result)
end

time_with_simd = median(times)

@printf("Julia with simd    %.3f µs\n", time_with_simd * 1e6)

@printf("ratio (without simd)/(with simd): %.2f\n", time_without_simd / time_with_simd)


text:
@testset "Density Util" begin
    Random.seed!(0)
    d, l = hcat(randn(5, 95), randn(5, 5) .+ 2), fill(:inlier, 100)
    l[end-4:end] .= :outlier
    gamma = 0.5
    c = KDECache(d, gamma)

    @testset "Threshold Strategies" begin
        @testset "FixedThresholdStrategy" begin
            n = 3
            s = FixedThresholdStrategy(n)
            @test n == calculate_threshold(s, d, l, gamma)
            @test n == calculate_threshold(s, c, l)
        end

        for x in [NumLabelTresholdStrategy, MaxDensityThresholdStrategy, OutlierPercentageThresholdStrategy]
            @testset "$x" begin
                @test_throws ArgumentError x(-1)
                @test_throws ArgumentError x(2)
                s = x(0.1)
                t, i, o = split_masks(s, c, l)
                @test t > 0
                @test all(xor.(i, o))
                @test any(o)
            end
        end

        for x in [OutlierPercentageThresholdStrategy, GroundTruthThresholdStrategy]
            @testset "$x" begin
                s = x()
                t, i, o = split_masks(s, c, l)
                @test t > 0
                @test all(xor.(i, o))
                @test any(o)
            end
        end
    end

    @testset "PWC" begin
        pwc = PWC(d, gamma)
        @test_throws ArgumentError predict(pwc, d)
        calculate_threshold!(pwc, l)
        pred = predict(pwc, d)
        @test length(pred) == size(d, 2)
        @test !all(pred .== :inlier)
    end
end


text:
struct Affine <: AbstractPredictorCorrectorAlgorithm{false} end

"""
    predict(alg::Affine, H, J_H, ∂H∂t, x, t, Δt)

Tangent predictor. `dx_xt` and `dt_xt` are the evaluations of partial derivatives from our homotopy at x and t, i.e.
``D_xH(x,t)`` and ``D_tH(x,t)``.
It is based on the explicit euler method:
```math
x(t)−∆t·D_xH^{-1}(x,t)·D_tH(x,t)
```

We have ``v = D_xH^{-1}(x,t)·D_tH(x,t)``,
therefore ``D_xH(x,t)·v = D_tH(x,t)``
and thus ``v = \\(D_xH(x,t), D_tH(x,t)``
"""
function predict(alg::Affine, H::AbstractHomotopy, J_H, ∂H∂t, x, t, Δt)
    x .- Δt .* \(J_H(x,t), ∂H∂t(x,t))
end

function correct!(
    u::Vector{T},
    alg::Affine,
    H::AbstractHomotopy,
    J_H,
    x::Vector{T},
    t,
    tol::Float64,
    max_iterations::Int
) where {T}
    N = length(x)
    res = zeros(x)
    dx_xt = zeros(T, N, N)
    u .= x
    for _ in 1:max_iterations
        res .= evaluate(H, u, t)
        if norm(res) < tol
            return true
        end
        dx_xt .= J_H(u, t)
        u .= u .- \(dx_xt, res)
    end

    return norm(res) < tol
end


text:
struct DiffuseLight <: Material
    emit::Texture
end

function scatter(light::DiffuseLight, ray::Ray, hp::HitPoint)::Tuple{Bool, Union{Nothing, AbstractVec}, Union{Nothing, Ray}}
    false, nothing, nothing
end

function emit(light::DiffuseLight, u::AbstractFloat, v::AbstractFloat, point::AbstractVec)::AbstractVec
    value(light.emit, u, v, point)
end

max_stars_repo_name: goooxu/JuliaRayTracing

max_issues_repo_name: goooxu/JuliaRayTracing

max_forks_repo_name: goooxu/JuliaRayTracing


text:
struct ZIBurrExpert{T<:Real} <: ZIContinuousExpert
    p::T
    k::T
    c::T
    λ::T
    ZIBurrExpert{T}(p::T, k::T, c::T, λ::T) where {T<:Real} = new{T}(p, k, c, λ)
end

function ZIBurrExpert(p::T, k::T, c::T, λ::T; check_args=true) where {T <: Real}
    check_args && @check_args(ZIBurrExpert, 0 <= p <= 1 && k >= zero(k) && c >= zero(c) && λ > zero(λ))
    return ZIBurrExpert{T}(p, k, c, λ)
end

#### Outer constructors
ZIBurrExpert(p::Real, k::Real, c::Real, λ::Real) = ZIBurrExpert(promote(p, k, c, λ)...)
ZIBurrExpert(p::Integer, k::Integer, c::Integer, λ::Integer) = ZIBurrExpert(float(p), float(k), float(c), float(λ))
ZIBurrExpert() = ZIBurrExpert(0.50, 1.0, 1.0, 1.0)

## Conversion
function convert(::Type{ZIBurrExpert{T}}, p::S, k::S, θ::S) where {T <: Real, S <: Real}
    ZIBurrExpert(T(p), T(k), T(c), T(λ))
end
function convert(::Type{ZIBurrExpert{T}}, d::ZIBurrExpert{S}) where {T <: Real, S <: Real}
    ZIBurrExpert(T(d.p), T(d.k), T(d.c), T(d.λ), check_args=false)
end
copy(d::ZIBurrExpert) = ZIBurrExpert(d.p, d.k, d.c, d.λ, check_args=false)

## Loglikelihood of Expoert
logpdf(d::ZIBurrExpert, x...) = Distributions.logpdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)
pdf(d::ZIBurrExpert, x...) = Distributions.pdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)
logcdf(d::ZIBurrExpert, x...) = Distributions.logcdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)
cdf(d::ZIBurrExpert, x...) = Distributions.cdf.(LRMoE.Burr(d.k, d.c, d.λ), x...)

## expert_ll, etc
expert_ll_exact(d::ZIBurrExpert, x::Real) = (x == 0.) ? log(p_zero(d)) :  log(1-p_zero(d)) + LRMoE.logpdf(d, x)
function expert_ll(d::ZIBurrExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_ll_pos = LRMoE.expert_ll(LRMoE.BurrExpert(d.k, d.c, d.λ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_ll = (yl == 0.) ? log.(p0 + (1-p0)*exp.(expert_ll_pos)) : log.(0.0 + (1-p0)*exp.(expert_ll_pos))
    expert_ll = (tu == 0.) ? log.(p0) : expert_ll
    return expert_ll
end
function expert_tn(d::ZIBurrExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_pos = LRMoE.expert_tn(LRMoE.BurrExpert(d.k, d.c, d.λ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_tn = (tl == 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_pos))
    expert_tn = (tu == 0.) ? log.(p0) : expert_tn
    return expert_tn
end
function expert_tn_bar(d::ZIBurrExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    expert_tn_bar_pos = LRMoE.expert_tn_bar(LRMoE.BurrExpert(d.k, d.c, d.λ), tl, yl, yu, tu)
    # Deal with zero inflation
    p0 = p_zero(d)
    expert_tn_bar = (tl > 0.) ? log.(p0 + (1-p0)*exp.(expert_tn_bar_pos)) : log.(0.0 + (1-p0)*exp.(expert_tn_bar_pos))
    return expert_tn_bar
end

exposurize_expert(d::ZIBurrExpert; exposure = 1) = d

## Parameters
params(d::ZIBurrExpert) = (d.p, d.k, d.c, d.λ)
p_zero(d::ZIBurrExpert) = d.p
function params_init(y, d::ZIBurrExpert)
    p_init = sum(y .== 0.0) / sum(y .>= 0.0)
    pos_idx = (y .> 0.0)

    c_init, k_init, θ_init = params(params_init(y[pos_idx], BurrExpert()))
    try
        return ZIBurrExpert(p_init, c_init, k_init, θ_init)
    catch;
        ZIBurrExpert()
    end
end

## KS stats for parameter initialization
function ks_distance(y, d::ZIBurrExpert)
    p_zero = sum(y .== 0.0) / sum(y .>= 0.0)
    return max(abs(p_zero-d.p), (1-d.p)*HypothesisTests.ksstats(y[y .> 0.0], LRMoE.Burr(d.k, d.c, d.λ))[2])
end

## Simululation
sim_expert(d::ZIBurrExpert) = (1 .- Distributions.rand(Distributions.Bernoulli(d.p), 1)[1]) .* Distributions.rand(LRMoE.Burr(d.k, d.c, d.λ), 1)[1]

## penalty
penalty_init(d::ZIBurrExpert) = [2.0 10.0 2.0 10.0 2.0 10.0]
no_penalty_init(d::ZIBurrExpert) = [1.0 Inf 1.0 Inf 1.0 Inf]
penalize(d::ZIBurrExpert, p) = (p[1]-1)*log(d.k) - d.k/p[2] + (p[3]-1)*log(d.c) - d.c/p[4] + (p[5]-1)*log(d.λ) - d.λ/p[6]

## statistics
mean(d::ZIBurrExpert) = (1-d.p)*mean(LRMoE.Burr(d.k, d.c, d.λ))
var(d::ZIBurrExpert) = (1-d.p)*var(LRMoE.Burr(d.k, d.c, d.λ)) + d.p*(1-d.p)*(mean(LRMoE.Burr(d.k, d.c, d.λ)))^2
quantile(d::ZIBurrExpert, p) = p <= d.p ? 0.0 : quantile(LRMoE.Burr(d.k, d.c, d.λ), p-d.p)
lev(d::ZIBurrExpert, u) = (1-d.p)*lev(BurrExpert(d.k, d.c, d.λ), u)
excess(d::ZIBurrExpert, u) = mean(d) - lev(d, u)

## EM: M-Step
function EM_M_expert(d::ZIBurrExpert,
                     tl, yl, yu, tu,
                     exposure,
                    #  expert_ll_pos,
                    #  expert_tn_pos,
                    #  expert_tn_bar_pos,
                     z_e_obs, z_e_lat, k_e;
                     penalty = true, pen_pararms_jk = [1.0 Inf 1.0 Inf])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Update zero probability
    expert_ll_pos = expert_ll.(LRMoE.BurrExpert(d.k, d.c, d.λ), tl, yl, yu, tu)
    expert_tn_bar_pos = expert_tn_bar.(LRMoE.BurrExpert(d.k, d.c, d.λ), tl, yl, yu, tu)

    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(yl, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    z_zero_e_lat = z_e_lat .* EM_E_z_zero_lat(tl, p_old, expert_tn_bar_pos)
    z_pos_e_lat = z_e_lat .- z_zero_e_lat
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, z_zero_e_lat, z_pos_e_lat, k_e)

    # Update parameters: call its positive part
    tmp_exp = BurrExpert(d.k, d.c, d.λ)
    tmp_update = EM_M_expert(tmp_exp,
                            tl, yl, yu, tu,
                            exposure,
                            # expert_ll_pos,
                            # expert_tn_pos,
                            # expert_tn_bar_pos,
                            # z_e_obs, z_e_lat, k_e,
                            z_pos_e_obs, z_pos_e_lat, k_e,
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIBurrExpert(p_new, tmp_update.k, tmp_update.c, tmp_update.λ)
end

## EM: M-Step, exact observations
function EM_M_expert_exact(d::ZIBurrExpert,
                     ye, exposure,
                     z_e_obs;
                     penalty = true, pen_pararms_jk = [Inf 1.0 Inf])

    # Old parameters
    p_old = p_zero(d)

    if p_old > 0.999999
        return d
    end

    # Update zero probability
    expert_ll_pos = expert_ll_exact.(LRMoE.BurrExpert(d.k, d.c, d.λ), ye)

    z_zero_e_obs = z_e_obs .* EM_E_z_zero_obs(ye, p_old, expert_ll_pos)
    z_pos_e_obs = z_e_obs .- z_zero_e_obs
    p_new = EM_M_zero(z_zero_e_obs, z_pos_e_obs, 0.0, 0.0, 0.0)

    # Update parameters: call its positive part
    tmp_exp = BurrExpert(d.k, d.c, d.λ)
    tmp_update = EM_M_expert_exact(tmp_exp,
                            ye, exposure,
                            z_pos_e_obs;
                            penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    return ZIBurrExpert(p_new, tmp_update.k, tmp_update.c, tmp_update.λ)
end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
# Autogenerated wrapper script for startin_jll for x86_64-w64-mingw32
export libstartin

JLLWrappers.@generate_wrapper_header("startin")
JLLWrappers.@declare_library_product(libstartin, "libstartin.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libstartin,
        "bin\\libstartin.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# # Test QCDNUM evolution options

# We can use [`QCDNUM.evolfg`](@ref) to evolve the PDF set for all flavours,
# or [`QCDNUM.evsgns`](@ref) to evolve an arbitrary set of singlet/non-singlet pdfs.
# The latter may offer potential for external parallelisation, but first we must verify
# that we can use it to recreate expected results.
#
# We do this by recreating the `example.f` demo using [`QCDNUM.evsgns`](@ref).
# This is based on the `testsgns.f` test job.

using QCDNUM
using Printf

# Let's start by defining some inputs

xmin = Float64.([1.0e-5, 0.2, 0.4, 0.6, 0.75])
iwt = Int32.([1, 2, 4, 8, 16])
nxsubg = 5
nxin = 100
iosp = 3

qq = Float64.([2e0, 1e4])
wt = Float64.([1e0, 1e0])
nqin = 60
itype = 1
as0 = 0.364
r20 = 2.0

q2c = 3.0
q2b = 25.0
q2t = 1e11
q0 = 2.0

def = Float64.([0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,
                0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,
                0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
                0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
                0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                0.,-1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,
                0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
                -1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. ])

isns =  Int32.([1,  2,  2,  2,  2,  2, -1, -2, -2, -2, -2, -2])

nfix = 6
aar = 1.0
bbr = 0.0;

# Now we can state our input PDF function:

function func1(ipdf, x)::Float64
    i = ipdf[]
    xb = x[]
    adbar = 0.1939875
    f = 0
    if (i == 0)
        ag = 1.7
        f = ag * xb^-0.1 * (1.0-xb)^5.0
    end
    if (i == 1)
        ad = 3.064320
        f = ad * xb^0.8 * (1.0-xb)^4.0
    end
    if (i == 2)
        au = 5.107200
        f = au * xb^0.8 * (1.0-xb)^3.0
    end
    if (i == 3)
        f = 0.0
    end
    if (i == 4)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0
    end
    if (i == 5)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
    end
    if (i == 6)
        xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
        xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
        f = 0.2 * (xdbar + xubar)
    end
    if (i == 7)
        f = 0.0
    end
    if (i == 8)
        xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
        xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
        f = 0.2 * (xdbar + xubar)
    end
    if (i == 9)
        f = 0.0
    end
    if (i == 10)
        xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
        xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
        f = 0.2 * (xdbar + xubar)
    end
    if (i == 11)
        f = 0.0
    end
    if (i == 12)
        xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
        xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
        f = 0.2 * (xdbar + xubar)
    end
    return f
end

func1_c = @cfunction(func1, Float64, (Ref{Int32}, Ref{Float64}));

# We can also make a function that passes an input PDF that is
# an output of another PDF. This will be useful later on.

function func2(ipdf, x)::Float64
    i = ipdf[]
    xb = x[]

    f = 0.0

    if (i == -1)
        iset = Int(QCDNUM.qstore("read", 1))
        iq0 = Int(QCDNUM.qstore("read", 2))
    else
        iset = Int(QCDNUM.qstore("read", 1))
        iq0 = Int(QCDNUM.qstore("read", 2))
        ix = QCDNUM.ixfrmx(xb)
        f = QCDNUM.bvalij(iset, i, ix, iq0, 1)
    end

    return f
end

func2_c = @cfunction(func2, Float64, (Ref{Int32}, Ref{Float64}));

# Now we can set up QCDNUM and print a summary of our inputs.

QCDNUM.qcinit(-6, " ")
nx = QCDNUM.gxmake(xmin, iwt, nxsubg, nxin, iosp) # make x grid
nq = QCDNUM.gqmake(qq, wt, 2, nqin) # make qq grid
nw = QCDNUM.fillwt(itype) # fill weight tables
QCDNUM.setord(3) # NLO in pQCD
QCDNUM.setalf(as0, r20) # alpha and scale
iqc = QCDNUM.iqfrmq(q2c)
iqb = QCDNUM.iqfrmq(q2b)
iqt = QCDNUM.iqfrmq(q2t)
QCDNUM.setcbt(nfix, iqc, iqb, iqt) # FFNS and thresholds
QCDNUM.setabr(aar, bbr) # Relationship between mu_F^2 and mu_R^2
iq0 = QCDNUM.iqfrmq(q0)

# Print summary
@printf("xmin = %0.2e, \t subgrids = %i\n", xmin[1], nxsubg)
@printf("qmin = %0.2e, \t qmax = %0.2e\n", qq[1], qq[2])
@printf("nfix = %i, \t\t q2c, b, t = %0.2e, %0.2e, %0.2e\n", nfix, q2c, q2b, q2t)
@printf("alpha = %0.2e, \t r20 = %0.2e\n", as0, r20)
@printf("aar = %0.2e, \t bbr = %0.2e\n", aar, bbr)
@printf("q20 = %0.2e\n", q0)

# Evolve with [`QCDNUM.evolfg`](@ref)
iset1 = 1
jtype = 10*iset1+itype
eps = QCDNUM.evolfg(jtype, func1_c, def, iq0)
QCDNUM.qstore("write", 1, float(iset1))
QCDNUM.qstore("write", 2, float(iq0))

# Evolve with [`QCDNUM.evsgns`](@ref)
iset2 = 2
jtype = 10*iset2+itype
QCDNUM.setint("edbg", 0) # turn off debug
eps = QCDNUM.evsgns(jtype, func2_c, isns, 12, iq0);

# Now we can compare the results from each case,
# by first defining a useful function

function compare(jset1, jset2, id, jq, ixmin)

    nx, xmi, xma, nq, qmi, qma, iord = QCDNUM.grpars()

    if (jq == 0)
        iq1 = 1
        iq2 = nq
        jxmin = 0
    else
        iq1 = jq
        iq2 = jq
        jxmin = ixmin
    end

    dif = 0.0
    for iq in [iq1, iq2]
        for ix in 1:nx
            val1 = QCDNUM.bvalij(jset1, id, ix, iq, 1)
            val2 = QCDNUM.bvalij(jset2, id, ix, iq, 1)
            difi = val1-val2
            dif = max(difi, abs(difi))
            if (ix < jxmin)
                println("(2I3,2F10.4,E13.3,2X,3I3)")
            end
        end
    end

    if (jxmin == 0)
        @printf("id = %i, \t dif = %0.5e\n", id, dif)
    end

end

for id in 0:12
    compare(iset1, iset2, id, 0, 1)
end


text:
using MusicXML
cd(@__DIR__)
@MX begin
# Reads musicxml file and then extracts the data, builds all the types and stores them in proper format.
scorepartwise = readmusicxml( "musescore.musicxml")

# Example 1:
# Prints Each instrument name and then the pitches

# Extracting each instrument information
scoreparts = scorepartwise.partlist.scoreparts

# Extracting parts
parts = scorepartwise.parts

# Extracting each part
for part in parts

    ind = Base.findfirst(x -> part.id == x.id, scoreparts) # returns the index of scorepart that matches the id of part

    # printing the instrument name
    println(scoreparts[ind].name)

    # Extracting each measure of the part
    iMsr=1
    for measure in part.measures
        # Extracting notes of each measure
        for note in measure.notes
            if !isnothing(note.pitch)

                println("Measure no. $iMsr")    # print measure number
                println(note.pitch)     # print pitch of the note
                println(note.duration)  # print duration of the note
            elseif !isnothing(note.unpitched)
                println("Measure no. $iMsr")    # print measure number
                println(note.unpitched) # print unpitched of the note
                println(note.duration)  # print duration of the note
            end

        end

        iMsr+=1
    end
end

end


text:
module collection

export
    KeyStore, connect!, destroy!, watch, unwatch

using Dates
using Printf

import JSON
import MsgPack

using ..session
using ..api
using ..api._storage

function _serialize_bytes(x)
    io = IOBuffer()
    serialize(io, x)
    take!(io)
end
_deserialize_bytes(x) = deserialize(IOBuffer(x))

# key serialiser/deserialiser pairs
const KEY_FORMAT_MAP = Dict{Symbol, Tuple{Function, Function}}(
    :json => (JSON.json, JSON.parse),
    :string => (string, identity)
)

# value serialiser/deserialiser pairs
const VAL_FORMAT_MAP = Dict{Symbol, Tuple{Function, Function}}(
    :json => (JSON.json, JSON.parse),
    :string => (string, identity),
    :data => (identity, identity),
    :julia => (_serialize_bytes, _deserialize_bytes),
    :msgpack => (MsgPack.pack, MsgPack.unpack)
)

"""
High-level container wrapping a Google Storage bucket
"""
struct KeyStore{K, V} <: AbstractDict{K, V}
    bucket_name::String
    session::GoogleSession
    key_decoder::Function
    key_encoder::Function
    reader::Function
    writer::Function
    gzip::Bool
    channel::Dict{Symbol, Any}
end

function KeyStore{K,V}(bucket_name::AbstractString; session::GoogleSession=get_session(storage),
                location::AbstractString="US", empty::Bool=false, gzip::Bool=true,
                key_format::Union{Symbol, AbstractString}=K <: String ? :string : :json,
                val_format::Union{Symbol, AbstractString}=V <: String ? :string : :json,
                debug=false) where {K,V}

    key_encoder, key_decoder = try KEY_FORMAT_MAP[Symbol(key_format)] catch
        error("Unknown key format: $key_format")
    end
    writer, reader = try VAL_FORMAT_MAP[Symbol(val_format)] catch
        error("Unknown value format: $val_format")
    end

    store = KeyStore{K,V}(bucket_name, session,
        (x) -> convert(K, key_decoder(x)), key_encoder,
        reader, writer,
        gzip, Dict{Symbol, Any}()
    )
    # establish availability of bucket
    connect!(store; location=location, empty=empty, debug=debug)
    store
end

function connect!(store::KeyStore; location::AbstractString="US",
                  empty::Bool=false, debug=false)
    response = storage(:Bucket, :get, store.bucket_name;
                       session=store.session, fields="", debug=debug)
    if iserror(response)
        code = response[:error][:code]
        if code == 404  # not found (available)
            response = storage(:Bucket, :insert; session=store.session,
                               data=Dict(:name => store.bucket_name, :location => location),
                               fields="", debug=debug)
            if iserror(response)
                error("Unable to create bucket: $(response[:error][:message])")
            end
        elseif code == 403  # forbidden (not available)
            error("Authorization failure or bucket name already taken: ",
                  response[:error][:message])
        else
            error("Error checking bucket: $(response[:error][:message])")
        end
    elseif empty
        empty!(store)
    end
    store
end

function destroy!(store::KeyStore{K, V}) where {K, V}
    response = storage(:Bucket, :delete, store.bucket_name; session=store.session, fields="")
    if iserror(response)
        error("Unable to delete bucket: $(response[:error][:message])")
    end
    nothing
end

function Base.print(io::IO, store::KeyStore{K, V}) where {K, V}
    print(io, @sprintf("""KeyStore{%s, %s}("%s")""", K, V, store.bucket_name))
end
Base.show(io::IO, store::KeyStore) = print(io, store)

function Base.setindex!(store::KeyStore{K, V}, val::V, key::K) where {K, V}
    name = store.key_encoder(key)
    data = store.writer(val)
    response = storage(:Object, :insert, store.bucket_name; session=store.session,
        name=name, data=data, gzip=store.gzip, content_type="application/octet-stream", fields=""
    )
    if iserror(response)
        error("Unable to set key '$key': $(response[:error][:message])")
    end
    store
end

function Base.getindex(store::KeyStore{K, V}, key::K) where {K, V}
    name = store.key_encoder(key)
    data = storage(:Object, :get, store.bucket_name, name; session=store.session)
    if iserror(data)
        throw(KeyError(key))
    end
    val = store.reader(data)
    convert(V, val)
end

function Base.get(store::KeyStore{K, V}, key::K, default) where {K, V}
    try
        return getindex(store, key)
    catch e
        if isa(e, KeyError)
            return default
        else
            throw(e)
        end
    end
end

function Base.haskey(store::KeyStore{K, V}, key::K) where {K, V}
    name = store.key_encoder(key)
    response = storage(:Object, :get, store.bucket_name, name; session=store.session, alt="", fields="")
    !iserror(response)
end

# WARNING: potential for race condition. don't zip keys with values... use collect instead
function Base.keys(store::KeyStore{K, V}) where {K, V}
    result = K[]
    for response in storage(:Object, :list, store.bucket_name; session=store.session, fields="items(name)")
        name = response[:name]
        key = store.key_decoder(name)
        if !isa(key, K)
            throw(TypeError(:keys, "$(store.bucket_name):$name", K, key))
        end
        push!(result, key)
    end
    result
end

# avoiding race condition where values might have been deleted after keys were generated
@inline function Base.values(store::KeyStore{K, V}) where {K, V}
    (x for x in (get(store, key, Nothing) for key in keys(store)) if x !== Voide)
end

function Base.delete!(store::KeyStore{K, V}, key::K) where {K, V}
    name = store.key_encoder(key)
    response = storage(:Object, :delete, store.bucket_name, name; session=store.session, fields="")
    if iserror(response)
        error("Unable to delete object: $(response[:error][:message])")
    end
    store
end

function Base.pop!(store::KeyStore{K, V}, key::K) where {K, V}
    val = getindex(store, key)
    delete!(store, key)
    val
end

function Base.pop!(store::KeyStore{K, V}, key::K, default) where {K, V}
    val = get(store, key, default)
    delete!(store, key)
    val
end

function Base.merge!(store::KeyStore{K, V}, d::AbstractDict{K, V}) where {K, V}
    for (k, v) in d
        store[k] = v
    end
end

Base.length(store::KeyStore) = length(collect(keys(store)))

Base.isempty(store::KeyStore) = length(store) == 0

function Base.empty!(store::KeyStore{K, V}) where {K, V}
    for object in storage(:Object, :list, store.bucket_name; session=store.session, fields="items(name)")
        response = storage(:Object, :delete, store.bucket_name, object[:name]; session=store.session, fields="")
        if iserror(response)
            error("Unable to delete object: $(response[:error][:message])")
        end
    end
    store
end

"""Skip over missing keys if any deleted since key list was geenrated"""
function fast_forward(store::KeyStore{K, V}, key_list) where {K, V}
    while !isempty(key_list)
        key = pop!(key_list)
        val = get(store, key, Nothing)
        if val !== Nothing
            return Pair{K, V}(key, val)
        end
    end
    nothing
end

#function Base.iterate(store::KeyStore{K, V}) where {K, V}
#
#end
#function Base.start(store::KeyStore{K, V}) where {K, V}
#    key_list = collect(keys(store))
#    return (fast_forward(store, key_list), key_list)
#end
#
#function Base.next(store::KeyStore{K, V}, state) where {K, V}
#    pair, key_list = state
#    return pair, (fast_forward(store, key_list), key_list)
#end
#
#function Base.done(store::KeyStore{K, V}, state) where {K, V}
#    pair, key_list = state
#    pair === nothing
#end
#
@inline function Base.IteratorSize(::Type{KeyStore{K, V}}) where {K, V}
    Base.SizeUnknown()
end

# notifications
function watch(store::KeyStore{K, V}, channel_id::AbstractString,
               address::AbstractString) where {K, V}
    if !isempty(store.channel)
        error("Already watching: $store.channel")
    end
    channel = storage(:Object, :watchAll, store.bucket_name;
        data=Dict(:type => "WEBHOOK", :address => address, :id => channel_id),
        session=store.session
    )
    if iserror(channel)
        error("Unable to watch bucket: $(channel[:error][:message])")
    end
    store.channel = channel
end

function unwatch(store::KeyStore{K, V}) where {K, V}
    response = storage(:Channel, :stop; data=store.channel, session=store.session)
    if iserror(response)
        error("Unable to unwatch bucket: $(response[:error][:message])")
    end
    empty!(store.channel)
    store
end

end


text:
module NLmodel

using JuMP
#using AmplNLWriter,
using Ipopt
#using CoinOptServices

function runModel(nodes, measuredNodeStateFull, LB, UB, expression, verbose)
    model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))

    weightRoot = 500
    weightMeasured = 10000000
    weightHard = 10000

    nodesList = collect(keys(nodes))
    measuredNodeState = filter((args)->indexin([args[1]], nodesList) != [0], measuredNodeStateFull)

    @variable(model, LB <= x[1:length(nodesList)] <= UB, start = 1)
    @variable(model, LB <= x_bar[1:length(nodesList)] <= UB, start = 1)
    @variable(model, LB <= p[1:length(nodesList)] <= UB)
    @variable(model, LB <= n[1:length(nodesList)] <= UB)

    numberMeasuredNodes = length(keys(measuredNodeState))

    if numberMeasuredNodes > 0
       @variable(model, m[1:numberMeasuredNodes])
    end

    rootIdxs = []
    variableIdxs = []
    measuredIdxs = indexin(collect(keys(measuredNodeState)), nodesList)
    for nodeName in keys(nodes)
        nodeIdxs = indexin([nodeName], nodesList)
        nodeIndex = nodeIdxs[1]

        measured = false
        j = 0
        for node in collect(keys(measuredNodeState))
            j = j + 1
            if measuredIdxs[j] == nodeIndex && nodes[nodeName].relation == "ROOT"
                measured = true
                rhs = measuredNodeState[node] < LB ? LB : measuredNodeState[node]
                @constraint(model, m[j] == rhs)
                break
            end
        end

        if nodes[nodeName].relation == "ROOT"
            if measured
                @constraint(model, p[nodeIndex] >= x[nodeIndex] - m[j])
                @constraint(model, n[nodeIndex] >= m[j] - x[nodeIndex])
            else
                @constraint(model, p[nodeIndex] >= x[nodeIndex] - 1)
                @constraint(model, n[nodeIndex] >= 1 - x[nodeIndex])
                push!(rootIdxs, nodeIndex)
            end
        else
            if measured
                @constraint(model, p[nodeIndex] >= x[nodeIndex] - m[j])
                @constraint(model, n[nodeIndex] >= m[j] - x[nodeIndex])
            else
                @constraint(model, p[nodeIndex] >= x[nodeIndex] - x_bar[nodeIndex])
                @constraint(model, n[nodeIndex] >= x_bar[nodeIndex] - x[nodeIndex])
                push!(variableIdxs, nodeIndex)
            end

            if nodes[nodeName].relation == "AND"
                parentIndexes = indexin(nodes[nodeName].parents, nodesList)
                if length(parentIndexes) == 1
                    @constraint(model, x[parentIndexes[1]] == x_bar[nodeIndex])
                else
                    @NLconstraint(model, x[parentIndexes[1]] * x[parentIndexes[2]] == x_bar[nodeIndex])
                end
            elseif nodes[nodeName].relation == "NEG"
                parentIndexes = indexin(nodes[nodeName].parents, nodesList)
                @NLconstraint(model, 1 / x[parentIndexes[1]] == x_bar[nodeIndex])
            elseif nodes[nodeName].relation == "ANDNEG"
                parentIndexes = indexin(nodes[nodeName].parents, nodesList)
                @NLconstraint(model, x[parentIndexes[1]] / x[parentIndexes[2]] == x_bar[nodeIndex])
            elseif nodes[nodeName].relation == "OR"
                posParentIdxs = indexin(nodes[nodeName].posParents, nodesList)
                count_expression = 0
                total_expression = 0
                for parent in nodes[nodeName].posParents
                    if haskey(expression, parent) && expression[parent] != ""
                        total_expression += parse(Float64, expression[parent])
                        count_expression += 1
                        ev = expression[parent]
                        if verbose
                            println("Parent node: $parent\tExpression value: $ev")
                        end
                    end
                end
                if verbose
                    println("")
                end

                if count_expression == 0
                    @constraint(model, sum(x[posParentIdxs[a]] for a = 1:length(posParentIdxs)) / length(posParentIdxs) ==  x_bar[nodeIndex])
                else
                    average_expression = total_expression / count_expression
                    evs = []
                    for parent in nodes[nodeName].posParents
                        if haskey(expression, parent) && expression[parent] != ""
                            push!(evs, parse(Float64, expression[parent]))
                        else
                            push!(evs, average_expression);
                        end
                    end

                    total_expression = sum(evs)

                    @constraint(model, sum(evs[a] * x[posParentIdxs[a]] for a = 1:length(posParentIdxs)) / total_expression ==  x_bar[nodeIndex])
                end
           end
        end
    end

    if !isempty(measuredIdxs)
        measuredIdxs = measuredIdxs[measuredIdxs .!= nothing]
    end

    @NLobjective(model,
                 Min,
                 weightHard * sum(p[variableIdxs[i]] + n[variableIdxs[i]] for i = 1:length(variableIdxs))^2
                 + weightMeasured * sum(p[measuredIdxs[j]] + n[measuredIdxs[j]] for j = 1:length(measuredIdxs))^2
                 + weightRoot * sum(p[rootIdxs[k]] + n[rootIdxs[k]] for k = 1:length(rootIdxs))^2)

    if verbose
        println("Solving Model")
        println(model)
    end

    optimize!(model)

    if verbose
        println("Objective value: ", JuMP.objective_value(model))
    end

    results = Dict()
    for i in eachindex(nodesList)
         results[nodesList[i]] = JuMP.value(x[i])
    end

    x_values = Dict()
    x_bar_values = Dict()
    for i in eachindex(nodesList)
        x_values[nodesList[i]] = JuMP.value(x[i])
        x_bar_values[nodesList[i]] = JuMP.value(x_bar[i])
    end

    return [results, x_values, x_bar_values]
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
function twopi_minus(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(6.28318530717958600, hi)
    t2, t3 = sub_2(2.4492935982947064e-16, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 5.9895396194366790e-33
    zlo += t7
    return DoubleFloat{Float64}(zhi, zlo)
end

twopi_minus(x::DoubleFloat{Float64}) = twopi_minus(HI(x), LO(x))

function pi_minus(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(3.14159265358979300, hi)
    t2, t3 = sub_2(1.2246467991473532e-16, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 2.9947698097183397e-33
    zlo += t7
    return DoubleFloat{Float64}(zhi, zlo)
end

pi_minus(x::DoubleFloat{Float64}) = pi_minus(HI(x), LO(x))

function halfpi_minus(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(1.57079632679489660, hi)
    t2, t3 = sub_2(6.1232339957367660e-17, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 1.4973849048591698e-33
    zlo += t7
    return DoubleFloat{Float64}(zhi, zlo)
end

halfpi_minus(x::DoubleFloat{Float64}) = halfpi_minus(HI(x), LO(x))

function qrtrpi_minus(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(0.78539816339744830, hi)
    t2, t3 = sub_2(3.0616169978683830e-17, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 7.4869245242958490e-34
    zlo += t7
    return DoubleFloat{Float64}(zhi, zlo)
end

qrtrpi_minus(x::DoubleFloat{Float64}) = qrtrpi_minus(HI(x), LO(x))

function minus_twopi(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(6.28318530717958600, hi)
    t2, t3 = sub_2(2.4492935982947064e-16, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 5.9895396194366790e-33
    zlo += t7
    return DoubleFloat{Float64}(-zhi, -zlo)
end

minus_twopi(x::DoubleFloat{Float64}) = minus_twopi(HI(x), LO(x))

function minus_pi(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(3.14159265358979300, hi)
    t2, t3 = sub_2(1.2246467991473532e-16, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 2.9947698097183397e-33
    zlo += t7
    return DoubleFloat{Float64}(-zhi, -zlo)
end

minus_pi(x::DoubleFloat{Float64}) = minus_pi(HI(x), LO(x))

function minus_halfpi(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(1.57079632679489660, hi)
    t2, t3 = sub_2(6.1232339957367660e-17, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 1.4973849048591698e-33
    zlo += t7
    return DoubleFloat{Float64}(-zhi, -zlo)
end

minus_halfpi(x::DoubleFloat{Float64}) = minus_halfpi(HI(x), LO(x))

function minus_qrtrpi(hi::T, lo::T) where {T<:Float64}
    zhi, t1 = sub_2(0.78539816339744830, hi)
    t2, t3 = sub_2(3.0616169978683830e-17, lo)
    t7, t4 = add_2(t1, t2)
    t5 = t3 + t4
    zlo = t5 - 7.4869245242958490e-34
    zlo += t7
    return DoubleFloat{Float64}(-zhi, -zlo)
end

minus_qrtrpi(x::DoubleFloat{Float64}) = minus_qrtrpi(HI(x), LO(x))


text:
"""
    x,w,T = nfieldtrans(M,n;basis="Hermite",ω=1.0,α=0.0)

Construct transforms and associated arrays for numerical quadrature evaluation of n-field integrals,
starting from a representation of the quantum state with respect to a particular basis of eigenstates.

# Arguments
`M`: number of modes in the c-field.

`n`: order of the field product.

### keyword arguments

`basis`: name of orthonormal eigenstate basis. Currently `"Hermite"` is implemented.

`ω`: frequency of the oscillator states in given direction, relative to chosen reference frequency.

`α`: extra parameter for Laguerre basis.

### Outputs
`T`: linear transformation matrix that takes spectral state coefficienets to a quadrature grid. In 1D:

```math
Tc = \\psi(x) \\equiv \\sum_{j=1}^{M}c_j\\phi_j(x)
```

for a state represented by `M` coefficients, the number of modes in the c-field.

 * `x` is the quadrature grid onto which ``\\psi(x)`` is mapped
 * `w` are weights such that the an exact integral may be carreid out.

The integral must be a product of order `n` in the field ``\\psi``.

The integral is performed by

1. Transforming to the quadrature grid using `T`

2. Constructing the product and then evaluating the sum:
```math
\\sum_jw_j*\\psi(x_j)^n = sum(w.*\\psi.^n)
```

## Examples

### C-field population
Compute the number of particles in the C-field, for a state of `M` modes:

```julia
M = 30
c = randn(M)+im*randn(M)
x,w,T = nfieldtrans(M,2)
ψ = T*c
N = sum(w.*abs2.(ψ))
73.24196674113007
```

Computes the integral ``N = \\int dx |\\psi(x)|^2`` as may be checked by direct summation:

```julia
sum(abs2.(c))
73.24196675017353
```

Relative error of numerical integrals will usually be smaller than ~1e-10.

### Interaction energy
The most common integral of this type involves a four-field product, `n=4`.
The PGPE interaction energy is of this form, as is the nonlinear term in the PGPE;
the exact propagation of the PGPE requires repeated use of this 4-field transformation:

```julia
x,w,T = nfieldtrans(M,4)
ψ = T*c
Uint = sum(w.*abs.(ψ).^4)
552.9762736751692
```
computes the integral ``U_{\\rm int}\\equiv \\int dx|\\psi(x)|^4`` to accuracy very close to working precision.

"""

function nfieldtrans(M,n;basis="Hermite",ω=1.0,α=0.0)
    iseven(n*M) ? K=Int(n*M/2) : K=Int((n*M+1)/2)
    if basis=="Hermite"
    x, w = gausshermite(K)
    w    = @. exp(log(w)+x^2)/sqrt(n*ω/2)
    T    = eigmat(M,x/sqrt(n*ω/2),ω=ω,basis=basis)
    elseif basis=="Laguerre"
      error(basis," ","basis not implemented yet.")
#==
      #sort out stable recursion for Laguerre
      M>26 && error("recursion unstable for M >26")
      #needs proper testing
      x, w = gausslaguerre(n)
      w    = w.*exp.(x)/(K*ω/2)./x.^α #careful wieght check needed for order K
      T    = eigmat(M,x/(K*ω/2),ω=ω,basis=basis,α=0.0)
      ==#
    else
        error(basis," ","basis not implemeted yet.")
    end
    return x,w,T
end


text:
include("../imports.jl")

@testset ExtendedTestSet "`ProjectiveTransforms`" begin
    @testset ExtendedTestSet "image" begin
        encoding = ProjectiveTransforms((32, 32))
        image = rand(RGB, 64, 96)
        block = FastAI.Image{2}()

        ## We run `ProjectiveTransforms` in the different [`Context`]s:
        imagetrain = encode(encoding, Training(), block, image)
        @test size(imagetrain) == (32, 32)

        imagevalid = encode(encoding, Validation(), block, image)
        @test size(imagevalid) == (32, 32)

        imageinference = encode(encoding, Inference(), block, image)
        @test size(imageinference) == (32, 48)

        ## During inference, the aspect ratio should stay the same
        @test size(image, 1) / size(image, 2) == size(imageinference, 1) / size(imageinference, 2)
    end

    @testset ExtendedTestSet "keypoints" begin
        encoding = ProjectiveTransforms((32, 48))
        ks = [SVector(0., 0), SVector(64, 96)]
        block = FastAI.Keypoints{2}(10)
        bounds = DataAugmentation.Bounds((1:32, 1:48))
        r = DataAugmentation.getrandstate(encoding.tfms[Training()])
        kstrain = encode(encoding, Training(), block, ks; state = (bounds, r))
        ksvalid = encode(encoding, Validation(), block, ks; state = (bounds, r))
        ksinference = encode(encoding, Inference(), block, ks; state = (bounds, r))

        @test kstrain == ksvalid == ksinference
    end

    @testset ExtendedTestSet "image and keypoints" begin
        encoding = ProjectiveTransforms((32, 32))
        image = rand(RGB, 64, 96)
        ks = [SVector(0., 0), SVector(64, 96)]
        blocks = (FastAI.Image{2}(), FastAI.Keypoints{2}(10))

        @test_nowarn encode(encoding, Training(), blocks, (image, ks))
        @test_nowarn encode(encoding, Validation(), blocks, (image, ks))
        @test_nowarn encode(encoding, Inference(), blocks, (image, ks))
    end

    #= depends on buffered interface
    @testset ExtendedTestSet "`run!`" begin
        encoding = ProjectiveTransforms((32, 32))
        image1 = rand(RGB, 64, 96)
        image2 = rand(RGB, 64, 96)
        buf = encode(encoding, Validation(), image1)
        cbuf = copy(buf)
        encode!(buf, encoding, Validation(), image1)
        # run and run! should return the same result
        @test buf ≈ cbuf

        encode!(buf, encoding, Validation(), image2)
        # buf should be modified on different input
        @test !(buf ≈ cbuf)
    end
    =#
end


text:
module Test

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

using Compat
using Compat.Test

include("config.jl")

include("modellike.jl")

include("contlinear.jl")
include("contconic.jl")
include("contquadratic.jl")

include("intlinear.jl")
include("intconic.jl")

include("nlp.jl")

include("UnitTests/unit_tests.jl")

end # module


text:
"iterator type over simplicial multi-indices"
struct MultiIndexSet
    d::Int # dimension
    N::Int # degree
end

Base.length(a::MultiIndexSet) = prod(ntuple(i->(a.N+i),a.d)) ÷ factorial(a.d)

function Base.iterate(set::MultiIndexSet)
    α = ntuple(x->zero(Int),set.d)
    return α,α
end
function Base.iterate(set::MultiIndexSet,α::NTuple)
    if last(α) == set.N
        return nothing
    else
        αnext = increment(α,set,1)
        return αnext,αnext
    end
end

"function increment(I::NTuple,set::MultiIndexSet,i::Int):
    increments a multi-index of arbitrary dimension and degree.
    Intended to be used within the interface for Base.iterate"
function increment(I::NTuple,set::MultiIndexSet,i::Int)
    if length(I) != set.d
        error(DimensionMismatch("Length of I = $(length(I)) while set.d = $(set.d)"))
    end
    if last(I)==set.N
        return (I[1:end-1]...,last(I)+1) # set "incorrect" state to end iteration
    end
    if sum(I) < set.N
        return (I[1:i-1]...,I[i]+1,I[i+1:end]...)
    else
        Inext = (zero.(I[1:i])...,I[i+1:end]...)
        increment(Inext,set,i+1) # recursive
    end
end


text:
using NestedEnvironments
using DifferentialEquations

struct IntegralRLEnv <: AbstractEnv
    env::InputAffineQuadraticCostEnv
    ∫r::IntegralRLEnv
end


text:
# ==============================
#  ======== Tested code ========
using Eirene
using Plots
using StatsPlots

#%%
"""
    get_bettis(results_eirene::Dict, max_dim::Integer; min_dim::Int=1)

Calls Eirene.betticurve for 'dim' in range from `min_dim` up to 'max_dim' and
stack the resulting Arrays into a vector.

The returned value is a Vector of Arrays{Float64,2}. Each array is of size
(n,2), where n is the maximal number of steps taken to compute Betti curve of dimensions
ranging form `min_dim` to `max_dim`. First column of each array contains numbered steps.
Second column are the Betti curve values for corresponding step.

Arrays in returned vector correspond to Betti curve dimensions form range
`min_dim` up to 'max_dim'.

"""
function get_bettis(results_eirene::Dict, max_dim::Integer; min_dim::Int = 1)
    bettis = Matrix{Float64}[]
    for d = min_dim:max_dim
        result = betticurve(results_eirene, dim = d)
        if isempty(result) && d > 1
            result = zeros(size(bettis[d-1]))
        end
        push!(bettis, result)
    end
    return bettis
end
# TODO add get_bettis_from_matrix, to wrap C= eirene...; get bettis

#%%

"""
  normalise_bettis(bettis::Vector)
    normalise_bettis(bettis::Array)

Normalise the number of steps for Betti curves. 'bettis' can be either vector of
arrays (each array contain Betti curve of different dimension) or an array
containing Betti curve of a single dimension.
"""
function normalise_bettis(bettis::Vector)

    @debug "Vector version"
    norm_bettis = copy(bettis)
    @debug "norm_bettis size :" size(norm_bettis)[1][1]

    max_dim = size(norm_bettis)[1]
    @debug "typeof(max_dim) :" typeof(max_dim[1])

    for d = 1:(max_dim)
        if !isempty(norm_bettis[d])
            norm_bettis[d][:, 1] /= findmax(norm_bettis[d][:, 1])[1]
        end
    end
    return norm_bettis

end

#%%
function normalise_bettis(bettis::Array)
    @debug "Array version"
    norm_bettis = copy(bettis)
    @debug "norm_bettis size :" size(norm_bettis)

    if !isempty(norm_bettis)
        norm_bettis[:, 1] /= findmax(norm_bettis[:, 1])[1]
    end
    return norm_bettis
end

#%%
# function vectorize_bettis(betti_curves::Array{Matrix{Float64,2}})
"""
    vectorize_bettis(betti_curves::Matrix{Float64})

Reshapes the 'betti_curves' from type Array{Matrices{Float64,2}} into
Matrix{Float64}.

The resulting matrix size is (n, k), where 'n' is equal to the number of
rows in each matrix, 'k' is equal to the number of matrices.

TODO: Change the name- it takse vector and returns a matrix.
TODO: get bettis could have an arguent betti_type which would determine resulting type
"""
function vectorize_bettis(betti_curves::Vector{Array{Float64,2}})

    first_betti = 1
    last_betti = size(betti_curves,1)
    return hcat([betti_curves[k][:, 2] for k = first_betti:last_betti]...)
end

#%%
@deprecate vectorize_bettis(eirene_results::Dict, maxdim::Integer, mindim::Integer) vectorize_bettis(betti_curves)

# ===
#%%
"""
  get_vectorized_bettis(results_eirene::Dict, max_dim::Integer; min_dim::Int=1)

Takes the eirene result and computes Betti curves for dimensions in range
'mindim:maxdim'. Every Betti curve is stored in successive column of the
resulting array.
TODO: this should return a matrix, where first col are indices and rest are B values (1st col is missing now)
"""
function get_vectorized_bettis(results_eirene::Dict, max_dim::Integer; min_dim::Int = 1)

    all_bettis = get_bettis(results_eirene, max_dim, min_dim = min_dim)
    bettis_vector = vectorize_bettis(all_bettis)

    return bettis_vector
end

# ==
#%%
"""
  plot_bettis(bettis; min_dim::Integer=1, betti_labels::Bool=true, default_labels::Bool=true kwargs...)

Creates a plot for set of betti numbers stored in `bettis` and return the
handler to the plot.

'kwargs' are plot parameters

Some of the possible 'kwargs' are:
  - title::String
  - legend:Bool
  - size::Tuple{T, T} where {T::Number}
  - lw::Integer or linewidth:Integer
(for more, see plots documentation):
TODO: min_dim is not included in all_dims variable
TODO: add change of x label based on x values- so it is either edge density for 0:1 range values or Filtration step otherwise
"""
function plot_bettis(bettis::Vector;
                        min_dim::Integer = 1,
                        use_edge_density::Bool=true,
                        betti_labels::Bool = true,
                        default_labels::Bool = true,
                        kwargs...)#; plot_size = (width=1200, height=800),
    max_dim = size(bettis, 1)
    all_dims = 1:max_dim

    if min_dim > max_dim
        throw(DomainError(
            min_dim,
            "\'min_dim\' must be greater that maximal dimension.",
        ))
    end

    lw_pos = findfirst(x -> x == :lw || x == :linewidth, keys(kwargs))
    if !isnothing(lw_pos)
        lw = kwargs[lw_pos]
    else
        lw = 2
    end

    # Create iterator for all loops
    all_iterations = 1:(max_dim) #TODO ths can not be starting from min_dim, because it may be 0

    if use_edge_density
        for p = all_iterations
            max_step = findmax(bettis[p][:, 1])[1]
            bettis[p][:, 1] ./= max_step
        end
    end

    colors_set = get_bettis_color_palete(min_dim=min_dim)
    plot_ref = plot(; kwargs...)
    # for p = min_dim:(max_dim) #TODO ths can not be starting from min_dim, because it may be 0
    # for p = all_iterations
    for (index, p) in enumerate(min_dim:max_dim)
        args = (lc = colors_set[index], linewidth = lw)
        if betti_labels
            args = (args..., label = "β$(p)")
        end
        plot!(bettis[index][:, 1], bettis[index][:, 2]; args...)
    end

    legend_pos = findfirst(x -> x == :legend, keys(kwargs))
    if !isnothing(legend_pos)
        plot!(legend = kwargs[legend_pos])
    else
        plot!(legend = betti_labels)
    end

    x_pos = findfirst(x -> x == :xlabel, keys(kwargs))
    y_pos = findfirst(x -> x == :ylabel, keys(kwargs))
    if !isnothing(x_pos)
        xlabel!(kwargs[x_pos])
    elseif default_labels
        xlabel!("Edge density")
    end
    if !isnothing(y_pos)
        ylabel!(kwargs[y_pos])
    elseif default_labels
        ylabel!("Number of cycles")
    end

    # set tlims to integer values
    max_ylim = findmax(ceil.(Int, ylims(plot_ref)))[1]
    if max_ylim <=3
        ylims!((0, 3))
    end

    if use_edge_density
        xlims!((0, 1))
    end

    return plot_ref
end

"""
  plot_bettis(bettis::Array; min_dim::Integer=1, betti_labels::Bool=true, default_labels::Bool=true kwargs...)

Creates a plot for set of betti numbers stored in `bettis` and return the
handler to the plot.

'kwargs' are plot parameters

Some of the possible 'kwargs' are:
  - title::String
  - legend:Bool
  - size::Tuple{T, T} where {T::Number}
  - lw::Integer or linewidth:Integer
(for more, see plots documentation):
TODO: add change of x label based on x values- so it is either edge density for 0:1 range values or Filtration step otherwise
"""
function plot_bettis(bettis::Array;
                        # min_dim::Integer = 1,
                        dims_range=1:size(bettis,2),
                        use_edge_density::Bool=true,
                        betti_labels::Bool = true,
                        default_labels::Bool = true,
                        normalised=true,
                        kwargs...)#; plot_size = (width=1200, height=800),
    max_dim = dims_range[end]
    min_dim = dims_range[1]

    if min_dim > max_dim
        throw(DomainError(
            min_dim,
            "\'min_dim\' must be greater that maximal dimension in \'bettis\'",
        ))
    end

    total_steps = size(bettis, 1)
    if normalised
        x_vals = range(0, stop=1, length=total_steps)
    else
        x_vals = range(0, stop=total_steps)
    end

    lw_pos = findfirst(x -> x == :lw || x == :linewidth, keys(kwargs))
    if !isnothing(lw_pos)
        lw = kwargs[lw_pos]
    else
        lw = 2
    end

    # if use_edge_density
    #     # for p = 1:(max_dim) #TODO ths can not be starting from min_dim, because it may be 0
    #     for (index, p) in enumerate(min_dim:max_dim)
    #         max_step = findmax(bettis[:, 1])[1]
    #         bettis[p][:, 1] ./=max_step
    #     end
    # end

    colors_set = get_bettis_color_palete(min_dim=min_dim)
    plot_ref = plot(; kwargs...)
    # for p = min_dim:(max_dim) #TODO ths can not be starting from min_dim, because it may be 0
    # for p = 1:(max_dim) #TODO ths can not be starting from min_dim, because it may be 0
    for (index, p) in enumerate(min_dim:max_dim)
        args = (lc = colors_set[index], linewidth = lw)
        if betti_labels
            args = (args..., label = "β$(p)")
        end
        plot!(x_vals, bettis[:, index]; args...)
    end

    legend_pos = findfirst(x -> x == :legend, keys(kwargs))
    if !isnothing(legend_pos)
        plot!(legend = kwargs[legend_pos])
    else
        plot!(legend = betti_labels)
    end

    x_pos = findfirst(x -> x == :xlabel, keys(kwargs))
    y_pos = findfirst(x -> x == :ylabel, keys(kwargs))
    if !isnothing(x_pos)
        xlabel!(kwargs[x_pos])
    elseif default_labels
        xlabel!("Edge density")
    end
    if !isnothing(y_pos)
        ylabel!(kwargs[y_pos])
    elseif default_labels
        ylabel!("Number of cycles")
    end

    # set tlims to integer values
    max_ylim = findmax(ceil.(Int, ylims(plot_ref)))[1]
    if max_ylim <=3
        ylims!((0, 3))
    end

    if use_edge_density
        xlims!((0, 1))
    end

    return plot_ref
end
#  ======= Untested code
# TODO add default kwargs paring function -> parse_kwargs()

"""
    plot_all_bettis ...
"""
function plot_all_bettis(bettis_collection;
                        min_dim::Integer = 1,
                        betti_labels::Bool = true,
                        default_labels::Bool = true,
                        normalised=true,
                        kwargs...)#; plot_size = (width=1200, height=800),
    total_dims = size(bettis_collection[1],2)

    lw_pos = findfirst(x -> x == :lw || x == :linewidth, keys(kwargs))
    if !isnothing(lw_pos)
        lw = kwargs[lw_pos]
    else
        lw = 2
    end

    colors_set = get_bettis_color_palete(min_dim=min_dim)
    max_y_val = find_max_betti(bettis_collection)

    plot_ref = plot(; kwargs...)
    for b = 1:total_dims
        args = (lc = colors_set[b], linewidth = lw, alpha=0.12,label=false, ylims=(0,max_y_val))
        for bettis = bettis_collection
            betti_vals = bettis[:,b]

            total_steps = size(bettis, 1)
            x_vals = range(0, stop=1, length=total_steps)

            plot!(x_vals, betti_vals; args...)
        end
        # my_label = "β$(b)"
        # betti_vals = results_d["bettis_collection"][:hc][end]
        # x_vals = range(0, stop=1, length=size(betti_vals, 1))
        # plot!(x_vals, betti_vals; lc = colors_set[b], linewidth = 1, alpha=0.1,label=my_label, ylims=(0,max_y_val))
    end
    plot!(legend=true)

    legend_pos = findfirst(x -> x == :legend, keys(kwargs))
    if !isnothing(legend_pos)
        plot!(legend = kwargs[legend_pos])
    else
        plot!(legend = betti_labels)
    end

    x_pos = findfirst(x -> x == :xlabel, keys(kwargs))
    y_pos = findfirst(x -> x == :ylabel, keys(kwargs))
    if !isnothing(x_pos)
        xlabel!(kwargs[x_pos])
    elseif default_labels
        xlabel!("Edge density")
    end
    if !isnothing(y_pos)
        ylabel!(kwargs[y_pos])
    elseif default_labels
        ylabel!("Number of cycles")
    end

    return plot_ref
end

"""
    find_max_betti(bettis_collection::Array)

Returns the highest Betti curve value from all dimensions.
"""
function find_max_betti(bettis_collection::Array)
    if typeof(bettis_collection) == Vector
        bettis_collection = vectorize_bettis(bettis_collection)
    end

    max_y_val = 0
    for betti_set in bettis_collection
        local_max = findmax(betti_set)[1]
        if local_max > max_y_val
            max_y_val = local_max
        end
    end
    return max_y_val
end

#  ======= Untested code == end
#%%

"""
  printready_plot_bettis(kwargs)

Creates a plot using 'plot_bettis' with arguments which were tested to be very
good for using them in prints. Used arguments are:

"""
function printready_plot_bettis(kwargs)
    return nothing
end

#%%
"""
  function get_bettis_color_palete()

Generates vector with colours used for Betti plots. Designed for Betti plots consistency.
"""
function get_bettis_color_palete(; min_dim = 1, use_set::Integer = 1)
    # TODO what does the number in the function below is used for?

    if use_set == 1
        cur_colors = [Gray(bw) for bw = 0.0:0.025:0.5]
        if min_dim == 0
            colors_set = [RGB(87 / 256, 158 / 256, 0 / 256)]
        else
            colors_set = []
        end
        max_RGB = 256
        colors_set = vcat(
            colors_set,
            [
                RGB(255 / max_RGB, 206 / max_RGB, 0 / max_RGB),
                RGB(248 / max_RGB, 23 / max_RGB, 0 / max_RGB),
                RGB(97 / max_RGB, 169 / max_RGB, 255 / max_RGB),
                RGB(163 / max_RGB, 0 / max_RGB, 185 / max_RGB),
                RGB(33 / max_RGB, 96 / max_RGB, 45 / max_RGB),
                RGB(4 / max_RGB, 0 / max_RGB, 199 / max_RGB),
                RGB(135 / max_RGB, 88 / max_RGB, 0 / max_RGB),
            ],
            cur_colors,
        )
    else
        use_set == 2
        cur_colors = get_color_palette(:auto, 1)
        cur_colors3 = get_color_palette(:lightrainbow, 1)
        cur_colors2 = get_color_palette(:cyclic1, 1)
        if min_dim == 0
            # colors_set =  [cur_colors[3], cur_colors[5], [:red], cur_colors[1]] #cur_colors[7],
            colors_set = [cur_colors3[3], cur_colors[5], cur_colors3[end], cur_colors[1]] #cur_colors[7],
        else
            colors_set = [cur_colors[5], cur_colors3[end], cur_colors[1]] #cur_colors[7],
            # colors_set =  [cur_colors[5], [:red], cur_colors[1], cur_colors[14]]
        end
        # for c =  [collect(11:25);]
        #     push!(colors_set, cur_colors2[c])
        # end
        colors_set = vcat(colors_set, [cur_colors2[c] for c in [collect(11:25);]])
    end

    return colors_set
end

# ==============================
#  ======= Untested code =======

# using Measures
# using Plots.PlotMeasures

#
# # Source: https://github.com/JuliaPlots/Plots.jl/issues/897
# function setdefaultplottingparams(;upscale=2)
#     #8x upscaling in resolution
#     fntsm = Plots.font("sans-serif", pointsize=round(12.0*upscale))
#     fntlg = Plots.font("sans-serif", pointsize=round(18.0*upscale))
#     default(titlefont=fntlg, guidefont=fntlg, tickfont=fntsm, legendfont=fntsm)
#     default(size=(800*upscale,600*upscale)) #Plot canvas size
#     default(dpi=500) #Only for PyPlot - presently broken
# end

#%%
"""
  plot_bettis_collection(bettis_collection, bett_num; step=1, show_plt=true, R=0., G=0.4, B=1.0)

PLots collection of Betti curves of rank 'bett-num'. Every successive plot has
lower opacity than predecessor.  'step' defines step between collection elements
that are ploted. By default, plot is displayed after carteation. This can be
disabled by setting 'show_plt' to false.

Color of the plot can be set with 'R', 'G', 'B' parameters.
"""
function plot_bettis_collection(bettis_collection,
                                bett_num,
                                max_rank;
                                step = 1,
                                show_plt = true,
                                R = 0.0,
                                G = 0.4,
                                B = 1.0)
    step > 0 || error("Steps should be natural number!")
    bettis_total = size(bettis_collection, 1)
    colors_set = zeros(Float64, bettis_total, 4)
    colors_set[:, 1] .= R
    colors_set[:, 2] .= G
    colors_set[:, 3] .= B
    max_betti = get_max_betti_from_collection(bettis_collection)
    @info "max_betti" max_betti

    x = 0
    y = bettis_total * 0.1
    va_range = collect(range(bettis_total + x, y, length = bettis_total))

    colors_set[:, 4] .= va_range / findmax(va_range)[1]
    rgba_set = RGBA[]
    for k = 1:size(colors_set, 1)
        push!(
            rgba_set,
            RGBA(colors_set[k, 1], colors_set[k, 2], colors_set[k, 3], colors_set[k, 4]),
        )
    end

    plt_reference = plot(1, title = "Betti curves collection, rank $(bett_num)", label = "")
    for b = 1:step:bettis_total
        betti = bettis_collection[b]
        x_vals_1 = (1:size(betti[:, bett_num], 1)) / size(betti[:, bett_num], 1)
        plot!(x_vals_1, betti[:, bett_num], lc = rgba_set[b], label = "rank=$(max_rank-b)")
        plot!(ylim = (0, max_betti))
    end
    xlabel!("Normalised steps")
    ylabel!("Number of cycles")
    plot!(legend = true)

    show_plt && display(plt_reference)
    return plt_reference
end

#%%

"""
    get_max_bettis(bettis)

Returns the maximal bettis of Betti curves for all dimensions.
"""
function get_max_bettis(bettis)

    all_max_bettis = findmax(bettis, dims=1)[1]

    return all_max_bettis
end

# TODO change name
# TODO check what for dim is used, change to min dim
function get_max_betti_from_collection(bettis_collection; dim = 1)
    max_betti = 0
    for betti in bettis_collection
        # global max_betti
        local_max = findmax(betti)[1]
        if (local_max > max_betti)
            max_betti = local_max
        end
    end
    return max_betti
end

#%%
"""
  plot_and_save_bettis(eirene_results, plot_title::String,
                results_path::String; extension = ".png",
                data_size::String="", do_save=true,
                extend_title=true, do_normalise=true, max_dim=3,
                legend_on=true)

Plot Betti curves from 0 up to `max_dim` using `eirene_results` from Eirene library and
returns handler for figure. Optionally, if `do_save` is set, saves the figure
or if `do_normalise` is set, sets the steps range to be normalised to the
horizontal axis maximal value.
"""
function plot_and_save_bettis(bettis,
                                plot_title::String,
                                results_path::String;
                                file_name = "",
                                extension = ".png",
                                do_save = true,
                                do_normalise = true,
                                min_dim = 0,
                                max_dim = 3,
                                legend_on = true,
                                kwargs...)
    bettis = get_bettis(eirene_results, max_dim)
    if do_normalise
        bettis = normalise_bettis(bettis)
    end
    plot_ref =
        plot_bettis(bettis, plot_title, legend_on = legend_on, min_dim = min_dim, kwargs...)

    if do_save
        if isempty(file_name)
            file_name = plot_title * extension
        elseif isempty(findall(x -> x == extension[2:end], split(file_name, ".")))
            #check for the extension in file name
            file_name *= extension
        end

        save_figure_with_params(
            plot_ref,
            results_path;
            extension = extension,
            prefix = split(file_name, ".")[1],
        )

    end
    return plot_ref
end

# TODO merge functions for getting betti curves
# Original function returns 2 different types of betti curves. If no default
# value parameters is given, it returns vector of matrices. If num of steps is
# given, then it return matrix maxdim x numsteps.
# """
# bettis_eirene(matr, maxdim; mintime=-Inf, maxtime=Inf, numofsteps=Inf, mindim=1)
#
# Takes the `matr` and computes Betti curves up to `maxdim`. Return matrix only
# with betti curve values
#
#
# Function taken from: https://github.com/alexyarosh/hyperbolic
# """
#%%

@deprecate bettis_eirene(matr, maxdim; mintime = -Inf, maxtime = Inf, numofsteps = Inf, mindim = 1) get_bettis(results_eirene, max_dim; min_dim = 1)
@deprecate get_bettis_from_image(img_name, plot_params; file_path = "", plot_heatmaps = true, save_heatmaps = false, plot_betti_figrues = true) get_bettis(results_eirene, max_dim; min_dim = 1)
@deprecate get_bettis_from_image2(img_name;file_path = "",plot_heatmaps = true, save_heatmaps = false, plot_betti_figrues = true) get_bettis_from_image(img_name, plot_params; file_path = "", plot_heatmaps = true, save_heatmaps = false, plot_betti_figrues = true)
@deprecate plot_and_save_bettis2(eirene_results, plot_title::String, results_path::String; file_name = "", extension = ".png", data_size::String = "", do_save = true, extend_title = true, do_normalise = true, min_dim = 0, max_dim = 3, legend_on = true) plot_and_save_bettis(bettis, plot_title::String, results_path::String; file_name = "", extension = ".png", do_save = true, do_normalise = true, min_dim = 0, max_dim = 3, legend_on = true, kwargs...)
@deprecate get_and_plot_bettis(eirene_results; max_dim = 3, min_dim = 1, plot_title = "", legend_on = false) get_bettis(results_eirene, max_dim; min_dim = 1)

#%%
"""
  lower_ordmat_resolution(ordered_matrix::Array, total_bins::Int)

Takes ordered matrix 'input_matrix' and reduces the resolution of values in the
matrix into 'total_bins' bins.
"""
function lower_ordmat_resolution(ordered_matrix::Array, total_bins::Int)
    new_ordered_matrix = zeros(size(ordered_matrix))
    max_val = findmax(ordered_matrix)[1]
    min_val = findmin(ordered_matrix)[1]

    bin_step = max_val ÷ total_bins
    old_bins = min_val:bin_step:max_val

    for bin = 1:total_bins
        @debug "First step threshold is $(old_bins[bin])"
        indices = findall(x -> (x >= old_bins[bin]), ordered_matrix)
        new_ordered_matrix[indices] .= bin - 1
    end

    @debug "Max_val in new matrix is " findmax(new_ordered_matrix)
    @debug "And should be " total_bins - 1

    return new_ordered_matrix
end

#%%
"""
  average_bettis(bettis_matrix; up_factor=8)

Takes the average values of betti curves stored in 'bettis_matrix'.

'bettis_matrix' consist of different simulations(first index of the matrix),
different ranks (third index of the matrix). Second index of the matrices
(saples) may vary accross many simulations and for this reason, all betti curves
are upsampled by a factor of 'upsample_factor' and then the average for every
dimension is computed.
"""
function average_bettis(bettis_matrix::Matrix; up_factor = 8)

    bettis_matrix_backup = copy(bettis_matrix)

    simulations = size(bettis_matrix, 1)
    dimensions = size(bettis_matrix[1], 1)

    max_samples = 0
    for k = 1:simulations
        # global max_samples
        current_len = length(bettis_matrix[k][1][:, 1])
        if max_samples < current_len
            max_samples = current_len
        end
    end

    bettis_size = size(bettis_matrix)

    total_upsamples = (max_samples - 1) * up_factor + 1
    x_resampled = range(0, 1, step = total_upsamples)

    avg_bettis = zeros(total_upsamples, dimensions)
    std_bettis = copy(avg_bettis)
    resampled_bettis = zeros(simulations, total_upsamples, dimensions)

    # resample betti curves
    for simulation = 1:simulations, betti = 1:dimensions
        resampled_bettis[simulation, :, betti] =
            upsample_vector2(bettis_matrix[simulation][betti][:, 2], total_upsamples)
    end

    # average and std Betti
    for dimension = 1:dimensions
        avg_bettis[:, dimension] = mean(resampled_bettis[:, :, dimension], dims = 1)
        std_bettis[:, dimension] = mean(resampled_bettis[:, :, dimension], dims = 1)
    end

    return avg_bettis, std_bettis
end

#%%
function upsample_vector2(input_vector, total_upsamples)
    total_orig_samples = size(input_vector, 1) - 1

    x_vals = range(0, 1, length = total_orig_samples + 1)
    spl = Spline1D(x_vals, input_vector)

    x_upsampled = range(0, 1, length = total_upsamples)
    y_upsampled = spl(x_upsampled)

    # ref = plot(range(0, 1, length=total_orig_samples), input_vector);
    # plot!(x_vals, y_upsampled);
    # display(ref)

    return y_upsampled
end

#%%
"""
  upsample_vector(input_vector; upsample_factor::Int=8)

Takes an 'input_vector' and returns a vector which has 'upsample_factor' many
times more samples. New samples are interpolated with 'spl' function from
'Dierckx' package.

"""
function upsample_vector(input_vector; upsample_factor::Int = 8)
    total_orig_samples = size(input_vector, 1) - 1
    total_samples = upsample_factor * total_orig_samples + 1

    x_vals = range(0, 1, length = total_orig_samples + 1)
    spl = Spline1D(x_vals, input_vector)

    x_upsampled = range(0, 1, length = total_samples)
    y_upsampled = spl(x_upsampled)

    # ref = plot(range(0, 1, length=total_orig_samples), input_vector);
    # plot!(x_vals, y_upsampled);
    # display(ref)

    return y_upsampled
end

# =========--=======-========-==========-=======-
# From bettis areas
# Area under Betti curve functions
#%%
"""
    get_area_under_betti_curve(betti_curves, min_dim, max_dim)

Computes the area under Betti curves stored in 'betti_curves', where each row is
a Betti curve and each column is a value.
"""
function get_area_under_betti_curve(betti_curves::Union{Matrix{Float64}, Array{Array{Float64,2}}};do_normalised::Bool=false)
    #TODO check this part
    if size(betti_curves,2) < 2
        bettis_vector = vectorize_bettis(betti_curves)
    else
        bettis_vector = betti_curves
    end
    # @info sum(bettis_vector, dims=1)
    bettis_area = sum(bettis_vector, dims=1)

    if do_normalised
        total_steps = size(bettis_vector,1)
        bettis_area ./= total_steps
    end
    # @info bettis_area
    return bettis_area
end

@deprecate get_area_under_betti_curve(C, min_dim, max_dim) get_area_under_betti_curve(betti_curves; do_normalised=false)

#%%
"""
    get_dataset_bettis_areas(dataset; min_dim::Integer=1, max_dim::Integer=3, return_matrix::Bool=true)

Computes topology of every matrix in dataset, computes Betti curves for dimensions
min_dim up to max_dim and returns vector (or matrix) of areas under Betti curves.
"""
function get_dataset_bettis_areas(dataset; min_dim::Integer=1, max_dim::Integer=3, return_matrix::Bool=true)
    areas_vector = Array[]
    for data = dataset
        @info "Computing topology."
        C = eirene(data, maxdim=max_dim,)
        matrix_bettis = get_bettis(C,max_dim, min_dim=min_dim)
        push!(areas_vector, get_area_under_betti_curve(matrix_bettis))
    end
    if return_matrix
        return vcat([areas_vector[k] for k=1:10]...)
    else
        return areas_vector
    end
end

#%%
"""
    get_area_boxes(areas_matrix, min_dim::Integer, max_dim::Integer)

Plots the boxplot of area under betti curves.
"""
function get_area_boxes(areas_matrix, min_dim::Integer, max_dim::Integer)
    bplot = StatsPlots.boxplot()

    data_colors = get_bettis_color_palete()

    for (index, value) in enumerate(min_dim:max_dim)
        StatsPlots.boxplot!(bplot, areas_matrix[:,index], labels="β$(value)", color=data_colors[value])
    end

    return bplot
end

# function get_bettis_collection_from_matrices(ordered_matrices_collection; max_dim::Int=3, min_dim::Int=1)
#     bettis_collection = Array[]
#
#     for matrix = ordered_matrices_collection
# 		@debug "Computing Bettis..."
# 		eirene_geom = eirene(matrix,maxdim=max_B_dim,model="vr")
#
# 		bettis = reshape_bettis(get_bettis(eirene_geom, max_B_dim))
# 		push!(bettis_collection, bettis)
#     end
#
#     return bettis_collection
# end
#
# #%%
# TODO find what are the alternative functions for the functions below
# @deprecate get_dataset_topology(dataset; min_dim::Integer=1, max_dim::Integer=3, get_curves::Bool=true, get_areas::Bool=true, get_persistence_diagrams::Bool=true, do_normalise::Bool=true)
# @deprecate get_bettis_collection(ordered_matrices_collection; max_B_dim=3)
# @deprecate reshape_bettis(bettis)
# @deprecate print_hmap_with_bettis(ordered_matrices_collection, bettis_collection, plot_data::PlottingData)
# @deprecate make_hm_and_betti_plot(ordered_geom_gr, bettis, title_hmap, title_bettis, max_betti)
# @deprecate matrix_analysis(test_data::PlottingData;generation_function=get_geom_matrix)
# @deprecate multiscale_matrix_testing(sample_space_dims = 3, maxsim = 5, min_B_dim = 1, max_B_dim = 3, size_start = 10, size_step = 5, size_stop = 50; do_random = true, control_saving = false, perform_eavl = false)
# @deprecate plot_betti_numbers(betti_numbers, edge_density, title="Geometric  matrix"; stop=0.6)