text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## RandomDevice

const BoolBitIntegerType = Union{Type{Bool},Base.BitIntegerType}
const BoolBitIntegerArray = Union{Array{Bool},Base.BitIntegerArray}

if Sys.iswindows()
    struct RandomDevice <: AbstractRNG
        buffer::Vector{UInt128}

        RandomDevice() = new(Vector{UInt128}(1))
    end

    function rand(rd::RandomDevice, T::BoolBitIntegerType)
        rand!(rd, rd.buffer)
        @inbounds return rd.buffer[1] % T
    end

    function rand!(rd::RandomDevice, A::BoolBitIntegerArray)
        ccall((:SystemFunction036, :Advapi32), stdcall, UInt8, (Ptr{Void}, UInt32),
              A, sizeof(A))
        A
    end
else # !windows
    struct RandomDevice <: AbstractRNG
        file::IOStream
        unlimited::Bool

        RandomDevice(unlimited::Bool=true) =
            new(open(unlimited ? "/dev/urandom" : "/dev/random"), unlimited)
    end

    rand(rd::RandomDevice, T::BoolBitIntegerType)   = read( rd.file, T)
    rand!(rd::RandomDevice, A::BoolBitIntegerArray) = read!(rd.file, A)
end # os-test

"""
    RandomDevice()

Create a `RandomDevice` RNG object.
Two such objects will always generate different streams of random numbers.
"""
RandomDevice

### generation of floats

rand(rng::RandomDevice, ::Type{Close1Open2}) =
    reinterpret(Float64, 0x3ff0000000000000 | rand(rng, UInt64) & 0x000fffffffffffff)

rand(rng::RandomDevice, ::Type{CloseOpen}) = rand(rng, Close1Open2) - 1.0

@inline rand(r::RandomDevice, T::BitFloatType) = rand_generic(r, T)

## MersenneTwister

const MTCacheLength = dsfmt_get_min_array_size()

mutable struct MersenneTwister <: AbstractRNG
    seed::Vector{UInt32}
    state::DSFMT_state
    vals::Vector{Float64}
    idx::Int

    function MersenneTwister(seed, state, vals, idx)
        length(vals) == MTCacheLength && 0 <= idx <= MTCacheLength ||
            throw(DomainError((length(vals), idx),
                      "`length(vals)` and `idx` must be consistent with $MTCacheLength"))
        new(seed, state, vals, idx)
    end
end

MersenneTwister(seed::Vector{UInt32}, state::DSFMT_state) =
    MersenneTwister(seed, state, zeros(Float64, MTCacheLength), MTCacheLength)

"""
    MersenneTwister(seed)

Create a `MersenneTwister` RNG object. Different RNG objects can have
their own seeds, which may be useful for generating different streams
of random numbers.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> x1 = rand(rng, 2)
2-element Array{Float64,1}:
 0.590845
 0.766797

julia> rng = MersenneTwister(1234);

julia> x2 = rand(rng, 2)
2-element Array{Float64,1}:
 0.590845
 0.766797

julia> x1 == x2
true
```
"""
MersenneTwister(seed) = srand(MersenneTwister(Vector{UInt32}(), DSFMT_state()), seed)

function copy!(dst::MersenneTwister, src::MersenneTwister)
    copy!(resize!(dst.seed, length(src.seed)), src.seed)
    copy!(dst.state, src.state)
    copy!(dst.vals, src.vals)
    dst.idx = src.idx
    dst
end

copy(src::MersenneTwister) =
    MersenneTwister(copy(src.seed), copy(src.state), copy(src.vals), src.idx)

==(r1::MersenneTwister, r2::MersenneTwister) =
    r1.seed == r2.seed && r1.state == r2.state && isequal(r1.vals, r2.vals) &&
    r1.idx == r2.idx

hash(r::MersenneTwister, h::UInt) = foldr(hash, h, (r.seed, r.state, r.vals, r.idx))

### low level API

@inline mt_avail(r::MersenneTwister) = MTCacheLength - r.idx
@inline mt_empty(r::MersenneTwister) = r.idx == MTCacheLength
@inline mt_setfull!(r::MersenneTwister) = r.idx = 0
@inline mt_setempty!(r::MersenneTwister) = r.idx = MTCacheLength
@inline mt_pop!(r::MersenneTwister) = @inbounds return r.vals[r.idx+=1]

function gen_rand(r::MersenneTwister)
    dsfmt_fill_array_close1_open2!(r.state, pointer(r.vals), length(r.vals))
    mt_setfull!(r)
end

@inline reserve_1(r::MersenneTwister) = (mt_empty(r) && gen_rand(r); nothing)
# `reserve` allows one to call `rand_inbounds` n times
# precondition: n <= MTCacheLength
@inline reserve(r::MersenneTwister, n::Int) = (mt_avail(r) < n && gen_rand(r); nothing)

### seeding

#### make_seed()

# make_seed produces values of type Vector{UInt32}, suitable for MersenneTwister seeding
function make_seed()
    try
        return rand(RandomDevice(), UInt32, 4)
    catch
        println(STDERR,
                "Entropy pool not available to seed RNG; using ad-hoc entropy sources.")
        seed = reinterpret(UInt64, time())
        seed = hash(seed, UInt64(getpid()))
        try
            seed = hash(seed, parse(UInt64,
                                    read(pipeline(`ifconfig`, `sha1sum`), String)[1:40],
                                    16))
        end
        return make_seed(seed)
    end
end

function make_seed(n::Integer)
    n < 0 && throw(DomainError(n, "`n` must be non-negative."))
    seed = UInt32[]
    while true
        push!(seed, n & 0xffffffff)
        n >>= 32
        if n == 0
            return seed
        end
    end
end

#### srand()

function srand(r::MersenneTwister, seed::Vector{UInt32})
    copy!(resize!(r.seed, length(seed)), seed)
    dsfmt_init_by_array(r.state, r.seed)
    mt_setempty!(r)
    return r
end

srand(r::MersenneTwister=GLOBAL_RNG) = srand(r, make_seed())
srand(r::MersenneTwister, n::Integer) = srand(r, make_seed(n))
srand(seed::Union{Integer,Vector{UInt32}}) = srand(GLOBAL_RNG, seed)

### Global RNG (must be defined after srand)

const GLOBAL_RNG = MersenneTwister(0)

### generation

#### helper functions

# precondition: !mt_empty(r)
@inline rand_inbounds(r::MersenneTwister, ::Type{Close1Open2}) = mt_pop!(r)
@inline rand_inbounds(r::MersenneTwister, ::Type{CloseOpen}) =
    rand_inbounds(r, Close1Open2) - 1.0
@inline rand_inbounds(r::MersenneTwister) = rand_inbounds(r, CloseOpen)

@inline rand_ui52_raw_inbounds(r::MersenneTwister) =
    reinterpret(UInt64, rand_inbounds(r, Close1Open2))
@inline rand_ui52_raw(r::MersenneTwister) = (reserve_1(r); rand_ui52_raw_inbounds(r))

@inline function rand_ui2x52_raw(r::MersenneTwister)
    reserve(r, 2)
    rand_ui52_raw_inbounds(r) % UInt128 << 64 | rand_ui52_raw_inbounds(r)
end

@inline function rand_ui104_raw(r::MersenneTwister)
    reserve(r, 2)
    rand_ui52_raw_inbounds(r) % UInt128 << 52 ⊻ rand_ui52_raw_inbounds(r)
end

rand_ui10_raw(r::MersenneTwister) = rand_ui52_raw(r)
rand_ui23_raw(r::MersenneTwister) = rand_ui52_raw(r)

#### floats

@inline rand(r::MersenneTwister, ::Type{I}) where {I<:FloatInterval} =
    (reserve_1(r); rand_inbounds(r, I))

@inline rand(r::MersenneTwister, T::BitFloatType) = rand_generic(r, T)

#### integers

@inline rand(r::MersenneTwister,
             ::Type{T}) where {T<:Union{Bool,Int8,UInt8,Int16,UInt16,Int32,UInt32}} =
                 rand_ui52_raw(r) % T

function rand(r::MersenneTwister, ::Type{UInt64})
    reserve(r, 2)
    rand_ui52_raw_inbounds(r) << 32 ⊻ rand_ui52_raw_inbounds(r)
end

function rand(r::MersenneTwister, ::Type{UInt128})
    reserve(r, 3)
    xor(rand_ui52_raw_inbounds(r) % UInt128 << 96,
        rand_ui52_raw_inbounds(r) % UInt128 << 48,
        rand_ui52_raw_inbounds(r))
end

rand(r::MersenneTwister, ::Type{Int64})  = reinterpret(Int64,  rand(r, UInt64))
rand(r::MersenneTwister, ::Type{Int128}) = reinterpret(Int128, rand(r, UInt128))

#### arrays of floats

function rand_AbstractArray_Float64!(r::MersenneTwister, A::AbstractArray{Float64},
                                     n=length(A),
                                     ::Type{I}=CloseOpen) where I<:FloatInterval
    # what follows is equivalent to this simple loop but more efficient:
    # for i=1:n
    #     @inbounds A[i] = rand(r, I)
    # end
    m = 0
    while m < n
        s = mt_avail(r)
        if s == 0
            gen_rand(r)
            s = mt_avail(r)
        end
        m2 = min(n, m+s)
        for i=m+1:m2
            @inbounds A[i] = rand_inbounds(r, I)
        end
        m = m2
    end
    A
end

rand!(r::MersenneTwister, A::AbstractArray{Float64}) = rand_AbstractArray_Float64!(r, A)

fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::Type{CloseOpen}) =
    dsfmt_fill_array_close_open!(s, A, n)

fill_array!(s::DSFMT_state, A::Ptr{Float64}, n::Int, ::Type{Close1Open2}) =
    dsfmt_fill_array_close1_open2!(s, A, n)

function rand!(r::MersenneTwister, A::Array{Float64}, n::Int=length(A),
               ::Type{I}=CloseOpen) where I<:FloatInterval
    # depending on the alignment of A, the data written by fill_array! may have
    # to be left-shifted by up to 15 bytes (cf. unsafe_copy! below) for
    # reproducibility purposes;
    # so, even for well aligned arrays, fill_array! is used to generate only
    # the n-2 first values (or n-3 if n is odd), and the remaining values are
    # generated by the scalar version of rand
    if n > length(A)
        throw(BoundsError(A,n))
    end
    n2 = (n-2) ÷ 2 * 2
    if n2 < dsfmt_get_min_array_size()
        rand_AbstractArray_Float64!(r, A, n, I)
    else
        pA = pointer(A)
        align = Csize_t(pA) % 16
        if align > 0
            pA2 = pA + 16 - align
            fill_array!(r.state, pA2, n2, I) # generate the data in-place, but shifted
            unsafe_copy!(pA, pA2, n2) # move the data to the beginning of the array
        else
            fill_array!(r.state, pA, n2, I)
        end
        for i=n2+1:n
            @inbounds A[i] = rand(r, I)
        end
    end
    A
end

@inline mask128(u::UInt128, ::Type{Float16}) =
    (u & 0x03ff03ff03ff03ff03ff03ff03ff03ff) | 0x3c003c003c003c003c003c003c003c00

@inline mask128(u::UInt128, ::Type{Float32}) =
    (u & 0x007fffff007fffff007fffff007fffff) | 0x3f8000003f8000003f8000003f800000

function rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}},
               ::Type{Close1Open2})
    T = eltype(A)
    n = length(A)
    n128 = n * sizeof(T) ÷ 16
    rand!(r, unsafe_wrap(Array, convert(Ptr{Float64}, pointer(A)), 2*n128),
          2*n128, Close1Open2)
    A128 = unsafe_wrap(Array, convert(Ptr{UInt128}, pointer(A)), n128)
    @inbounds for i in 1:n128
        u = A128[i]
        u ⊻= u << 26
        # at this point, the 64 low bits of u, "k" being the k-th bit of A128[i] and "+"
        # the bit xor, are:
        # [..., 58+32,..., 53+27, 52+26, ..., 33+7, 32+6, ..., 27+1, 26, ..., 1]
        # the bits needing to be random are
        # [1:10, 17:26, 33:42, 49:58] (for Float16)
        # [1:23, 33:55] (for Float32)
        # this is obviously satisfied on the 32 low bits side, and on the high side,
        # the entropy comes from bits 33:52 of A128[i] and then from bits 27:32
        # (which are discarded on the low side)
        # this is similar for the 64 high bits of u
        A128[i] = mask128(u, T)
    end
    for i in 16*n128÷sizeof(T)+1:n
        @inbounds A[i] = rand(r, T) + oneunit(T)
    end
    A
end

function rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}},
               ::Type{CloseOpen})
    rand!(r, A, Close1Open2)
    I32 = one(Float32)
    for i in eachindex(A)
        @inbounds A[i] = Float32(A[i])-I32 # faster than "A[i] -= one(T)" for T==Float16
    end
    A
end

rand!(r::MersenneTwister, A::Union{Array{Float16},Array{Float32}}) = rand!(r, A, CloseOpen)

#### arrays of integers

function rand!(r::MersenneTwister, A::Array{UInt128}, n::Int=length(A))
    if n > length(A)
        throw(BoundsError(A,n))
    end
    Af = unsafe_wrap(Array, convert(Ptr{Float64}, pointer(A)), 2n)
    i = n
    while true
        rand!(r, Af, 2i, Close1Open2)
        n < 5 && break
        i = 0
        @inbounds while n-i >= 5
            u = A[i+=1]
            A[n]    ⊻= u << 48
            A[n-=1] ⊻= u << 36
            A[n-=1] ⊻= u << 24
            A[n-=1] ⊻= u << 12
            n-=1
        end
    end
    if n > 0
        u = rand_ui2x52_raw(r)
        for i = 1:n
            @inbounds A[i] ⊻= u << (12*i)
        end
    end
    A
end

# A::Array{UInt128} will match the specialized method above
function rand!(r::MersenneTwister, A::Base.BitIntegerArray)
    n = length(A)
    T = eltype(A)
    n128 = n * sizeof(T) ÷ 16
    rand!(r, unsafe_wrap(Array, convert(Ptr{UInt128}, pointer(A)), n128))
    for i = 16*n128÷sizeof(T)+1:n
        @inbounds A[i] = rand(r, T)
    end
    A
end

#### from a range

@inline function rand_lteq(r::AbstractRNG, randfun, u::U, mask::U) where U<:Integer
    while true
        x = randfun(r) & mask
        x <= u && return x
    end
end

function rand(rng::MersenneTwister, r::UnitRange{T}) where T<:Union{Base.BitInteger64,Bool}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    m = last(r) % UInt64 - first(r) % UInt64
    bw = (64 - leading_zeros(m)) % UInt # bit-width
    mask = (1 % UInt64 << bw) - (1 % UInt64)
    x = bw <= 52 ? rand_lteq(rng, rand_ui52_raw, m, mask) :
                   rand_lteq(rng, rng->rand(rng, UInt64), m, mask)
    (x + first(r) % UInt64) % T
end

function rand(rng::MersenneTwister, r::UnitRange{T}) where T<:Union{Int128,UInt128}
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    m = (last(r)-first(r)) % UInt128
    bw = (128 - leading_zeros(m)) % UInt # bit-width
    mask = (1 % UInt128 << bw) - (1 % UInt128)
    x = bw <= 52  ? rand_lteq(rng, rand_ui52_raw, m % UInt64, mask % UInt64) % UInt128 :
        bw <= 104 ? rand_lteq(rng, rand_ui104_raw, m, mask) :
                    rand_lteq(rng, rng->rand(rng, UInt128), m, mask)
    x % T + first(r)
end

### randjump

"""
    randjump(r::MersenneTwister, jumps::Integer,
             [jumppoly::AbstractString=dSFMT.JPOLY1e21]) -> Vector{MersenneTwister}

Create an array of the size `jumps` of initialized `MersenneTwister` RNG objects. The
first RNG object given as a parameter and following `MersenneTwister` RNGs in the array are
initialized such that a state of the RNG object in the array would be moved forward (without
generating numbers) from a previous RNG object array element on a particular number of steps
encoded by the jump polynomial `jumppoly`.

Default jump polynomial moves forward `MersenneTwister` RNG state by `10^20` steps.
"""
function randjump(mt::MersenneTwister, jumps::Integer, jumppoly::AbstractString)
    mts = MersenneTwister[]
    push!(mts, mt)
    for i in 1:jumps-1
        cmt = mts[end]
        push!(mts, MersenneTwister(copy(cmt.seed), dSFMT.dsfmt_jump(cmt.state, jumppoly)))
    end
    return mts
end

randjump(r::MersenneTwister, jumps::Integer) = randjump(r, jumps, dSFMT.JPOLY1e21)

max_stars_repo_name: nextjournal/julia-lang

max_issues_repo_name: nextjournal/julia-lang

max_forks_repo_name: nextjournal/julia-lang


text:
using Makie, MakieLayout
using AbstractPlotting: px
using GLFW; GLFW.WindowHint(GLFW.FLOATING, 1)
using KernelDensity
using StatsBase
using DataFrames

function defaultaxis(T::Type, args...)
    error("No default axis type defined for plot type $T.")
end

defaultaxis(T::Type{Lines}, args...) = LAxis

struct CreatesLayout end
struct CreatesNoLayout end

struct MarginDensityScatter end

createslayout(x::T) where T = error("""createslayout is not defined for type $T but this trait
    dispatch function is necessary for all recipe types.""")
createslayout(::Type{Lines}) = CreatesNoLayout()
createslayout(::Type{MarginDensityScatter}) = CreatesLayout()

"""Non-mutating general signature before layout trait dispatch"""
function myplot(T::Type, args...; kwargs...)
    myplot(createslayout(T), T, args...; kwargs...)
end

"""Mutating general signature before layout trait dispatch"""
function myplot!(T::Type, args...; kwargs...)
    myplot!(createslayout(T), T, args...; kwargs...)
end

"""
Recipes that don't create a layout, without scene or layout.
Returns created scene, layout, axis and plot object.
"""
function myplot(::CreatesNoLayout, T::Type, args...; kwargs...)
    scene, layout = layoutscene()
    ax = defaultaxis(T, args...)(scene)
    layout[1, 1] = ax
    plotobj = plot!(ax, T, args...)
    (scene = scene, layout = layout, axis = ax, plot = plotobj)
end

"""
Recipes that don't create a layout, scene and layout position.
Returns created axis and plot object.
"""
function myplot!(::CreatesNoLayout, T::Type, scene::Scene, gp::GridPosition, args...; kwargs...)
    ax = defaultaxis(T, args...)(scene)
    gp[] = ax
    plotobj = plot!(ax, T, args...)
    (axis = ax, plot = plotobj)
end

"""
Recipes that don't create a layout, mutating into axis.
Returns created plot object.
"""
function myplot!(::CreatesNoLayout, T::Type, ax::LAxis, args...; kwargs...)
    plotobj = plot!(ax, T, args...)
    plotobj
end

"""
Recipes that create a layout, without scene or layout.
Returns created scene, top layout, plot layout, axes and plot objects.
"""
function myplot(::CreatesLayout, T::Type, args...; kwargs...)
    scene, layout = layoutscene()
    result = myplot!(CreatesLayout(), T, scene, layout[1, 1], args...; kwargs...)
    # (scene = scene, layout = layout, axis)
    (scene = scene, toplayout = layout, layout = result.layout, axes = result.axes, plots = result.plots)
end

"""
Recipes that create a layout, with scene and layout position.
"""
function myplot!(::CreatesLayout, T::Type, scene::Scene, gp::GridPosition, args...; kwargs...)

    axs, layout = create_axes_and_layout(scene, T(), args...)
    gp[] = layout

    plots = myplot!(T, axs, args...; kwargs...)
    # (scene = scene, layout = layout, axis)
    (layout = layout, axes = axs, plots = plots)
end

function Base.display(nt::NamedTuple{(:scene, :layout, :axis, :plot), Tuple{S,G,A,P}}) where {S<:SceneLike, G, A, P}
    display(nt.scene)
end
function Base.display(nt::NamedTuple{(:scene, :toplayout, :layout, :axes, :plots), Tuple{S,TG,G,A,P}}) where {S<:SceneLike, TG,G, A, P}
    display(nt.scene)
end
function Base.display(nt::NamedTuple{(:scene, :axis, :plot), Tuple{S,A,P}}) where {S<:SceneLike, A, P}
    display(nt.scene)
end

function create_axes_and_layout(scene, type::T, args...) where T
    error("create_axes_and_layout is not defined for type $T but that function is necessary for CreatesLayout recipes.")
end

function create_axes_and_layout(scene, ::MarginDensityScatter, args...; xlabel = " ", ylabel = " ")
    mainax = LAxis(scene, xlabel = xlabel, ylabel = ylabel)
    topax = LAxis(scene; xlabelvisible = false, xticklabelsvisible = false, xticksvisible = false)
    rightax = LAxis(scene; ylabelvisible = false, yticklabelsvisible = false, yticksvisible = false)

    linkxaxes!(mainax, topax)
    linkyaxes!(mainax, rightax)

    tightlimits!(topax, Bottom())
    tightlimits!(rightax, Left())

    layout = GridLayout(2, 2; rowsizes = [Relative(1/4), Auto()], colsizes = [Auto(), Relative(1/4)])

    layout[2, 1] = mainax
    layout[1, 1] = topax
    layout[2, 2] = rightax

    (axes = (main = mainax, top = topax, right = rightax), layout = layout)
end

function create_axes_and_layout(scene, ::MarginDensityScatter, df::DataFrame)
    if size(df, 2) != 2
        error("2 columns needed")
    end

    create_axes_and_layout(scene, MarginDensityScatter();
        xlabel = string(names(df)[1]), ylabel = string(names(df)[2]))
end

"""

"""
function myplot!(::Type{MarginDensityScatter}, axs, x, y; kwargs...)

    xkde = kde(x; npoints = 300)
    ykde = kde(y; npoints = 300)

    scat = scatter!(axs.main, x, y; markersize = 10px, color = :red, strokecolor = :black, strokewidth = 1)
    kde1 = poly!(axs.top, Point2f0.(xkde.x, xkde.density), color = (:red, 0.5))
    kde2 = poly!(axs.right, Point2f0.(ykde.density, ykde.x), color = (:red, 0.5))

    autolimits!(axs.main)

    (scatter = scat, xkde = kde1, ykde = kde2)
end

function myplot!(::Type{MarginDensityScatter}, axs, df::DataFrame; kwargs...)
    myplot!(MarginDensityScatter, axs, df[!, 1], df[!, 2])
end

##

scene, toplayout, layout, axs, plotobjs = myplot(MarginDensityScatter, randn(200), randn(200))
myplot!(MarginDensityScatter, scene, toplayout[1, 2], randn(200), randn(200))
ax, line  = myplot!(Lines, scene, toplayout[2, 1:2], randn(100))
myplot!(Lines, ax, randn(100) .+ 5)
scene

##
df = DataFrame(xvalues = randn(100), yvalues = randn(100))
scene, topl, layout, axs, plotobjs = myplot(MarginDensityScatter, df)
myplot!(MarginDensityScatter, axs, randn(100) .+ 3, randn(100) .+ 3)
myplot!(MarginDensityScatter, axs, randn(100) .+ 5, randn(100) .+ 5)
scene

myplot!(MarginDensityScatter, scene, topl[1, 2], randn(100), randn(100));

ext: jl


text:
__precompile__()
module OLPK

using Dierckx
using TwoFAST
using OffsetArrays
using Memoization #Need to go through and memoize many fftlog functions, in progress (Seems to throw an error look into)

const N = 16*1024
const kmino = exp(-25)
const kmaxo = exp(25)
const rmini = exp(-11)#exp(-14.555)
const rmaxi = exp(11)#exp(11)
const r0o = 1/kmaxo
const r0i = 1/rmaxi

export BiasDict13, BiasDictlhd, BiasDict22
export extend_func, eval_extend_func
export Plhd, Pllhd, Plhdmultipole
export ξ
export best_q, henry_q, show_best_q
export Legendre

include("GeneralFunctions/BiasDictConstructors.jl")
include("GeneralFunctions/ExtendedFunctions.jl")
include("GeneralFunctions/Plhd.jl")
include("GeneralFunctions/xi.jl")
include("GeneralFunctions/Bestq.jl")
include("GeneralFunctions/Legendre.jl")

export generateCln
export M

include("P13/Cln.jl")
include("P13/M13.jl")

export generateMnli, generateMnlismall, generateMnliFast

include("P22/Mnli.jl")
include("P22/MnliFast.jl")

export GrowthPowerDict13, P_13ln, I1, I2, I3, I4, I5, P13ggs, P13lggs, P13ggsmultipole

include("P13/P13.jl")

export P22lggs, P22ggsmultipole, intMnlitesting

include("P22/P22.jl")

export BiasDictAll

function BiasDictAll(;b1, bη, bK², bδη, bη², bKKpara, btd, bδΠ2para, bηΠ2para, bΠ2Kpara, bΠ2para, bΠ3para,
    P0_ϵ, b∇2δ, β∇2v, β∂2parav, P2_ϵ, P2_ϵε_η, b2)
    bd13 = BiasDict13(b1=b1, bη=bη, bK²=bK², bδη=bδη, bη²=bη², bKKpara=bKKpara, btd=btd, bδΠ2para=bδΠ2para,
        bηΠ2para=bηΠ2para, bΠ2Kpara=bΠ2Kpara, bΠ2para=bΠ2para, bΠ3para=bΠ3para)
    bd22 = BiasDict22(b1=b1, bη=bη, b2=b2, bK²=bK², bδη=bδη, bη²=bη², bKKpara=bKKpara, bΠ2para=bΠ2para)
    bdlhd = BiasDictlhd(b1=b1, bη=bη, P0_ϵ=P0_ϵ, b∇2δ=b∇2δ, β∇2v=β∇2v, β∂2parav=β∂2parav, P2_ϵ=P2_ϵ, P2_ϵε_η=P2_ϵε_η)

    return merge(bd13, bd22, bdlhd)
end

export PNLOggs

function PNLOggs(Pk, μs; GeneralBiasDict, f)
    k, P22 = P22ggsmultipole(Pk, μs, BiasOperatorDict22=GeneralBiasDict, f=f)
    _, P13 = P13ggs(Pk, μs, BiasOperatorDict13=GeneralBiasDict, f=f)
    Pl = Plhd(Pk, k, μs, BiasOperatorDictlhd=GeneralBiasDict, f=f)
    PNLO = @. Pl + P22 + 2P13

    return k, PNLO
end

export PNLOlggs

function PNLOlggs(Pk, l; GeneralBiasDict, f)
    r, Mnli = generateMnliFast(Pk, GeneralBiasDict, f)
    k, P22l = P22lggs(Pk, l, BiasOperatorDict22=GeneralBiasDict, f=f, Mnli=Mnli, r=r)

    Cln = generateCln(GeneralBiasDict, f)
    _, P13l = P13lggs(Pk, l, BiasOperatorDict13=GeneralBiasDict, f=f, Cln=Cln)

    Pll = Pllhd(Pk, k, l, BiasOperatorDictlhd=GeneralBiasDict, f=f)

    PNLOl = @. Pll + P22l + 2P13l

    return k, PNLOl
end

end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
include("header.jl")

#=
# The Blasius and Falkner-Skan boundary layers

This notebook focuses on the solution of the boundary layer equations when the
external velocity outside the boundary layer has the form of flow past a wedge,
as in the diagram below:

<img src="https://raw.githubusercontent.com/UCLAMAEThreads/MAE150A/master/notebook/WedgeDiagram.svg" alt="wedge diagram" width="600" align="left"/>
=#

#=
For the flow past a wedge with half-angle $\beta/2$, potential flow theory predicts
that the external velocity has the form,
$$U_e(x) = A x^{m}$$
where $A$ is a positive constant and
$$ m = \frac{\beta}{2-\beta}$$

Our focus here is on the boundary layer that develops along the wall of the wedge.
The boundary-layer flow has a **self-similar form** as it develops along the $x$ direction. Just as the Rayleigh flow develops with the same profile shape in time, with a time-varying thickness $\delta(t) = (\nu t)^{1/2}$, the boundary layer's profile has the same shape at every $x$ location, but with **spatially-varying thickness** proportional to
$$\delta(x) = \left( \frac{\nu x}{U_e(x)}\right)^{1/2}$$

This means that the velocity tangent to the wall can be written in the form
$$ u(x,y) = U_e(x) F'(\eta) $$
where
$$\eta = \dfrac{y}{\delta(x)}$$

With this form, the boundary layer equations reduce to an ordinary differential
equation for $F$ in this case,

$$ F''' + \frac{1}{2} (m+1) F F'' + m (1-F'^2)  = 0$$

This is called the **Falkner-Skan equation**. The equation is non-linear, however,
and must be solved computationally. The tools to solve it are provided for you here.
=#

#=
#### Boundary conditions
The boundary conditions on the flow are no-slip ($u = 0$ at $y = 0$), no-flow-through
($v = 0$ at $y = 0$). Also, the flow must approach the external velocity outside
the boundary layer. This means $u \rightarrow U_e$ as $y \rightarrow \infty$.
=#

#=
#### Special cases
The case $\beta = 0$ (so that $m = 0$) corresponds to a flat plate, parallel to
free stream, and $U_e(x) = U_\infty$. This is called a **Blasius boundary layer**.

When $\beta > 0$, the external flow is accelerating, and the pressure gradient
is said to be **favorable**.

When $\beta < 0$, the external flow is decelerating, and the pressure gradient
is adverse. The lowest admissible value of $\beta$ is approximately $-0.1999$.
Below this, the boundary layer equations do not hold.

#### Other boundary-layer quantities
*Every* measure of thickness of the boundary layer is proportional to $\delta(x)$.
However, we need the solution of the Falkner-Skan equation to determine what
the constant of proportionality is.

For example, the 99% thickness is

$$\delta_{99}(x) = C_{1}\delta(x)$$

the displacement thickness is

$$\delta^{*}(x) = C_{2} \delta(x)$$

and the momentum thickness is

$$\theta(x) = C_{3} \delta(x)$$

These constants $C_1$, $C_2$ and $C_3$ depend on $\beta$, and they can only be
found by solving the Falkner-Skan problem.

Another quantity of interest is the **skin friction coefficient**,

$$c_{f} = \dfrac{\tau_{w}}{\frac{1}{2} \rho U_e^2} $$

It can be shown that the skin friction coefficient of a Falkner-Skan boundary layer must be of the form

$$c_{f} = \dfrac{C_4}{(U_e(x) x/\nu)^{1/2}} $$

As with the other constants, $C_4$ can only be found by solving the Falkner-Skan
equation.
=#

# ### Set up the module
using MAE250B
#-
using Plots

#=
### Generate the velocity profile and other solution features
Here, we will compute the self-similar velocity profile and compute other
boundary-layer quantities. The only required input for this `falknerskan` function
is $\beta$. It returns the all of the aspects of the solution described above.
=#

β = 0.0
sol = falknerskan(β);

#=
For example, the proportionality constant on the 99 percent thickness is
=#
d99(sol)

#=
which means that, for $\beta = 0$,

$$\delta_{99}(x) = 4.91 \sqrt{\dfrac{\nu x}{U_e}} $$

or, in other words,

$$\dfrac{\delta_{99}(x)}{x} = \dfrac{4.91}{\sqrt{Re_x}} $$

Similarly, the momentum thickness:
=#
theta(sol)

#=
so that

$$\dfrac{\theta(x)}{x} = \dfrac{0.664}{\sqrt{Re_x}} $$
=#

#=
#### Let's plot the profile
=#
plot(u(sol),η(sol),xlim=(0,2),ylim=(0,Inf),xlabel=L"F'(\eta) = u/U_e",ylabel=L"\eta = y/\delta(x)")

#=
The vertical velocity profile is
=#
plot(v(sol),η(sol),xlim=(0,2),ylim=(0,Inf),xlabel=L"v/U_e/Re_x^{1/2}",ylabel=L"\eta = y/\delta(x)")

#=
#### Drag force
When $\beta = 0$ (Blasius), the drag coefficient on a plate of length $L$ is equal to

$$C_D = \dfrac{4\theta(L)}{L} = \dfrac{2.656}{\sqrt{Re_L}}$$
=#

#=
If you wish to specify a wall velocity, you can use the optional argument `Vw = `,
e.g.,
=#
sol = falknerskan(β,Vw=0.2)
plot(u(sol),η(sol),xlim=(0,2),ylim=(0,Inf),xlabel=L"F'(\eta) = u/U_e",ylabel=L"\eta = y/\delta(x)")


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# This file contains markdown extensions designed to make documenting
# Julia easy peasy.
#
# We start by borrowing GitHub's `fencedcode` extension – more to follow.

include("interp.jl")

@flavor julia [blocktex, blockinterp, hashheader, list, indentcode, fencedcode,
               blockquote, admonition, footnote, github_table, horizontalrule, setextheader, paragraph,

               linebreak, escapes, tex, interp, en_dash, inline_code,
               asterisk_bold, underscore_bold, asterisk_italic, underscore_italic, image, footnote_link, link, autolink]

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
using JuMP, GLPK, Cbc

function graph_coloring(;fast=true, benchmark=false)
    #m = Model(with_optimizer(GLPK.Optimizer, msg_lev = 3))
    m = Model(with_optimizer(Cbc.Optimizer, logLevel=0))

    @variable(m, max_color, Int)
    @variable(m, c[1:49,1:4], Bin)

    @objective(m, Min, max_color)

    for s=1:49
        @constraint(m, sum(c[s,1:4]) == 1)
    end

    if fast
        for s=1:49
            @constraint(m, sum(c[s,i]*i for i=1:4) <= max_color)
        end
    else
        for s=1:49, i=1:4
            @constraint(m, max_color >= i*c[s,i])
        end
    end

    for i=1:4
        @constraint(m, c[1,i] + c[11,i] <= 1)
        @constraint(m, c[1,i] + c[8,i] <= 1)
        @constraint(m, c[11,i] + c[8,i] <= 1)
        @constraint(m, c[11,i] + c[22,i] <= 1)
        @constraint(m, c[11,i] + c[24,i] <= 1)
        @constraint(m, c[24,i] + c[22,i] <= 1)
        @constraint(m, c[24,i] + c[36,i] <= 1)
        @constraint(m, c[22,i] + c[8,i] <= 1)
        @constraint(m, c[22,i] + c[23,i] <= 1)
        @constraint(m, c[22,i] + c[36,i] <= 1)
        @constraint(m, c[8,i] + c[2,i] <= 1)
        @constraint(m, c[8,i] + c[6,i] <= 1)
        @constraint(m, c[8,i] + c[23,i] <= 1)
        @constraint(m, c[23,i] + c[6,i] <= 1)
        @constraint(m, c[23,i] + c[31,i] <= 1)
        @constraint(m, c[23,i] + c[41,i] <= 1)
        @constraint(m, c[23,i] + c[36,i] <= 1)
        @constraint(m, c[36,i] + c[31,i] <= 1)
        @constraint(m, c[36,i] + c[41,i] <= 1)
        @constraint(m, c[2,i] + c[4,i] <= 1)
        @constraint(m, c[2,i] + c[5,i] <= 1)
        @constraint(m, c[2,i] + c[6,i] <= 1)
        @constraint(m, c[6,i] + c[5,i] <= 1)
        @constraint(m, c[6,i] + c[15,i] <= 1)
        @constraint(m, c[6,i] + c[31,i] <= 1)
        @constraint(m, c[31,i] + c[15,i] <= 1)
        @constraint(m, c[31,i] + c[33,i] <= 1)
        @constraint(m, c[31,i] + c[37,i] <= 1)
        @constraint(m, c[31,i] + c[41,i] <= 1)
        @constraint(m, c[41,i] + c[37,i] <= 1)
        @constraint(m, c[41,i] + c[40,i] <= 1)
        @constraint(m, c[4,i] + c[10,i] <= 1)
        @constraint(m, c[4,i] + c[5,i] <= 1)
        @constraint(m, c[5,i] + c[10,i] <= 1)
        @constraint(m, c[5,i] + c[13,i] <= 1)
        @constraint(m, c[5,i] + c[15,i] <= 1)
        @constraint(m, c[15,i] + c[13,i] <= 1)
        @constraint(m, c[15,i] + c[35,i] <= 1)
        @constraint(m, c[15,i] + c[33,i] <= 1)
        @constraint(m, c[33,i] + c[35,i] <= 1)
        @constraint(m, c[33,i] + c[37,i] <= 1)
        @constraint(m, c[37,i] + c[46,i] <= 1)
        @constraint(m, c[37,i] + c[40,i] <= 1)
        @constraint(m, c[40,i] + c[46,i] <= 1)
        @constraint(m, c[40,i] + c[47,i] <= 1)
        @constraint(m, c[10,i] + c[7,i] <= 1)
        @constraint(m, c[10,i] + c[13,i] <= 1)
        @constraint(m, c[13,i] + c[7,i] <= 1)
        @constraint(m, c[13,i] + c[26,i] <= 1)
        @constraint(m, c[13,i] + c[35,i] <= 1)
        @constraint(m, c[35,i] + c[26,i] <= 1)
        @constraint(m, c[35,i] + c[32,i] <= 1)
        @constraint(m, c[35,i] + c[39,i] <= 1)
        @constraint(m, c[35,i] + c[46,i] <= 1)
        @constraint(m, c[46,i] + c[39,i] <= 1)
        @constraint(m, c[46,i] + c[43,i] <= 1)
        @constraint(m, c[46,i] + c[47,i] <= 1)
        @constraint(m, c[47,i] + c[43,i] <= 1)
        @constraint(m, c[7,i] + c[26,i] <= 1)
        @constraint(m, c[7,i] + c[49,i] <= 1)
        @constraint(m, c[26,i] + c[21,i] <= 1)
        @constraint(m, c[26,i] + c[32,i] <= 1)
        @constraint(m, c[32,i] + c[21,i] <= 1)
        @constraint(m, c[32,i] + c[25,i] <= 1)
        @constraint(m, c[32,i] + c[29,i] <= 1)
        @constraint(m, c[32,i] + c[34,i] <= 1)
        @constraint(m, c[32,i] + c[39,i] <= 1)
        @constraint(m, c[39,i] + c[34,i] <= 1)
        @constraint(m, c[39,i] + c[38,i] <= 1)
        @constraint(m, c[39,i] + c[44,i] <= 1)
        @constraint(m, c[39,i] + c[42,i] <= 1)
        @constraint(m, c[39,i] + c[43,i] <= 1)
        @constraint(m, c[43,i] + c[42,i] <= 1)
        @constraint(m, c[49,i] + c[21,i] <= 1)
        @constraint(m, c[49,i] + c[25,i] <= 1)
        @constraint(m, c[21,i] + c[25,i] <= 1)
        @constraint(m, c[42,i] + c[44,i] <= 1)
        @constraint(m, c[42,i] + c[48,i] <= 1)
        @constraint(m, c[25,i] + c[17,i] <= 1)
        @constraint(m, c[25,i] + c[29,i] <= 1)
        @constraint(m, c[3,i] + c[12,i] <= 1)
        @constraint(m, c[12,i] + c[9,i] <= 1)
        @constraint(m, c[12,i] + c[14,i] <= 1)
        @constraint(m, c[9,i] + c[14,i] <= 1)
        @constraint(m, c[9,i] + c[16,i] <= 1)
        @constraint(m, c[14,i] + c[16,i] <= 1)
        @constraint(m, c[14,i] + c[19,i] <= 1)
        @constraint(m, c[14,i] + c[18,i] <= 1)
        @constraint(m, c[18,i] + c[19,i] <= 1)
        @constraint(m, c[16,i] + c[18,i] <= 1)
        @constraint(m, c[16,i] + c[17,i] <= 1)
        @constraint(m, c[16,i] + c[20,i] <= 1)
        @constraint(m, c[17,i] + c[20,i] <= 1)
        @constraint(m, c[17,i] + c[28,i] <= 1)
        @constraint(m, c[17,i] + c[30,i] <= 1)
        @constraint(m, c[17,i] + c[29,i] <= 1)
        @constraint(m, c[20,i] + c[28,i] <= 1)
        @constraint(m, c[30,i] + c[27,i] <= 1)
        @constraint(m, c[30,i] + c[29,i] <= 1)
        @constraint(m, c[30,i] + c[34,i] <= 1)
        @constraint(m, c[27,i] + c[34,i] <= 1)
        @constraint(m, c[29,i] + c[34,i] <= 1)
        @constraint(m, c[34,i] + c[38,i] <= 1)
        @constraint(m, c[38,i] + c[45,i] <= 1)
        @constraint(m, c[45,i] + c[44,i] <= 1)
        @constraint(m, c[44,i] + c[48,i] <= 1)
    end

    optimize!(m)
    if !benchmark
        println(JuMP.objective_value(m))
        println(JuMP.termination_status(m))
        _, color_washington = findmax(JuMP.value.(c[1,:]))
        println("color_washington: ", color_washington)
    end
end


text:
using GR
x = 0:0.1:100
y = sin.(x)
xlim([0, 120])
ylim([-1.5, 1.5])
xlabel("time")
ylabel("sin-value")
title("sin-plot")
plot(x, y)

max_stars_repo_name: gassive/texJuliaCopyMatlab

max_issues_repo_name: gassive/texJuliaCopyMatlab

max_forks_repo_name: gassive/texJuliaCopyMatlab


text:
function PiecewiseBicubicSplineInterpolatingFunction(arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return PiecewiseBicubicSplineInterpolatingFunction((Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end

function is_valid_point(obj::PiecewiseBicubicSplineInterpolatingFunction, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "isValidPoint", jboolean, (jdouble, jdouble), arg0, arg1)
end

function value(obj::PiecewiseBicubicSplineInterpolatingFunction, arg0::RealFieldElement, arg1::RealFieldElement)
    return jcall(obj, "value", RealFieldElement, (RealFieldElement, RealFieldElement), arg0, arg1)
end

function value(obj::PiecewiseBicubicSplineInterpolatingFunction, arg0::jdouble, arg1::jdouble)
    return jcall(obj, "value", jdouble, (jdouble, jdouble), arg0, arg1)
end


text:
abstract type ModelStatistics end

"""
    BetaStatistics(numsuccesses, numtrials)

Sample statistics of data needed to update a Beta distribution.
"""
struct BetaStatistics <: ModelStatistics
    s::Int
    n::Int
    function BetaStatistics(; s::Int, n::Int)
        # TODO: validate input data
        s <= n ||
            error("Number of trials should be equal or larger than number of successes.")
        return new(s, n)
    end

    function BetaStatistics(data::Vector{T}) where {T<:Real}
        return new(sum(data), length(data))
    end
end

"""
    GammaStatistics <: ModelStatistics

Sample statistics of data needed to update a Gamma distribution.
"""
struct GammaStatistics <: ModelStatistics
    n::Int
    x̄::Real

    function GammaStatistics(; n, x̄)
        # TODO: validate input data
        return new(n, x̄)
    end

    function GammaStatistics(data::Vector{T}) where {T<:Real}
        return new(length(data), mean(data))
    end
end

"""
    NormalStatistics <: ModelStatistics

Sample statistics of data related to Normal distribution.

## Methods

    update!(stats_old, stats_new)

Batch update for Normal statistics.

    merge(twostats)

Convert statistics of two independent samples into one NormalStatistics.

- Mean is the difference between group 1 and group 2, minus the difference in
  the null hypothesis.
- The standard deviation is calculated with pooled standard deviation.
- The sample size is the effective sample size.

    effect(stats)

Calculate effect size for one NormalStatistics or TwoNormalStatistics.

## References

- [Batch updates for simple statistics](
https://notmatthancock.github.io/2017/03/23/simple-batch-stat-updates.html)
"""
struct NormalStatistics <: ModelStatistics
    n::Real # effective sample size can be non-integer
    meanx::Real
    sdx::Real

    function NormalStatistics(; n, meanx, sdx)
        return new(n, meanx, sdx)
    end

    function NormalStatistics(data::Vector{T}) where {T<:Real}
        n = length(data)
        meanx = mean(data)
        sdx = std(data)
        return new(n, meanx, sdx)
    end
end

TwoNormalStatistics = SVector{2, NormalStatistics}

function update!(stats_old::T, stats_new::T) where {T<:NormalStatistics}
    n1 = stats_old.n
    m1 = stats_old.meanx
    s1 = stats_old.sdx
    n2 = stats_new.n
    m2 = stats_new.meanx
    s2 = stats_new.sdx

    new_n = n1 + n2
    n1_ratio = n1 / (n1 + n2)
    n2_ratio = n2 / (n1 + n2)
    new_meanx = n1_ratio * m1 + n2_ratio * m2
    new_sdx = sqrt(n1_ratio * s1^2 + n2_ratio * s2^2 + n1 * n2 / (n1 + n2)^2 * (m1 - m2)^2)

    return NormalStatistics(; n=new_n, meanx=new_meanx, sdx=new_sdx)
end

function update!(
    stats_old::TwoNormalStatistics, stats_new::TwoNormalStatistics)
    stats1 = update!(stats_old[1], stats_new[1])
    stats2 = update!(stats_new[1], stats_new[2])
    return TwoNormalStatistics(stats1, stats2)
end

function Base.merge(twostats::TwoNormalStatistics)
    stats1 = twostats[1]
    stats2 = twostats[2]
    n1 = stats1.n
    n2 = stats2.n
    eff_n = effsamplesize(n1, n2)
    sd1 = stats1.sdx
    sd2 = stats2.sdx

    Δ = (stats1.meanx - stats2.meanx)
    σ = pooledsd(sd1, sd2, n1, n2)
    return NormalStatistics(meanx=Δ, sdx=σ, n=eff_n)
end

function effectsize(stats::NormalStatistics; μ0=0.0)
    (stats.n > 0 && stats.sdx > 0) || error("invalid statistics $(stats)")
    return (stats.meanx - μ0)/stats.sdx
end

function effectsize(twosamplestats::TwoNormalStatistics; μ0=0.0)
    stats = merge(twosamplestats)
    return stats.meanx / stats.sdx
end

"""
    LogNormalStatistics <: ModelStatistics

Sample statistics of data related to Log-Normal distribution.
"""
struct LogNormalStatistics <: ModelStatistics
    n::Int
    meanlogx::Real
    sdlogx::Real

    function LogNormalStatistics(; n, meanlogx, sdlogx)
        return new(n, meanlogx, sdlogx)
    end

    function LogNormalStatistics(data::Vector{T}) where {T<:Real}
        data > zero(data) || error("data must be positive for lognormal data")
        n = length(data)
        logdata = log.(data)
        meanlogx = mean(logdata)
        sdlogx = std(logdata)
        return new(n, meanlogx, sdlogx)
    end
end

getall(stats::NormalStatistics) = (stats.n, stats.meanx, stats.sdx)
getall(stats::LogNormalStatistics) = (stats.n, stats.meanlogx, stats.sdlogx)

# tstat functions
function tstat(stats::NormalStatistics)
    return stats.meanx / stats.sdx * sqrt(stats.n)
end

function tstatpooled(twostats::TwoNormalStatistics)
    n1 = twostats[1].n
    n2 = twostats[2].n
    es = effectsize(twostats)
    ess = effsamplesize(n1, n2)
    return es*sqrt(ess)
end

function tstatwelch(twostats::TwoNormalStatistics)
    stats1 = twostats[1]
    stats2 = twostats[2]
    mu1 = stats1.meanx
    mu2 = stats2.meanx
    sd1 = stats1.sdx
    sd2 = stats2.sdx
    n1 = stats1.n
    n2 = stats2.n
    se = sqrt(sd1^2/n1 + sd2^2/n2)
    return (mu1 - mu2) / se
end

"""
    StudentTStatistics <: ModelStatistics

Sample statistics of data related to Student's t distribution.
"""
struct StudentTStatistics <: ModelStatistics
    t::Real
    dof::Real
    n::Real
end

StudentTStatistics(;t, dof, n) = StudentTStatistics(t, dof, n)

function StudentTStatistics(stats::NormalStatistics)
    t = tstat(stats)
    dof = stats.n - 1
    n = stats.n
    return StudentTStatistics(t, dof, n)
end

function StudentTStatistics(twostats::TwoNormalStatistics; pooled=true)
    n1 = twostats[1].n
    n2 = twostats[2].n
    sd1 = twostats[1].sdx
    sd2 = twostats[2].sdx

    if pooled == true
        t = tstatpooled(twostats)
        dof = dofpooled(n1, n2)
        n = effsamplesize(n1, n2)
    else
        t = tstatwelch(twostats)
        dof = dofwelch(sd1, sd2, n1, n2)
        n = effsamplesize(n1, n2)
    end

    return StudentTStatistics(t, dof, n)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiRbacV1beta1Role <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiRbacV1beta1PolicyRule} } # spec name: rules

    function IoK8sApiRbacV1beta1Role(;apiVersion=nothing, kind=nothing, metadata=nothing, rules=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1beta1Role, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiRbacV1beta1Role, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiRbacV1beta1Role, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiRbacV1beta1Role, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        o
    end
end # type IoK8sApiRbacV1beta1Role

const _property_map_IoK8sApiRbacV1beta1Role = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("rules")=>Symbol("rules"))
const _property_types_IoK8sApiRbacV1beta1Role = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("rules")=>"Vector{IoK8sApiRbacV1beta1PolicyRule}")
Base.propertynames(::Type{ IoK8sApiRbacV1beta1Role }) = collect(keys(_property_map_IoK8sApiRbacV1beta1Role))
Swagger.property_type(::Type{ IoK8sApiRbacV1beta1Role }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiRbacV1beta1Role[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1beta1Role }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1beta1Role[property_name]

function check_required(o::IoK8sApiRbacV1beta1Role)
    true
end

function validate_property(::Type{ IoK8sApiRbacV1beta1Role }, name::Symbol, val)
end


text:
# n is number of robot joints
if !(@isdefined LinkFrameCache)
    @computed struct LinkFrameCache{n,S}
        link::RigidBody{S}
        path::RigidBodyDynamics.Graphs.TreePath
        transform_jacobian::CacheElement{Vector{SMatrix{12,n,S,n*12}}}
        transform_jacobian_dot::CacheElement{Vector{SMatrix{12,n,S,n*12}}}
    end
end

function LinkFrameCache(link, state::MechanismState{X}) where X
    n = num_positions(state)
    link_path = path(state.mechanism, root_body(state.mechanism), link)
    transform_jacobian = CacheElement([@SMatrix zeros(12,n)])
    transform_jacobian_dot = CacheElement([@SMatrix zeros(12,n)])
    LinkFrameCache{n,X}(link, link_path, transform_jacobian, transform_jacobian_dot)
end

function setdirty!(link_cache::LinkFrameCache)
    CustomCollections.setdirty!(link_cache.transform_jacobian)
    CustomCollections.setdirty!(link_cache.transform_jacobian_dot)
    nothing
end

function link_frame_transform(xm, state_cache, link)
    state = state_cache[eltype(xm)]
    set_configuration!(state, xm)
    T = transform_to_root(state, link)
    R, t = SVector{9,eltype(xm)}(reshape(rotation(T),9)), SVector{3,eltype(xm)}(translation(T))
    [R; t]
end

@inline function update_link_transform_jacobian!(transform_jacobian::CacheElement, xm,
        state_cache, link)
    isdirty(transform_jacobian) && _update_link_transform_jacobian!(transform_jacobian, xm,
        state_cache, link)
    nothing
end

function _update_link_transform_jacobian!(transform_jacobian::CacheElement, xm, state_cache, link)
    f = (xm) -> link_frame_transform(xm, state_cache, link)
    transform_jacobian.data[1] = ForwardDiff.jacobian(f, xm)
    transform_jacobian.dirty = false
    nothing
end

@inline function update_link_transform_jacobian_dot!(transform_jacobian_dot::CacheElement,
        xm, vm, state_cache, link)
    isdirty(transform_jacobian_dot) && _update_link_transform_jacobian_dot!(transform_jacobian_dot,
        xm, vm, state_cache, link)
    nothing
end

function _update_link_transform_jacobian_dot!(transform_jacobian_dot::CacheElement,
        xm, vm, state_cache, link)
    n, m = size(transform_jacobian_dot.data[1])
    f = (xm) -> link_frame_transform(xm, state_cache, link)
    Jf = (xm) -> ForwardDiff.jacobian(f, xm)
    ∇xf = ForwardDiff.jacobian(xm -> reshape(Jf(xm), m*n), xm)
    Jfdot = SMatrix{n,m,eltype(xm)}(reshape(∇xf*vm, n, m))::SArray{Tuple{n,m},eltype(xm),2,m*n}

    transform_jacobian_dot.data[1] = Jfdot
    transform_jacobian_dot.dirty = false
    nothing
end


text:
#=
--------------------------------------------------------------------------------------------------------------------------------
Types and constructors
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    DFMSettings(...)

Define an immutable structure used to initialise the estimation routine for a DFM.

# Arguments
- `Y`: observed measurements (`nxT`)
- `n`: Number of series
- `T`: Number of observations
- `lags`: Order of the autoregressive polynomial of each common cycle
- `n_trends`: Number of trends
- `n_drifts`: Number of drifts
- `n_cycles`: Number of cycles
- `n_non_stationary`: n_trends + n_drifts
- `m`: n_non_stationary + n_cycles*lags + n
- `trends_skeleton`: The basic structure for the trends loadings (or nothing)
- `cycles_skeleton`: The basic structure for the cycles loadings (or nothing)
- `drifts_selection`: BitArray{1} identifying which trend has a drift to estimate (or nothing)
- `trends_free_params`: BitArray{2} identifying the trend loadings to estimate (or nothing)
- `cycles_free_params`: BitArray{2} identifying the cycle loadings to estimate (or nothing)
- `λ`: overall shrinkage hyper-parameter for the elastic-net penalty
- `α`: weight associated to the LASSO component of the elastic-net penalty
- `β`: additional shrinkage for distant lags (p>1)
- `Γ`: Diagonal matrix used to input the hyperparameters in the penalty computation for the common cycles
- `Γ_idio`: Diagonal matrix used to input the hyperparameters in the penalty computation for the idiosyncratic cycles
- `Γ_extended`: Diagonal matrix used to input the hyperparameters in the estimation
- `ε`: Small number (default: 1e-4)
- `tol`: tolerance used to check convergence (default: 1e-4)
- `max_iter`: maximum number of iterations for the estimation algorithm (default: 1000)
- `prerun`: number of iterations prior the actual estimation algorithm (default: 2)
- `check_quantile`: check the quantile of the relative change of the parameters for convergence purposes (default: false)
- `verb`: Verbose output (default: true)
"""
struct DFMSettings <: EstimSettings
    Y::Union{FloatMatrix, JMatrix{Float64}}
    n::Int64
    T::Int64
    lags::Int64
    n_trends::Int64
    n_drifts::Int64
    n_cycles::Int64
    n_non_stationary::Int64
    m::Int64
    trends_skeleton::Union{FloatMatrix, Nothing}
    cycles_skeleton::Union{FloatMatrix, Nothing}
    drifts_selection::Union{BitVector, Nothing}
    trends_free_params::Union{BitMatrix, Nothing}
    cycles_free_params::Union{BitMatrix, Nothing}
    λ::Float64
    α::Float64
    β::Float64
    Γ::DiagMatrix
    Γ_idio::Float64
    Γ_extended::DiagMatrix
    ε::Float64
    tol::Float64
    max_iter::Int64
    prerun::Int64
    check_quantile::Bool
    verb::Bool
end

#=
DFMSettings constructor (no input provided for basic skeletons and free parameters)
- This constructor returns the default implementation for a stationary DFM
=#
function DFMSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, lags::Int64, n_cycles::Int64, λ::Float64, α::Float64, β::Float64; ε::Float64=1e-4, tol::Float64=1e-4, max_iter::Int64=1000, prerun::Int64=2, check_quantile::Bool=false, verb::Bool=true)

    # Make sure there is at least one common cycle to estimate
    check_bounds(n_cycles, 0);

    # Compute missing dimensions
    n, T = size(Y);
    n_trends = 0;
    n_drifts = 0;
    n_non_stationary = 0;
    m = n_non_stationary + n_cycles*lags + n;

    # Compute missing input for estimation
    cycles_skeleton = hcat([[zeros(i-1); 1; 2*ones(n-i)] for i=1:n_cycles]...);
    cycles_free_params = cycles_skeleton .> 1;

    # Compute penalty data
    Γ = build_Γ(1, lags, λ, β);
    Γ_idio = build_Γ(1, 1, λ, β)[1];
    Γ_extended = cat(dims=[1,2], [Γ_idio for i=1:n]..., [Γ for i=1:n_cycles]...) |> Array |> Diagonal;

    # Return DFMSettings
    return DFMSettings(Y, n, T, lags, n_trends, n_drifts, n_cycles, n_non_stationary, m, nothing, cycles_skeleton, nothing, nothing, cycles_free_params, λ, α, β, Γ, Γ_idio, Γ_extended, ε, tol, max_iter, prerun, check_quantile, verb);
end

#=
DFMSettings constructor (input provided only for the cycles skeleton and free parameters)
- This constructor returns the default custom implementation for a stationary DFM
=#
function DFMSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, lags::Int64, cycles_skeleton::FloatMatrix, cycles_free_params::BitMatrix, λ::Float64, α::Float64, β::Float64; ε::Float64=1e-4, tol::Float64=1e-4, max_iter::Int64=1000, prerun::Int64=2, check_quantile::Bool=false, verb::Bool=true)

    # Compute missing dimensions
    n, T = size(Y);
    n_trends = 0;
    n_drifts = 0;
    n_cycles = size(cycles_skeleton, 2);
    n_non_stationary = 0;
    m = n_non_stationary + n_cycles*lags + n;

    # Compute penalty data
    Γ = build_Γ(1, lags, λ, β);
    Γ_idio = build_Γ(1, 1, λ, β)[1];
    Γ_extended = cat(dims=[1,2], [Γ_idio for i=1:n]..., [Γ for i=1:n_cycles]...) |> Array |> Diagonal;

    # Return DFMSettings
    return DFMSettings(Y, n, T, lags, n_trends, n_drifts, n_cycles, n_non_stationary, m, nothing, cycles_skeleton, nothing, nothing, cycles_free_params, λ, α, β, Γ, Γ_idio, Γ_extended, ε, tol, max_iter, prerun, check_quantile, verb);
end

#=
DFMSettings constructor (input provided for the cycles and trends skeleton / free parameters)
- This constructor returns the default custom implementation for a non-stationary DFM
=#
function DFMSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, lags::Int64, trends_skeleton::FloatMatrix, cycles_skeleton::FloatMatrix, drifts_selection::BitVector, trends_free_params::BitMatrix, cycles_free_params::BitMatrix, λ::Float64, α::Float64, β::Float64; ε::Float64=1e-4, tol::Float64=1e-4, max_iter::Int64=1000, prerun::Int64=2, check_quantile::Bool=false, verb::Bool=true)

    # Compute missing dimensions
    n, T = size(Y);
    n_trends = size(trends_skeleton, 2);
    n_drifts = sum(drifts_selection);
    n_cycles = size(cycles_skeleton, 2);
    n_non_stationary = n_trends + n_drifts;
    m = n_non_stationary + n_cycles*lags + n;

    # Compute penalty data
    Γ = build_Γ(1, lags, λ, β);
    Γ_idio = build_Γ(1, 1, λ, β)[1];
    Γ_extended = cat(dims=[1,2], Diagonal(zeros(n_non_stationary, n_non_stationary)), [Γ_idio for i=1:n]..., [Γ for i=1:n_cycles]...) |> Array |> Diagonal;

    # Return DFMSettings
    return DFMSettings(Y, n, T, lags, n_trends, n_drifts, n_cycles, n_non_stationary, m, trends_skeleton, cycles_skeleton, drifts_selection, trends_free_params, cycles_free_params, λ, α, β, Γ, Γ_idio, Γ_extended, ε, tol, max_iter, prerun, check_quantile, verb);
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Initialisation
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    check_model_bounds(estim::DFMSettings)

Check whether the arguments for the estimation of a DFM model are correct.
"""
function check_model_bounds(estim::DFMSettings)

    # Simple bounds
    check_bounds(estim.lags, 0);
    check_bounds(estim.n_cycles, 0);
    check_bounds(estim.λ, 0);
    check_bounds(estim.α, 0, 1);
    check_bounds(estim.β, 1);
    check_bounds(estim.max_iter, 3);
    check_bounds(estim.max_iter, estim.prerun);
    check_bounds(estim.n, 2);

    # Advanced checks

    if eltype(estim.trends_skeleton) <: FloatMatrix
        if size(estim.trends_skeleton, 2) != size(estim.trends_free_params, 2)
            throw(DomainError);
        end
    end

    if eltype(estim.cycles_skeleton) <: FloatMatrix
        if size(estim.cycles_skeleton, 2) != size(estim.cycles_free_params, 2)
            throw(DomainError);
        end
    end

    if eltype(estim.drifts_selection) <: BitVector
        if size(estim.trends_skeleton, 2) != length(estim.drifts_selection)
            throw(DomainError);
        end
    end

    if maximum(diff(sum(estim.cycles_skeleton .!= 0, dims=2), dims=1)) .> 1
        throw(error("The structure of `cycles_skeleton` is incorrect. Cycles should be added incrementally."));
    end

    if maximum(diff(sum(estim.cycles_free_params .!= 0, dims=2), dims=1)) .> 1
        throw(error("The structure of `cycles_free_params` breaks `cycles_skeleton`."));
    end
end

"""
    initialise_trends(estim::DFMSettings)

Return initialised parameters for the trends.
"""
function initialise_trends(estim::DFMSettings, common_trends::Union{FloatMatrix, JMatrix{Float64}})

    # TBD: Check with one trend only!!!!!

    # Estimate drift where necessary
    drift_trends = mean_skipmissing(diff(common_trends, dims=2))[:];
    variance_trends = std_skipmissing(diff(common_trends, dims=2))[:].^2;

    #=
    Memory pre-allocation for the state-space parameters linked to the non-stationary components
    Note: similarly to Watson and Engle (1983), estim.n_non_stationary lags of the main non-stationary components are added to compute PPs in the ECM algorithm
    =#

    B_trends = zeros(estim.n, 2*estim.n_non_stationary);
    C_trends = zeros(2*estim.n_non_stationary, 2*estim.n_non_stationary);
    D_trends = zeros(2*estim.n_non_stationary, estim.n_non_stationary);
    Q_trends = zeros(estim.n_non_stationary, estim.n_non_stationary);
    X0_trends = zeros(2*estim.n_non_stationary);
    P0_trends = zeros(2*estim.n_non_stationary, 2*estim.n_non_stationary);

    # Initialise counter
    i = 0;

    # Loop over the trends
    for j=1:estim.n_trends

        B_trends[:, 1+i] = estim.trends_skeleton[:,j];
        n_previous_drifts = ifelse(j>1, sum(estim.drifts_selection[1:j-1]), 0)::Int64;

        # Trends with drift
        if estim.drifts_selection[j]
            C_trends[1+i:4+i, 1+i:4+i] = [1 0 1 0; 1 0 0 0; 0 0 1 0; 0 0 1 0];
            D_trends[1+i, j+n_previous_drifts] = 1;
            D_trends[3+i, j+1+n_previous_drifts] = 1;
            Q_trends[j+n_previous_drifts, j+n_previous_drifts] = variance_trends[j];
            Q_trends[j+1+n_previous_drifts, j+1+n_previous_drifts] = estim.ε;
            X0_trends[3+i] = drift_trends[j];
            P0_trends[1+i:4+i, 1+i:4+i] = (1/estim.ε)*Matrix(I,4,4);
            i += 4;

        # Driftless trends
        else
            C_trends[1+i:2+i, 1+i:2+i] = [1 0; 1 0];
            D_trends[1+i, j+n_previous_drifts] = 1;
            Q_trends[j+n_previous_drifts, j+n_previous_drifts] = variance_trends[j];
            P0_trends[1+i:2+i, 1+i:2+i] = (1/estim.ε)*Matrix(I,2,2);
            i += 2;
        end
    end

    # Return output
    return B_trends, C_trends, D_trends, Q_trends, X0_trends, P0_trends;
end

"""
    initialise_common_cycle(estim::DFMSettings, residual_data::FloatMatrix, coordinates_current_block::IntVector)

Initialise current common cycle via PCA.
"""
function initialise_common_cycle(estim::DFMSettings, residual_data::FloatMatrix, coordinates_current_block::IntVector)

    # Convenient shortcuts
    data_current_block = residual_data[coordinates_current_block, :];
    data_current_block_standardised = standardise(data_current_block);

    # Compute PCA loadings
    eigen_val, eigen_vect = eigen(Symmetric(cov(data_current_block_standardised, dims=2)));
    loadings = eigen_vect[:, sortperm(-abs.(eigen_val))[1]];

    # Compute PCA factor
    pc1 = permutedims(loadings)*data_current_block_standardised;

    # Rescale PCA loadings to match the original scale
    loadings .*= std(data_current_block, dims=2)[:];

    # Rescale PC1 wrt the first series
    pc1 .*= loadings[1];
    loadings /= loadings[1];

    # Estimate ridge autoregressive coefficients
    pc1_y, pc1_x = lag(pc1, estim.lags);
    ar_coeff = pc1_y*pc1_x'/Symmetric(pc1_x*pc1_x' + estim.Γ);
    enforce_causality_and_invertibility!(ar_coeff);

    # Estimate var-cov matrix of the residuals
    ar_residuals = pc1_y - ar_coeff*pc1_x;
    ar_residuals_variance = (ar_residuals*ar_residuals')/length(ar_residuals);

    if estim.lags > 1

        # Initialise ridge loadings associated to the factor and its lags
        complete_loadings = zeros(length(loadings), estim.lags);
        complete_loadings[1,1] = 1.0;

        # Regular VAR configuration
        B = [1.0 zeros(1, estim.lags-1)];
        R = Symmetric(estim.ε * ones(1,1));
        C = companion_form(ar_coeff, extended=false);
        D = [1.0; zeros(estim.lags-1, 1)];
        Q = Symmetric(ar_residuals_variance);

        # Add estim.lags missings before pc1
        sspace = KalmanSettings(hcat(missing .* ones(1, estim.lags-1), pc1), B, R, C, D, Q, compute_loglik=false);
        status = kfilter_full_sample(sspace);
        smoothed_states, _ = ksmoother(sspace, status);

        # Backcast factor up to t=-estim.lags
        pc1_extended = mapreduce(Xt -> sspace.B*Xt, hcat, smoothed_states);
        pc1_extended_x = vcat([pc1_extended[:, estim.lags-j:end-j] for j=0:estim.lags-1]...);

        # Compute ridge loadings associated to the factor and its lags
        complete_loadings[2:end,:] = @views data_current_block[2:end,:]*pc1_extended_x'/Symmetric(pc1_extended_x*pc1_extended_x' + estim.Γ);

        # Explained data
        explained_data = complete_loadings*pc1_extended_x;

    else
        complete_loadings = loadings;
        explained_data = complete_loadings*pc1;
    end

    # Return output
    return complete_loadings, explained_data, ar_coeff, ar_residuals_variance;
end

"""
    initialise_cycles(estim::DFMSettings, data::FloatMatrix)

Return initialised parameters for the cycles.
"""
function initialise_cycles(estim::DFMSettings, data::FloatMatrix)

    #=
    Memory pre-allocation for the state-space parameters referring to the cycles (both common and idiosyncratic)
    Note: similarly to Watson and Engle (1983), n_cycles+n lags of the main stationary components are added to compute PPs in the ECM algorithm
    =#

    B_cycles = zeros(estim.n, estim.n_cycles*estim.lags + estim.n_cycles + 2*estim.n);
    C_cycles = zeros(estim.n_cycles*estim.lags + estim.n_cycles + 2*estim.n, estim.n_cycles*estim.lags + estim.n_cycles + 2*estim.n);
    D_cycles = zeros(estim.n_cycles*estim.lags + estim.n_cycles + 2*estim.n, estim.n_cycles + estim.n);
    Q_cycles = zeros(estim.n_cycles+estim.n, estim.n_cycles+estim.n);
    X0_cycles = zeros(size(B_cycles,2));
    P0_cycles = zeros(estim.n_cycles*estim.lags + estim.n_cycles + 2*estim.n, estim.n_cycles*estim.lags + estim.n_cycles + 2*estim.n);

    # Set-up loadings for idiosyncratic components
    B_cycles[:, 1:2*estim.n] = cat(dims=[1,2], [[1 0] for i=1:estim.n]...);

    # Determine which variables can load on each cycle
    boolean_coordinates_blocks = (estim.cycles_skeleton .!= 0) .| estim.cycles_free_params;
    coordinates_blocks = [findall(boolean_coordinates_blocks[:,i]) for i=1:estim.n_cycles];

    # Initialise common cycles iteratively via PCA
    residual_data = copy(data);

    for i=1:estim.n_cycles

        # Pointers
        coordinates_current_block = coordinates_blocks[i];

        # Initialise i-th common cycle
        complete_loadings, explained_data, ar_coeff, ar_residuals_variance = initialise_common_cycle(estim, residual_data, coordinates_current_block);

        #=
        Update state-space representation
        Note: the idiosyncratic components are placed before the common ones, as in the reference paper
        =#

        # Position of the i-th common cycle in the state-space representation
        coordinates_pc1 = 1+(i-1)*(estim.lags+1)+2*estim.n;

        # Convenient short cut for `C_cycles` and `P0_cycles`
        current_companion = companion_form(ar_coeff, extended=true);
        current_cov_companion = Symmetric(cat(dims=[1,2], ar_residuals_variance[1], zeros(estim.lags, estim.lags)));

        # Update coefficients
        B_cycles[coordinates_current_block, coordinates_pc1:coordinates_pc1+estim.lags-1] = complete_loadings;
        C_cycles[coordinates_pc1:coordinates_pc1+estim.lags, coordinates_pc1:coordinates_pc1+estim.lags] = current_companion;
        D_cycles[coordinates_pc1, i+estim.n] = 1;
        Q_cycles[i+estim.n, i+estim.n] = ar_residuals_variance[1];

        # Update `P0_cycles` in a block-wise fashion
        P0_cycles[coordinates_pc1:coordinates_pc1+estim.lags, coordinates_pc1:coordinates_pc1+estim.lags] = Array(solve_discrete_lyapunov(current_companion, current_cov_companion));

        # Recompute residual data
        residual_data[coordinates_current_block, :] .-= explained_data;
    end

    # Initialise idiosyncratic cycles as white noises
    for i=1:estim.n

        # Estimate var-cov matrix of the residuals
        idio_resid = permutedims(diff(residual_data[i,:]));
        idio_resid_var = (idio_resid*idio_resid')/length(idio_resid);

        #=
        Update state-space representation
        Note: the idiosyncratic components are placed before the common ones, as in the reference paper
        =#

        # Position of the i-th idiosyncratic component
        coordinates_idio = 1+(i-1)*2; # place the idiosyncratic components before the common ones as in the reference paper

        # Convenient short cut for `C_cycles` and `P0_cycles`
        current_companion = companion_form(zeros(1,1), extended=true);
        current_cov_companion = Symmetric(cat(dims=[1,2], idio_resid_var[1], zeros(1,1)));

        # Update coefficients
        C_cycles[coordinates_idio:coordinates_idio+1, coordinates_idio:coordinates_idio+1] = current_companion;
        D_cycles[coordinates_idio, i] = 1;
        Q_cycles[i, i] = idio_resid_var[1];

        # Update `P0_cycles` in a block-wise fashion
        P0_cycles[coordinates_idio:coordinates_idio+1, coordinates_idio:coordinates_idio+1] = Array(solve_discrete_lyapunov(current_companion, current_cov_companion));
    end

    # Return output
    return B_cycles, C_cycles, D_cycles, Q_cycles, X0_cycles, P0_cycles;
end

"""
    initialise(estim::DFMSettings, trends_skeleton::Nothing)

Initialise the ECM algorithm for stationary DFM models.

    initialise(estim::DFMSettings, trends_skeleton::FloatMatrix)

Initialise the ECM algorithm for non-stationary DFM models.

    initialise(estim::DFMSettings)

Initialise the ECM algorithm for DFM models.

# Arguments
- `estim`: settings used for the estimation

# References
Pellegrino (2022)
"""
function initialise(estim::DFMSettings, trends_skeleton::Nothing)

    # Interpolate data
    data = interpolate_series(estim.Y, estim.n, estim.T);

    # Build state-space parameters
    B, C, D, Q, X0, P0 = initialise_cycles(estim, data);
    Q = Symmetric(Q);
    R = Symmetric(Matrix(estim.ε * I, estim.n, estim.n));

    # Generate sspace
    sspace = KalmanSettings(estim.Y, B, R, C, D, Q, X0, P0, compute_loglik=false);

    # TBA coordinates etc
    error("This is not finalised yet!");

    # Return output
    return sspace;
end

function initialise(estim::DFMSettings, trends_skeleton::FloatMatrix)

    # Compute individual trends
    sqrt_T = floor(sqrt(estim.T)) |> Int64;
    trends = centred_moving_average(forward_backwards_rw_interpolation(estim.Y, estim.n, estim.T), estim.n, estim.T, 2*sqrt_T+1);

    # Removing trailing and leading missings
    first_ind = findfirst(sum(ismissing.(trends), dims=1) .< estim.n)[2];
    last_ind = findlast(sum(ismissing.(trends), dims=1) .< estim.n)[2];
    trends = trends[:, first_ind:last_ind] |> FloatMatrix;

    # Compute common trends. `common_trends` is equivalent to `trends` if there aren't common trends to compute.
    common_trends = ones(estim.n_trends, estim.T-2*sqrt_T);
    for i=1:estim.n_trends
        coordinates_current_block = findall(estim.trends_skeleton[:, i] .!= 0);
        common_trends[i,:] = median(trends[coordinates_current_block, :] ./ estim.trends_skeleton[coordinates_current_block, i], dims=1);
    end

    # Interpolate detrended data
    detrended_data = interpolate_series(estim.Y[:, first_ind:last_ind] - estim.trends_skeleton*common_trends, estim.n, estim.T-2*sqrt_T);

    # Build state-space parameters
    B_trends, C_trends, D_trends, Q_trends, X0_trends, P0_trends = initialise_trends(estim, common_trends);
    B_cycles, C_cycles, D_cycles, Q_cycles, X0_cycles, P0_cycles = initialise_cycles(estim, detrended_data);
    B = hcat(B_trends, B_cycles);
    R = Symmetric(Matrix(estim.ε * I, estim.n, estim.n));
    C = cat(dims=[1,2], C_trends, C_cycles);
    D = cat(dims=[1,2], D_trends, D_cycles);
    Q = Symmetric(cat(dims=[1,2], Q_trends, Q_cycles));

    # Initial conditions
    X0 = vcat(X0_trends, X0_cycles);
    P0 = Symmetric(cat(dims=[1,2], P0_trends, P0_cycles));

    # Generate sspace
    sspace = KalmanSettings(estim.Y, B, R, C, D, Q, X0, P0, compute_loglik=false);

    # `coordinates_transition_current` identifies the states for which there is an associated variance that is allowed to differ from zero.
    coordinates_transition_current = findall(sum(D, dims=2)[:] .== 1);

    # Coordinates per type
    coordinates_transition_non_stationary = coordinates_transition_current[coordinates_transition_current .<= 2*estim.n_non_stationary];
    coordinates_transition_stationary = coordinates_transition_current[coordinates_transition_current .> 2*estim.n_non_stationary];
    coordinates_transition_idio_cycles = coordinates_transition_stationary[1:end-estim.n_cycles];
    coordinates_transition_common_cycles = coordinates_transition_stationary[end-estim.n_cycles+1:end];

    # coordinates_transition_lagged and *_PPs
    coordinates_transition_lagged = sort(vcat(coordinates_transition_non_stationary, coordinates_transition_idio_cycles, [coordinates_transition_common_cycles .+ i for i=0:estim.lags-1]...));
    coordinates_transition_PPs = coordinates_transition_lagged .+ 1;

    # `coordinates_transition_P0` identifies the entry in P0 that the cm-step should recompute
    coordinates_transition_P0 = findall(P0[:] .!= 0.0);

    # Comment this out if you want to keep the original diffuse initialisation for the non-stationary components
    # filter!(coordinate->coordinate.I[1] > size(P0_trends,1), coordinates_transition_P0);

    #=
    `coordinates_measurement_states` identifies the states for which the associated loadings are allowed to differ from zero.
    Note: the most compact version for this representation would be setdiff(coordinates_transition_lagged, coordinates_transition_drifts).
          However, this implementation of the DFM also includes the columns for the drifts, which always have zero loadings.
          The main advantage of this choice is that estim.Γ_extended can be also used for the loadings cm step.
    =#

    coordinates_measurement_states = copy(coordinates_transition_lagged);

    # Convenient views for using sspace.B in the expected logliklihood and cm steps calculations
    B_star = @view sspace.B[:, coordinates_measurement_states];

    # Coordinates free parameters (B)
    cartesian_B = CartesianIndices(B_star);
    free_params_B_trends = hcat([ifelse(estim.drifts_selection[i], hcat(estim.trends_free_params[:,i], ones(estim.n).==0), estim.trends_free_params[:,i]) for i=1:estim.n_trends]...);
    free_params_B_idio_cycles = zeros(estim.n, estim.n) .== 1;
    free_params_B_common_cycles = hcat(permutedims([estim.cycles_free_params[:,i] for i=1:estim.n_cycles, j=1:estim.lags])...);
    coordinates_free_params_B = cartesian_B[hcat(free_params_B_trends, free_params_B_idio_cycles, free_params_B_common_cycles)];

    # Convenient views for using sspace.C in the expected logliklihood and cm steps calculations
    C_star = @view sspace.C[coordinates_transition_current, coordinates_transition_lagged];

    # Coordinates free parameters (C)
    cartesian_C = CartesianIndices(C_star);
    free_params_C_trends = zeros(estim.n_non_stationary, estim.n_non_stationary) .== 1;
    free_params_C_idio_cycles = Matrix(I, estim.n, estim.n);
    free_params_C_common_cycles = cat(dims=[1,2], [ones(1, estim.lags) for i=1:estim.n_cycles]...) .== 1;
    coordinates_free_params_C = cartesian_C[cat(dims=[1,2], free_params_C_trends, free_params_C_idio_cycles, free_params_C_common_cycles)];

    # View on Q from DQD
    Q_view = @view sspace.DQD.data[coordinates_transition_current, coordinates_transition_current];

    # Return output
    return sspace, B_star, C_star, Q_view, coordinates_measurement_states, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0, coordinates_free_params_B, coordinates_free_params_C;
end

function initialise(estim::DFMSettings)
    return initialise(estim, estim.trends_skeleton);
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Convergence check
--------------------------------------------------------------------------------------------------------------------------------
=#

function vec_ecm_params(estim::DFMSettings, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

    # Setup and memory pre-allocation
    counter = 1;
    output = zeros(length(coordinates_free_params_B) + length(coordinates_free_params_C) + size(Q_view,1));

    # Add parameters in B_star to output
    for ij in coordinates_free_params_B
        output[counter] = B_star[ij];
        counter += 1;
    end

    # Add parameters in C_star to output
    for ij in coordinates_free_params_C
        output[counter] = C_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for i in axes(Q_view,1)
        output[counter] = Q_view[i,i];
        counter += 1;
    end

    # Return output
    return output;
end

function vec_ecm_params!(estim::DFMSettings, output::FloatVector, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

    # Setup and memory pre-allocation
    counter = 1;

    # Add parameters in B_star to output
    for ij in coordinates_free_params_B
        output[counter] = B_star[ij];
        counter += 1;
    end

    # Add parameters in C_star to output
    for ij in coordinates_free_params_C
        output[counter] = C_star[ij];
        counter += 1;
    end

    # Add parameters in Q_view to output
    for i in axes(Q_view,1)
        output[counter] = Q_view[i,i];
        counter += 1;
    end
end

#=
--------------------------------------------------------------------------------------------------------------------------------
CM-step
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    cm_step!(estim::DFMSettings, sspace::KalmanSettings, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, F_raw::FloatMatrix, G::FloatMatrix, H_raw::FloatMatrix, M::FloatArray, N::Vector{SparseMatrixCSC{Float64, Int64}}, O::Array{VectorsArray{Float64},1}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

Run the CM-step for elastic-net DFMs.
"""
function cm_step!(estim::DFMSettings, sspace::KalmanSettings, B_star::SubArray{Float64}, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, F_raw::FloatMatrix, G::FloatMatrix, H_raw::FloatMatrix, M::FloatArray, N::Vector{SparseMatrixCSC{Float64, Int64}}, O::Array{VectorsArray{Float64},1}, coordinates_free_params_B::CoordinatesVector, coordinates_free_params_C::CoordinatesVector)

    # Make sure `F_raw` and `H_raw` are symmetric
    F = Symmetric(F_raw)::SymMatrix;
    H = Symmetric(H_raw)::SymMatrix;

    # Loop over the free parameters in B
    for ij in coordinates_free_params_B

        # Compute remaining sufficient statistics
        sum_numerator, sum_denominator, i, j = cm_step_time_loop(sspace, B_star, ij, N, O);

        # Compute new value for B_star[ij]
        B_star[ij] = soft_thresholding(M[ij] - sum_numerator, 0.5*estim.α*estim.ε*estim.Γ_extended[j, j]);
        B_star[ij] /= sum_denominator + (1-estim.α)*estim.ε*estim.Γ_extended[j, j];
    end

    # Copy the sspace.C components referring to the common cycles
    C_cycles_old = sspace.C[2*(estim.n_non_stationary+estim.n)+1:end, 2*(estim.n_non_stationary+estim.n)+1:end];
    C_cycles_new = @view sspace.C[2*(estim.n_non_stationary+estim.n)+1:end, 2*(estim.n_non_stationary+estim.n)+1:end];

    # Compute inverse of Q
    inv_Q = inv(Symmetric(Q_view));

    # Counter for common cycles and base coordinates
    n_previous_common_cycles = 0;
    base_coordinates_cycle = collect(1:estim.lags);

    # Loop over the free parameters
    for (position, ij) in enumerate(coordinates_free_params_C)

        # Coordinates
        i,j = ij.I;

        # Compute new value for sspace.C[ij]
        C_star[ij] = soft_thresholding(turbo_dot(inv_Q[:,i], G[:,j] - C_star*H[:,j]) + inv_Q[i,i]*C_star[ij]*H[j,j], 0.5*estim.α*estim.Γ_extended[j, j]);
        C_star[ij] /= inv_Q[i,i]*H[j,j] + (1-estim.α)*estim.Γ_extended[j, j];

        # Adjust autoregressive coefficient of the idiosyncratic cycles to enforce causality
        if position <= estim.n
            if abs(C_star[ij]) > 0.98
                C_star[ij] = sign(C_star[ij])*0.98;
            end

        # Adjust autoregressive coefficients of the common cycles to enforce causality
        elseif mod(position-estim.n, estim.lags) == 0

            # Position of the current cycle in `C_cycles_new` and `C_cycles_old`
            coordinates_current_cycle = base_coordinates_cycle .+ n_previous_common_cycles*(estim.lags+1);

            # Convenient views
            C_current_cycle_new = @view C_cycles_new[coordinates_current_cycle, coordinates_current_cycle];
            C_current_cycle_old = @view C_cycles_old[coordinates_current_cycle, coordinates_current_cycle];

            # Adjust coefficients (if needed)
            enforce_causality_and_invertibility!(C_current_cycle_new, C_current_cycle_old);

            # Update counter
            n_previous_common_cycles += 1;
        end
    end

    # CM-step for Q
    Q_cm_step = (F-G*C_star'-C_star*G'+C_star*H*C_star')/estim.T; # this matrix product could be improved, since it also considers the off-diagonal elements (which are not used in this cm-step)
    for i=1:estim.n_non_stationary + estim.n + estim.n_cycles
        Q_view[i,i] = Q_cm_step[i,i];
    end
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Model validation
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    rescale_estim_params!(coordinates_params_rescaling::VectorsArray{Int64}, estim::DFMSettings, std_presample::FloatMatrix)

Rescale loadings associated to the trends.
"""
function rescale_estim_params!(coordinates_params_rescaling::VectorsArray{Int64}, estim::DFMSettings, std_presample::FloatMatrix)
    for (coordinate_std, coordinate_skeleton) in coordinates_params_rescaling
        estim.trends_skeleton[coordinate_skeleton] = 1/std_presample[coordinate_std];
    end
end

"""
    DFMSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::Nothing, p::Int64, λ::Float64, α::Float64, β::Float64)
    DFMSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::NamedTuple, p::Int64, λ::Float64, α::Float64, β::Float64)

Generate the relevant DFMSettings structure by using `model_inputs` and the candidate hyperparameters `p`, `λ`, `α`, `β`.
"""
DFMSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::Nothing, p::Int64, λ::Float64, α::Float64, β::Float64) = DFMSettings(Y, p, model_args..., λ, α, β, verb=false);
DFMSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::NamedTuple, p::Int64, λ::Float64, α::Float64, β::Float64) = DFMSettings(Y, p, model_args..., λ, α, β; model_kwargs...);


text:
# test to start debugging adding more than one hybrid
# claudia may 2015

if !isdefined(:individualtest) individualtest = false; end

if(individualtest)
    include("../src/types.jl")
    include("../src/functions.jl")
end

seed=1234
tree = "(((((((1,2),3),4),5),(6,7)),(8,9)),10);"
currT0 = readTopologyLevel1(tree);
#printEdges(currT0)
besttree = deepcopy(currT0);
srand(1234);
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
success || error("added hybrid not successful")
flag || error("added hybrid not successful")
flag2 || error("added hybrid not successful")
flag3 || error("added hybrid not successful")
!nocycle || error("added hybrid not successful")
hybrid.hybrid || error("added hybrid not successful")
hybrid.number == 11 || error("added hybrid not successful")
besttree.edge[12].inCycle == 11 || error("incycle not correct")
besttree.edge[13].inCycle == 11 || error("incycle not correct")
besttree.edge[18].inCycle == 11 || error("incycle not correct")
besttree.edge[20].inCycle == 11 || error("incycle not correct")
hybrid.k == 4 || error("not correct hybrid.k")
!hybrid.isBadDiamondI || error("thinks its bad diamond I")
!hybrid.isBadDiamondII || error("thinks its bad diamond II")
#printEdges(besttree)
#printNodes(besttree)
#writeTopology(besttree,true)
success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree); #will add a bad triangle
!success || error("added bad triangle and did not notice")
flag || error("added bad triangle and did not notice")
!flag2 || error("added bad triangle and did not notice")
flag3 || error("added bad triangle and did not notice")
!nocycle || error("added bad triangle and did not notice")
hybrid.k == 3 || error("added bad triangle and did not notice")
hybrid.isVeryBadTriangle || error("added bad triangle and did not notice")
ed = hybridEdges(hybrid)
ed[1].isMajor || error("in bad triangle, update major hybrid not working, and it should")
ed[1].gamma > 0.5 || error("in bad triangle, update major hybrid not working, and it should")
ed[1].hybrid || error("in bad triangle, update major hybrid not working, and it should")

# did not recognize as bad diamond II
tree = "(6,(5,#H7:0.0):9.970714072991349,(3,(((2,1):0.2950382234364404,4):0.036924483697671304)#H7:0.00926495670648208):1.1071489442240392);"
net = readTopologyLevel1(tree);
net.node[10].isBadDiamondII || error("does not recognize as bad diamond II")

#print("NO ERRORS!!")

# printEdges(besttree)
# printNodes(besttree)
# writeTopology(besttree,true)

#----add hyb by parts
# include("../src/types.jl")
# include("../src/functions.jl")
# seed=1234
# tree = "(((((((1,2),3),4),5),(6,7)),(8,9)),10);"
# currT0 = readTopologyUpdate(tree);
# printEdges(currT0)
# besttree = deepcopy(currT0);
# srand(1234);
# success,hybrid,flag,nocycle,flag2,flag3 = addHybridizationUpdate!(besttree);
# printEdges(besttree)
# ##srand(5678);
# hybrid = addHybridization!(besttree)
# printEdges(besttree)
# flag=false
# nocycle=true
# flag, nocycle, edgesInCycle, nodesInCycle = updateInCycle!(besttree,hybrid)
# printEdges(besttree)
# updateMajorHybrid!(besttree,hybrid)
# besttree.root=9
# plot(besttree)
# writeTopology(besttree)
# flag2, edgesGammaz = updateGammaz!(besttree,hybrid,false)
# flag3, edgesRoot = updateContainRoot!(besttree,hybrid)
# plot(besttree)


text:
module InfiniteOpt

# Import and export JuMP
import Reexport
Reexport.@reexport using JuMP

# Import the necessary packages.
import MathOptInterface
import Distributions
import DataStructures
import FastGaussQuadrature

# Make useful aliases (note we get MOI and MOIU from JuMP)
const JuMPC = JuMP.Containers
const MOIUC = MOIU.CleverDicts

# Import the Collections module
include("Collections/Collections.jl")

# Import core methods
include("datatypes.jl")
include("infinite_domains.jl")
include("scalar_parameters.jl")
include("array_parameters.jl")
include("variable_basics.jl")
include("infinite_variables.jl")
include("semi_infinite_variables.jl")
include("point_variables.jl")
include("finite_variables.jl")
include("expressions.jl")
include("measures.jl")

# Import and export MeasureToolbox
include("MeasureToolbox/MeasureToolbox.jl")
Reexport.@reexport using .MeasureToolbox

# import more core methods
include("derivatives.jl")
include("constraints.jl")
include("macros.jl")
include("objective.jl")
include("measure_expansions.jl")
include("derivative_evaluations.jl")
include("optimize.jl")
include("results.jl")
include("show.jl")
include("utilities.jl")
include("general_variables.jl")

# Import and export TranscriptionOpt
include("TranscriptionOpt/TranscriptionOpt.jl")
Reexport.@reexport using .TranscriptionOpt

# Deprecations introduced with v0.4.0
Base.@deprecate(ParameterBounds, DomainRestrictions)
Base.@deprecate(has_parameter_bounds, has_domain_restrictions)
Base.@deprecate(parameter_bounds, domain_restrictions)
Base.@deprecate(set_parameter_bounds, set_domain_restrictions)
Base.@deprecate(add_parameter_bounds, add_domain_restrictions)
Base.@deprecate(delete_parameter_bounds, delete_domain_restrictions)
Base.@deprecate(IntervalSet, IntervalDomain)
Base.@deprecate(UniDistributionSet, UniDistributionDomain)
Base.@deprecate(MultiDistributionSet, MultiDistributionDomain)
Base.@deprecate(CollectionSet, CollectionDomain)
Base.@deprecate(supports_in_set, supports_in_domain)
Base.@deprecate(collection_sets, collection_domains)
Base.@deprecate(infinite_set, infinite_domain)
Base.@deprecate(set_infinite_set, set_infinite_domain)
Base.@deprecate(InfiniteParameterFunction, ParameterFunction)

for op in (:has_parameter_bounds, :parameter_bounds, :set_parameter_bounds,
           :add_parameter_bounds, :delete_parameter_bounds)
    @eval begin
        function $op(::JuMP.AbstractVariableRef, args...; kwargs...)
            error("The use of parameter bounds (now called domain ",
                  "domain restrictions) has been discontinued for finite ",
                  "variables. The preferred syntax is to specify ",
                  "`DomainRestrictions` in connection with constraints.")
        end
    end
end

macro independent_parameter(model, args...)
    error("`@independent_parameter` has been deprecated in favor of ",
          "`@infinite_parameter`. ",
          "\n\nOld Syntax: `@independent_parameter(model, par[...] in set, ",
          "kwargs...)`\nNew Syntax: `@infinite_parameter(model, par[...] in ",
          "set, independent = true, kwargs...)`.")
end
macro dependent_parameters(model, args...)
    error("`@Dependent_parameters` has been deprecated in favor of ",
          "`@infinite_parameter`. ",
          "\n\nOld Syntax: `@dependent_parameters(model, par[...] in set, ",
          "kwargs...)`\nNew Syntax: `@infinite_parameter(model, par[...] in ",
          "set, kwargs...)`.")
end
macro infinite_variable(model, args...)
    error("`@infinite_variable` has been deprecated in favor of `@variable`. ",
          "\n\nOld Syntax: `@infinite_variable(model, var[idxs...](params...), ",
          "args..., kwargs...)`\nNew Syntax: `@variable(model, var[idxs...], ",
          "Infinite(params...), args..., kwargs...)`.")
end
macro point_variable(model, args...)
    error("`@point_variable` has been deprecated in favor of `@variable`. ",
          "\n\nOld Syntax: `@point_variable(model, ivref(param_values...), ",
          "var_expr, args..., kwargs...)`\nNew Syntax: `@variable(model, ",
          "var_expr, Point(ivref, param_values...), args..., kwargs...)`.")
end
macro hold_variable(model, args...)
    error("`@hold_variable` has been deprecated in favor of `@variable`.",
          "\n\nOld Syntax: `@hold_variable(model, var_expr, args..., ",
          "kwargs...)`\nNew Syntax: `@variable(model, var_expr, args..., ",
          "kwargs...)`.")
end
macro derivative_variable(model, args...)
    error("`@derivative_variable` has been deprecated in favor of `@variable`. ",
          "\n\nOld Syntax: `@derivative_variable(model, d(ivref)/d(pref), ",
          "var_expr, kwargs...)`\nNew Syntax: `@variable(model, var_expr, ",
          "Deriv(ivref, pref), kwargs...)`.")
end
macro BDconstraint(model, args...)
    error("`@BDconstraint` has been deprecated in favor of `@constraint`. ",
          "\n\nOld Syntax: `@BDconstraint(model, ref_expr(conditions...), ",
          "constr_expr, kwargs...)`\nNew Syntax: `@constraint(model, ref_expr, ",
          "constr_expr, DomainRestrictions(conditions...), kwargs...)`.")
end
macro set_parameter_bounds(args...)
    error("`@set_parameter_bounds` has been deprecated in favor of ",
          "`set_domain_restrictions`. Macro based modification is no longer ",
          "supported.")
end
macro add_parameter_bounds(args...)
    error("`@add_parameter_bounds` has been deprecated in favor of ",
          "`add_domain_restrictions`. Macro based modification is no longer ",
          "supported.")
end

# Define additional stuff that should not be exported
const _EXCLUDE_SYMBOLS = [Symbol(@__MODULE__), :eval, :include]

# Following JuMP, export everything that doesn't start with a _
for sym in names(@__MODULE__, all = true)
    sym_string = string(sym)
    if sym in _EXCLUDE_SYMBOLS || startswith(sym_string, "_") || startswith(sym_string, "@_")
        continue
    end
    if !(Base.isidentifier(sym) || (startswith(sym_string, "@") && Base.isidentifier(sym_string[2:end])))
        continue
    end
    @eval export $sym
end

end # module


text:
density_flux_code =
"""
function density_flux!(::Type{Val{:jac}},jac,P,Pmat,Tmat,Jx,Jy,V_x,V_y,dx,dy)
    ny = length(indices(Tmat)[1])-2
    nx = length(indices(Tmat)[2])-2
    diag_0_indices = diagind(jac,0)  # dPij.
    diag_m1_indices = [diagind(jac,nx*ny-1);diagind(jac,-1)]  # dPijm1.
    # diag_mnxm1_indices = diagind(jac,-nx)  # dPim1jm1.
    diag_mnx_indices = [diagind(jac,nx*ny-nx);diagind(jac,-nx)]  # dPim1j.
    # diag_mnxp1_indices = diagind(jac,-nx)  # dPim1jp1.
    diag_p1_indices = [diagind(jac,1);diagind(jac,-nx*ny+1)]  # dPijp1.
    diag_pnx_indices = [diagind(jac,nx);diagind(jac,-nx*ny+nx)]  # dPip1j
    diag_p1_indices = [diagind(jac,1);diagind(jac,-nx*ny+1)]  # dPijp1.
    diag_pnx_indices = [diagind(jac,nx);diagind(jac,-nx*ny+nx)]  # dPip1j
    for row in 1:nx*ny
        i,j = ind2sub((nx,ny),row)
        jac[diag_0_indices[row]] = $(dP[1])
        jac[diag_m1_indices[row]] = $(dP[2])
        jac[diag_mnx_indices[row]] = $(dP[3])
        jac[diag_p1_indices[row]] = $(dP[4])
        jac[diag_pnx_indices[row]] = $(dP[5])
    end
    jac[diagind(jac,nx-1)[1:nx:end-nx]] .= jac[diagind(jac,-1)[nx:nx:end]]
    jac[diagind(jac,-nx+1)[1:nx:end-nx]] .= jac[diagind(jac,1)[nx:nx:end]]
    jac[diagind(jac,-1)[nx:nx:end]] = 0
    jac[diagind(jac,1)[nx:nx:end]] = 0
    jac[diagind(jac,nx*ny-1)] = 0
    jac[diagind(jac,-nx*ny+1)] = 0
    jac[diagind(jac,nx*ny-nx)] = 0
    jac[diagind(jac,-nx*ny+nx)] = 0
    jac
end
"""
clipboard(density_flux_code)
##
dT = Array{String}(5)
for (index,var) in enumerate([T_22,T_12,T_21,T_32,T_23])
    var = diff(density_flux,var) |> simplify |> string
    for replacement in replacements
        var = replace(var,replacement...)
    end
    var = replace(var," ","")
    dT[index] = var
end
dT

temperature_flux_code =
"""
function temperature_flux!(::Type{Val{:jac}},jac,T,Pmat,Tmat,Jx,Jy,V_x,V_y,dx,dy)
    ny = length(indices(Tmat)[1])-2
    nx = length(indices(Tmat)[2])-2
    diag_0_indices = diagind(jac,0)  # dPij.
    diag_m1_indices = [diagind(jac,nx*ny-1);diagind(jac,-1)]  # dPijm1.
    # diag_mnxm1_indices = diagind(jac,-nx)  # dPim1jm1.
    diag_mnx_indices = [diagind(jac,nx*ny-nx);diagind(jac,-nx)]  # dPim1j.
    # diag_mnxp1_indices = diagind(jac,-nx)  # dPim1jp1.
    diag_p1_indices = [diagind(jac,1);diagind(jac,-nx*ny+1)]  # dPijp1.
    diag_pnx_indices = [diagind(jac,nx);diagind(jac,-nx*ny+nx)]  # dPip1j
    diag_p1_indices = [diagind(jac,1);diagind(jac,-nx*ny+1)]  # dPijp1.
    diag_pnx_indices = [diagind(jac,nx);diagind(jac,-nx*ny+nx)]  # dPip1j
    for row in 1:nx*ny
        i,j = ind2sub((nx,ny),row)
        jac[diag_0_indices[row]] = $(dT[1])
        jac[diag_m1_indices[row]] = $(dT[2])
        jac[diag_mnx_indices[row]] = $(dT[3])
        jac[diag_p1_indices[row]] = $(dT[4])
        jac[diag_pnx_indices[row]] = $(dT[5])
    end
    jac[diagind(jac,nx-1)[1:nx:end-nx]] .= jac[diagind(jac,-1)[nx:nx:end]]
    jac[diagind(jac,-nx+1)[1:nx:end-nx]] .= jac[diagind(jac,1)[nx:nx:end]]
    jac[diagind(jac,-1)[nx:nx:end]] = 0
    jac[diagind(jac,1)[nx:nx:end]] = 0
    jac[diagind(jac,nx*ny-1)] = 0
    jac[diagind(jac,-nx*ny+1)] = 0
    jac[diagind(jac,nx*ny-nx)] = 0
    jac[diagind(jac,-nx*ny+nx)] = 0
    jac
end
"""
clipboard(temperature_flux_code)

function stencil_indices!(inds, i, j, nn)::Array{Bool, 1}
    inds[:] .= true
    j < 1 && (inds[1:3] .= false)
    j > nn && (inds[7:9] .= false)
    # i <= 1 && (inds[[1,4]] .= false)
    # i >= nn && (inds[[6,9]] .= false)
    inds
end

function wrap_index(index,i,j,nx,ny)
    # 1 <= j <= ny && return index
    # i < 1 && (index = index+nx+1)
    # i > nx && (index = index-nx-1)
    1 <= index <= nx*ny ? index : mod(index,nx)+1
end

function array_index(row,i,j,nx,ny)
    ((1 <= i <= nx) && (1 <= j <= ny)) && return sub2ind((nx,ny),i,j)
    i == 0 && return sub2ind((nx,ny),nx,mod(j,ny)+1)
    i == nx+1 && return sub2ind((nx,ny),1,mod(j,ny)+1)
    j == 0 && return row
    j == ny+1 && return row
    # i == -1 &&
    # i < 1 && return mod(sub2ind((nx,ny),i,j),nx*ny)+1
    # i > nx && return mod(-sub2ind((nx,ny),i,j),nx)+1
    # 1 <= row <= nx*ny ? row : mod(row,nx)+1
end

# function density_flux!(::Type{Val{:jac}},jac,P,Pmat,Tmat,Jx,Jy,V_x,V_y,dx,dy)
#     # @show "hi"
#     ny = length(indices(Tmat)[1])-2
#     nx = length(indices(Tmat)[2])-2
#     inds = Array{Bool}(9)
#     arrayinds = Array{Int64}(9)
#     jac[:,:] = 0.0
#     for row in 1:nx*ny
#         i,j = ind2sub((nx,ny),row)
#         # arrayinds[:] .= Int64[wrap_index(sub2ind((nx,ny),ii,jj),ii,jj,nx,ny) for ii in (i-1):(i+1),
#         #                                                                    jj in (j-1):(j+1)][:]
#         arrayinds[:] .= Int64[array_index(row,ii,jj,nx,ny) for ii in (i-1):(i+1),
#                                                                jj in (j-1):(j+1)][:]
#         # stencil_indices!(inds,i,j,nx)  # TODO make this work for nx != ny.
#         # inds = [1 <= jj <= ny for ii in (i-1):(i+1), jj in (j-1):(j+1)][:]
#         inds = fill(true,9)
#         jac[row,arrayinds[inds]] .= [0.0;
#                                (Tmat[i,j-1]+Tmat[i,j]-V_y[i,j]*dy)/(2*dx*dy);
#                                0.0;
#                                (Tmat[i-1,j]+Tmat[i,j]-V_x[i,j]*dx)/(2*dx*dy);
#                                -(Tmat[i-1,j]+Tmat[i,j-1]+4*Tmat[i,j]+Tmat[i,j+1]+Tmat[i+1,j]+dx*(V_x[i,j]-V_x[i+1,j])+dy*(V_y[i,j]-V_y[i,j+1]))/(2*dx*dy);
#                                (Tmat[i,j]+Tmat[i+1,j]+V_x[i+1,j]*dx)/(2*dx*dy);
#                                0.0;
#                                (Tmat[i,j]+Tmat[i,j+1]+V_y[i,j+1]*dy)/(2*dx*dy);
#                                0.0][inds]
#     end
#     # jac[:,1:ny:nx*ny] = 0
#     jac[diagind(jac,1)[nx:nx:nx*ny-nx+1]] = 0
#     # jac[end,1] = 0
#     jac
# end
#
# density_flux!(::Type{Val{:jac}},jac,P,params,t) = density_flux!(Val{:jac},jac,P,params...)

# function flux_jacobian!(jac,du,u,A,coupling,T0,V,Vxshift,Vyshift,V_x,V_y,dx,dy)
#     # TODO Make T0 a vector.
#     ny,nx = size(V_x)[1],size(V_x)[2]
#     P = reshape(view(u,1:nx*ny),ny,nx)
#     T = reshape(view(u,(nx*ny+1):2nx*ny),ny,nx)
#
#     dP22 = fill(zero(eltype(V_x)),ny,nx)
#     for i in 2:nx, j in 2:ny
#         dP22[j,i] = (V_x[j,i]-V_x[j,i-1])/2dy + (V_x[j,i]-V_y[j-1,i])/2dx
#     end
#     jac[diagind(jac,0)[1:nx*ny]] .= dP22[:]
# end
#
#
# function get_flux(potential, gradpotential_x, gradpotential_y, coupling)
#     @syms _x _y real=true
#     @syms dx dy real=true positive=true
#     @syms P11 P12 P13 P21 P22 P23 P31 P32 P33 real=true positive=true
#     @syms T11 T12 T13 T21 T22 T23 T31 T32 T33 real=true positive=true
#
#     density_flux = (Jx22-Jx12)/dy + (Jy22-Jy21)/dx
#     energy_flux = A*(((Jx22*(V22+V12)/2 - Jx12*Vxshift[1:ny,2:nx+1]
#                                -coupling*(T_x[1:ny,1:nx]-T_x[1:ny,2:nx+1]))/dy
#                              +(Jy[1:ny,1:nx].*Vyshift[1:ny,1:nx]-Jy[2:ny+1,1:nx].*Vyshift[2:ny+1,1:nx]
#                                -coupling*(T_y[1:ny,1:nx]-T_y[2:ny+1,1:nx]))/dx))
#     temperature_flux = energy_flux-A*V22*density_flux
#
#     density_flux, temperature_flux
# end


text:
# I wish yo find an explanation for the infinite loop thing problem. Also motivated because this issue arrise when adding noise and in Tim's code this wasn't an Issue.
# I modified to noise model so it is normalized and we can actually speak of noise strength. I believe this to be the cause, as in Tim's code there was no normalization issue. I want to quantify this missmatch and maybe conclude that the infinite loop thing arise just because there is too much noise involved.
push!(LOAD_PATH, "../models")
push!(LOAD_PATH, "../")
push!(LOAD_PATH, "../SparseInverseProblems/src")
using SparseInverseProblemsMod
using SuperResModels
using TestCases
using Utils
using PyCall
@pyimport numpy as np
@pyimport matplotlib.pyplot as plt

## Parameters

# Static Parameters
f_c = 20
x_max = 1.0
num_x = 10*f_c
filter_x = ones(2*f_c+1)

# Dynamic Parameters
K  = 2
tau = 1.0/(K*1.0)
v_max = 0.5
num_v = 10

# Variability of weights
minWeights = 1
maxWeights = 1

# Variability of number of particles
min_number_part = 4
max_number_part = 10

# Number of generated examples
num_trials = 1000

# test case
test_case = () -> TestCases.cloud_1d_full(x_max, v_max, minWeights, maxWeights, K, tau, rand(min_number_part:max_number_part))
# WARNING: changing anything in the test_case will generate troubles if the rejection_sampling algorithm
# in Utils.js is not either updated for it, or disabled.

model_static = SuperResModels.Fourier1d(f_c, x_max, filter_x, num_x)
model_dynamic = SuperResModels.DynamicFourier1d(model_static, v_max, tau, K, num_v*K)

### Location of data folder
dataFolder = "data/1Dsimulations"

bins = np.load("../"*dataFolder*"/separationDistribBins.npy")
density = np.load("../"*dataFolder*"/separationDistribVal.npy")

    K = div(length(model_dynamic.times)-1,2)
    tau = maximum(model_dynamic.times)/K
    # Dynamic case
    # We generate the target measure
    results = zeros(num_trials,2)
    for i in 1:num_trials
	(thetas, weights) =  Utils.Rejection_sampling(test_case, bins, density, K, tau, model_static.x_max)
        target = Utils.generate_target(model_dynamic, thetas, weights)
        # We generate our model noise
        noise = randn(size(target))
        Normalizer = norm(target[:],2)/norm(noise[:],2)
        results[i,1] = Normalizer
	results[i,2] = length(weights)
    end

indexMax = indmax(results[:,1])
println(results)
println(" Maximum is : ", results[indexMax,1], ", with ", results[indexMax,2], " particles. ")