text:
##########################################################################################
# Uncertain values represented by theoretical distributions.
# Base stats can be estimated directly from the distributions, so no need for resampling.
# These methods definitions are just for compatibility with other uncertain values, which
# may not have analytically determined statistics associated with them.
##########################################################################################
StatsBase.mode(uv::TheoreticalDistributionScalarValue, n::Int) = mode(uv)
Statistics.mean(uv::TheoreticalDistributionScalarValue, n::Int) = mean(uv)
Statistics.median(uv::TheoreticalDistributionScalarValue, n::Int) = median(uv)
Statistics.middle(uv::TheoreticalDistributionScalarValue, n::Int) = middle(uv)
Statistics.quantile(uv::TheoreticalDistributionScalarValue, q, n::Int) = quantile(uv, q)
Statistics.std(uv::TheoreticalDistributionScalarValue, n::Int) = std(uv)
Statistics.var(uv::TheoreticalDistributionScalarValue, n::Int) = var(uv)

StatsBase.mode(uv::UncertainScalarTheoreticalFit, n::Int) = mode(uv)
Statistics.mean(uv::UncertainScalarTheoreticalFit, n::Int) = mean(uv)
Statistics.median(uv::UncertainScalarTheoreticalFit, n::Int) = median(uv)
Statistics.middle(uv::UncertainScalarTheoreticalFit, n::Int) = middle(uv)
Statistics.quantile(uv::UncertainScalarTheoreticalFit, q, n::Int) = quantile(uv, q)
Statistics.std(uv::UncertainScalarTheoreticalFit, n::Int) = std(uv)
Statistics.var(uv::UncertainScalarTheoreticalFit, n::Int) = var(uv)


text:
"""
    HMM{Tr<:DiscreteMarkovChain,Em}

Hidden Markov Model with arbitrary transition model (must be a discrete Markov chain) and emission distributions.

# Fields
- `transitions::Tr`: state evolution process.
- `emissions::Vector{Em}`: one emission distribution per state.
"""
struct HiddenMarkovModel{Tr<:DiscreteMarkovChain,Em} <: AbstractMeasure
    transitions::Tr
    emissions::Vector{Em}
end

const HMM = HiddenMarkovModel

## Access

transitions(hmm::HMM) = hmm.transitions
initial_distribution(hmm::HMM) = initial_distribution(transitions(hmm))
transition_matrix(hmm::HMM) = transition_matrix(transitions(hmm))

emissions(hmm::HMM) = hmm.emissions
emission(hmm::HMM, s::Int) = hmm.emissions[s]
nb_states(hmm::HMM) = length(emissions(hmm))

## Simulation

function Base.rand(rng::AbstractRNG, hmm::HMM, T::Integer)
    states = rand(rng, transitions(hmm), T)
    observations = [rand(rng, emission(hmm, states[t])) for t = 1:T]
    return states, observations
end

Base.rand(hmm::HMM, T::Integer) = rand(GLOBAL_RNG, hmm, T)

ext: jl


text:
#Solves Poisson or Helmoltz equation with separable tensorial operator.
# $T_{xy} X_{xy} = B_{xy}$
# Let $T_{xy} = Tx⊗Iy+Ix⊗Ty+λIx⊗Iy$
# via tensorial expression of $T^{-1}_{xy}
# Let $M_x$ and $M_y$ be the diagonalization matrices of $T_x$ and $T_y$ :
# $Tx=M_x I_x M_x^{-1}$ and $Tx=M_y I_x M_y^{-1}$
# ..
# See for example : Plagne, Laurent, and Jean-Yves Berthou.
# "Tensorial basis spline collocation method for Poisson's equation."
# Journal of Computational Physics 157.2 (2000): 419-440.

using LinearAlgebra
export TensorialOperator

struct TridiagOperator{A2D,A1D}
   m::A2D
   mi::A2D
   d::A1D
end

function TridiagOperator(A,T,n,h,cl,α,regularization=false)
   k1=α*K1(n,h,cl)
   if regularization
      k1[1,1]=3/2*k1[1,1];
   end
   tri=SymTridiagonal(collect(diag(k1,0)),collect(diag(k1,1)))
   e=eigen(tri)
   m=e.vectors
   d=e.values
   mi=inv(m)
   TridiagOperator{A{T,2},A{T,1}}(m,mi,d)
end

# Tx⊗Iy+Ix⊗Ty+λIx⊗Iy
struct TensorialOperator{A2D,A1D}
   Tx::TridiagOperator{A2D,A1D}
   Ty::TridiagOperator{A2D,A1D}
   λ::Float64
   nx::Int
   ny::Int
   txy::A2D
   tyx::A2D
   tyx2::A2D
   Dxy::A2D
end

function TensorialOperator(A,T,nx,ny,hx,hy,clx,cly,αx,αy,λ,regularization=false)
   Tx=TridiagOperator(A,T,nx,hx,clx,αx,regularization)
   Ty=TridiagOperator(A,T,ny,hy,cly,αy,regularization)
   Dxy=zeros(nx,ny)

   txd=Array(Tx.d)
   tyd=Array(Ty.d)

   Dxy .= T(1.0) ./ (T(λ) .+ kron(txd,ones(ny)') .+ kron(ones(nx),tyd'))

   TensorialOperator{A{T,2},A{T,1}}(Tx,Ty,λ,nx,ny,zeros(nx,ny),zeros(ny,nx),zeros(ny,nx),Dxy)
end

function solve!(Rxy,Bxy,Txy::TensorialOperator)

   Tx,Ty=Txy.Tx,Txy.Ty
   Mx,Mxm1,Dx=Tx.m,Tx.mi,Tx.d
   My,Mym1,Dy=Ty.m,Ty.mi,Ty.d
   txy,tyx,tyx2=Txy.txy,Txy.tyx,Txy.tyx2
   nx,ny,λ=Txy.nx,Txy.ny,Txy.λ

   LinearAlgebra.mul!(txy,Mxm1,Bxy)
   LinearAlgebra.mul!(tyx,Mym1,transpose(txy))

   tyx .*= (Txy.Dxy)'

   LinearAlgebra.mul!(tyx2,My,tyx)
   LinearAlgebra.mul!(Rxy,Mx,transpose(tyx2))

end

function speye(n)
   SparseMatrixCSC(1.0*SparseArrays.I, n, n)
end

function K1(n,h,a11)
# a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
   h2m1=1/(h^2)
   a=spdiagm(-1 => -ones(n-1), 0 => 2ones(n), 1 => -ones(n-1))
   a[1,1]=a11
   a[n,n]=a11
   a.*=h2m1
   a
end

function K2(n,h,a11,ann)
# a11: Neumann=1, Dirichlet=2, Dirichlet mid=3;
   h2m1=1/(h^2)
   a=spdiagm(-1 => -ones(n-1), 0 => 2ones(n), 1 => -ones(n-1))
   a[1,1]=a11
   a[n,n]=ann
   a.*=h2m1
   a
end


text:
include("./Coupling_coefficients.jl")
include("./Coupling_coefficients_interpolated.jl")
include("./Legendre_polynomials.jl")
include("./Spherical_harmonics.jl")
using Distributions
using Plots

mutable struct Cluster_initial
    #### This structure will hold all the data necessary to represent an astrophysical cluster ####
    #### But I will only use it to represent initial clusters i.e. clusters of 16 disks        ####
    # physical parameters
    G::Float64 # constant of graviation, typically 1.0
    M::Float64 # mass of the central black hole, typically 1.0

    # size parameters of the problem
    n_disks::Int64 # number of disks,typically 16
    N::Int64 # number of stars per disk, typically 512
    lmax::Int64 # number of modes considered in the multipole expansion

    # parameters of the distribution
    m_min::Float64 # minimal mass, typically 1.0
    m_max::Float64 # maximal mass, typically 100.0
    γₘ::Float64 # index of the power law for the distribution in mass, typically -2.0

    a_min::Float64 # minimal semi-major axis, typically 1.0
    a_max::Float64 # maximal semi-major axis, typically 100.0
    γₐ::Float64 # index of the power law for the distribution in semi-major axis, typically 0.0

    e_min::Float64 # minimal eccentricity, typically 0.0
    e_max::Float64 # maximal eccentricity, typically 0.3
    γₑ::Float64 # index of the power law for the distribution in eccentricity, typically 1.0

    σ::Float64 # thickness of the disks i.e. for any star in any disk, we have L_star . L_disk >= (1 - σ). σ is typically worth 0.006
    κ::Float64 # von Mises - Fisher parameter for the distribution of the disks, typically worth 0.0

    # description of the cluster
    K::Array{Tuple{Float64,Float64,Float64}, 1} # List of the orbital parameters of all stars in the cluster -- to be chosen to match the distribution
    L::Array{Tuple{Float64,Float64}, 1} # list of the initial orientation of the stellar annulii -- to be chosen to match the distribution

    # invariants of motion
    Etot::Float64 # total energy of the cluster -- invariant of motion -- to be computed once the orientations are chosen
    Ltot::Float64 # total energy of the cluster -- invariant of motion -- to be computed once the orientations are chosen
end

function Cluster_initial(n_disks::Int64, N::Int64, lmax::Int64;
                        G::Float64 = 1.0, M::Float64 = 1.0,
                        m_min::Float64 = 1.0, m_max::Float64 = 100.0, γₘ::Float64 = -2.0,
                        a_min::Float64 = 1.0, a_max::Float64 = 100.0, γₐ::Float64 = 0.0,
                        e_min::Float64 = 0.0, e_max::Float64 = 0.3, γₑ::Float64 = 1.0,
                        σ::Float64 = 0.006, κ::Float64 = 0.0)
    ############
    # This function creates a cluster structure demanding only the size parameters n_disks, and lmax.
    # To do that, it provides typical values for the the parameters of the distribution, from Szölgyen 2018 page 2 right column
    # Then it samples at random n_disks of N stars from the distribution
    # And finally it calls 2 functions that compute the total energy and total angular momentum
    ############

    # chose the orbital parameters
    K = Array{Tuple{Float64,Float64,Float64}, 1}(undef, n_disks * N)

    for i = 1:n_disks # i will always denote the number of the disk
        for j = 1:N # j will always denote the number of the star, within the disk i
            # draw m, a, e
            m = draw_from_power_law(m_min, m_max, γₘ)
            a = draw_from_power_law(a_min, a_max, γₐ)
            e = draw_from_power_law(e_min, e_max, γₑ)

            index = (i - 1) * N + j # index will always denote the absolute number of the star
            K[index] = (m, a, e) # I hope this is a tuple
        end
    end

    # chose the orbital orientations
    L = Array{Tuple{Float64,Float64}, 1}(undef, n_disks * N)

    # first let's choose the position of the disks' centers
    List_centers = Array{Tuple{Float64,Float64}, 1}(undef, n_disks)

    if (κ == 0)
        for i = 1:n_disks
            θ_disk = acos(2 * (rand() - 0.5)) # I want a uniform distribution on the sphere, not on θ
            ϕ_disk = 2 * π * rand()

            List_centers[i] = (θ_disk, ϕ_disk)
        end

    elseif (κ > 0)
        L_z = [0.0, 0.0, 1.0]
        vMF = VonMisesFisher(L_z, κ) # Defines the distribution we will use to draw the disks' centers.
                                     # With κ = 0, it is the uniform distrbution on the sphere.
        for i = 1:n_disks
            l_drawn = rand(vMF, 1) # this gives the cartesian coordinates

            θ_disk = acos(l_drawn[3]) # WHich I need to translate back to angular coordinates
            if (l_drawn[2] > 0)
                ϕ_disk = acos(l_drawn[1] / sin(θ_disk))
            elseif (l_drawn[2] < 0)
                ϕ_disk = - acos(l_drawn[1] / sin(θ_disk))
            elseif ((l_drawn[2] == 0) && (l_drawn[1] > 0))
                ϕ_disk = 0
            elseif ((l_drawn[2] == 0) && (l_drawn[1] < 0))
                ϕ_disk = π
            elseif ((l_drawn[2] == 0) && (l_drawn[2] == 0))
                ϕ_disk = 0 # I prefer to explicitely break the pole issue, but it should never happen anyways
            else
                println("There is something wrong : $(l_drawn[1]), $(l_drawn[2]), $(l_drawn[3])")
            end

            List_centers[i] = (θ_disk, ϕ_disk) # now I can save it
        end

    else
        println("κ should be positive : κ = $(κ)")
    end

    for i = 1:n_disks
        θ_disk, ϕ_disk = List_centers[i]

        for j = 1:N
            # first I choose point uniformly on the northern polar cap of the standard sphere
            θ_1 = acos(1 - σ * rand())
            ϕ_1 = 2 * π * rand()

            # then I describe this point in standard cartesian coordinates
            x_1 = sin(θ_1) * cos(ϕ_1)
            y_1 = sin(θ_1) * sin(ϕ_1)
            z_1 = cos(θ_1)

            # then I perfom the first rotation θ = 0 -> θ = θ_disk i.e. the rotation around axis y
            x_2 = x_1 * cos(θ_disk) + z_1 * sin(θ_disk)
            y_2 = y_1
            z_2 = - x_1 * sin(θ_disk) + z_1 * cos(θ_disk)

            # then I perform the second rotation ϕ = 0 -> ϕ = ϕ_disk i.e. the rotation around the new axis z
            x_3 = x_2 * cos(ϕ_disk) - y_2 * sin(ϕ_disk)
            y_3 = x_2 * sin(ϕ_disk) + y_2 * cos(ϕ_disk)
            z_3 = z_2

            # then I translate back to angular coordinates
            θ = acos(z_3)

            if (y_3 > 0)
                ϕ = acos(x_3 / sin(θ))
            elseif (y_3 < 0)
                ϕ = 2 * π - acos(x_3 / sin(θ))
            elseif ((x_3 > 0) && (y_3 == 0))
                ϕ = 0.0
            elseif ((x_3 < 0) && (y_3 == 0))
                ϕ = π
            elseif ((x_3 == 0) && (y_3 == 0))
                ϕ = 0.0 # I prefer to explicitely break the pole issue, but it should never happen anyways
            else
                println("There is something wrong : $(x_3), $(y_3), $(z_3)")
            end

            index = (i - 1) * N + j
            L[index] = (θ, ϕ)
        end
    end

    # create a Cluster_initial object with all these infos, by calling the standard initializer of the structure
    Etot_t = Inf # t is for temporary -- This is a mock energy, the correct one will be computed just afterwards
    Ltot_t = Inf # t is for temporary -- This is a mock momentum, the correct one will be computed just afterwards

    clust = Cluster_initial(G, M,
                        n_disks, N, lmax,
                        m_min, m_max, γₘ,
                        a_min, a_max, γₐ,
                        e_min, e_max, γₑ,
                        σ, κ,
                        K, L,
                        Etot_t, Ltot_t)

    # compute the total energy and total momentum of the cluster
    update_Ltot!(clust)
    update_Etot_v3!(clust) # v3 is the fastest of our three algorithms for computing the energy of a cluster

    return(clust)
end

function draw_from_power_law(y_min::Float64, y_max::Float64, γ::Float64)
    #### Returns a Float64. It is a number drawn at random from a power law of index γ between y_min and y_max ####
    #### It assumes that y_min is greater than y_max                                                           ####
    # first draw at random uniformly between 0 and 1
    r = rand()

    # then transform that random variable R into random variable Y that follows the right power law
    y = 0.0 # This is just to help the compiler with scopes
    if (γ > - 1.0)
        x = y_min^(γ + 1) + r * (y_max^(γ + 1) - y_min^(γ + 1))
        y = x^(1 / (γ + 1))
    elseif (γ == -1.0)
        x = log(y_min) + r * log(y_max / y_min)
        y = exp(x)
    else
        x = y_max^(γ + 1) + r * (y_min^(γ + 1) - y_max^(γ + 1))
        y = x^(1 / (γ + 1))
    end

    return(y)
end

function update_Ltot!(clust::Cluster_initial)
    #### Updates the total angular momentum (in norm) of a cluster ####
    # First let's compute the vector total angular momentum
    Ltot_vect = [0.0 0.0 0.0] # Will store the vector total angular momentum

    for i = 1:clust.n_disks
        for j = 1:clust.N
            index = (i - 1) * clust.N + j

            m, a, e = clust.K[index]
            θ, ϕ = clust.L[index]

            x = sin(θ) * cos(ϕ)
            y = sin(θ) * sin(ϕ)
            z = cos(θ)
            l = m * sqrt(clust.G * clust.M * a * (1 - e^(2)))

            Ltot_vect[1] += l * x
            Ltot_vect[2] += l * y
            Ltot_vect[3] += l * z
        end
    end

    # Then let's imagine that we rotate the sphere so that the vector total angular momentum is positively aligned with z
    # in these conditions, we only need to remember the norm of the total angular momentum
    Ltot = sqrt(Ltot_vect[1]^(2) + Ltot_vect[2]^(2) + Ltot_vect[3]^(2))
    clust.Ltot = Ltot
end

function update_Etot!(clust::Cluster_initial)
    #### Updates the total energy of a cluster using the formula Etot = - 0.5 ∑ᵢ ∑ⱼ ∑ₗ (2l + 1) / (4 π) × Jₗ[Kᵢ, Kⱼ] Pₗ(Lᵢ . Lⱼ) ####
    #### and complete coupling coefficients.                                                                                ####
    # Passing by reference the clsuter's fields we will need
    n_disks, N, lmax = clust.n_disks, clust.N, clust.lmax
    L, K = clust.L, clust.K

    # Getting the constants we will need
    lmax_half = div(lmax, 2)
    n_bins = 10^3

    # Preallocating some arrays we will need
    Pₗ = zeros(lmax_half) # Will hold the Legendre polynomials
    Pl = zeros(lmax_half) # Read the comments in coupling_coefficients_v2 for this one
    Cₗ = zeros(lmax_half) # Read the comments in coupling_coefficients_v2 for this one
    R₁ = zeros(n_bins) # Read the comments in coupling_coefficients_v2 for this one
    R₂ = zeros(n_bins) # Read the comments in coupling_coefficients_v2 for this one
    Iₗ = zeros(lmax_half) # Read the comments in coupling_coefficients_v2 for this one
    Jₗ = zeros(lmax_half) # Will hold the coupling coefficients

    # Let's do the triple sum
    Etot = 0.0

    for star_1 = 1:(n_disks * N)
        θ₁, ϕ₁ = L[star_1]
        x₁ = sin(θ₁) * cos(ϕ₁)
        y₁ = sin(θ₁) * sin(ϕ₁)
        z₁ = cos(θ₁)

        K₁ = K[star_1]

        for star_2 = 1:star_1 # we use the symetry of the dot product and of the coupling coefficients
            θ₂, ϕ₂ = L[star_2]
            x₂ = sin(θ₂) * cos(ϕ₂)
            y₂ = sin(θ₂) * sin(ϕ₂)
            z₂ = cos(θ₂)

            K₂ = K[star_2]

            x = x₁ * x₂ + y₁ * y₂ + z₁ * z₂ # This is L₁ . L₂

            Legendre_polynomials!(x, lmax, Pₗ = Pₗ) # Pₗ now holds the (even) Legendre polynomials in L₁ . L₂
            Coupling_coefficients_v2!(K₁, K₂, lmax, Pₗ = Pl, Cₗ = Cₗ, R₁ = R₁, R₂ = R₂, Iₗ = Iₗ, Jₗ = Jₗ) # Jₗ now holds the coupling coefficents between K₁ and K₂

            for l_half = 1:lmax_half
                l = 2 * l_half
                pref = -0.5 *  (2 * l + 1) / (4 * π)

                if (star_1 == star_2) # because we use the symmetry, we have to account twice for each non-diagonal term and once for the diagonal terms
                    Etot += pref * Jₗ[l_half] * Pₗ[l_half]
                else
                    Etot += 2 * pref * Jₗ[l_half] * Pₗ[l_half]
                end
            end
        end
    end

    clust.Etot = Etot
end

function update_Etot_v2!(clust::Cluster_initial)
    #### Updates the total energy of a cluster using the formula                                                                               ####
    #### Etot = - 0.5 ∑ᵢ,ⱼ ∑ₗ,ₘ Yₗₘ(Lᵢ) Yₗₘ(Lⱼ) × (G mᵢ mⱼ 4 π Pₗ(0)² / (2 l + 1)) × ∫dEᵢ / π ∫dEⱼ / π μ(rᵢ) μ(rⱼ) Min(rᵢ, rⱼ)ˡ / Max(rᵢ, rⱼ)ˡ⁺¹  ####
    #### and the multipole algorithm. Should theoretically be the fastest, but with this implementation it isn't                               ####
    #### Please feel free to improve this implementation and use it                                                                            ####

    # Passing by reference the cluster's fields we will need
    n_disks, N, lmax, G = clust.n_disks, clust.N, clust.lmax, clust.G
    K, L = clust.K, clust.L

    # getting the constants we will need
    lmax_half = div(lmax, 2) # number of interesting modes
    n_stars = n_disks * N # number of stars
    n_bins = 10^3 # Number of Riemann nodes in the integrals
    n_nodes = n_stars * n_bins # Number of nodes in the multipole algorithm

    # preallocating the arrays we will need
    Pₗ = zeros(lmax_half) # Will store the Legendre polynomials
    γₗ = zeros(lmax_half) # Will store the prefactors γₗ = 4 π Pₗ(0)² / (2 l + 1)
    ordered_n = Array{Tuple{Float64,Int64}, 1}(undef, n_nodes) # Will store the nodes (r, star) of the multipole algorithm
    Cos_E = zeros(n_bins) # Will store the nodes of the Riemann sums
    P = zeros(2 * lmax + 1, lmax_half) # Will hold the partial sums P of the multipole algorithm
    Q = zeros(2 * lmax + 1, lmax_half) # Will hold the partial sums Q of the multipole algorithm
    δP = zeros(2 * lmax + 1, lmax_half) # Will hold the steps δP of the multipole algorithm
    δQ = zeros(2 * lmax + 1, lmax_half) # Will hold the steps δQ of the multipole algorithm
    S = zeros(2 * lmax + 1, lmax_half) # Will hold the final sums of the multipole algorithm
    Yₗₘ = zeros(2 * lmax + 1, lmax_half) # Will hold the spherical harmonics
    Pₗₘ = zeros(lmax + 1, lmax_half) # Read the comments in spherical_harmonics! for this one

    # Let's compute the prefactors
    prefactor = G / (n_bins^2)

    Legendre_polynomials!(0.0, lmax, Pₗ = Pₗ) # Pₗ now hold the (even) Legendre polynomials in 0
    for l_half = 1:lmax_half
        l = 2 * l_half
        γₗ[l_half] = -0.5 * 4 * π * Pₗ[l_half]^(2) / (2 * l + 1)
    end

    # Let's compute the nodes of the Riemann sums
    for k = 1:n_bins
        Cos_E[k] = cos(π * (2 * k - 1) / (2 * n_bins))  # we use the mid-point rule for faster convergence
    end

    # Let's create the unordered list of nodes (r, star)
    for star = 1:n_stars
        m, a, e = K[star]

        for k = 1:n_bins
            index = (star - 1) * n_bins + k
            cos_E = Cos_E[k]

            r = a * (1 - e * cos_E)

            ordered_n[index] = (r, star)
        end
    end

    # finally let's order the list of nodes -- I think it uses the quicksort or Pivot algorithm
    sort!(ordered_n, by = x -> x[1])

    ### Let's compute the energy

    # let's start with P
    r_min, mass_min, μ_min = 0.0, 0.0, 0.0 # this is just to help with scope
    Yₗₘ_min = zeros(2 * lmax + 1, lmax_half) # this is just to help with scope

    for n = 1:n_nodes
        # acess the parameters of this node
        r, star = ordered_n[n]
        mass, a, e = K[star]
        θ, ϕ = L[star]

        # get μ for this node
        μ = r / a

        # get the spherical harmonics for this node
        Spherical_harmonics!(θ, ϕ, lmax, Yₗₘ = Yₗₘ, Pₗₘ = Pₗₘ)

        for l_half = 1:lmax_half
            l = 2 * l_half

            for m = (-l):l
                Y = Yₗₘ[lmax + 1 + m, l_half]
                Y_min = Yₗₘ_min[lmax + 1 + m, l_half]

                if (n > 1)
                    P[lmax + 1 + m, l_half] *= Y * μ * mass / (Y_min * μ_min * mass_min)
                    P[lmax + 1 + m, l_half] *= (r_min / r)^(l + 1)
                end

                δP[lmax + 1 + m, l_half] = (Y * μ * mass)^(2) / r
                P[lmax + 1 + m, l_half] += δP[lmax + 1 + m, l_half]
                S[lmax + 1 + m, l_half] += P[lmax + 1 + m, l_half]
            end
        end

        r_min, mass_min, μ_min = r, mass, μ
        Yₗₘ_min[:, :] = Yₗₘ
    end

    # let's do the same for Q
    r_sup, mass_sup, μ_sup = 0.0, 0.0, 0.0 # this is just to help with scope
    Yₗₘ_sup = zeros(2 * lmax + 1, lmax_half) # this is just to help with scope

    for n = n_nodes:(-1):1
        # acess the parameters of this node
        r, star = ordered_n[n]
        mass, a, e = K[star]
        θ, ϕ = L[star]

        # get μ for this node
        μ = r / a

        # get the spherical harmonics for this node -- 5 ms
        Spherical_harmonics!(θ, ϕ, lmax, Yₗₘ = Yₗₘ, Pₗₘ = Pₗₘ)

        for l_half = 1:lmax_half
            l = 2 * l_half

            for m = (-l):l
                Y = Yₗₘ[lmax + 1 + m, l_half]
                Y_sup = Yₗₘ_sup[lmax + 1 + m, l_half]

                if (n < n_nodes)
                    Q[lmax + 1 + m, l_half] *= Y * μ * mass / (Y_sup * μ_sup * mass_sup)
                    Q[lmax + 1 + m, l_half] *= (r / r_sup)^(l)

                    δQ[lmax + 1 + m, l_half] = Y * μ * mass * Y_sup * μ_sup * mass_sup * (r / r_sup)^(l) * (1 / r_sup)
                end

                Q[lmax + 1 + m, l_half] += δQ[lmax + 1 + m, l_half]
                S[lmax + 1 + m, l_half] += Q[lmax + 1 + m, l_half]
            end
        end

        r_sup, mass_sup, μ_sup = r, mass, μ
        Yₗₘ_sup[:, :] .= Yₗₘ
    end

    # now we have all we need to compute the energy
    Etot = 0.0

    for l_half = 1:lmax_half
        l = 2 * l_half
        γ_l = γₗ[l_half]

        for m = (-l):l
            Etot += prefactor * γ_l * S[lmax + 1 + m, l_half]
        end
    end

    # Let's update the energy attribute of the cluster (a Float cannot be passed by reference)
    clust.Etot = Etot
end

function update_Etot_v3!(clust::Cluster_initial)
    #### Updates the total energy of a cluster using the formula Etot = - 0.5 ∑ᵢ ∑ⱼ ∑ₗ (2l + 1) / (4 π) × Jₗ[Kᵢ, Kⱼ] Pₗ(Lᵢ . Lⱼ) ####
    #### and interpolated coupling coefficients. This is the fastest of my three methods (for now).                            ####
    # Passing by reference the cluster's fields that we will need
    n_disks, N, lmax = clust.n_disks, clust.N, clust.lmax
    L, K = clust.L, clust.K

    # getting the constants we will need
    lmax_half = div(lmax, 2)
    n_stars = n_disks * N

    # Preallocating an array we will need
    Pₗ = zeros(lmax_half) # will hold the Legendre polynomials

    # Let's do the triple sum
    Etot = 0.0

    for star_1 = 1:n_stars
        K₁ = K[star_1]

        θ₁, ϕ₁ = L[star_1]
        x₁ = sin(θ₁) * cos(ϕ₁)
        y₁ = sin(θ₁) * sin(ϕ₁)
        z₁ = cos(θ₁)

        for star_2 = 1:star_1 # we use the symetry of the dot product and of the coupling coefficients
            K₂ = K[star_2]

            θ₂, ϕ₂ = L[star_2]
            x₂ = sin(θ₂) * cos(ϕ₂)
            y₂ = sin(θ₂) * sin(ϕ₂)
            z₂ = cos(θ₂)

            x = x₁ * x₂ + y₁ * y₂ + z₁ * z₂ # This is L₁ . L₂

            Legendre_polynomials!(x, lmax, Pₗ = Pₗ) # Pₗ now holds the (even) Legendre polynomials at point L₁ . L₂

            for l_half = 1:lmax_half
                l = 2 * l_half

                pref = -0.5 * (2 * l + 1) / (4 * π)
                Jl = Coupling_coefficient_interpolated_v2(l, K₁, K₂) # this is 3 times slower than theory, I don't know why
                Pl = Pₗ[l_half]

                if (star_1 == star_2) # because we use the symmetry, we have to account twice for each non-diagonal term and once for the diagonal terms
                    Etot += pref * Jl * Pl
                else
                    Etot += 2 * pref * Jl * Pl
                end
            end
        end
    end

    # let's update Etot
    clust.Etot = Etot
end

function get_ES_from_clust(clust::Cluster_initial)
    #### Returns a tuple of 2 Float64. They are the normalized total energy and the total angular momentum of the cluster ####
    # getting the constants we will need
    n_stars = length(clust.K)

    # Getting Etot_norm
    J_typ = clust.G * clust.m_min^2 / clust.a_min # typical value of J (circular orbits and twice the same population) -- To be improved
    Etot_norm = clust.Etot / (J_typ * n_stars^2) # Our normalization is Etot_norm = Etot / (J_typ N^2)

    # Getting s (= Ltot_norm)
    Ltot_max = 0.0
    for k = 1:n_stars
        mₖ, aₖ, eₖ = clust.K[k]
        lₖ = mₖ * sqrt(clust.G * clust.M * aₖ * (1 - eₖ^2))
        Ltot_max += lₖ
    end
    s = clust.Ltot / Ltot_max # our normalization is s = Ltot / Ltot_max

    # Returning the tuple
    return((Etot_norm, s))
end

function plot_that_cluster(clust::Cluster_initial)
    #### Returns a plot object. It represents the positions of the stars on the unit sphere, the color being their momentum ####
    plotlyjs()

    ### plot the sphere
    X(z,θ) = sqrt(1 - z^2) * cos(θ)
    Y(z,θ) = sqrt(1 - z^2) * sin(θ)
    Z(z,θ) = z
    C(z, θ) = 0.0 # to impose a uniform color

    zs = range(-1, 1, length = 50)
    θs = range(0, 2 * pi, length = 50)

    my_plot = surface(X.(zs', θs), Y.(zs', θs), Z.(zs', θs),
                        color = "grey", alpha = 0.4,
                        showaxis = false, grid = false)

    ### Plot the stars

    # first we need to compute their cartesian position, and momentums (for the colors)
    x = zeros(clust.n_disks * clust.N)
    y = zeros(clust.n_disks * clust.N)
    z = zeros(clust.n_disks * clust.N)
    Momentums = zeros(clust.n_disks * clust.N)

    for i = 1:clust.n_disks
        for j = 1:clust.N
            index = (i - 1) * clust.N + j

            m, a, e = clust.K[index]
            θ, ϕ = clust.L[index]

            x[index] = sin(θ) * cos(ϕ)
            y[index] = sin(θ) * sin(ϕ)
            z[index] = cos(θ)
            Momentums[index] = log10(m * sqrt(clust.G * clust.M * a * (1 - e^(2))))
        end
    end

    # now we can plot the stars as a scatter plot
    L_min = log10(clust.m_min * sqrt(clust.G * clust.M * clust.a_min * (1 - clust.e_max^(2))))
    L_max = log10(clust.m_max * sqrt(clust.G * clust.M * clust.a_max * (1 - clust.e_min^(2))))

    my_plot = scatter!(x, y, z,
                    xlims = (-1.0, 1.0), ylims = (-1.0, 1.0), zlims = (-1.0, 1.0),
                    zcolor = Momentums, color = :thermal, clims = (L_min, L_max),
                    markersize = 1.0, markerstrokewidth = 0,
                    cbar = false, leg = false)

    # now let's return the figure
    return(my_plot)
end

max_stars_repo_name: NathanMagnan/VRROOMe

max_issues_repo_name: NathanMagnan/VRROOMe

max_forks_repo_name: NathanMagnan/VRROOMe


text:
ENV["MPLBACKEND"]="agg"
using Documenter, ExtendableGrids, Literate, GridVisualize,Pluto
using GridVisualize.FlippableLayout
import PyPlot
using Test

plotting=joinpath(@__DIR__,"..","examples","plotting.jl")
include(plotting)

function rendernotebook(name)
    println("rendernotebook($(name))")
    input=joinpath(@__DIR__,"..","examples",name*".jl")
    output=joinpath(@__DIR__,"src",name*".html")
    session = Pluto.ServerSession();
    notebook = Pluto.SessionActions.open(session, input; run_async=false)
    html_contents = Pluto.generate_html(notebook)
    write(output, html_contents)
end

include("makeplots.jl")

example_md_dir  = joinpath(@__DIR__,"src","examples")

function mkdocs()

    Literate.markdown(plotting, example_md_dir, documenter=false,info=false)
    rendernotebook("plutovista")
    makeplots(example_md_dir, Plotter=PyPlot)
    generated_examples=joinpath.("examples",filter(x->endswith(x, ".md"),readdir(example_md_dir)))
    makedocs(sitename="GridVisualize.jl",
             modules = [GridVisualize],
             doctest = false,
             clean = false,
             authors = "J. Fuhrmann",
             repo="https://github.com/j-fu/GridVisualize.jl",
             pages=[
                 "Home"=>"index.md",
                 "Public API"=> "api.md",
                 "Private API"=> "privapi.md",
                 "Examples" => generated_examples,
             ])
end

mkdocs()

deploydocs(repo = "github.com/j-fu/GridVisualize.jl.git",devbranch = "main")


text:
using naumer_Dimensionality_2022
using Random: MersenneTwister
using LinearAlgebra: I, tr, norm, svd
using CSV
using BenchmarkTools

dim_lst = 1:50
const μ = 1.0
const σ² = 1.0
const τ = 0.05
const stepcount = 1000
const T = 10

const stepsize = 1e4

timing = zeros(length(dim_lst))

function FunctionToTime(dim,μ,σ²,τ,stepcount,T,stepsize)
    # Initialize system
    system = VanDerPolSystem_expanded(μ,dim)
    crlb = I(dimension(system))

    state = ones(dimension(system))

    # Run 1D approximate optimal sampling procedure
    for j in 1:stepcount

        _jacobian = flowJacobian(state, T, system)
        ~, ~, V = svd(_jacobian)
        limitvector = Array(V[:,1])

        # Decide action
        action = limitvector + .01 * randn(length(limitvector))
        action ./= norm(action)
        for step in 1:100
            action .= actiongradientstep_1DApprox(action, limitvector, crlb, σ², stepsize*dimension(system) / step^(2.0/3.0))
        end

        # Update System
        jacobian = flowJacobian(state, τ, system)
        crlb = updateCRLB(crlb, action, jacobian, σ²)
        state = flow(state, τ, system)
    end
    return tr(crlb)
end

FunctionToTime(dim_lst[1],μ,σ²,τ,stepcount,T,stepsize)

for i in 1:length(dim_lst)
    global timing
    timing[i] = @belapsed FunctionToTime($(dim_lst[i]),μ,σ²,τ,stepcount,T,stepsize) samples=10
end

CSV.write("DimensionTiming.csv", ( timing = timing, junk = rand(length(timing))))


text:
module AttitudeWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const AbstractGNSSAttitudeProvider = @jimport org.orekit.gnss.attitude.AbstractGNSSAttitudeProvider
const Attitude = @jimport org.orekit.attitudes.Attitude
const AttitudeProvider = @jimport org.orekit.attitudes.AttitudeProvider
const BeidouGeo = @jimport org.orekit.gnss.attitude.BeidouGeo
const BeidouIGSO = @jimport org.orekit.gnss.attitude.BeidouIGSO
const BeidouMeo = @jimport org.orekit.gnss.attitude.BeidouMeo
const Class = @jimport java.lang.Class
const DerivativeStructure = @jimport org.hipparchus.analysis.differentiation.DerivativeStructure
const ExtendedPVCoordinatesProvider = @jimport org.orekit.utils.ExtendedPVCoordinatesProvider
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldAttitude = @jimport org.orekit.attitudes.FieldAttitude
const FieldDerivativeStructure = @jimport org.hipparchus.analysis.differentiation.FieldDerivativeStructure
const FieldPVCoordinatesProvider = @jimport org.orekit.utils.FieldPVCoordinatesProvider
const FieldTurnSpan = @jimport org.orekit.gnss.attitude.FieldTurnSpan
const Frame = @jimport org.orekit.frames.Frame
const GNSSAttitudeContext = @jimport org.orekit.gnss.attitude.GNSSAttitudeContext
const GNSSAttitudeProvider = @jimport org.orekit.gnss.attitude.GNSSAttitudeProvider
const GNSSFieldAttitudeContext = @jimport org.orekit.gnss.attitude.GNSSFieldAttitudeContext
const GPSBlockIIA = @jimport org.orekit.gnss.attitude.GPSBlockIIA
const GPSBlockIIF = @jimport org.orekit.gnss.attitude.GPSBlockIIF
const GPSBlockIIR = @jimport org.orekit.gnss.attitude.GPSBlockIIR
const Galileo = @jimport org.orekit.gnss.attitude.Galileo
const GenericGNSS = @jimport org.orekit.gnss.attitude.GenericGNSS
const Glonass = @jimport org.orekit.gnss.attitude.Glonass
const JString = @jimport java.lang.String
const Object = @jimport java.lang.Object
const PVCoordinatesProvider = @jimport org.orekit.utils.PVCoordinatesProvider
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const TimeStampedAngularCoordinates = @jimport org.orekit.utils.TimeStampedAngularCoordinates
const TimeStampedFieldAngularCoordinates = @jimport org.orekit.utils.TimeStampedFieldAngularCoordinates
const TurnSpan = @jimport org.orekit.gnss.attitude.TurnSpan

export AbstractGNSSAttitudeProvider
export BeidouGeo
export BeidouIGSO
export BeidouMeo
export FieldTurnSpan
export GNSSAttitudeContext
export GNSSAttitudeProvider
export GNSSFieldAttitudeContext
export GPSBlockIIA
export GPSBlockIIF
export GPSBlockIIR
export Galileo
export GenericGNSS
export Glonass
export TurnSpan
export beta
export beta_ds
export compute_phi
export equals
export get_attitude
export get_class
export get_date
export get_default_yaw_rate
export get_delta_ds
export get_mu_rate
export get_orbit_angle_since_midnight
export get_secured_beta
export get_sv_bcos
export get_turn_duration
export get_turn_end_date
export get_turn_span
export get_turn_start_date
export get_yaw_end
export get_yaw_start
export hash_code
export in_orbit_plane_absolute_angle
export in_sun_side
export in_turn_time_range
export linear_model_still_active
export nominal_yaw
export notify
export notify_all
export orbit_normal_yaw
export set_half_span
export set_up_turn_region
export time_since_turn_start
export to_string
export turn_corrected_attitude
export update_end
export update_start
export validity_end
export validity_start
export wait
export yaw_rate

include("abstract_gnss_attitude_provider.jl")
include("beidou_geo.jl")
include("beidou_igs_o.jl")
include("beidou_meo.jl")
include("field_turn_span.jl")
include("galileo.jl")
include("generic_gns_s.jl")
include("glonass.jl")
include("gnss_attitude_context.jl")
include("gnss_attitude_provider.jl")
include("gnss_field_attitude_context.jl")
include("gps_block_ii_a.jl")
include("gps_block_ii_f.jl")
include("gps_block_ii_r.jl")
include("turn_span.jl")

end


text:
include("LL.jl")
#= include("CLL.jl") =#

function Simple_Smoother(x :: Array{Float64,1},
                  left :: Float64 = 0.25,
                  mid  :: Float64 = 0.5,
                  right :: Float64 = (1- left - mid))

  smo = fill(NaN, length(x))
  smo[2:end-1] = left* x[1:end-2] + mid * x[2:end-1] + right * x[3:end]
  smo[1] = (left * x[1] + right * x[2])/ (left + right)
  smo[end] = (left * x[end-1] + right * x[end])/ (left + right)
  return smo
end
export Simple_Smoother


text:
using Currencies
using FixedPointDecimals

using Test
using Dates

# Get currencies for tests
@usingcurrencies USD, CAD, EUR, GBP, JPY, AUD, INR
@usingcurrencies CNY  # test one-currency version

# Data tests
include("data.jl")

# Basic functionality tests
include("monetary.jl")
include("basket.jl")
include("mixed.jl")

# Currencies
include("data-access.jl")

# Display tests
include("display.jl")

# Custom currencies
include("custom.jl")

# Computations tests
include("valuation.jl")

# README & Doc examples
include("examples.jl")


text:
# multimodal_gmm_langevin_floor
#
# Considers the stein discrepancy for a 2 mixture GMM with gap between modes
# gap > 0. It compares a bad distribution [which only samples from one mode]
# and the actual iid sequence.

include("experiment_utils.jl")

# this controls the seed used for the approximating the wasserstein metric
@parseintcli seed "s" "seed" 7
# the gap between the two modes
@parseintcli gap "g" "gap" 3
# should we use multiple cores?
@parseintcli numcores "n" "numcores" 1
# the dimension of the space
@parseintcli d "d" "dimension" 1
# the gap between the two modes
@parsestringcli discrepancytype "t" "discrepancytype" "kernel"
# the gap between the two modes
@parsestringcli sampler "q" "sampler" "iid"
# parallelize it
if (numcores > 1)
    if (numcores > CPU_CORES)
        error(@sprintf("Requested %d cores, only %d available.", numcores, CPU_CORES))
    end
    addprocs(numcores - 1)
    @assert (numcores == nprocs())
end

using Distributions: MvNormal

using SteinDistributions:
    SteinGMM,
    gradlogdensity,
    cdf
using SteinDiscrepancy:
    SteinInverseMultiquadricKernel,
    ksd,
    gsd,
    wasserstein1d

# select largest n
maxn = 30000
# Select an optimization problem solver
if discrepancytype == "graph"
#    solver = "clp"
    import Gurobi
    solver = Gurobi.GurobiSolver(Threads=numcores)
end
# the number of components
K = 2
# distribution sampled from
distname = "$K-gmm"
# target distribution
gmmmus = gap .* [-0.5, 0.5]
gmmmus = hcat(gmmmus, zeros(2, d-1))
gmmsigmas = eye(d)
gmmweights = [0.5, 0.5]
target = SteinGMM(gmmmus, gmmsigmas)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# Sample sizes at which optimization problem will be solved
ns = vcat([5],
          10:10:min(100,maxn),
          200:100:min(1000,maxn),
          2000:1000:min(10000,maxn),
          20000:10000:min(100000,maxn))

if sampler == "iid"
    X = @setseed rand(target, maxn)
    if discrepancytype == "kernel"
        # this is so the g-h plot in paper are in same range
        X = -X
    end
elseif sampler == "unimodal"
    wrongtarget = MvNormal(vec(gmmmus[1,:]), gmmsigmas)
    X = @setseed rand(wrongtarget, maxn)'
end

println("Beginning optimization!")
for i in ns
    # initialize all the variables
    objectivevalue = nothing
    edgetime = nothing
    solvetime = nothing
    points = X[1:i,:]
    g = nothing
    gradg = nothing
    operatorg = nothing
    if discrepancytype == "graph"
        res = gsd(points=X[1:i,:],
                  gradlogdensity=gradlogp,
                  solver=solver)
        objectivevalue = res.objectivevalue
        edgetime = res.edgetime
        solvetime = res.solvetime
        points = res.points
        g = res.g
        gradg = res.gradg
        operatorg = res.operatorg
    elseif discrepancytype == "kernel"
        kernel = SteinInverseMultiquadricKernel()
        res = ksd(points=X[1:i,:],
                  gradlogdensity=gradlogp,
                  kernel=kernel)
        objectivevalue = sqrt(res.discrepancy2)
        solvetime = res.solvetime
    elseif discrepancytype == "wasserstein"
        @assert (d == 1)
        (objectivevalue, error) = wasserstein1d(X[1:i,:],
                                                targetcdf=x::Float64 -> cdf(target, x))
    end

    println("\tn = $(i), objective = $(objectivevalue)")
    # Package and save results
    instance_data = Dict{Any, Any}(
        "seed" => seed,
        "mus" => gmmmus,
        "sigmas" => gmmsigmas,
        "distweights" => gmmweights,
        "sampler" => sampler,
        "gap" => gap,
        "distname" => distname,
        "X" => points,
        "g" => g,
        "gradg" => gradg,
        "operatorg" => operatorg,
        "n" => i,
        "d" => d,
        "objectivevalue" => objectivevalue,
        "discrepancytype" => discrepancytype,
        "edgetime" => edgetime,
        "solvetime" => solvetime,
        "ncores" => nprocs(),
    )

    save_json(
        instance_data;
        dir="multimodal_gmm_langevin_floor",
        distname=distname,
        sampler=sampler,
        discrepancytype=discrepancytype,
        n=i,
        d=d,
        gap=gap,
        seed=seed,
        numcores=nprocs(),
    )
end

println("COMPLETE!")


text:
# This file was generated, do not modify it. # hide
X = copy(dfX1)
y = copy(dfY1)

coerce!(X, autotype(X, :string_to_multiclass))
yv = Vector(y[:, 1])

LinearRegressorPipe = @pipeline(Standardizer(),
                                OneHotEncoder(drop_last = true),
                                LinearRegressor())

LinearModel = machine(LinearRegressorPipe, X, yv)
fit!(LinearModel)
fp = fitted_params(LinearModel)


text:
using Documenter
using FGenerators

makedocs(
    sitename = "FGenerators",
    format = Documenter.HTML(),
    modules = [FGenerators]
)

deploydocs(;
    repo = "github.com/JuliaFolds/FGenerators.jl",
    push_preview = true,
)


text:
# Command: Find Element
# Refer to src/01_structs/06_Element.jl


text:
import DrWatson: quickactivate, savename
quickactivate(@__DIR__, "Chemostat_Heerden2013")

import SparseArrays

# -------------------------------------------------------------------
# Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in the Julia Pkg REPL to install the
# package, then you must activate the package enviroment (see README)
import Chemostat_Heerden2013
const Hd  = Chemostat_Heerden2013.HeerdenData;
const BD  = Chemostat_Heerden2013.BegData;
const iJR = Chemostat_Heerden2013.iJR904

# -------------------------------------------------------------------
# run add "https://github.com/josePereiro/Chemostat" in the
# julia Pkg REPL for installing the package
import Chemostat
import Chemostat.LP.MathProgBase
const ChU = Chemostat.Utils
const ChSS = Chemostat.SteadyState
const ChLP = Chemostat.LP
const ChEP = Chemostat.MaxEntEP
const ChSU = Chemostat.SimulationUtils

using Plots
using Plots.Measures
using Serialization
using Statistics
using Printf

## -------------------------------------------------------------------
var_order = 49
fname = savename("exploration_", (;var_order), "jld")
fpath = joinpath(iJR.MODEL_PROCESSED_DATA_DIR, fname)
D = deserialize(fpath)
@info "Data info" var_order length(D["epouts"]) D["error"]

model = D["model"]
M, N = size(model)
obj_ider = iJR.BIOMASS_IDER
obj_idx = ChU.rxnindex(model, obj_ider)
glc_ider = "EX_glc_LPAREN_e_RPAREN_"
glc_idx = ChU.rxnindex(model, glc_ider)
epouts = D["epouts"]
sepouts = sort(collect(epouts), by = first)
betas = first.(sepouts)
fbaout = ChLP.fba(model, obj_idx)
ex_glc = fbaout.v[glc_idx]
@assert ex_glc != 0.0
EPS = 1e-8

## -------------------------------------------------------------------
# Meta
const fileid = replace(basename(@__FILE__), ".jl" => "")

## -------------------------------------------------------------------
# Tools
mysavefig(p, fname::String) = (savefig(p, joinpath(iJR.MODEL_FIGURES_DIR, fname)); p)
mysavefig(fname::String, p) = mysavefig(p, fname)
sci(n) = @sprintf("%0.0e", n)
setfontsize!(p, s) = plot!(p; xtickfontsize = s, titlefontsize = s,
    ytickfontsize = s, xguidefontsize = s, yguidefontsize = s, legendfontsize = s)

## -------------------------------------------------------------------
# abs(fba - ep) / ex_glc
let
    p = plot(title = iJR.PROJ_IDER,
        xlabel = "log10( beta )",
        ylabel = "log10( abs( fba - ep ) / ex_glc )"
    )
    f(x) = log10(x)
    mat = zeros(length(sepouts), N)
    for (i, (beta, epout)) in enumerate(sepouts)
        mat[i, :] = (abs.(fbaout.v .- epout.av) ./ abs(ex_glc)) .+ EPS
    end
    plot!(p, first.(sepouts), f.(mat);
        label = "", color = :black, alpha = 0.1
    )
    plot!(p, first.(sepouts), f.(mat[:, obj_idx] .+ EPS);
        label = "", color = :blue, ls = :dash, alpha = 0.8, lw = 3
    )

    # saving
    mysavefig(p, "$(fileid)_fba_ep_diff_vs_beta.png")
end

## -------------------------------------------------------------------
# normalized abs(fba - ep) / max_diff
let
    p1 = plot(title = iJR.PROJ_IDER,
        xlabel = "beta",
        ylabel = "abs( fba - ep ) / max_diff"
    )
    f(x) = x
    mat = zeros(length(sepouts), N)
    for (i, (beta, epout)) in enumerate(sepouts)
        mat[i, :] = (abs.(fbaout.v .- epout.av) ./ abs(ex_glc)) .+ EPS
    end

    norm = ones(length(sepouts)) * maximum.(eachcol(mat))'
    mat ./= norm

    plot!(p1, first.(sepouts), f.(mat);
        label = "", color = :black, alpha = 0.1
    )
    plot!(p1, first.(sepouts), f.(mat[:, obj_idx] .+ EPS);
         label = "", color = :blue, ls = :dash, alpha = 0.8, lw = 3
    )

    # histograms
    kwargs = (;
        normalize = :pdf,
        bins = 200,
        xlabel = "prob. density",
        ylabel = "abs( fba - ep ) / max_diff",
        orientation = :h,
        color = :black,
    )
    p2 = histogram(mat[begin, :];
        label = string("beta: ", round(Int, betas[begin])),
        kwargs...
    )
    p3 = histogram(mat[end, :];
        label = string("beta: ", sci(betas[end])),
        kwargs...
    )
    setfontsize!(p2, 8)
    setfontsize!(p3, 8)

    kwargs = (;
        layout = @layout([a [b ;c]]),
        font = 10
    )
    p = plot(p1, p2, p3; kwargs...)

    # saving
    mysavefig(p, "$(fileid)_normalized_fba_ep_diff_vs_beta.png")
end

## -------------------------------------------------------------------
# stoi err
let
    p = plot(title = iJR.PROJ_IDER,
        xlabel = "beta", ylabel = "log10( stoi_err / ex_glc )")
    f(x) = log10(x)
    mat = zeros(length(sepouts), M)
    for (i, (beta, epout)) in enumerate(sepouts)
        mat[i, :] = (abs.(ChU.stoi_err(model, epout)) ./ abs(ex_glc)) .+ EPS
    end
    plot!(p, first.(sepouts), f.(mat);
        label = "", color = :black, alpha = 0.1
    )
    plot!(p, first.(sepouts), f.(mean.(eachrow(mat)));
        label = "mean", color = :red, alpha = 0.8,
        lw = 5, ls = :dash
    )

    # saving
    mysavefig(p, "$(fileid)_normalized_stoi_err_vs_beta.png")
end


text:
struct UserStrategy <: Strategy end

function BoardGames.getvarsnames(::UserStrategy)
    return String[]
end

BoardGames.name(::UserStrategy) = "User"

Base.copy(s::UserStrategy) = s

BoardGames.getvalues(s::UserStrategy) = ()


text:
using Test
using ITensors

using ITensors: ∑, ∏

@testset "Simple trotterization" begin
  H = ∑{Op}() + ("X", 1) + ("Y", 1)

  s = siteinds("Qubit", 1)

  for nsteps in [10, 100, 1000]
    expHᵉˣᵃᶜᵗ = ITensor(exp(H), s)
    @test expHᵉˣᵃᶜᵗ ≈ ITensor(exp(H; alg=Trotter{1}(nsteps)), s) rtol = 1 / nsteps
    @test expHᵉˣᵃᶜᵗ ≈ ITensor(exp(H; alg=Trotter{2}(nsteps)), s) rtol = (1 / nsteps)^2
    @test expHᵉˣᵃᶜᵗ ≈ ITensor(exp(H; alg=Trotter{4}(nsteps)), s) rtol = (1 / nsteps)^2
    @test expHᵉˣᵃᶜᵗ ≈ ITensor(exp(H; alg=Trotter{8}(nsteps)), s) rtol = (1 / nsteps)^2

    # Convert to ITensors
    t = 1.0
    Uᵉˣᵃᶜᵗ = ∏([ITensor(exp(im * t * H), s)])
    U = ∏{ITensor}(exp(im * t * H; alg=Trotter{2}(nsteps)), s)
    ψ₀ = onehot(s .=> "0")
    Uᵉˣᵃᶜᵗψ₀ = Uᵉˣᵃᶜᵗ(ψ₀)
    Uψ₀ = U(ψ₀)
    @test Uᵉˣᵃᶜᵗψ₀ ≈ Uψ₀ rtol = (1 / nsteps)^2
  end
end

function heisenberg(N)
  os = Sum{Op}()
  for j in 1:(N - 1)
    os += "Sz", j, "Sz", j + 1
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
  end
  return os
end

@testset "Heisenberg Trotter" begin
  N = 4
  ℋ = heisenberg(N)
  s = siteinds("S=1/2", N)
  ψ₀ = MPS(s, n -> isodd(n) ? "↑" : "↓")
  t = 1.0
  for nsteps in [10, 100]
    for order in [1, 2, 4]
      𝒰 = exp(im * t * ℋ; alg=Trotter{order}(nsteps))
      U = ∏{ITensor}(𝒰, s)
      H = ITensor(ℋ, s)
      Uʳᵉᶠψ₀ = replaceprime(exp(im * t * H) * prod(ψ₀), 1 => 0)
      atol = max(1e-6, 1 / nsteps^order)
      @test prod(U(ψ₀)) ≈ Uʳᵉᶠψ₀ atol = atol
    end
  end
end


text:
# Autogenerated wrapper script for ATK_jll for x86_64-apple-darwin
export libatk

using Glib_jll
JLLWrappers.@generate_wrapper_header("ATK")
JLLWrappers.@declare_library_product(libatk, "@rpath/libatk-1.0.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(Glib_jll)
    JLLWrappers.@init_library_product(
        libatk,
        "lib/libatk-1.0.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Formulas for representing and working with linear-model-type expressions
# Original by Harlan D. Harris.  Later modifications by John Myles White
# and Douglas M. Bates.

## Formulas are written as expressions and parsed by the Julia parser.
## For example :(y ~ a + b + log(c))
## In Julia the & operator is used for an interaction.  What would be written
## in R as y ~ a + b + a:b is written :(y ~ a + b + a&b) in Julia.
## The equivalent R expression, y ~ a*b, is the same in Julia

## The lhs of a one-sided formula is 'nothing'
## The rhs of a formula can be 1

type Formula
    lhs::Union(Symbol, Expr, Nothing)
    rhs::Union(Symbol, Expr, Integer)
end

macro ~(lhs, rhs)
    ex = Expr(:call,
              :Formula,
              Base.Meta.quot(lhs),
              Base.Meta.quot(rhs))
    return ex
end

type Terms
    terms::Vector
    eterms::Vector        # evaluation terms
    factors::Matrix{Int8} # maps terms to evaluation terms
    order::Vector{Int}    # orders of rhs terms
    response::Bool        # indicator of a response, which is eterms[1] if present
    intercept::Bool       # is there an intercept column in the model matrix?
end

type ModelFrame
    df::AbstractDataFrame
    terms::Terms
    msng::BitArray
end

type ModelMatrix{T <: Union(Float32, Float64)}
    m::Matrix{T}
    assign::Vector{Int}
end

Base.size(mm::ModelMatrix) = size(mm.m)
Base.size(mm::ModelMatrix, dim...) = size(mm.m, dim...)

function Base.show(io::IO, f::Formula)
    print(io,
          string("Formula: ",
                 f.lhs == nothing ? "" : f.lhs, " ~ ", f.rhs))
end

## Return, as a vector of symbols, the names of all the variables in
## an expression or a formula
function allvars(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    [[allvars(a) for a in ex.args[2:end]]...]
end
allvars(f::Formula) = unique(vcat(allvars(f.rhs), allvars(f.lhs)))
allvars(sym::Symbol) = [sym]
allvars(v::Any) = Array(Symbol, 0)

# special operators in formulas
const specials = Set([:+, :-, :*, :/, :&, :|, :^])

function dospecials(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in specials) return ex end
    excp = copy(ex)
    excp.args = vcat(a1,map(dospecials, ex.args[2:end]))
    if a1 != :* return excp end
    aa = excp.args
    a2 = aa[2]
    a3 = aa[3]
    if length(aa) > 3
        excp.args = vcat(a1, aa[3:end])
        a3 = dospecials(excp)
    end
    :($a2 + $a3 + $a2 & $a3)
end
dospecials(a::Any) = a

const associative = Set([:+,:*,:&])       # associative special operators

## If the expression is a call to the function s return its arguments
## Otherwise return the expression
function ex_or_args(ex::Expr,s::Symbol)
    if ex.head != :call error("Non-call expression encountered") end
    excp = copy(ex)
    a1 = ex.args[1]
    a2 = map(condense, ex.args[2:end])
    if a1 == s return a2 end
    excp.args = vcat(a1, a2)
    excp
end
ex_or_args(a,s::Symbol) = a

## Condense calls like :(+(a,+(b,c))) to :(+(a,b,c))
## Also need to work out how to distribute & over +
function condense(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    a1 = ex.args[1]
    if !(a1 in associative) return ex end
    excp = copy(ex)
    excp.args = vcat(a1, map(x->ex_or_args(x,a1), ex.args[2:end])...)
    excp
end
condense(a::Any) = a

getterms(ex::Expr) = (ex.head == :call && ex.args[1] == :+) ? ex.args[2:end] : ex
getterms(a::Any) = a

ord(ex::Expr) = (ex.head == :call && ex.args[1] == :&) ? length(ex.args)-1 : 1
ord(a::Any) = 1

const nonevaluation = Set([:&,:|])        # operators constructed from other evaluations
## evaluation terms - the (filtered) arguments for :& and :|, otherwise the term itself
function evt(ex::Expr)
    if ex.head != :call error("Non-call expression encountered") end
    if !(ex.args[1] in nonevaluation) return ex end
    filter(x->!isa(x,Number), vcat(map(getterms, ex.args[2:end])...))
end
evt(a) = {a}

function Terms(f::Formula)
    rhs = condense(dospecials(f.rhs))
    tt = getterms(rhs)
    if !isa(tt,AbstractArray) tt = [tt] end
    tt = tt[!(tt .== 1)]             # drop any explicit 1's
    noint = (tt .== 0) | (tt .== -1) # should also handle :(-(expr,1))
    tt = tt[!noint]
    oo = int(map(ord, tt))           # orders of interaction terms
    if !issorted(oo)                 # sort terms by increasing order
        pp = sortperm(oo)
        tt = tt[pp]
        oo = oo[pp]
    end
    etrms = map(evt, tt)
    haslhs = f.lhs != nothing
    if haslhs
        unshift!(etrms, {f.lhs})
        unshift!(oo, 1)
    end
    ev = unique(vcat(etrms...))
    facs = int8(hcat(map(x->(s=Set(x);map(t->int8(t in s), ev)),etrms)...))
    Terms(tt, ev, facs, oo, haslhs, !any(noint))
end

## Default NA handler.  Others can be added as keyword arguments
function na_omit(df::DataFrame)
    cc = complete_cases(df)
    df[cc,:], cc
end

## Trim the pool field of da to only those levels that occur in the refs
function dropUnusedLevels!(da::PooledDataArray)
    rr = da.refs
    uu = unique(rr)
    if length(uu) == length(da.pool) return da end
    T = eltype(rr)
    su = sort!(uu)
    dict = Dict(su, one(T):convert(T,length(uu)))
    da.refs = map(x->dict[x], rr)
    da.pool = da.pool[uu]
    da
end
dropUnusedLevels!(x) = x

function ModelFrame(f::Formula, d::AbstractDataFrame)
    trms = Terms(f)
    df, msng = na_omit(DataFrame(map(x -> d[x], trms.eterms)))
    names!(df, convert(Vector{Symbol}, map(string, trms.eterms)))
    for c in eachcol(df) dropUnusedLevels!(c[2]) end
    ModelFrame(df, trms, msng)
end
ModelFrame(ex::Expr, d::AbstractDataFrame) = ModelFrame(Formula(ex), d)

function StatsBase.model_response(mf::ModelFrame)
    mf.terms.response || error("Model formula one-sided")
    convert(Array, mf.df[bool(mf.terms.factors[:,1])][:,1])
end

function contr_treatment(n::Integer, contrasts::Bool, sparse::Bool, base::Integer)
    if n < 2 error("not enought degrees of freedom to define contrasts") end
    contr = sparse ? speye(n) : eye(n) .== 1.
    if !contrasts return contr end
    if !(1 <= base <= n) error("base = $base is not allowed for n = $n") end
    contr[:,vcat(1:(base-1),(base+1):end)]
end
contr_treatment(n::Integer,contrasts::Bool,sparse::Bool) = contr_treatment(n,contrasts,sparse,1)
contr_treatment(n::Integer,contrasts::Bool) = contr_treatment(n,contrasts,false,1)
contr_treatment(n::Integer) = contr_treatment(n,true,false,1)
cols(v::PooledDataVector) = contr_treatment(length(v.pool))[v.refs,:]
cols(v::DataVector) = reshape(float64(v.data), (length(v),1))
cols(v::Vector) = float64(v)

function isfe(ex::Expr)                 # true for fixed-effects terms
    if ex.head != :call error("Non-call expression encountered") end
    ex.args[1] != :|
end
isfe(a) = true

## Expand the columns in an interaction term
function expandcols(trm::Vector)
    if length(trm) == 1 return float64(trm[1]) end
    if length(trm) == 2
        a = float64(trm[1])
        b = float64(trm[2])
        nca = size(a,2)
        ncb = size(b,2)
        return hcat([a[:,i].*b[:,j] for i in 1:nca, j in 1:ncb]...)
    end
    error("code for 3rd and higher order interactions not yet written")
end

nc(trm::Vector) = *([size(x,2) for x in trm]...)

function ModelMatrix(mf::ModelFrame)
    trms = mf.terms
    aa = {{ones(size(mf.df,1),int(trms.intercept))}}
    asgn = zeros(Int, (int(trms.intercept)))
    fetrms = bool(map(isfe, trms.terms))
    if trms.response unshift!(fetrms,false) end
    ff = trms.factors[:,fetrms]
    ## need to be cautious here to avoid evaluating cols for a factor with many levels
    ## if the factor doesn't occur in the fetrms
    rows = vec(bool(sum(ff,[2])))
    ff = ff[rows,:]
    cc = [cols(x[2]) for x in eachcol(mf.df[:,rows])]
    for j in 1:size(ff,2)
        trm = cc[bool(ff[:,j])]
        push!(aa, trm)
        asgn = vcat(asgn, fill(j, nc(trm)))
    end
    ModelMatrix{Float64}(hcat([expandcols(t) for t in aa]...), asgn)
end

termnames(term::Symbol, col) = [string(term)]
function termnames(term::Symbol, col::PooledDataArray)
    levs = levels(col)
    [string(term, " - ", levs[i]) for i in 2:length(levs)]
end

function coefnames(fr::ModelFrame)
    if fr.terms.intercept
        vnames = UTF8String["(Intercept)"]
    else
        vnames = UTF8String[]
    end
    # Need to only include active levels
    for term in fr.terms.terms
        if isa(term, Expr)
            if term.head == :call && term.args[1] == :&
                a = term.args[2]
                b = term.args[3]
                for lev1 in termnames(a, fr.df[a]), lev2 in termnames(b, fr.df[b])
                    push!(vnames, string(lev1, " & ", lev2))
                end
            else
                error("unrecognized term $term")
            end
        else
            append!(vnames, termnames(term, fr.df[term]))
        end
    end
    return vnames
end


text:
abstract type AbstractSuperPixelAlgorithm end

struct SimpleLinearIterativeClustering <: AbstractSuperPixelAlgorithm
    m::Int
    S::Int
end
const SLIC = SimpleLinearIterativeClustering
SLIC() = SLIC(10)

(alg::SLIC)(img::AbstractArray{<:Number, 2}) = alg(Lab.(Gray.(img)))
(alg::SLIC)(img::AbstractArray{<:Colorant, 2}) = alg(Lab.(img))
function (alg::SLIC)(img::AbstractArray{Lab, 2})
    throw("Not Implemented Error")
end

function _dist(x::Pixel, y::Pixel, alg::SLIC)
    m, S = alg.m, alg.S
    d_lab = Euclidean()(x.color, y.color)
    d_xy = sqrt(sum(@. abs2(x.pos.I - y.pos.I)))
    return d_lab + m/S * d_xy
end


text:
@testset "1016.binary-string-with-substrings-representing-1-to-n.jl" begin
    @test query_string("0110", 3) == true
    @test query_string("0110", 4) == false
end


text:
export get_weights, get_values, get_reference_data

function get_weights(R::Real, past::Day, future::Day, k_gen::Day)
    [ R^(7 / Dates.value(k_gen)) * ones(Dates.value(past)); ones(1 + Dates.value(future)) ]
end

get_weights(R::Real, past::Int, future::Int, k_gen::Int) = get_weights(R, Day(past), Day(future), Day(k_gen))

function get_values(df::DataFrame, present, past::Day, future::Day; column::Symbol = :cases)
    # ToDo: check type of present
    df[ present - past .<= df.days .<= present + future , column]
end

get_values(df::DataFrame, present, past::Int, future::Int; column::Symbol = :cases) = get_values(df, present, Day(past), Day(future); column = column)

function get_reference_data(df::DataFrame; days_col::String, data_col::String, kind::String)
    @assert kind ∈ ( "R", "cases" ) "kind $kind not supported"
    days = df[!, days_col]
    data = df[!, data_col]
    inds = .!ismissing.(data)
    DataFrame( Dict("days" => Vector{Dates.Date}(days[inds]), "$(kind)" => Vector(data[inds]) ) )
end

function get_reference_data(d::Dict; days_col::String, data_col::String, kind::String)
    d_ = deepcopy(d)
    for (name, df) in d
        d_[name] = get_reference_data(df, days_col=days_col, data_col=data_col, kind=kind)
    end
    d_
end


text:
using MRIReco
using Test
using LinearAlgebra
using HTTP
using Random
using StatsBase

include("testTrajectories.jl")
include("testIO.jl")
include("testSimulation.jl")
include("testOperators.jl")
include("testReconstruction.jl")
include("testBrukerFile.jl")
include("testISMRMRD.jl")
include("testCoilsens.jl")


text:
type SharedMemoryModel <: Model
    U::ModelFactor
    P::ModelFactor
    nfactors::Int
    lambda::Float64
    lambdaI::Nullable{ModelFactor}
    Pinv::Nullable{ModelFactor}
end

nusers(model::SharedMemoryModel) = size(model.U, 2)
nitems(model::SharedMemoryModel) = size(model.P, 2)
nfactors(model::SharedMemoryModel) = model.nfactors

function share!(model::SharedMemoryModel)
    isa(model.U, SharedArray) || (model.U = share(model.U))
    isa(model.P, SharedArray) || (model.P = share(model.P))

    lambdaI = get(model.lambdaI)
    isa(lambdaI, SharedArray) || (model.lambdaI = share(lambdaI))

    if !isnull(model.Pinv)
        Pinv = get(model.Pinv)
        isa(Pinv, SharedArray) || (model.Pinv = share(Pinv))
    end
    nothing
end

function localize!(model::SharedMemoryModel)
    isa(model.U, SharedArray) && (model.U = copy(model.U))
    isa(model.P, SharedArray) && (model.P = copy(model.P))
    nothing
end

function clear(model::SharedMemoryModel)
    model.lambdaI = nothing
    model.Pinv = nothing
end

function ensure_loaded(model::SharedMemoryModel)
    model.lambdaI = model.lambda * eye(model.nfactors)
    nothing
end

function pinv(model::SharedMemoryModel)
    if isnull(model.Pinv)
        # since: I = (P * Pt) * inv(P * Pt)
        # Pinv = Pt * inv(P * Pt)
        P = model.P
        PT = P'
        model.Pinv = PT * inv(P * PT)
    end
    get(model.Pinv)
end
mul_pinv(model::SharedMemoryModel, v) = v * pinv(model)

function prep{TI<:Inputs}(inp::TI, nfacts::Int, lambda::Float64)
    ensure_loaded(inp)
    t1 = time()
    @logmsg("preparing inputs...")

    nu = nusers(inp)
    ni = nitems(inp)

    U = zeros(nfacts, nu)
    P = rand(nfacts, ni)
    for idx in 1:ni
        P[1,idx] = mean(all_user_ratings(inp, idx))
    end

    lambdaI = lambda * eye(nfacts)
    model = SharedMemoryModel(U, P, nfacts, lambda, lambdaI, nothing)

    t2 = time()
    @logmsg("prep time: $(t2-t1)")
    model
end


text:
include("euler/euler.jl")
using .BigIntegers: BigInteger
using .Numbers: digits_sum
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int)::Int = digits_sum(prod(map(BigInteger, 1:n)))

compute(10)

compute(100)

@benchmark compute(100)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
function add_project_toml() :: Int64
    touch("Project.toml")

    open("Project.toml", "w") do io_project
        open("../sample/Project.toml") do io_sample
            write(io_project, io_sample)
        end
    end
end

#=function compare_http_header(headers::Array, key::String, value::String) :: Nothing
    for item in headers
        if item[1]==key
            @test item[2]==value
        end
    end
end

function delete_docker_container_image_and_project(;container_id::String, image_id::String) :: Nothing
    run(`docker container stop $container_id`)
    run(`docker image rmi $image_id -f`)
    cd("..")
    rm("demo", recursive=true)
end=#

#=function extract_html_body_content(html_body::Array{UInt8,1}) :: String
    return split(
        split(String(html_body), "
")[2],
        "
"
    )[1]
end=#

function project_settings() :: Nothing
    cd("settings")

    touch("dev.jl")
    open("dev.jl", "w") do io
        write(io, ":dev = true \n")
        write(io, ":server_host = 0.0.0.0\n")
    end

    open("Global.jl", "a+") do io
        write(io, "include(\"dev.jl\")")
    end

    cd("..")
end


text:
RpcApi.SimpleOrganism[RpcApi.SimpleOrganism(0x0000000000005998,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_3(var_2::Int16=-475)
            return var_2
        end
        local var_1::Int8 = 118
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        if var_1 !== var_1
            Creature.stepLeft(o)
            Creature.stepLeft(o)
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        return true
    end),[110,300,98,0,1,1,1],0,1000,1,364621,0x00000036,Dict{Int16,Int16}(),[178,153],12),RpcApi.SimpleOrganism(0x00000000000065db,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_3(var_2::Int16=-475)
            Creature.stepLeft(o)
            return var_2
        end
        local var_1::Int8 = 118
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        var_1 = var_1 * var_1
        if var_1 !== var_1
            Creature.stepUp(o)
            Creature.stepDown(o)
        end
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        return true
    end),[110,300,98,0,1,1,1],0,1000,1,301713,0x0000003c,Dict{Int16,Int16}(),[238,157],12),RpcApi.SimpleOrganism(0x0000000000005f73,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 118
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.stepRight(o)
        Creature.stepRight(o)
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        if var_1 !== var_1
            Creature.stepLeft(o)
        end
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.stepUp(o)
        return true
    end),[110,300,98,0,1,1,1],0,1000,14,234083,0x0000003f,Dict{Int16,Int16}(),[195,172],11),RpcApi.SimpleOrganism(0x0000000000006117,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_1::Int8 = 118
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            local var_6::Int8 = 0
            local var_5::Int64 = 5157083918494790572
            Creature.eatUp(c,o,Int(var_6))
            Creature.stepRight(o)
            Creature.stepLeft(o)
            Creature.stepUp(o)
        end
        Creature.stepDown(o)
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        var_1 = var_1 + var_1
        Creature.stepDown(o)
        return true
    end),[110,300,98,0,1,1,1],0,1000,1,227042,0x0000004b,Dict{Int16,Int16}(),[31,92],21),RpcApi.SimpleOrganism(0x000000000000659b,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_8::Int64 = 804755539837399077
        function func_6(var_5::Bool=false)
            return var_5
        end
        local var_1::Int8 = 118
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        if var_1 <= var_1
        end
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.stepRight(o)
        return true
    end),[110,300,98,0,1,1,1],0,1000,0,223926,0x00000041,Dict{Int16,Int16}(),[142,4],9),RpcApi.SimpleOrganism(0x00000000000063b8,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_3(var_2::String="coZMo3wf")
            Creature.stepRight(o)
            return var_2
        end
        local var_1::Int8 = 118
        function func_6(var_5::String="SYGCsvm3")
            return var_5
        end
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            Creature.stepRight(o)
        end
        Creature.eatRight(c,o,Int(var_1))
        Creature.stepDown(o)
        return true
    end),[110,300,98,0,1,1,1],0,1000,1,223547,0x0000003e,Dict{Int16,Int16}(),[130,21],13),RpcApi.SimpleOrganism(0x00000000000057d9,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_6::String = "ZyhtgXbe"
        function func_5(var_4::Int8=-70)
            Creature.stepLeft(o)
            return var_4
        end
        local var_1::Int8 = 118
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
            Creature.stepUp(o)
        end
        return true
    end),[110,300,98,0,1,1,1],4,1000,1,211404,0x00000041,Dict{Int16,Int16}(),[83,177],11),RpcApi.SimpleOrganism(0x00000000000060ac,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_6::String = "ZyhtgXbe"
        function func_5(var_4::Int8=-70)
            for i::Int8 = 1:div(var_4,_LOOP_STEPS_DIVIDER)
            end
            Creature.eatUp(c,o,Int(var_4))
            return var_4
        end
        local var_1::Int8 = 118
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        Creature.eatDown(c,o,Int(var_1))
        return true
    end),[110,300,98,0,1,1,1],0,1000,1,206309,0x00000041,Dict{Int16,Int16}(),[84,177],11),RpcApi.SimpleOrganism(0x000000000000667a,:(function (c::Config.ConfigData,o::Creature.Organism)
        local var_9::Int64 = 695849987980146412
        function func_5(var_4::Int16=-30639)
            local var_8::Bool = true
            local var_7::Int8 = 64
            var_4 = var_4 + var_4
            return var_4
        end
        local var_1::Int8 = 118
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        Creature.eatLeft(c,o,Int(var_1))
        Creature.eatUp(c,o,Int(var_1))
        return true
    end),[110,300,98,0,1,1,1],0,1000,1,203438,0x00000046,Dict{Int16,Int16}(),[132,182],12),RpcApi.SimpleOrganism(0x00000000000065b2,:(function (c::Config.ConfigData,o::Creature.Organism)
        function func_9(var_7::Int16=2352,var_8::Int64=9133158049229628223)
            if var_7 < var_7
                Creature.stepUp(o)
            end
            return var_7
        end
        local var_1::Int8 = 118
        Creature.eatDown(c,o,Int(var_1))
        Creature.eatRight(c,o,Int(var_1))
        Creature.eatDown(c,o,Int(var_1))
        for i::Int8 = 1:div(var_1,_LOOP_STEPS_DIVIDER)
        end
        Creature.stepRight(o)
        return true
    end),[110,300,98,0,1,1,1],22,1000,1,202868,0x00000050,Dict{Int16,Int16}(),[12,177],9)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# Correlation

#######################################
#
#   Helper functions
#
#######################################

default_laglen(lx::Int) = min(lx-1, int(10log10(lx)))
check_lags(lx::Int, lags::AbstractVector) = (maximum(lags) < lx || error("lags must be less than the sample length."))

function demean_col!{T<:RealFP}(z::Vector{T}, x::Matrix{T}, j::Int, demean::Bool)
    m = size(x, 1)
    @assert m == length(z)
    b = m * (j-1)
    if demean
        s = zero(T)
        for i = 1 : m
            s += x[b + i]
        end
        mv = s / m
        for i = 1 : m
            z[i] = x[b + i] - mv
        end
    else
        copy!(z, 1, x, b+1, m)
    end
    z
end

#######################################
#
#   Auto-correlations
#
#######################################

default_autolags(lx::Int) = 0 : default_laglen(lx)

_autodot{T<:RealFP}(x::Vector{T}, lx::Int, l::Int) = dot(x, 1:lx-l, x, 1+l:lx)

## autocov

function autocov!{T<:RealFP}(r::RealVector, x::Vector{T}, lags::IntegerVector; demean::Bool=true)
    lx = length(x)
    m = length(lags)
    length(r) == m || raise_dimerror()
    check_lags(lx, lags)

    z::Vector{T} = demean ? x - mean(x) : x
    for k = 1 : m  # foreach lag value
        r[k] = _autodot(z, lx, lags[k]) / lx
    end
    return r
end

function autocov!{T<:RealFP}(r::RealMatrix, x::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    size(r) == (m, ns) || raise_dimerror()
    check_lags(lx, lags)

    z = Array(T, lx)
    for j = 1 : ns
        demean_col!(z, x, j, demean)
        for k = 1 : m
            r[k,j] = _autodot(z, lx, lags[k]) / lx
        end
    end
    return r
end

function autocov{T<:Real}(x::Vector{T}, lags::IntegerVector; demean::Bool=true)
    autocov!(Array(fptype(T), length(lags)), float(x), lags; demean=demean)
end

function autocov{T<:Real}(x::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    autocov!(Array(fptype(T), length(lags), size(x,2)), float(x), lags; demean=demean)
end

autocov{T<:Real}(x::VecOrMat{T}; demean::Bool=true) = autocov(x, default_autolags(size(x,1)); demean=demean)

## autocor

function autocor!{T<:RealFP}(r::RealVector, x::Vector{T}, lags::IntegerVector; demean::Bool=true)
    lx = length(x)
    m = length(lags)
    length(r) == m || raise_dimerror()
    check_lags(lx, lags)

    z::Vector{T} = demean ? x - mean(x) : x
    zz = dot(z, z)
    for k = 1 : m  # foreach lag value
        r[k] = _autodot(z, lx, lags[k]) / zz
    end
    return r
end

function autocor!{T<:RealFP}(r::RealMatrix, x::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    size(r) == (m, ns) || raise_dimerror()
    check_lags(lx, lags)

    z = Array(T, lx)
    for j = 1 : ns
        demean_col!(z, x, j, demean)
        zz = dot(z, z)
        for k = 1 : m
            r[k,j] = _autodot(z, lx, lags[k]) / zz
        end
    end
    return r
end

function autocor{T<:Real}(x::Vector{T}, lags::IntegerVector; demean::Bool=true)
    autocor!(Array(fptype(T), length(lags)), float(x), lags; demean=demean)
end

function autocor{T<:Real}(x::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    autocor!(Array(fptype(T), length(lags), size(x,2)), float(x), lags; demean=demean)
end

autocor{T<:Real}(x::VecOrMat{T}; demean::Bool=true) = autocor(x, default_autolags(size(x,1)); demean=demean)

#######################################
#
#   Cross-correlations
#
#######################################

default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l)

_crossdot{T<:RealFP}(x::Vector{T}, y::Vector{T}, lx::Int, l::Int) = (l >= 0 ? dot(x, 1:lx-l, y, 1+l:lx) : dot(x, 1-l:lx, y, 1:lx+l))

## crosscov

function crosscov!{T<:RealFP}(r::RealVector, x::Vector{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    lx = length(x)
    m = length(lags)
    (length(y) == lx && length(r) == m) || raise_dimerror()
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x - mean(x) : x
    zy::Vector{T} = demean ? y - mean(y) : y
    for k = 1 : m  # foreach lag value
        r[k] = _crossdot(zx, zy, lx, lags[k]) / lx
    end
    return r
end

function crosscov!{T<:RealFP}(r::RealMatrix, x::Matrix{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    (length(y) == lx && size(r) == (m, ns)) || raise_dimerror()
    check_lags(lx, lags)

    zx = Array(T, lx)
    zy::Vector{T} = demean ? y - mean(y) : y
    for j = 1 : ns
        demean_col!(zx, x, j, demean)
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / lx
        end
    end
    return r
end

function crosscov!{T<:RealFP}(r::RealMatrix, x::Vector{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    lx = length(x)
    ns = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, ns)) || raise_dimerror()
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x - mean(x) : x
    zy = Array(T, lx)
    for j = 1 : ns
        demean_col!(zy, y, j, demean)
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / lx
        end
    end
    return r
end

function crosscov!{T<:RealFP}(r::AbstractArray{T,3}, x::Matrix{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    lx = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, nx, ny)) || raise_dimerror()
    check_lags(lx, lags)

    # cached (centered) columns of x
    zxs = Array(Vector{T}, 0)
    sizehint(zxs, nx)
    for j = 1 : nx
        xj = x[:,j]
        if demean
            mv = mean(xj)
            for i = 1 : lx
                xj[i] -= mv
            end
        end
        push!(zxs, xj)
    end

    zx = Array(T, lx)
    zy = Array(T, lx)
    for j = 1 : ny
        demean_col!(zy, y, j, demean)
        for i = 1 : nx
            zx = zxs[i]
            for k = 1 : m
                r[k,i,j] = _crossdot(zx, zy, lx, lags[k]) / lx
            end
        end
    end
    return r
end

function crosscov{T<:Real}(x::Vector{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    crosscov!(Array(fptype(T), length(lags)), float(x), float(y), lags; demean=demean)
end

function crosscov{T<:Real}(x::Matrix{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    crosscov!(Array(fptype(T), length(lags), size(x,2)), float(x), float(y), lags; demean=demean)
end

function crosscov{T<:Real}(x::Vector{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    crosscov!(Array(fptype(T), length(lags), size(y,2)), float(x), float(y), lags; demean=demean)
end

function crosscov{T<:Real}(x::Matrix{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    crosscov!(Array(fptype(T), length(lags), size(x,2), size(y,2)), float(x), float(y), lags; demean=demean)
end

crosscov{T<:Real}(x::VecOrMat{T}, y::VecOrMat{T}; demean::Bool=true) = crosscov(x, y, default_crosslags(size(x,1)); demean=demean)

## crosscor

function crosscor!{T<:RealFP}(r::RealVector, x::Vector{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    lx = length(x)
    m = length(lags)
    (length(y) == lx && length(r) == m) || raise_dimerror()
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x - mean(x) : x
    zy::Vector{T} = demean ? y - mean(y) : y
    sc = sqrt(dot(zx, zx) * dot(zy, zy))
    for k = 1 : m  # foreach lag value
        r[k] = _crossdot(zx, zy, lx, lags[k]) / sc
    end
    return r
end

function crosscor!{T<:RealFP}(r::RealMatrix, x::Matrix{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    lx = size(x, 1)
    ns = size(x, 2)
    m = length(lags)
    (length(y) == lx && size(r) == (m, ns)) || raise_dimerror()
    check_lags(lx, lags)

    zx = Array(T, lx)
    zy::Vector{T} = demean ? y - mean(y) : y
    yy = dot(zy, zy)
    for j = 1 : ns
        demean_col!(zx, x, j, demean)
        sc = sqrt(dot(zx, zx) * yy)
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / sc
        end
    end
    return r
end

function crosscor!{T<:RealFP}(r::RealMatrix, x::Vector{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    lx = length(x)
    ns = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, ns)) || raise_dimerror()
    check_lags(lx, lags)

    zx::Vector{T} = demean ? x - mean(x) : x
    zy = Array(T, lx)
    xx = dot(zx, zx)
    for j = 1 : ns
        demean_col!(zy, y, j, demean)
        sc = sqrt(xx * dot(zy, zy))
        for k = 1 : m
            r[k,j] = _crossdot(zx, zy, lx, lags[k]) / sc
        end
    end
    return r
end

function crosscor!{T<:RealFP}(r::AbstractArray{T,3}, x::Matrix{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    lx = size(x, 1)
    nx = size(x, 2)
    ny = size(y, 2)
    m = length(lags)
    (size(y, 1) == lx && size(r) == (m, nx, ny)) || raise_dimerror()
    check_lags(lx, lags)

    # cached (centered) columns of x
    zxs = Array(Vector{T}, 0)
    sizehint(zxs, nx)
    xxs = Array(T, nx)

    for j = 1 : nx
        xj = x[:,j]
        if demean
            mv = mean(xj)
            for i = 1 : lx
                xj[i] -= mv
            end
        end
        push!(zxs, xj)
        xxs[j] = dot(xj, xj)
    end

    zx = Array(T, lx)
    zy = Array(T, lx)
    for j = 1 : ny
        demean_col!(zy, y, j, demean)
        yy = dot(zy, zy)
        for i = 1 : nx
            zx = zxs[i]
            sc = sqrt(xxs[i] * yy)
            for k = 1 : m
                r[k,i,j] = _crossdot(zx, zy, lx, lags[k]) / sc
            end
        end
    end
    return r
end

function crosscor{T<:Real}(x::Vector{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    crosscor!(Array(fptype(T), length(lags)), float(x), float(y), lags; demean=demean)
end

function crosscor{T<:Real}(x::Matrix{T}, y::Vector{T}, lags::IntegerVector; demean::Bool=true)
    crosscor!(Array(fptype(T), length(lags), size(x,2)), float(x), float(y), lags; demean=demean)
end

function crosscor{T<:Real}(x::Vector{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    crosscor!(Array(fptype(T), length(lags), size(y,2)), float(x), float(y), lags; demean=demean)
end

function crosscor{T<:Real}(x::Matrix{T}, y::Matrix{T}, lags::IntegerVector; demean::Bool=true)
    crosscor!(Array(fptype(T), length(lags), size(x,2), size(y,2)), float(x), float(y), lags; demean=demean)
end

crosscor{T<:Real}(x::VecOrMat{T}, y::VecOrMat{T}; demean::Bool=true) = crosscor(x, y, default_crosslags(size(x,1)); demean=demean)

#######################################
#
#   Partial auto-correlations
#
#   TODO: the codes below need cleanup.
#
#######################################

function pacf_regress!{T<:RealFP}(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector, mk::Integer)
    lx = size(X, 1)
    tmpX = ones(T, lx, mk + 1)
    for j = 1 : size(X,2)
        for l = 1 : mk
            for i = 1+l:lx
                tmpX[i,l+1] = X[i-l,j]
            end
        end
        for i = 1 : length(lags)
            l = lags[i]
            sX = sub(tmpX, 1+l:lx, 1:l+1)
            r[i,j] = (cholfact!(sX'sX)\(sX'sub(X, 1+l:lx, j)))[end]
        end
    end
    r
end

function pacf_yulewalker!{T<:RealFP}(r::RealMatrix, X::Matrix{T}, lags::IntegerVector, mk::Integer)
    tmp = Array(T, mk)
    for j = 1 : size(X,2)
        acfs = autocor(X[:,j], 1:mk)
        for i = 1 : length(lags)
            l = lags[i]
            r[i,j] = l == 0 ? one(T) : l == 1 ? acfs[i] : -durbin!(sub(acfs, 1:l), tmp)[l]
        end
    end
end

function pacf!{T<:RealFP}(r::RealMatrix, X::Matrix{T}, lags::IntegerVector; method::Symbol=:regression)
    lx = size(X, 1)
    m = length(lags)
    minlag, maxlag = minmax(lags)
    (0 <= minlag && 2maxlag < lx) || error("Invalid lag value.")
    size(r) == (m, size(X,2)) || raise_dimerror()

    if method == :regression
        pacf_regress!(r, X, lags, maxlag)
    elseif method == :yulewalker
        pacf_yulewalker!(r, X, lags, maxlag)
    else
        error("Invalid method: $method")
    end
    return r
end

function pacf{T<:Real}(X::Matrix{T}, lags::IntegerVector; method::Symbol=:regression)
    pacf!(Array(fptype(T), length(lags), size(X,2)), float(X), lags; method=method)
end

function pacf{T<:Real}(x::Vector{T}, lags::IntegerVector; method::Symbol=:regression)
    vec(pacf(reshape(x, length(x), 1), lags, method=method))
end


text:
function Vector2DFormat()
    return Vector2DFormat(())
end

function Vector2DFormat(arg0::JString, arg1::JString, arg2::JString)
    return Vector2DFormat((JString, JString, JString), arg0, arg1, arg2)
end

function Vector2DFormat(arg0::JString, arg1::JString, arg2::JString, arg3::NumberFormat)
    return Vector2DFormat((JString, JString, JString, NumberFormat), arg0, arg1, arg2, arg3)
end

function Vector2DFormat(arg0::NumberFormat)
    return Vector2DFormat((NumberFormat,), arg0)
end

function format(obj::Vector2DFormat, arg0::JVector, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (JVector, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::VectorFormat, arg0::JVector)
    return jcall(obj, "format", JString, (JVector,), arg0)
end

function get_available_locales(::Type{VectorFormat})
    return jcall(VectorFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_format(obj::VectorFormat)
    return jcall(obj, "getFormat", NumberFormat, ())
end

function get_prefix(obj::VectorFormat)
    return jcall(obj, "getPrefix", JString, ())
end

function get_separator(obj::VectorFormat)
    return jcall(obj, "getSeparator", JString, ())
end

function get_suffix(obj::VectorFormat)
    return jcall(obj, "getSuffix", JString, ())
end

function get_vector2d_format(::Type{Vector2DFormat})
    return jcall(Vector2DFormat, "getVector2DFormat", Vector2DFormat, ())
end

function get_vector2d_format(::Type{Vector2DFormat}, arg0::Locale)
    return jcall(Vector2DFormat, "getVector2DFormat", Vector2DFormat, (Locale,), arg0)
end

function parse(obj::Vector2DFormat, arg0::JString)
    return jcall(obj, "parse", Vector2D, (JString,), arg0)
end

function parse(obj::Vector2DFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", Vector2D, (JString, ParsePosition), arg0, arg1)
end


text:
# 'cstrt_graham_ngon' provides the vertices coordinates (a,b) of a small
# n-gon for n = 2m and m >= 3
function cstrt_graham_ngon(n)
if mod(n,2) == 0 && n>=4
    a = zeros(n-1); b = zeros(n-1);
    for i = 1:Int(n/2-1)
        a[i] = sin(2*pi*i/(n-1))/(2*cos(pi/(2*n-2))); a[n-i] = -a[i];
        b[i] = (1 - cos(2*pi*i/(n-1)))/(2*cos(pi/(2*n-2))); b[n-i] = b[i];
    end
    a[Int(n/2)] = 0; b[Int(n/2)] = 1;
    return(a,b)
end
end

max_stars_repo_name: cbingane/optigon

max_issues_repo_name: cbingane/optigon

max_forks_repo_name: cbingane/optigon


text:
using Redis: RedisConnectionBase, execute_command, flatten

mutable struct Cache
    labels::Vector{String}
    relationshipTypes::Vector{String}
    propertyKeys::Vector{String}
    Cache() = new(Vector{String}(), Vector{String}(), Vector{String}())
end

struct Graph
    id::String
    redis_conn::RedisConnectionBase
    nodes::Dict
    edges::Vector{Edge}
    _cache::Cache
    Graph(id, redis_conn) = new(id, redis_conn, Dict(), Vector{Edge}(), Cache())
end

function addnode!(g::Graph, node::Node)
    g.nodes[node.label] = node
end

function getnode(g::Graph, nodeid::Integer)
    res = query(g, "MATCH (node) WHERE ID(node) = $nodeid RETURN node")
    @assert !isempty(res) "There is no node with id=$nodeid"
    return res.results[1]
end

function call_procedure(g::Graph, procedure::AbstractString, args...; kwargs...)
    params = join(args, ",")
    query_str = "CALL $procedure($params)"

    yield_values = get(kwargs, :y, Vector())
    if length(yield_values) != 0
        values = join(yield_values, ",")
        query_str = query_str * " YIELD $values"
    end

    return query(g, query_str)
end

function lazy_update_cache!(cache::Cache, g::Graph, attr_name::String, required_idx::Int)
    if required_idx > length(getproperty(cache, Symbol(attr_name)))
        updated_value = call_procedure(g, "db.$attr_name").results
        setproperty!(cache, Symbol(attr_name), updated_value)
    end
end

function getlabel(g::Graph, idx::Int)
    required_idx = idx + 1
    lazy_update_cache!(g._cache, g, "labels", required_idx)
    return g._cache.labels[required_idx]
end

function getrelation(g::Graph, idx::Int)
    required_idx = idx + 1
    lazy_update_cache!(g._cache, g, "relationshipTypes", required_idx)
    return g._cache.relationshipTypes[required_idx]
end

function getprop(g::Graph, idx::Int)
    required_idx = idx + 1
    lazy_update_cache!(g._cache, g, "propertyKeys", required_idx)
    return g._cache.propertyKeys[required_idx]
end

function addedge!(g::Graph, edge::Edge)
    @assert haskey(g.nodes, edge.src_node.label) && haskey(g.nodes, edge.dest_node.label)
    push!(g.edges, edge)
end

function commit(g::Graph)
    if (length(g.nodes) == 0 && length(g.edges) == 0)
        return
    end

    items = vcat(string.(values(g.nodes)), string.(g.edges))

    query_str = "CREATE " * join(items, ",")
    return query(g, query_str)
end

function query(g::Graph, q::String)
    response = execute_command(g.redis_conn, flatten(["GRAPH.QUERY", g.id, q, "--compact"]))
    return QueryResult(g, response)
end

function flush!(g::Graph)
    commit(g)
    g.nodes = Dict()
    g.edges = Array{Edge, 1}()
end

function delete(g::Graph)
    execute_command(g.redis_conn, flatten(["GRAPH.DELETE", g.id]))
end

function merge(g::Graph, pattern::String)
    query_str = "MERGE $pattern"
    return query(g, query_str)
end