text:
struct EnumerableDrop{T,S} <: Enumerable
    source::S
    n::Int
end

function drop(source::Enumerable, n::Integer)
    T = eltype(source)
    S = typeof(source)
    return EnumerableDrop{T,S}(source, Int(n))
end

Base.IteratorSize(::Type{EnumerableDrop{T,S}}) where {T,S} = (Base.IteratorSize(S) isa Base.HasLength || Base.IteratorSize(S) isa Base.HasShape) ? Base.HasLength() : Base.SizeUnknown()

Base.eltype(::Type{EnumerableDrop{T,S}}) where {T,S} = T

Base.length(iter::EnumerableDrop{T,S}) where {T,S} = max(length(iter.source)-iter.n,0)

function Base.iterate(iter::EnumerableDrop{T,S}) where {T,S}
    ret = iterate(iter.source)
    for i in 1:iter.n
        if ret===nothing
            return nothing
        else
            ret = iterate(iter.source, ret[2])
        end
    end
    return ret
end

function Base.iterate(iter::EnumerableDrop{T,S}, state) where {T,S}
    return iterate(iter.source, state)
end


text:
export bisection, hybridsolver, hybridinverse!

"Implement the classical Bisection algorithm on the segment `[xm, xp]` of the real function `f`."
function bisection(x, fx, xm, fm, xp, f)
@assert xm < xp "Error in the order of the argument"
    if sign(fx) == sign(f)
    xp = x
    f = fx
    else
    xm = x
    fm = fx
    end
    # More stable than 0.5*(xp +xm)
    x = xm + 0.5*(xp - xm)
    return x, fx, xm, fm, xp, f
end

# Combine Bisection and Newton method, this method has guaranteed convergence.
# Convergence order should be between linear and quadratic
# http://www.m-hikari.com/ams/ams-2017/ams-53-56-2017/p/hahmAMS53-56-2017.pdf
function hybridsolver(f, g, out, a, b; ϵx = 1e-7, ϵf = 1e-7, niter = 100)
    dxold = abs(b-a)
    dx = dxold
    fout = f(out)
    gout = g(out)

    fa = f(a)
    fb = f(b)

    @inbounds for j=1:niter
        # Bisect if Newton out of range, or not decreasing fast enough.
        if ((out - b)*gout - fout)*((out - a)*gout - fout) > 0.0 || abs(2.0*fout) >  abs(dxold * gout)
            dxold = dx
            dx = 0.5*(b-a)
            out = a + dx
            if isapprox(a, out, atol = ϵx)
                return out
            end
        else #Newton step is acceptable
            dxold = dx
            dx    = fout/gout
            tmp   = out
            out = out - dx
            if isapprox(tmp, out, atol = ϵx)
                return out
            end
        end
        # Convergence criterion
        if abs(dx) < ϵx || abs(fout) < ϵf
            return out
        end
        # The one new function evaluation per iteration
        fout = f(out)
        gout = g(out)
        # Maintain the bracket on the root
        if fout<0.0
            a = out
        else
            b = out
        end
    end

    return out
end

"""
$(TYPEDSIGNATURES)

Invert the ensemble matrix `X` such that R(X) = F
"""
function hybridinverse!(X, F, R::IntegratedFunction, S::Storage; niter= 1000, ϵx = 1e-8, ϵf = 1e-8, P::Parallel = serial)
    Nψ = R.Nψ
    Nx = R.Nx
    NxX, Ne = size(X)

    @assert NxX == R.Nx "Wrong dimension of the sample X"

    cache  = zeros(Ne, Nψ)
    cache_vander = zeros(Ne, maximum(R.f.idx[:,Nx])+1)
    fout = zeros(Ne)

    # Remove f(x_{1:k-1},0) from the output F
    @avx for i=1:Ne
        f0i = zero(Float64)
        for j=1:Nψ
            f0i += (S.ψoffψd0[i,j])*R.f.coeff[j]
        end
        F[i] -= f0i
    end

    # lower and upper brackets of the ensemble members
    xk = view(X, Nx, :)
    xa = copy(xk)
    xb = copy(xk)
    σ = std(xa)
    xa .-= 2.0*σ
    xb .+= 2.0*σ
    fa = zeros(Ne)
    fb = zeros(Ne)

    ##  Find a bracket for the different samples
    factor = 1.6
    bracketed = false
    @inbounds for j=1:niter
        functionalf!(fa, xa, cache, cache_vander, S.ψoff, F, R)
        functionalf!(fb, xb, cache, cache_vander, S.ψoff, F, R)
        # We know that the function is strictly increasing
        if all(fa .< 0.0) && all(fb .> 0.0)
            bracketed = true
            break
        end
        @inbounds for i=1:Ne
            if fa[i]*fb[i] > 0.0
                Δ = factor*(xb[i] - xa[i])
                if abs(fa[i]) < abs(fb[i])
                    xa[i] -= Δ
                else
                    xb[i] += Δ
                end
            end
        end
    end

    @assert bracketed == true "Maximal number of iterations reached, without bracket"

    dx = zeros(Ne)
    gout = zeros(Ne)
    # Initial guess: mid point of the bracket
    @avx for i=1:Ne
        xai = xa[i]
        xbi = xb[i]
        dx[i] = xbi - xai
        xk[i] = xai + 0.5*(dx[i])
    end

    dxold = copy(dx)

    functionalf!(fout, xk, cache, cache_vander, S.ψoff, F, R)
    functionalg1D!(gout, xk, cache, cache_vander, S.ψoff, F, R)
    convergence = false

    @inbounds for j=1:niter
        @inbounds for i=1:Ne
            # Bisect if Newton out of range, or not decreasing fast enough.
            if ((xk[i] - xb[i])*gout[i] - fout[i])*((xk[i] - xa[i])*gout[i] - fout[i]) > 0.0 ||
                abs(2.0*fout[i]) >  abs(dxold[i] * gout[i])
                dxold[i] = dx[i]
                dx[i] = 0.5*(xb[i]-xa[i])
                xk[i] = xa[i] + dx[i]
                # if isapprox(xa[i], xk[i], atol = ϵx)
                # end
            else #Newton step is acceptable
                dxold[i] = dx[i]
                dx[i]    = fout[i]/gout[i]
                xk[i] -=  dx[i]
                # if isapprox(dx[i], 0.0, atol = ϵx)
                # # if isapprox(tmp, xk[i], atol = ϵx)
                #     continue
                # end
            end

            # if abs(dx[i]) < ϵx || abs(fout[i]) < ϵf
            #     continue
            # end
        end
        # Convergence criterion
        # if norm(xa - xk, Inf) < ϵx || norm(dx, Inf) < ϵx || norm(fout, Inf) < ϵf
        # if norm(dx, Inf) < ϵx || norm(fout, Inf) < ϵf
        if norm(fout) < ϵf
            convergence = true
            break
        end

        # Evaluate the function and its gradient for the different samples
        functionalf!(fout, xk, cache, cache_vander, S.ψoff, F, R)
        functionalg1D!(gout, xk, cache, cache_vander, S.ψoff, F, R)

        # Maintain the bracket on the root
        @inbounds for i=1:Ne
            if fout[i]<0.0
                xa[i] = xk[i]
            else
                xb[i] = xk[i]
            end
        end
    end

    @assert convergence == true "Inversion did not converge"
end

"""
$(TYPEDSIGNATURES)

Invert the ensemble matrix `X` such that C(X) = F
"""
hybridinverse!(X, F, C::HermiteMapComponent, S::Storage; P::Parallel = serial) = hybridinverse!(X, F, C.I, S; P = P)

"""
$(TYPEDSIGNATURES)

Invert the ensemble matrix `X` such that L.C(L.L*X) = F
"""
function hybridinverse!(X::Array{Float64,2}, F, L::LinHermiteMapComponent, S::Storage; P::Parallel = serial)
    # Pay attention that S is computed in the renormalized space for improve performance !!!
    transform!(L.L, X)
    hybridinverse!(X, F, L.C, S; P = P)
    itransform!(L.L, X)
end

# function hybridinverse!(X, F, R::IntegratedFunction, S::Storage; niter= 100, ϵx = 1e-7, ϵf = 1e-7, P::Parallel = serial)
#     Nψ = R.Nψ
#     Nx = R.Nx
#     NxX, Ne = size(X)
#
#     @assert NxX == R.Nx "Wrong dimension of the sample X"
#
#     cache  = zeros(Ne, Nψ)
#     cache_vander = zeros(Ne, maximum(R.f.idx[:,Nx])+1)
#     fout = zeros(Ne)
#
#     # Remove f(x_{1:k-1},0) from the output F
#     @avx for i=1:Ne
#         f0i = zero(Float64)
#         for j=1:Nψ
#             f0i += (S.ψoffψd0[i,j])*R.f.coeff[j]
#         end
#         F[i] -= f0i
#     end
#
#     # We can tighter the brackets, as we know that g(dxf)>0 for every dxf value
#     # so either the upper or the lower bound is zero for the different smaples
#
#     # lower and upper brackets of the ensemble members
#     xk = view(X, Nx, :)
#     # xbound = copy(xk)
#     # fbound = zeros(Ne)
#     #
#     σ = std(xk)
#     # # We can reduce by two the number of evaluations, if we use the sign of F
#     xbound = xk +2σ*sign.(F)
#     # @inbounds for i=1:Ne
#     #     if F[i]>0.0
#     #         # implies the root is >0
#     #         # xa[i] = 0.0
#     #         # xb[i] = xk[i] + 2σ
#     #         xbound[i] += 2σ
#     #     else
#     #         # xa[i] = xk[i] - 2σ
#     #         # xb[i] = 0.0
#     #         xbound[i] -= 2σ
#     #     end
#     # end
#     # Evaluate at 0.0 f
#     # We can rewrite the bracketing with a single function evaluation for xbound
#     factor = 1.6
#     bracketed = false
#     for j=1:niter
#         functionalf!(fbound, xbound, cache, cache_vander, S.ψoff, F, R)
#         # We know that the function is strictly increasing
#         if all(sign.(fbound) == sign.(xbound))
#             bracketed = true
#             break
#         end
#         @inbounds for i=1:Ne
#             if sign(fbound[i]) != sign(xbound[i])
#                 Δ = factor*(xbound[i] - 0.0)
#                 if xbound[i]<0.0
#                     xbound[i] -= Δ
#                 else
#                     xbound[i] += Δ
#                 end
#             end
#         end
#     end
#     #
#     @assert bracketed == true "Maximal number of iterations reached, without bracket"
#
#     # Let create xa, xb, fa, fb based on xbound
#     xa = zeros(Ne)
#     xb = zeros(Ne)
#     fa = zeros(Ne)
#     fb = zeros(Ne)
#     #
#     @inbounds for i=1:Ne
#         if xbound[i]<0.0
#             xa[i] = xbound[i]
#             # fa[i] = fbound[i]
#             # fb[i] = F[i]
#         else
#             xb[i] = xbound[i]
#             # fa[i] = F[i]
#             # fb[i] = fbound[i]
#         end
#     end
#     @assert all(xa .)
#     functionalf!(fa, xa, cache, cache_vander, S.ψoff, F, R)
#     functionalf!(fb, xb, cache, cache_vander, S.ψoff, F, R)
#
#     xk = view(X, Nx, :)
#     xa = copy(xk)
#     xb = copy(xk)
#     σ = std(xa)
#     xa .-= 2.0*σ
#     xb .+= 2.0*σ
#
#     fa = zeros(Ne)
#     fb = zeros(Ne)
#
#     #  Find a bracket for the different samples
#     factor = 1.6
#     bracketed = false
#     @inbounds for j=1:niter
#         functionalf!(fa, xa, cache, cache_vander, S.ψoff, F, R)
#         functionalf!(fb, xb, cache, cache_vander, S.ψoff, F, R)
#         # We know that the function is strictly increasing
#         if all(fa .< 0.0) && all(fb .> 0.0)
#             bracketed = true
#             break
#         end
#         @inbounds for i=1:Ne
#             if fa[i]*fb[i] > 0.0
#                 Δ = factor*(xb[i] - xa[i])
#                 if abs(fa[i]) < abs(fb[i])
#                     xa[i] -= Δ
#                 else
#                     xb[i] += Δ
#                 end
#             end
#         end
#     end
#     @assert bracketed == true "Maximal number of iterations reached, without bracket"
#
#     dx = zeros(Ne)
#     gout = zeros(Ne)
#     # Initial guess: mid point of the bracket
#     @avx for i=1:Ne
#         xai = xa[i]
#         xbi = xb[i]
#         dx[i] = xbi - xai
#         xk[i] = xai + 0.5*(dx[i])
#     end
#
#     dxold = copy(dx)
#
#     functionalf!(fout, xk, cache, cache_vander, S.ψoff, F, R)
#     functionalg1D!(gout, xk, cache, cache_vander, S.ψoff, F, R)
#     convergence = false
#     @inbounds for j=1:niter
#         @inbounds for i=1:Ne
#             # Bisect if Newton out of range, or not decreasing fast enough.
#             if ((xk[i] - xb[i])*gout[i] - fout[i])*((xk[i] - xa[i])*gout[i] - fout[i]) > 0.0 ||
#                 abs(2.0*fout[i]) >  abs(dxold[i] * gout[i])
#                 dxold[i] = dx[i]
#                 dx[i] = 0.5*(xb[i]-xa[i])
#                 xk[i] = xa[i] + dx[i]
#                 # if isapprox(xa[i], xk[i], atol = ϵx)
#                 # end
#             else #Newton step is acceptable
#                 dxold[i] = dx[i]
#                 dx[i]    = fout[i]/gout[i]
#                 xk[i] -=  dx[i]
#                 # if isapprox(dx[i], 0.0, atol = ϵx)
#                 # # if isapprox(tmp, xk[i], atol = ϵx)
#                 #     continue
#                 # end
#             end
#
#             # if abs(dx[i]) < ϵx || abs(fout[i]) < ϵf
#             #     continue
#             # end
#         end
#         # Convergence criterion
#         # @show min(norm(xa - xk, Inf), norm(xb - xk, Inf)), norm(dx, Inf), norm(fout, Inf)
#         if min(norm(xa - xk, Inf), norm(xb - xk, Inf))< ϵx || norm(dx, Inf) < ϵx || norm(fout, Inf) < ϵf
#             convergence = true
#             break
#         end
#
#         # Evaluate the function and its gradient for the different samples
#         functionalf!(fout, xk, cache, cache_vander, S.ψoff, F, R)
#         functionalg1D!(gout, xk, cache, cache_vander, S.ψoff, F, R)
#
#         # Maintain the bracket on the root
#         @inbounds for i=1:Ne
#             if fout[i]<0.0
#                 xa[i] = xk[i]
#             else
#                 xb[i] = xk[i]
#             end
#         end
#     end
#     @assert convergence == true "Inversion did not converge"
# end

ext: jl


text:
# %% Simulate progression of wealth over many rounds with different fees
# using CairoMakie
using Distributions
using DataFrames
using Gadfly
import Cairo, Fontconfig
using Printf

# %% Functions
g = Geometric(0.5)
toss() = rand(g)+1

function symlog(x)
    return sign(x)*log10(abs(x)+1)
end

function calcWinnings(N::Int,P::Int, F::Float64=log2(N)+2)
    winnings = zeros(Float64, N, P)
    for p in 1:P
        for i in 2:N
            winnings[i,p] = winnings[i-1,p] + 2^toss() - F
        end
    end
    return winnings
end

# %% Figure setup
theme = Theme(major_label_font="PragmataPro Mono",
              minor_label_font="PragmataPro Mono",
              key_title_font="PragmataPro Mono Bold",
              key_position=:none)
Gadfly.push_theme(theme)

nbins = 1000

function plotWinnings(winnings, fee)
    fig = plot(
        winnings,
        x=Row.index,
        y=Col.value,
        Geom.hexbin(xbincount=nbins, ybincount=nbins),
        Guide.xlabel("ROUNDS"),
        Guide.ylabel("CUMULATIVE PAYOUT (SYMLOG)"),
        Guide.xticks(ticks=[0 2500 5000 7500 10000]),
        Guide.yticks(label=false),
        Guide.title(@sprintf("Fee Per Round = %s", fee)),
        Coord.cartesian(xmin=0, xmax=1e4, ymin=-8, ymax=8))
    return fig
end

# %% Plot figures
wN = calcWinnings(Int(1e4), Int(1e3)) .|> symlog
draw(PDF("../../figures/winnings-N.pdf"), plotWinnings(wN, "log2(10000)+2"))

# %%
df2_wide = DataFrame(w2, :auto)
cols = names(df2_wide)
df2_wide.x = 1:(size(df2_wide)[1])
df2_wide.f .= "P=\$2"
df2 = stack(df2_wide, cols)

df100_wide = DataFrame(w100, :auto)
cols = names(df100_wide)
df100_wide.x = 1:(size(df100_wide)[1])
df100_wide.f .= "P=\$100"
df100 = stack(df100_wide, cols)

df = vcat(df2, df100)

# %%
p = plot(df,
         x=:x,
         y=:value,
         ygroup=:f,
         Geom.subplot_grid(Geom.hexbin(xbincount=nbins, ybincount=nbins),
                           Coord.cartesian(xmin=0, xmax=1e4, ymin=-8, ymax=8),
                           # Guide.yticks(ticks=nothing), Does not work :(
                           Guide.xticks(ticks=[0 2500 5000 7500 10000]),
                           ),
         Guide.xlabel("ROUNDS"),
         Guide.ylabel("CUMULATIVE PAYOUT (SYMLOG)"),
         Guide.title("Distribution of Payouts over Time"))

draw(PDF("../../figures/winnings-2-100.pdf"), p)
# %%

Guide.yticks

max_stars_repo_name: akoen/petersburg

max_issues_repo_name: akoen/petersburg

max_forks_repo_name: akoen/petersburg


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# ╔═╡ 9fd2d94c-4e6c-11ec-040b-31d105963e37
# Activate the shared project environment
begin
    import Pkg
    Pkg.activate(Base.current_project())

	using IntrinsicDecoherence
end

# ╔═╡ f1827968-cc5f-4f76-8846-69cd02e21bd3
# Calling plotting suite
begin
	using PyPlot
	using LaTeXStrings
	ioff()
	rc("text", usetex = true)
end

# ╔═╡ 677404c2-d0df-4e00-9cc2-dfb0976fe9b1
md"
# Intrinsic decoherence dynamics in the three-coupled harmonic oscillators interaction

Authors: **_Alejandro R. Urzúa & Héctor M. Moya Cessa_**

Instituto Nacional de Astrofísica, Óptica y Electrónica

([arurz@inaoep.mx](mailto:arurz@inaoep.mx), [hmmc@inaoep.mx](mailto:hmmc@inaoep.mx))

Developer: AR. Urzúa (2021)

 [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
"

# ╔═╡ 6562f5c9-6b21-47a2-a85e-dfbc6435c32a
md"
## Initial problem

The principal idea is to solve the mean value evolution

$$\begin{equation}
	\braket{\hat{A}}  = e^{-\gamma t}\sum_{k=0}^{\infty}\frac{ (\gamma t)^k}{k!} \left\langle \psi_{k} | \hat{A} | \psi_{k}\right\rangle,
\end{equation}$$

giving that the $k$-th element of the wavefunction $\psi$ is given by

$$\begin{equation}
    \rho_{k}=\ket{\psi_{k}}\bra{\psi_k}, \qquad \ket{\psi_{k}}=e^{-\mathrm{i} k\frac{\hat{H}}{\gamma}}\ket{\psi(0)},
\end{equation}$$

with the Hamiltonian of the system defined as

$$\newcommand{\op}[1]{\hat{#1}}\begin{equation}
\begin{aligned}
    \op{H} = \omega\left(\op{a}_{1}^{\dagger}\op{a}_{1} + \op{a}_{2}^{\dagger}\op{a}_{2} + \op{a}_{3}^{\dagger}\op{a}_{3}\right) + \lambda\left(\op{a}_{1}^{\dagger}\op{a}_{2} + \op{a}_{2}^{\dagger}\op{a}_{1}\right)
    + g\left[\op{a}_{3}\left(\op{a}_{1}^{\dagger}+\op{a}_{2}^{\dagger}\right)+\op{a}_{3}^{\dagger}(\op{a}_{1}+\op{a}_{2})\right].
\end{aligned}
\end{equation}$$
"

# ╔═╡ 09a75736-ae84-4889-8e58-b1a8e16d8231
md"
## Solution

The solution we find for the $k$-th element of the wavefunction $\psi$ is given by

$$\newcommand{\op}[1]{\hat{#1}}\newcommand{\mi}{\mathrm{i}}\begin{align*}
    \ket{\psi_{k}} &= \op{R}_{12}^{\dagger}\op{R}_{2}^{\dagger}\,e^{-\mi k\frac{\op{H}_{2}}{\gamma}}\ket{\alpha\cos\theta}_{1}\ket{\alpha\cos\phi\sin\theta}_{2}\ket{-\alpha\sin\phi\sin\theta}_{3}\\
    &= \op{R}_{12}^{\dagger}\op{R}_{2}^{\dagger}\ket{e^{-\frac{\mi k}{\gamma}\Omega_{2}}\alpha\cos\theta}_{1}\ket{e^{-\frac{\mi k}{\gamma}\omega_{-}}\alpha\cos\phi\sin\theta}_{2}\ket{-e^{-\frac{\mi k}{\gamma}\Omega}\alpha\sin\phi\sin\theta}_{3},
\end{align*}$$

with $\gamma$ the decay parameter, $\theta = \pi/4$, and
$$\begin{equation*}
    2\phi = \arctan\left(\frac{2\sqrt{2}g}{\omega_{+}-\omega}\right),
\end{equation*}$$ where g is the coupling strength, and $\omega_{\pm} = \omega \pm \lambda$.
"

# ╔═╡ fa63f69e-d586-485c-89c0-1556addb738d
# Set the constant parameters
begin
	const α = 4
	const ω = 4
	const λ = 0.5
end;

# ╔═╡ 6a7d48b3-beff-453c-99d3-7d4b9184565a
md"

## Case 1. $\gamma = 10$

"

# ╔═╡ d845eb59-ace3-4280-b5ab-de280ec11787
# For γ = 10
begin
	t_list = range(0, stop = 70, length = 1000)

	texv1_a = [exv_a(ω, λ, 0.1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv1_a1 = [exv_a₁(ω, λ, 0.1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv1_a2 = [exv_a₂(ω, λ, 0.1, 10, α, t_list[τ]) for τ in 1:length(t_list)]

	texv2_a = [exv_a(ω, λ, 0.5, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv2_a1 = [exv_a₁(ω, λ, 0.5, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv2_a2 = [exv_a₂(ω, λ, 0.5, 10, α, t_list[τ]) for τ in 1:length(t_list)]

	texv3_a = [exv_a(ω, λ, 1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv3_a1 = [exv_a₁(ω, λ, 1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
	texv3_a2 = [exv_a₂(ω, λ, 1, 10, α, t_list[τ]) for τ in 1:length(t_list)]
end;

# ╔═╡ 568f94f7-716e-4322-9d3c-265a2b82a4a5
begin
	figs, (axs1, axs2, axs3) = subplots(3,1, figsize = (6,10), sharex = true)

	axs1.plot(t_list, real(texv1_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axs1.plot(t_list, real(texv1_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axs1.plot(t_list, real(texv1_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axs2.plot(t_list, real(texv2_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axs2.plot(t_list, real(texv2_a2), "b", ms = 4, lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axs2.plot(t_list, real(texv2_a), "g", linestyle = (0, (5, 5)), lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axs3.plot(t_list, real(texv3_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axs3.plot(t_list, real(texv3_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axs3.plot(t_list, real(texv3_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axs1.text(69, 0, L"(a)", fontsize=16)
	axs2.text(69, 0, L"(b)", fontsize=16)
	axs3.text(69, 0, L"(c)", fontsize=16)

	axs1.set_yticks([0, 4, 8, 12, 16])
	axs2.set_yticks([0, 4, 8, 12, 16])
	axs3.set_yticks([0, 4, 8, 12, 16])

	axins = axs3.inset_axes([0.7, 0.6, 0.25, 0.25])
	axins.plot(t_list, real(texv3_a1), "r")
	axins.plot(t_list, real(texv3_a2), "b")
	axins.plot(t_list, real(texv3_a), "g")
	axins.set_xlim([0,8])
	axins.set_title(L"Close view $0 \leq t \leq 8$", fontsize = 12)
	setp(axins.get_xticklabels(), fontsize = 12)
	setp(axins.get_yticklabels(), fontsize = 12)
	axins.set_yticks([0,8,16])
	axins.set_xticks([0,4,8])

	axs1.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axs1.get_xticklabels(), fontsize = 18)
	setp(axs1.get_yticklabels(), fontsize = 18)
	#axs1.legend(fontsize = 13)
	axs1.annotate(L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle", xy=(12, 12), xytext=(20, 15), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axs1.annotate(L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle", xy=(21, 10.6), xytext=(29, 13), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axs1.annotate(L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle", xy=(27.1, 1.6), xytext=(33, 4), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)

	axs2.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axs2.get_xticklabels(), fontsize = 18)
	setp(axs2.get_yticklabels(), fontsize = 18)
	#axs2.legend(fontsize = 13)
	axs2.annotate(L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle", xy=(8.5, 11.5), xytext=(20, 15), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axs2.annotate(L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle", xy=(23, 3.8), xytext=(27, 6.5), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axs2.annotate(L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle", xy=(21, 3.2), xytext=(25, 0), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=-0.2",
                                  "fc"=>"w")
	)

	axs3.set_xlabel(L"t\; [s]", fontsize = 18)
	axs3.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axs3.get_xticklabels(), fontsize = 18)
	setp(axs3.get_yticklabels(), fontsize = 18)
	#axs3.legend(fontsize = 13)
	axs3.annotate(L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle", xy=(8.7, 9.15), xytext=(15, 12), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axs3.annotate(L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle", xy=(22.8, 7.8), xytext=(29, 10.5), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axs3.annotate(L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle", xy=(21, 3.92), xytext=(27, 0.5), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=-0.2",
                                  "fc"=>"w")
	)

	tight_layout()

	#savefig("γ10_sp.pdf", dpi = 300)
	#savefig("γ10_sp.png", dpi = 300)
end

# ╔═╡ dfd3f4be-2fae-4a48-8e78-c40d7252b711
figs

# ╔═╡ 7f491f8e-a79e-4e6e-ae2d-5d67af390512
md"
## Case 2. $\gamma = 100$
"

# ╔═╡ 49f28ae0-b7ed-4570-aa2b-80cbeeeccb81
# For γ = 100
begin
	t_lisT = range(0, stop = 7, length = 500)

	texV1_a = [exv_a(ω, λ, 0.1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV1_a1 = [exv_a₁(ω, λ, 0.1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV1_a2 = [exv_a₂(ω, λ, 0.1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]

	texV2_a = [exv_a(ω, λ, 0.5, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV2_a1 = [exv_a₁(ω, λ, 0.5, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV2_a2 = [exv_a₂(ω, λ, 0.5, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]

	texV3_a = [exv_a(ω, λ, 1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV3_a1 = [exv_a₁(ω, λ, 1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
	texV3_a2 = [exv_a₂(ω, λ, 1, 100, α, t_lisT[τ]) for τ in 1:length(t_lisT)]
end;

# ╔═╡ 6bfbc529-17b9-4e4e-bcdb-c7e3c72cec04
begin
	figS, (axS1, axS2, axS3) = subplots(3,1, figsize = (6,10), sharex = true)

	axS1.plot(t_lisT, real(texV1_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axS1.plot(t_lisT, real(texV1_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axS1.plot(t_lisT, real(texV1_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axS2.plot(t_lisT, real(texV2_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axS2.plot(t_lisT, real(texV2_a2), "b", ms = 3, lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axS2.plot(t_lisT, real(texV2_a), "g", linestyle = (0, (5, 5)), lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axS3.plot(t_lisT, real(texV3_a1), "r", lw = 3, label = L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle")
	axS3.plot(t_lisT, real(texV3_a2), "b", lw = 3, label = L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle")
	axS3.plot(t_lisT, real(texV3_a), "g", lw = 3, label = L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle")

	axS1.text(6.9, 0, L"(a)", fontsize=16)
	axS2.text(6.9, 0, L"(b)", fontsize=16)
	axS3.text(6.9, 0, L"(c)", fontsize=16)

	axS1.set_yticks([0, 4, 8, 12, 16])
	axS2.set_yticks([0, 4, 8, 12, 16])
	axS3.set_yticks([0, 4, 8, 12, 16])

	axS1.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axS1.get_xticklabels(), fontsize = 18)
	setp(axS1.get_yticklabels(), fontsize = 18)
	#axS1.legend(fontsize = 13)
	axS1.annotate(L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle", xy=(0.55, 14.5), xytext=(1.6, 15), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axS1.annotate(L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle", xy=(3.55, 13.5), xytext=(4.5, 14), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axS1.annotate(L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle", xy=(3, 1), xytext=(3, 5), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)

	axS2.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axS2.get_xticklabels(), fontsize = 18)
	setp(axS2.get_yticklabels(), fontsize = 18)
	#axS2.legend(fontsize = 13)
	axS2.annotate(L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle", xy=(0.4, 14.5), xytext=(1.6, 15), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axS2.annotate(L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle", xy=(2.06, 7), xytext=(2.06, 10), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axS2.annotate(L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle", xy=(3.2, 3.13), xytext=(4, 5), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)

	axS3.set_xlabel(L"t\; [s]", fontsize = 18)
	axS3.set_ylabel(L"\textrm{Number of photons}", fontsize = 20)
	setp(axS3.get_xticklabels(), fontsize = 18)
	setp(axS3.get_yticklabels(), fontsize = 18)
	#axS3.legend(fontsize = 13)
	axS3.annotate(L"\langle \hat{a}_{3}^{\dagger}\hat{a}_{3} \rangle", xy=(0.3, 14.2), xytext=(1.6, 15), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axS3.annotate(L"\langle \hat{a}_{2}^{\dagger}\hat{a}_{2} \rangle", xy=(4.6, 13.5), xytext=(6, 14), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=0.2",
                                  "fc"=>"w")
	)
	axS3.annotate(L"\langle \hat{a}_{1}^{\dagger}\hat{a}_{1} \rangle", xy=(3.9, 3.14), xytext=(4, 7), size= 12,
		bbox=Dict("boxstyle"=>"round4", "fc"=>"w"),
            arrowprops=Dict("arrowstyle"=>"-|>",
                        "connectionstyle"=>"arc3,rad=-0.2",
                                  "fc"=>"w")
	)

	tight_layout()

	#savefig("γ100_sp.pdf", dpi = 300)
	#savefig("γ100_sp.png", dpi = 300)
end

# ╔═╡ 7d3ad796-1855-44e1-9505-9543587df7f8
figS

# ╔═╡ Cell order:
# ╟─677404c2-d0df-4e00-9cc2-dfb0976fe9b1
# ╟─9fd2d94c-4e6c-11ec-040b-31d105963e37
# ╟─f1827968-cc5f-4f76-8846-69cd02e21bd3
# ╟─6562f5c9-6b21-47a2-a85e-dfbc6435c32a
# ╟─09a75736-ae84-4889-8e58-b1a8e16d8231
# ╟─fa63f69e-d586-485c-89c0-1556addb738d
# ╟─6a7d48b3-beff-453c-99d3-7d4b9184565a
# ╟─d845eb59-ace3-4280-b5ab-de280ec11787
# ╟─568f94f7-716e-4322-9d3c-265a2b82a4a5
# ╟─dfd3f4be-2fae-4a48-8e78-c40d7252b711
# ╟─7f491f8e-a79e-4e6e-ae2d-5d67af390512
# ╟─49f28ae0-b7ed-4570-aa2b-80cbeeeccb81
# ╟─6bfbc529-17b9-4e4e-bcdb-c7e3c72cec04
# ╟─7d3ad796-1855-44e1-9505-9543587df7f8

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


text:
using Distributions, DataFrames, Optim, DataFramesMeta, StatsModels, Plots, StatsPlots, GLM

### Funciones ###

# Especificación de Cₜ
function Cₜ(α::Real, τ::Real)
    return (-log(1-τ))^(1/α)
end

# Especificación de función acumulada de la nueva reparametricación
function F_Wᵣ(Y::Real,Qₜ::Real,α::Real,τ::Real)
    β = Qₜ/Cₜ(α,τ)
    W = Weibull(α,β)
    return cdf(W,Y)
end

# Especificación del parámetro Qₜ en base a los datos y la nueva función de enlace
function Qₜ(β::AbstractArray, Database::AbstractArray)
    return exp.(Database * β)
end

# Especificación de la función generadora de números aleatorios
function rand_Wᵣ(n::Int,Qₜ::Real,α::Real,τ::Real)
    β = Qₜ/Cₜ(α, τ)
    W = Weibull(α,β)
    random = rand(W,n)
    return random
end

# Espeficiación de la base de datos de covariables simuladas
function DF_simulation(n::Int)
    X₁ = Normal(2,0.25)
    X₂ = Beta(2,3)
    X₃ = Gamma(1,0.5)
    DF = DataFrame(X₁ = rand(X₁,n),
                   X₂= rand(X₂,n),
                   X₃ = rand(X₃,n))
    return DF
end

# Especificación de la base de datos final, que contiene la variable censurada
function DF_F_Simulation(df::AbstractDataFrame,βs::AbstractArray,τ::Real)
    M = size(df,1)
    des_matrix = convert(Matrix,hcat(ones(M),df))
    Qₜᵢ = Qₜ(βs,des_matrix)
    Y =  first.(rand_Wᵣ.(1,Qₜᵢ,1,0.5))
    min_Y = minimum(Y); Q₉_Y = quantile(Y,0.8); interval = (Q₉_Y-min_Y)/6
    # Creación del intervalo de censura
    range = [min_Y:interval:Q₉_Y;Inf]
    # Inicialización de variable censurada: límites superiores (ϕₛ) e inferiores (ϕᵢ:)
    ϕᵢ = zeros(M)
    ϕₛ = zeros(M)
    # Censura de la variable Y: Definición del límite superior ϕₛ
    for i in 1:length(Y)
        for j in 1:length(range)
            if Y[i] < range[j]
                ϕₛ[i] = range[j]
                break
            end
        end
    end
    # Censura de la variable Y: Definición del límite inferior ϕᵢ
    for i in 1:length(Y)
        for j in 1:length(range)
            if Y[i] > reverse(range)[j]
                ϕᵢ[i] = reverse(range)[j]
                break
            end
        end
    end
    # Creación de la base de datos finalj
    F_df = hcat(DataFrame(ϕᵢ = ϕᵢ, ϕₛ = ϕₛ),df)
    return F_df
end

sim = DF_simulation(1000)
β_sim = [0.9,0.4,0.6,0.9]
t_sim = 0.5
final_db = DF_F_Simulation(sim,β_sim,t_sim)

## Recordatorio: Los límites inferiores y superiores requieren tener la notación indicada en el paper (ϕᵢ, ϕₛ)'
function reg_Wᵣ(df::AbstractDataFrame,Lᵢ::Int, Lₛ::Int,param::AbstractVector)
    Lₛ = 2
    covariates = select(final_db,Not(Lₛ))
    ## Creación de la matriz de diseño
    names_cov = Symbol.(names(covariates))
    f = Term(:ϕᵢ) ~ sum(term.(names_cov[2:length(names_cov)]))
    design_matrix = modelmatrix(ModelFrame(f,covariates))
    ## Creación de los betas relacionados
    M₀ = lm(f,covariates)
    first.(coef(M₀))

    βs = [first.(log.(first.([first.(coef(lm(f,covariates))),2])))]

end

function log_Lₖ(param::AbstractVector,df::AbstractDataFrame,τ::real)
    n = length(param) - 1
    for i in 1:n
        βs[i] = param[i]
    end
    α = param[n+1]
    Qₜ = Qₜ(βs,sim)
    Lₖ =

    end

end

function reg_Wᵣ(DF::AbstractDataFrame,ϕᵢ::Int,ϕᵢ₊₁::Int,t::Real)
    Parameters =[1.0,1.0,1.0,1.0,2.0]
    lower = [0.01, 0.01, 0.01, 0.01,0]

    func = TwiceDifferentiable(Lₖ(Parameters,DF,1,2,0.5))
    optimize(vars -> Lₖ(vars,DF,1,2,0.5),lower, Parameters)
end

reg_Wᵣ(DF,1,2,0.5)

max_stars_repo_name: jamanrique/MSc-Tesis-PUCP

max_issues_repo_name: jamanrique/MSc-Tesis-PUCP

max_forks_repo_name: jamanrique/MSc-Tesis-PUCP


text:
function boardtest(lib::Lib) where Lib<:PolyhedraLibrary
    A1 = -eye(Int, 9) # x >= 0
    b1 = zeros(Int, 9)
    A2 = eye(Int, 9) # x <= 1
    b2 = ones(Int, 9)
    A3 = zeros(Int, 9, 9)
    b3 = 3 * ones(Int, 9)
    i = 1
    for a = 1:3
        for b = (a+1):3
            for c = 1:3
                for d = (c+1):3
                    ac = a + (c-1) * 3
                    ad = a + (d-1) * 3
                    bc = b + (c-1) * 3
                    bd = b + (d-1) * 3
                    A3[i, ac] = 1
                    A3[i, ad] = 1
                    A3[i, bc] = 1
                    A3[i, bd] = 1
                    i += 1
                end
            end
        end
    end
    A = [A1; A2; A3]
    b = [b1; b2; b3]
    ine = hrep(A, b)
    poly = polyhedron(ine, lib)
    @test nrays(poly) == 0
    @test npoints(poly) == 614
    @test !isempty(poly)
    ext  = MixedMatVRep(vrep(poly))
    target = ones(Int, 9) * (3 // 4)
    ok = false
    for i = 1:size(ext.V, 1)
        # In julia v0.4 [i,:] returns a row matrix and in v0.5 it is
        # a 1D vector hence the use of vec
        if vec(ext.V[i,:]) == target
            ok = true
        end
    end
    @test ok

    cutA = ones(Int, 1, 9)
    cutb = 6
    Acut = [cutA; A]
    bcut = [cutb; b]
    inecut = hrep(Acut, bcut)
    polycut = polyhedron(inecut, lib)
    @test !isempty(polycut)
    #(isredundant, certificate) = isredundantinequality(polycut, 1)
    #@test !isredundant
    #@test certificate == target
    @test !isredundant(polycut, first(eachindex(halfspaces(polycut))))
    # @test BitSet() == gethredundantindices(polycut) # TODO reactivate it when I figure out why it makes LRS tests fail
    #(issredundant, scertificate) = isstronglyredundantinequality(polycut, 1)
    #@test !issredundant
    #@test scertificate == target
    #@test BitSet([]) == getstronglyredundantinequalities(polycut)
end


text:
function count_evaluations(::Type{LeastSquaresFactory}, arg0::LeastSquaresProblem, arg1::Incrementor)
    return jcall(LeastSquaresFactory, "countEvaluations", LeastSquaresProblem, (LeastSquaresProblem, Incrementor), arg0, arg1)
end

function create(::Type{LeastSquaresFactory}, arg0::MultivariateJacobianFunction, arg1::RealVector, arg2::RealVector, arg3::ConvergenceChecker, arg4::jint, arg5::jint)
    return jcall(LeastSquaresFactory, "create", LeastSquaresProblem, (MultivariateJacobianFunction, RealVector, RealVector, ConvergenceChecker, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5)
end

function create(::Type{LeastSquaresFactory}, arg0::MultivariateJacobianFunction, arg1::RealVector, arg2::RealVector, arg3::RealMatrix, arg4::ConvergenceChecker, arg5::jint, arg6::jint)
    return jcall(LeastSquaresFactory, "create", LeastSquaresProblem, (MultivariateJacobianFunction, RealVector, RealVector, RealMatrix, ConvergenceChecker, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5, arg6)
end

function create(::Type{LeastSquaresFactory}, arg0::MultivariateJacobianFunction, arg1::RealVector, arg2::RealVector, arg3::RealMatrix, arg4::ConvergenceChecker, arg5::jint, arg6::jint, arg7::jboolean, arg8::ParameterValidator)
    return jcall(LeastSquaresFactory, "create", LeastSquaresProblem, (MultivariateJacobianFunction, RealVector, RealVector, RealMatrix, ConvergenceChecker, jint, jint, jboolean, ParameterValidator), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7, arg8)
end

function create(::Type{LeastSquaresFactory}, arg0::MultivariateVectorFunction, arg1::MultivariateMatrixFunction, arg2::Vector{jdouble}, arg3::Vector{jdouble}, arg4::RealMatrix, arg5::ConvergenceChecker, arg6::jint, arg7::jint)
    return jcall(LeastSquaresFactory, "create", LeastSquaresProblem, (MultivariateVectorFunction, MultivariateMatrixFunction, Vector{jdouble}, Vector{jdouble}, RealMatrix, ConvergenceChecker, jint, jint), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function evaluation_checker(::Type{LeastSquaresFactory}, arg0::ConvergenceChecker)
    return jcall(LeastSquaresFactory, "evaluationChecker", ConvergenceChecker, (ConvergenceChecker,), arg0)
end

function model(::Type{LeastSquaresFactory}, arg0::MultivariateVectorFunction, arg1::MultivariateMatrixFunction)
    return jcall(LeastSquaresFactory, "model", MultivariateJacobianFunction, (MultivariateVectorFunction, MultivariateMatrixFunction), arg0, arg1)
end

function weight_diagonal(::Type{LeastSquaresFactory}, arg0::LeastSquaresProblem, arg1::RealVector)
    return jcall(LeastSquaresFactory, "weightDiagonal", LeastSquaresProblem, (LeastSquaresProblem, RealVector), arg0, arg1)
end

function weight_matrix(::Type{LeastSquaresFactory}, arg0::LeastSquaresProblem, arg1::RealMatrix)
    return jcall(LeastSquaresFactory, "weightMatrix", LeastSquaresProblem, (LeastSquaresProblem, RealMatrix), arg0, arg1)
end


text:
using Distributed
using ClusterManagers
# addprocs()
addprocs(LocalAffinityManager(np = parse(Int, ARGS[1]), mode = BALANCED, affinities = Int[]))
@everywhere push!(LOAD_PATH, pwd());
# println(ARGS)
println("=== Compiling ===\n");
@everywhere using main_module
println("=== Compiled ===");
# println(length(workers()))
# testraid()
main()
check_leaks()

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
#===============================================================================
                Extensions of Trajectories for GuidProp
===============================================================================#
function Trajectories.trajectory(
        P::GuidProp,
        v::Type=DD.default_type(P),
        w::Type=DD.default_wiener_type(P)
    )
    trajectory(time(P), P.P_target, v, w)
end

function Trajectories.trajectory(
        PP::AbstractArray{<:GuidProp},
        v::Type=DD.default_type(PP[1]),
        w::Type=DD.default_wiener_type(PP[1])
    )
    (
        process = [
            DD._process_traj(time(P), P.P_target, v) for P in PP
        ],
        wiener = [
            DD._wiener_traj(time(P), P.P_target, w) for P in PP
        ],
    )
end

#===============================================================================
                preconditioned Crank-Nicolson scheme
===============================================================================#
const Mutable = Array

function crank_nicolson!(y°, y, ρ) # For immutable types
    λ = sqrt(1-ρ^2)
    for i in 1:length(y)
        y°[i] = λ*y°[i] + ρ*y[i]
    end
end

function crank_nicolson!(y°::Vector{T}, y, ρ) where {T<:Mutable} #NOTE GPUs will need to be treated separately
    λ = sqrt(1-ρ^2)
    for i in 1:length(y)
        mul!(y°[i], y[i], true, ρ, λ)
    end
end

#===============================================================================
                Simple sampling over a single interval
===============================================================================#
"""
    Base.rand(
        [rng::Random.AbstractRNG], P::GuidProp, y1=zero(P); f=DD.__DEFAULT_F
    )

Sample a trajectory of a guided proposal `P` started from `y1`. Initialize
containers in the background and compute the functional `f` at the time of
sampling.
"""
function Base.rand(P::GuidProp, y1=zero(P); f=DD.__DEFAULT_F)
    rand(Random.GLOBAL_RNG, P, y1, DD.ismutable(y1); f=f)
end

function Base.rand(
        rng::Random.AbstractRNG, P::GuidProp, y1=zero(P); f=DD.__DEFAULT_F
    )
    rand(rng, P, y1, DD.ismutable(y1); f=f)
end

function Base.rand(
        rng::Random.AbstractRNG,
        P::GuidProp, y1::K, v::Val{false};
        f=DD.__DEFAULT_F
    ) where K
    w0 = (
        DD.default_wiener_type(P) <: Number ?
        zero(eltype(K)) :
        zero(similar_type(K, Size(DD.dim_wiener(P.P_target))))
    )
    Wnr = Wiener()
    X, W = trajectory(P, typeof(y1), typeof(w0))
    success, f_accum = false, nothing
    while !success
        rand!(rng, Wnr, W, w0)
        success, f_accum = DD.solve!(X, W, P, y1; f=f)
    end
    typeof(f) != DD._DEFAULT_F && return X, W, Wnr, f_accum
    X, W, Wnr
end

function Base.rand(
        rng::Random.AbstractRNG,
        P::GuidProp, y1::K, v::Val{true};
        f=DD.__DEFAULT_F
    ) where K
    error("in-place not implemented")
end

#===============================================================================
                in-place sampling over a single interval
===============================================================================#

#=---- Vanilla ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG],
        P::GuidProp,
        X, W, y1=zero(P);
        f=DD.__DEFAULT_F, Wnr=Wiener()
    )

Sample a trajectory of a guided proposal `P` started from `y1`. Use containers
`X` and `W` to save the results. Compute the functional `f` at the time of
sampling.
"""
function Random.rand!(
        P::GuidProp,
        X, W, y1=zero(P);
        f=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        Random.GLOBAL_RNG,
        P, X, W, y1, DD.ismutable(y1);
        f=f, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X, W, y1=zero(P);
        f=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        rng,
        P, X, W, y1, DD.ismutable(y1);
        f=f, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X, W, y1::K, v::Val{false};
        f=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    rand!(rng, Wnr, W)
    DD.solve!(X, W, P, y1; f=f)
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X, W, y1::K, v::Val{true};
        f=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    error("in-place not implemented")
end

#=---- with Crank-Nicolson scheme ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG],
        P::GuidProp,
        X°, W°, W, ρ, y1=zero(P);
        f=DD.__DEFAULT_F, Wnr=Wiener()
    )

Sample a trajectory of a guided proposal `P` started from `y1`. Use containers
`X°` and `W°` to save the results. Use a preconditioned Crank-Nicolson scheme
with memory parameter `ρ` and a previously sampled Wiener noise `W`. Compute the
functional `f` at the time of sampling.
"""
function Random.rand!(
        P::GuidProp,
        X°, W°, W, ρ, y1=zero(P);
        f=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        Random.GLOBAL_RNG,
        P, X°, W°, W, ρ, y1, DD.ismutable(y1);
        f=f, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X°, W°, W, ρ, y1=zero(P);
        f=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        rng,
        P, X°, W°, W, ρ, y1, DD.ismutable(y1);
        f=f, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X°, W°, W, ρ, y1::K, v::Val{false};
        f=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    rand!(rng, Wnr, W°)
    crank_nicolson!(W°.x, W.x, ρ)
    DD.solve!(X°, W°, P, y1; f=f)
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X°, W°, W, ρ, y1::K, v::Val{true};
        f=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    error("in-place not implemented")
end

#=---- with log-likelihood ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG]
        P::GuidProp,
        X, W, v::Val{:ll}, y1=zero(P);
        Wnr=Wiener()
    )

Sample a trajectory of a guided proposal `P` started from `y1`. Use containers
`X` and `W` to save the results. Compute log-likelihood (only path contribution)
along the way.
"""
function Random.rand!(
        P::GuidProp,
        X, W, v::Val{:ll}, y1=zero(P);
        Wnr=Wiener(), skip=0
    )
    rand!(
        Random.GLOBAL_RNG,
        P, X, W, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X, W, v::Val{:ll}, y1=zero(P);
        Wnr=Wiener(), skip=0
    )
    rand!(
        rng,
        P, X, W, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X, W, ::Val{:ll}, y1::K, ::Val{false};
        Wnr=Wiener(), skip=0
    ) where K
    rand!(rng, Wnr, W)
    solve_and_ll!(X, W, P, y1; skip=skip)
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X, W, ::Val{:ll}, y1::K, ::Val{true};
        Wnr=Wiener(), skip=0
    ) where K
    error("in-place not implemented")
end

#=---- with log-likelihood and Crank-Nicolson scheme ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG]
        P::GuidProp,
        X°, W°, W, ρ, v::Val{:ll}, y1=zero(P);
        Wnr=Wiener()
    )

Sample a trajectory of a guided proposal `P` started from `y1`. Use containers
`X°` and `W°` to save the results. Use a preconditioned Crank-Nicolson scheme
with memory parameter `ρ` and a previously sampled Wiener noise `W`. Compute
log-likelihood (only path contribution) along the way.
"""
function Random.rand!(
        P::GuidProp,
        X°, W°, W, ρ, v::Val{:ll}, y1=zero(P);
        Wnr=Wiener(), skip=0
    )
    rand!(
        Random.GLOBAL_RNG,
        P, X°, W°, W, ρ, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X°, W°, W, ρ, v::Val{:ll}, y1=zero(P);
        Wnr=Wiener(), skip=0
    )
    rand!(
        rng,
        P, X°, W°, W, ρ, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X°, W°, W, ρ, ::Val{:ll}, y1::K, ::Val{false};
        Wnr=Wiener(), skip=0
    ) where K
    rand!(rng, Wnr, W°)
    crank_nicolson!(W°.x, W.x, ρ)
    solve_and_ll!(X°, W°, P, y1; skip=skip)
end

function Random.rand!(
        rng::Random.AbstractRNG,
        P::GuidProp,
        X°, W°, W, ρ, ::Val{:ll}, y1::K, ::Val{true};
        Wnr=Wiener(),
    ) where K
    error("in-place not implemented")
end

#===============================================================================
                Simple sampling over multiple intervals
===============================================================================#
"""
    Base.rand(
        [rng::Random.AbstractRNG],
        PP::AbstractArray{<:GuidProp}, y1=zero(PP[1]); f=DD.__DEFAULT_F
    )

Sample a trajectory started from `y1`, defined for multiple guided proposals
`PP` that correspond to consecutive intervals. Initialize containers in the
background and compute the functionals `f` (one for each interval) at the time
of sampling.
"""
function Base.rand(
        PP::AbstractArray{<:GuidProp},
        y1=zero(PP[1]);
        f=DD.__DEFAULT_F
    )
    rand(Random.GLOBAL_RNG, PP, y1; f=f)
end

function Base.rand(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        y1::K;
        f=DD.__DEFAULT_F
    ) where K
    results = map(1:length(PP)) do i
        result = rand(rng, PP[i], y1, DD.ismutable(y1); f=f[i])
        y1 = result[1].x[end]
        result
    end
    XX = map(r->r[1], results)
    WW = map(r->r[2], results)
    Wnr = results[1][3]
    typeof(f) != DD._DEFAULT_F && return XX, WW, Wnr, map(r->r[4], results)
    XX, WW, Wnr
end

#===============================================================================
                out-of-place sampling over multiple intervals
===============================================================================#

#=---- Vanilla ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG]
        PP::AbstractArray{<:GuidProp},
        XX, WW, y1=zero(PP[1]);
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener()
    )

Sample a trajectory started from `y1` over multiple intervals for guided
proposals `PP` that correspond to consecutive intervals. Use containers `XX` and
`WW` to save the results. Compute the functionals `f` (one for each interval) at
the time of sampling and store the results in `f_out`.
"""
function Random.rand!(
        PP::AbstractArray{<:GuidProp},
        XX, WW, y1=zero(PP[1]);
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        Random.GLOBAL_RNG,
        PP, XX, WW, y1, DD.ismutable(y1);
        f=f, f_out=f_out, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX, WW, y1=zero(PP[1]);
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        rng,
        PP, XX, WW, y1, DD.ismutable(y1);
        f=f, f_out=f_out, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX, WW, y1::K, v::Val{false};
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    for i in eachindex(PP)
        rand!(rng, Wnr, WW[i])
        success, f_out[i] = rand!(
            rng, PP[i], XX[i], WW[i], y1, v; f=f[i], Wnr=Wnr
        )
        success || return false
        y1 = XX[i].x[end]
    end
    true
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX, WW, y1::K, v::Val{true};
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    error("in-place not implemented")
end

#=---- with Crank-Nicolson scheme ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG]
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρρ, y1=zero(PP[1]);
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener()
    )

Sample a trajectory started from `y1` over multiple intervals for guided
proposals `PP` that correspond to consecutive intervals. Use containers `XX°`
and `WW°` to save the results. Use a preconditioned Crank-Nicolson scheme
with memory parameter `ρ` (one for the entire interval) and a previously sampled
Wiener noise `WW`. Compute the functionals `f` (one for each interval) at
the time of sampling and store the results in `f_out`.
"""
function Random.rand!(
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, y1=zero(PP[1]);
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        Random.GLOBAL_RNG,
        PP, XX°, WW°, WW, ρ, y1, DD.ismutable(y1);
        f=f, f_out=f_out, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, y1=zero(PP[1]);
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener()
    )
    rand!(
        rng,
        PP, XX°, WW°, WW, ρ, y1, DD.ismutable(y1);
        f=f, f_out=f_out, Wnr=Wnr
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, y1::K, v::Val{false};
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    for i in eachindex(PP)
        success, f_out[i] = rand!(
            rng, PP[i], XX°[i], WW°[i], WW[i], ρ, y1, v; f=f[i], Wnr=Wnr
        )
        success || return false
        y1 = XX°[i].x[end]
    end
    true
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, y1::K, v::Val{true};
        f=DD.__DEFAULT_F, f_out=DD.__DEFAULT_F, Wnr=Wiener(),
    ) where K
    error("in-place not implemented")
end

#=---- with log-likelihood ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG]
        PP::AbstractArray{<:GuidProp},
        XX, WW, v::Val{:ll}, y1=zero(PP[1]);
        Wnr=Wiener()
    )

Sample a trajectory started from `y1` over multiple intervals for guided
proposals `PP` that correspond to consecutive intervals. Use containers `XX°`
and `WW°` to save the results. Compute log-likelihood (path contribution AND
end-points contribution) along the way.
"""
function Random.rand!(
        PP::AbstractArray{<:GuidProp},
        XX, WW, v::Val{:ll}, y1=zero(PP[1]);
        Wnr=Wiener(), skip=0
    )
    rand!(
        Random.GLOBAL_RNG,
        PP, XX, WW, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX, WW, v::Val{:ll}, y1=zero(PP[1]);
        Wnr=Wiener(), skip=0
    )
    rand!(
        rng,
        PP, XX, WW, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX, WW, v::Val{:ll}, y1::K, m::Val{false};
        Wnr=Wiener(), skip=0
    ) where K
    ll_tot = loglikhd_obs(PP[1], y1)
    for i in eachindex(PP)
        success, ll = rand!(
            rng, PP[i], XX[i], WW[i], v, y1, m; Wnr=Wnr, skip=skip
        )
        success || return false, ll
        ll_tot += ll
        y1 = XX[i].x[end]
    end
    true, ll_tot
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX, WW, ::Val{:ll}, y1::K, ::Val{true};
        Wnr=Wiener(), skip=0
    ) where K
    error("in-place not implemented")
end

#=---- with log-likelihood and Crank-Nicolson scheme ----=#
"""
    Random.rand!(
        [rng::Random.AbstractRNG]
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρρ, v::Val{:ll}, y1=zero(PP[1]);
        Wnr=Wiener()
    )

Sample a trajectory started from `y1` over multiple intervals for guided
proposals `PP` that correspond to consecutive intervals. Use containers `XX°`
and `WW°` to save the results. Use a preconditioned Crank-Nicolson scheme
with memory parameters `ρ` (one for the entire interval) and a previously
sampled Wiener noise `WW`. Compute log-likelihood (path contribution AND
end-points contribution) along the way.
"""
function Random.rand!(
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, v::Val{:ll}, y1=zero(PP[1]);
        Wnr=Wiener(), skip=0
    )
    rand!(
        Random.GLOBAL_RNG,
        PP, XX°, WW°, WW, ρ, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, v::Val{:ll}, y1=zero(PP[1]);
        Wnr=Wiener(), skip=0
    )
    rand!(
        rng,
        PP, XX°, WW°, WW, ρ, v, y1, DD.ismutable(y1);
        Wnr=Wnr, skip=skip
    )
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, v::Val{:ll}, y1::K, m::Val{false};
        Wnr=Wiener(), skip=0
    ) where K
    ll_tot = loglikhd_obs(PP[1], y1)
    for i in eachindex(PP)
        success, ll = rand!(
            rng, PP[i], XX°[i], WW°[i], WW[i], ρ, v, y1, m;
            Wnr=Wnr, skip=skip
        )
        success || return false, ll
        ll_tot += ll
        y1 = XX°[i].x[end]
    end
    true, ll_tot
end

function Random.rand!(
        rng::Random.AbstractRNG,
        PP::AbstractArray{<:GuidProp},
        XX°, WW°, WW, ρ, ::Val{:ll}, y1::K, ::Val{true};
        Wnr=Wiener(), skip=0
    ) where K
    error("in-place not implemented")
end

#===============================================================================
                                Forward guide
===============================================================================#
"""
alias to rand!
"""
const forward_guide! = rand!

#===============================================================================
                                Backward filter
===============================================================================#
"""
alias to `recompute_guiding_term!`
"""
const backward_filter! = recompute_guiding_term!


text:
#testbug_GB16
include("GB17.jl")
using DelimitedFiles
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV
using DataFrames

rng = Random.seed!(2019)

gbdt = GB16.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                 # Learning rate
  num_iterations = 100,                # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB16.GBLearner(gbdt, :regression)

x, y = GB16.test_data(500; p = 600, sigma = sqrt(1))
x_test, y_test = GB16.test_data(500; p = 600, sigma = sqrt(1))
mu = GB16.fried(x)
mu_test = GB16.fried(x_test)

# initialize the hypers
hypers = GB16.Hypers(x, y, 1, 1, 0.8, 2)
hypers.σ = sqrt(1)
hypers.τ = 3 / 10 * var(y) / gbdt.num_iterations
hypers.α = 0.95
hypers.γ = 2.0
hypers.s =  ones(size(x, 2)) / size(x, 2)
hypers

# initialize opt
opts = GB16.Opts()
opts.num_burn = 5

rmse2, Store_s = GB16.TBD_predict_and_rmse!(
    x,
    y,
    x_test,
    y_test,
    mu_test,
    gbdt,
    hypers,
    opts,
    rng)

hypers
sum(hypers.varcount) # 2765(burn 460)
sum(hypers.varcount[1:5]) # 1218(burn 460)
sum(hypers.s[1:5])

hypers.s

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
module ArrayTools

export
    ..,
    ArrayAxes,
    ArrayAxis,
    ArraySize,
    MaybeArrayAxes,
    MaybeArrayAxis,
    RubberIndex,
    all_match,
    allof,
    anyof,
    axis_limits,
    bcastcopy,
    bcastlazy,
    bcastsize,
    cartesian_indices,
    check_size,
    colons,
    common_indices,
    has_standard_indexing,
    noneof,
    promote_eltype,
    reversemap,
    all_indices,
    same_axes,
    same_size,
    same_standard_size,
    split_interval,
    standard_size,
    strictmap!,
    to_axes,
    to_axis,
    to_int,
    to_size,
    to_type,
    # storage trait
    StorageType,
    AnyStorage,
    FlatStorage,
    to_flat_array,
    is_flat_array,
    # Fast arrays and indexing trait:
    IndexingType,
    FastIndexing,
    AnyIndexing,
    to_fast_array,
    is_fast_array

using Base: OneTo, axes1, tail
import Base: dotview, getindex, setindex!, to_indices

include("traits.jl")
include("utils.jl")
include("indexing.jl")
include("broadcasting.jl")
include("rubberindex.jl")

include("PseudoArrays.jl")
using .PseudoArrays

include("AnnotatedArrays.jl")
using .AnnotatedArrays

end # module


text:
"""
$(SIGNATURES)
Sets the number of threads that will be used by BLAS library.
!Notes: Set the thread number with env JULIA_NUM_THREADS to your desired number of threads.
For example: `JULIA_NUM_THREADS=16 julia`

# Arguments
- `nthread` : desired number of threads

# Output:
Shows number of threads set.

"""
function set_blas_threads(nthread)
    LinearAlgebra.BLAS.set_num_threads(nthread)
    core_nums = ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
    println("Number of threads using: $core_nums")
end


text:
struct Burn
    orbit::Orbit
    Δv̄::SVector{3,Float64}
    time::Float64
end

function Burn(time, prograde, normal, radial, orb::Orbit)
    rot = orientation_rotation(time_to_true(time, orb), orb)
    Δv̄ = rot*@SVector([prograde, normal, -radial])
    return Burn(orb, Δv̄, time)
end

# methods for instantaneous orientation (prograde, radial, normal)

"""
    rotmat = orientation_rotation(pos, vel)
    rotmat = orientation_rotation(θ, orb)

Computes the rotation matrix which defines, the prograde, radial, and normal directions
for an orbital state vector `pos`, `vel`.

If a true anomaly `θ` and Orbit `orb` are provided, they are used to compute the state vector,
which will be used to compute the rotation matrix.
"""
orientation_rotation(pos, vel, h=cross(pos, vel)
    ) = basis_rotation(normalize(vel), normalize(h))

orientation_rotation(θ, orb::Orbit) = orientation_rotation(state_vector(θ, orb)...)

"""
    rotmat = time_orientation_rotation(t, orb)

Computes the rotation matrix which defines, the prograde, radial, and normal directions
for the orbital state vector of Orbit `orb` at time `t`.
"""
time_orientation_rotation(t, orb::Orbit) = orientation_rotation(time_to_true(t, orb), orb)

"""
    prograde, normal, radial = orientation_components(brn)

Computes the prograde, normal, and radial components of the Burn `brn`.
"""
function orientation_components(brn::Burn)
    vec = time_orientation_rotation(brn.time, brn.orbit) \ brn.Δv̄
    return vec[1], vec[2], -vec[3]
end

apply_burn(brn::Burn
    ) = Orbit(brn.time, time_orbital_position(brn.time, brn.orbit), time_orbital_velocity(brn.time, brn.orbit)+brn.Δv̄, brn.orbit.primary)

### descriptive display ###

function Base.show(io::IO, brn::Burn)
    o_Δv̄ = orientation_components(brn)
    println(io,
        "Instantaneous burn in orbit above $(brn.orbit.primary.name):\n",
        "  UT:\t\t$(brn.time) s\n",
        "  total Δv:\t$(norm(brn.Δv̄)) m/s\n",
        "  prograde:\t$(o_Δv̄[1]) m/s\n",
        "  normal:\t$(o_Δv̄[2]) m/s\n",
        "  radial:\t$(o_Δv̄[3]) m/s\n",
    )
    return
end


text:
"""
    edges(polygon)

Return iterator over the edges of a polygon.

```jldoctest
julia> collect(edges([1, 2, 3]))
3-element Array{Tuple{Int64,Int64},1}:
 (1, 2)
 (2, 3)
 (3, 1)
```
"""
function edges(polygon)
    zip(polygon, [polygon[2:end]..., polygon[1]])
end


text:
using Tree
using Base.Test

tests = [
          "Tree",
          "node"
        ]

for test in tests
    file_path = joinpath(dirname(@__FILE__), "test_$test.jl")
    println("Running tests for $file_path ...")
    include(file_path)
end


text:
using COSMO
using Plots
using PortfolioOpt
using PortfolioOpt.TestUtils:
    backtest_po, get_test_data, mean_variance,
    percentchange, timestamp, rename!

############ Read Prices #############
Prices = get_test_data()
numD, numA = size(Prices) # A: Assets    D: Days

############# Calculating returns #####################
returns_series = percentchange(Prices)
returns = values(returns_series)
# risk free asset
rf = fill(0.0, numD) # fill(3.2e-4, numD) # readcsv(".\\rf.csv")
rm = mean(returns; dims=2)

############ Efficient frontier ###################

DEFAULT_SOLVER = optimizer_with_attributes(
    COSMO.Optimizer, "verbose" => false, "max_iter" => 9000000
)

# General Parameters
t = 120
k_back = 60
Σ, r̄ = mean_variance(returns[(t - k_back):(t - 1), :])

# Parameter Strategies
δ = 0.028 * 6.05
Δ = std(returns[(end - k_back):end, :]; dims=1)'[:,1] / 5

formulation_markowitz = MeanVariance(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
)
formulation_soyster = RobustBertsimas(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
    uncertainty_delta = Δ,
    bertsimas_budjet = 6.0,
)
formulation_bertsimas = RobustBertsimas(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
    uncertainty_delta = Δ,
    bertsimas_budjet = 3.0,
)
formulation_bental = RobustBenTal(;
    predicted_mean = r̄,
    predicted_covariance = Σ,
    uncertainty_delta = δ
)

# Loop range and log (no robustness)
range_R = 0:0.00005:0.0033; # Used for `po_min_variance_limit_return`
len = size(range_R, 1)
R_noR = zeros(len);
E_noR = zeros(len);
σ_noR = zeros(len);
R_soy = zeros(len);
x_noR = zeros(len, size(r̄, 1));

# Loop no robust
iter = 0
for R_0 in range_R
    global iter += 1
    model = po_min_variance_limit_return(formulation_markowitz, R_0; rf = rf[t])
    x_noR[iter, :], R_noR[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])
    E_noR[iter] = sum(x_noR[iter, :]'r̄)
    σ_noR[iter] = sqrt(sum(x_noR[iter, :]'Σ * x_noR[iter, :]))
end

# Loop range and log
range_R = 0:0.00005:0.0015; # Used for `po_min_variance_limit_return`
len = size(range_R, 1)
R_ber = zeros(len);
E_ber = zeros(len);
σ_ber = zeros(len);
R_soy = zeros(len);
E_soy = zeros(len);
σ_soy = zeros(len);
R_ben_tal = zeros(len);
E_ben_tal = zeros(len);
σ_ben_tal = zeros(len);
x_soy = zeros(len, size(r̄, 1));
x_ber = zeros(len, size(r̄, 1));
x_ben_tal = zeros(len, size(r̄, 1));

# Loop robust
iter = 0
for R_0 in range_R
    global iter += 1
    model = po_min_variance_limit_return(formulation_soyster, R_0; rf = rf[t])
    x_soy[iter, :], R_soy[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])
    model = po_min_variance_limit_return(formulation_bertsimas, R_0; rf = rf[t])
    x_ber[iter, :], R_ber[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])
    model = po_min_variance_limit_return(formulation_bental, R_0; rf = rf[t])
    x_ben_tal[iter, :], R_ben_tal[iter] = compute_solution(model, DEFAULT_SOLVER), value.(model[:R])

    E_ber[iter] = sum(x_ber[iter, :]'r̄)
    E_soy[iter] = sum(x_soy[iter, :]'r̄)
    E_ben_tal[iter] = sum(x_ben_tal[iter, :]'r̄)

    σ_ber[iter] = sqrt(sum(x_ber[iter, :]'Σ * x_ber[iter, :]))
    σ_soy[iter] = sqrt(sum(x_soy[iter, :]'Σ * x_soy[iter, :]))
    σ_ben_tal[iter] = sqrt(sum(x_ben_tal[iter, :]'Σ * x_ben_tal[iter, :]))
end

# plot true frontier
plt = plot(
    σ_noR,
    E_noR;
    title="Efficient Frontier Robust Mean-Variance",
    xlabel="Standard Deviation",
    ylabel="Expected Return",
    label="Markowitz",
    legend=:outertopright,
);
plot!(plt, σ_ber, E_ber; label="Bertsimas");
plot!(plt, σ_soy, E_soy; label="Soyster");
plot!(plt, σ_ben_tal, E_ben_tal; label="BenTal");
plt

# plot believed frontier
plt = plot(
    σ_noR,
    R_noR;
    title="Efficient Frontier Robust Mean-Variance",
    xlabel="Standard Deviation",
    ylabel="Optimized Return (R)",
    label="Markowitz",
    legend=:outertopright,
);
plot!(plt, σ_ber, R_ber; label="Bertsimas");
plot!(plt, σ_soy, R_soy; label="Soyster");
plot!(plt, σ_ben_tal, R_ben_tal; label="BenTal");
plt


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

"""
    DensityRatioEstimator

A method for density ratio estimation.
"""
abstract type DensityRatioEstimator end

"""
    densratio(x_nu, x_de, dre; [optlib])

Estimate density ratio `p_nu(x_de) / p_de(x_de)` with estimator
`dre` and indexable collections of numerator and denominator
samples, `x_nu` and `x_de`.

Optionally choose an optimization library `optlib` from
the list below:

* `JuliaLib`  - Pure Julia implementation
* `OptimLib`  - Implementation with Optim.jl
* `ConvexLib` - Implementation with Convex.jl
* `JuMPLib`   - Implementation with JuMP.jl

See also [`densratiofunc`](@ref).
"""
densratio(x_nu, x_de, dre::DensityRatioEstimator;
          optlib=default_optlib(typeof(dre))) =
  _densratio(x_nu, x_de, dre, optlib)

"""
    densratiofunc(x_nu, x_de, dre; [optlib])

Similar to [`densratio`](@ref), but return a function
`r(x) = p_nu(x) / p_de(x)` that can be evaluated at a
new unseen sample `x`.

See also [`densratio`](@ref).

### Notes

Only some estimators define a ratio function that can
be evaluated outside `x_de`.
"""
densratiofunc(x_nu, x_de, dre::DensityRatioEstimator;
              optlib=default_optlib(typeof(dre))) =
  _densratiofunc(x_nu, x_de, dre, optlib)

"""
    default_optlib(dre)

Return default optimization library for density ratio
estimator `dre`. The function can also be called on the
type `typeof(dre)`.
"""
default_optlib(dre::DensityRatioEstimator) =
  default_optlib(typeof(dre))

"""
    available_optlib(dre)

Return list of implementations available via different
optimization frameworks.
"""
available_optlib(dre::DensityRatioEstimator) =
  available_optlib(typeof(dre))

###################################################
## functions to be implemented by new estimators ##
###################################################

_densratio(x_nu, x_de, dre::DensityRatioEstimator,
           optlib::Type{OptimizationLibrary}) =
  @error "not implemented"

_densratiofunc(x_nu, x_de, dre::DensityRatioEstimator,
               optlib::Type{OptimizationLibrary}) =
  @error "not implemented"

default_optlib(dre::Type{DensityRatioEstimator}) =
  @error "not implemented"

available_optlib(dre::Type{DensityRatioEstimator}) =
  @error "not implemented"