text:
# Invertible network layer from Putzky and Welling (2019): https://arxiv.org/abs/1911.10914
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export NetworkLoop

"""
    L = NetworkLoop(nx, ny, n_in, n_hidden, batchsize, maxiter, Ψ; k1=4, k2=3, p1=0, p2=1, s1=4, s2=1) (2D)

    L = NetworkLoop(nx, ny, nz, n_in, n_hidden, batchsize, maxiter, Ψ; k1=4, k2=3, p1=0, p2=1, s1=4, s2=1) (3D)

 Create an invertibel recurrent inference machine (i-RIM) consisting of an unrooled loop
 for a given number of iterations.

 *Input*:

 - `nx`, `ny`, `nz`, `n_in`, `batchsize`: spatial dimensions, number of channels and batchsize of input tensor

 - `n_hidden`: number of hidden units in residual blocks

 - `maxiter`: number unrolled loop iterations

 - `Ψ`: link function

 - `k1`, `k2`: stencil sizes for convolutions in the residual blocks. The first convolution
   uses a stencil of size and stride `k1`, thereby downsampling the input. The second
   convolutions uses a stencil of size `k2`. The last layer uses a stencil of size and stride `k1`,
   but performs the transpose operation of the first convolution, thus upsampling the output to
   the original input size.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`) in
   residual block

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`) in
   residual block

 *Output*:

 - `L`: invertible i-RIM network.

 *Usage:*

 - Forward mode: `η_out, s_out = L.forward(η_in, s_in, d, A)`

 - Inverse mode: `η_in, s_in = L.inverse(η_out, s_out, d, A)`

 - Backward mode: `Δη_in, Δs_in, η_in, s_in = L.backward(Δη_out, Δs_out, η_out, s_out, d, A)`

 *Trainable parameters:*

 - None in `L` itself

 - Trainable parameters in the invertible coupling layers `L.L[i]`, and actnorm layers
   `L.AN[i]`, where `i` ranges from `1` to the number of loop iterations.

 See also: [`CouplingLayerIRIM`](@ref), [`ResidualBlock`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
struct NetworkLoop <: InvertibleNetwork
    L::Union{AbstractArray{CouplingLayerIRIM, 1}, AbstractArray{CouplingLayerHINT}}
    AN::AbstractArray{ActNorm, 1}
    Ψ::Function
end

@Flux.functor NetworkLoop

# 2D Constructor
function NetworkLoop(nx, ny, n_in, n_hidden, batchsize, maxiter, Ψ; k1=4, k2=3, p1=0, p2=1, s1=4, s2=1, type="additive")

    if type == "additive"
        L = Array{CouplingLayerIRIM}(undef, maxiter)
    elseif type == "HINT"
        L = Array{CouplingLayerHINT}(undef, maxiter)
    end

    AN = Array{ActNorm}(undef, maxiter)
    for j=1:maxiter
        if type == "additive"
            L[j] = CouplingLayerIRIM(nx, ny, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)
        elseif type == "HINT"
            L[j] = CouplingLayerHINT(nx, ny, n_in, n_hidden, batchsize; logdet=false, permute="both", k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)
        end
        AN[j] = ActNorm(1)
    end

    return NetworkLoop(L, AN, Ψ)
end

# 3D Constructor
function NetworkLoop(nx, ny, nz, n_in, n_hidden, batchsize, maxiter, Ψ; k1=4, k2=3, p1=0, p2=1, s1=4, s2=1, type="additive")

    if type == "additive"
        L = Array{CouplingLayerIRIM}(undef, maxiter)
    elseif type == "HINT"
        L = Array{CouplingLayerHINT}(undef, maxiter)
    end
    AN = Array{ActNorm}(undef, maxiter)
    for j=1:maxiter
        if type == "additive"
            L[j] = CouplingLayerIRIM(nx, ny, nz, n_in, n_hidden, batchsize; k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)
        elseif type == "HINT"
            L[j] = CouplingLayerHINT(nx, ny, nz, n_in, n_hidden, batchsize; logdet=false,
                permute="both", k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2)
        end
        AN[j] = ActNorm(1)
    end

    return NetworkLoop(L, AN, Ψ)
end

# 2D Forward loop: Input (η, s), Output (η, s)
function forward(η::AbstractArray{Float32, 4}, s::AbstractArray{Float32, 4}, d, J, UL::NetworkLoop)

    # Dimensions
    nx, ny, n_s, batchsize = size(s)
    n_in = n_s + 1
    maxiter = length(UL.L)
    N = cuzeros(η, nx, ny, n_in-2, batchsize)

    for j=1:maxiter
        g = J'*(J*reshape(UL.Ψ(η), :, batchsize) - reshape(d, :, batchsize))
        g = reshape(g, nx, ny, 1, batchsize)
        gn = UL.AN[j].forward(g)   # normalize
        s_ = s + tensor_cat(gn, N)

        ηs = UL.L[j].forward(cat(η, s_; dims=3))
        η = ηs[:, :, 1:1, :]
        s = ηs[:, :, 2:end, :]
    end
    return η, s
end

# 3D Forward loop: Input (η, s), Output (η, s)
function forward(η::AbstractArray{Float32, 5}, s::AbstractArray{Float32, 5}, d, J, UL::NetworkLoop)

    # Dimensions
    nx, ny, nz, n_s, batchsize = size(s)
    n_in = n_s + 1
    maxiter = length(UL.L)
    N = cuzeros(η, nx, ny, nz, n_in-2, batchsize)

    for j=1:maxiter
        g = J'*(J*reshape(UL.Ψ(η), :, batchsize) - reshape(d, :, batchsize))
        g = reshape(g, nx, ny, nz, 1, batchsize)
        gn = UL.AN[j].forward(g)   # normalize
        s_ = s + tensor_cat(gn, N)

        ηs = UL.L[j].forward(tensor_cat(η, s_))
        η = ηs[:, :, :, 1:1, :]
        s = ηs[:, :, :, 2:end, :]
    end
    return η, s
end

# 2D Inverse loop: Input (η, s), Output (η, s)
function inverse(η::AbstractArray{Float32, 4}, s::AbstractArray{Float32, 4}, d, J, UL::NetworkLoop)

    # Dimensions
    nx, ny, n_s, batchsize = size(s)
    n_in = n_s + 1
    maxiter = length(UL.L)
    N = cuzeros(η, nx, ny, n_in-2, batchsize)

    for j=maxiter:-1:1
        ηs_ = UL.L[j].inverse(tensor_cat(η, s))
        η = ηs_[:, :, 1:1, :]
        s_ = ηs_[:, :, 2:end, :]

        g = J'*(J*reshape(UL.Ψ(η), :, batchsize) - reshape(d, :, batchsize))
        g = reshape(g, nx, ny, 1, batchsize)
        gn = UL.AN[j].forward(g)   # normalize
        s = s_ - tensor_cat(gn, N)
    end
    return η, s
end

# 3D Inverse loop: Input (η, s), Output (η, s)
function inverse(η::AbstractArray{Float32, 5}, s::AbstractArray{Float32, 5}, d, J, UL::NetworkLoop)

    # Dimensions
    nx, ny, nz, n_s, batchsize = size(s)
    n_in = n_s + 1
    maxiter = length(UL.L)
    N = cuzeros(η, nx, ny, nz, n_in-2, batchsize)

    for j=maxiter:-1:1
        ηs_ = UL.L[j].inverse(tensor_cat(η, s))
        η = ηs_[:, :, :, 1:1, :]
        s_ = ηs_[:, :, :, 2:end, :]

        g = J'*(J*reshape(UL.Ψ(η), :, batchsize) - reshape(d, :, batchsize))
        g = reshape(g, nx, ny, nz, 1, batchsize)
        gn = UL.AN[j].forward(g)   # normalize
        s = s_ - tensor_cat(gn, N)
    end
    return η, s
end

# 2D Backward loop: Input (Δη, Δs, η, s), Output (Δη, Δs, η, s)
function backward(Δη::AbstractArray{Float32, 4}, Δs::AbstractArray{Float32, 4},
    η::AbstractArray{Float32, 4}, s::AbstractArray{Float32, 4}, d, J, UL::NetworkLoop)

    # Dimensions
    nx, ny, n_s, batchsize = size(s)
    n_in = n_s + 1
    maxiter = length(UL.L)
    N = cuzeros(Δη, nx, ny, n_in-2, batchsize)
    typeof(Δs) == Float32 && (Δs = 0f0.*s)  # make Δs zero tensor

    for j=maxiter:-1:1
        Δηs_, ηs_ = UL.L[j].backward(tensor_cat(Δη, Δs), tensor_cat(η, s))

        # Inverse pass
        η = ηs_[:, :, 1:1, :]
        s_ = ηs_[:, :, 2:end, :]
        g = J'*(J*reshape(UL.Ψ(η), :, batchsize) - reshape(d, :, batchsize))
        g = reshape(g, nx, ny, 1, batchsize)
        gn = UL.AN[j].forward(g)   # normalize
        s = s_ - tensor_cat(gn, N)

        # Gradients
        Δs = Δηs_[:, :, 2:end, :]
        Δgn = Δs[:, :, 1:1, :]
        Δg = UL.AN[j].backward(Δgn, gn)[1]
        Δη = reshape(J'*J*reshape(Δg, :, batchsize), nx, ny, 1, batchsize) + Δηs_[:, :, 1:1, :]
    end
    return Δη, Δs, η, s
end

# 3D Backward loop: Input (Δη, Δs, η, s), Output (Δη, Δs, η, s)
function backward(Δη::AbstractArray{Float32, 5}, Δs::AbstractArray{Float32, 5},
    η::AbstractArray{Float32, 5}, s::AbstractArray{Float32, 5}, d, J, UL::NetworkLoop)

    # Dimensions
    nx, ny, nz, n_s, batchsize = size(s)
    n_in = n_s + 1
    maxiter = length(UL.L)
    N = cuzeros(Δη, nx, ny, nz, n_in-2, batchsize)
    typeof(Δs) == Float32 && (Δs = 0f0.*s)  # make Δs zero tensor

    for j=maxiter:-1:1
        Δηs_, ηs_ = UL.L[j].backward(tensor_cat(Δη, Δs), tensor_cat(η, s))

        # Inverse pass
        η = ηs_[:, :, :, 1:1, :]
        s_ = ηs_[:, :, :, 2:end, :]
        g = J'*(J*reshape(UL.Ψ(η), :, batchsize) - reshape(d, :, batchsize))
        g = reshape(g, nx, ny, nz, 1, batchsize)
        gn = UL.AN[j].forward(g)   # normalize
        s = s_ - tensor_cat(gn, N)

        # Gradients
        Δs = Δηs_[:, :, :, 2:end, :]
        Δgn = Δs[:, :, :, 1:1, :]
        Δg = UL.AN[j].backward(Δgn, gn)[1]
        Δη = reshape(J'*J*reshape(Δg, :, batchsize), nx, ny, nz, 1, batchsize) + Δηs_[:, :, :, 1:1, :]
    end
    return Δη, Δs, η, s
end

# Clear gradients
function clear_grad!(UL::NetworkLoop)
    maxiter = length(UL.L)
    for j=1:maxiter
        clear_grad!(UL.L[j])
        clear_grad!(UL.AN[j])
        UL.AN[j].s.data = nothing
        UL.AN[j].b.data = nothing
    end
end

# Get parameters (do not update actnorm weights)
function get_params(UL::NetworkLoop)
    maxiter = length(UL.L)
    p = get_params(UL.L[1])
    if maxiter > 1
        for j=2:maxiter
            p = cat(p, get_params(UL.L[j]); dims=1)
        end
    end
    return p
end

# Put parameters
function put_params!(UL::NetworkLoop, Params::Array{Any,1})
    maxiter = length(UL.L)

    if UL.L[1] isa CouplingLayerHINT
        idx_s = 1
        counter = 0
        nlayers_x = length(UL.L[1].CL)
        isnothing(UL.L[1].C) ? counter += 5*nlayers_x : counter += 5*nlayers_x + 3
        put_params!(UL.L[1], Params[idx_s:idx_s+counter-1])

        if maxiter > 1
            for j=2:maxiter
                idx_s += counter
                counter = 0
                nlayers_x = length(UL.L[j].CL)
                isnothing(UL.L[j].C) ? counter += 5*nlayers_x : counter += 5*nlayers_x + 3
                put_params!(UL.L[j], Params[idx_s:idx_s+counter-1])
            end
        end

    elseif UL.L[1] isa CouplingLayerIRIM
        idx_s = 1
        counter = 8
        put_params!(UL.L[1], Params[idx_s:idx_s+counter-1])

        if maxiter > 1
            for j=2:maxiter
                idx_s += counter
                put_params!(UL.L[j], Params[idx_s:idx_s+counter-1])
            end
        end
    end
end


text:
"""
Compute the expected count of each item categories.
For the dentype = :DN
"""
function expectedcount(post, pind, nind, iind, pars, U)
    N = zeros(Float64, size(pars, 1), size(post, 2))
    r = [zeros(Float64, length(unique(k)), size(post, 2)) for k in eachcol(U)]
    for (i, par) in enumerate(pars)
        group = par.group
        ind = vcat([pind[g] for g in group])[1]
        respind = [i ∈ iind[p] for p in ind]
        ind = ind[respind]
        x::Array{Int64, 1} = @views U[findall(respind), i] .+ one(eltype(U[:, i])) # response vector (original data are starts from 0.)
        for (s::Int64, p::Int64) in zip(x, ind)
            r[i][s, :] += @views post[p, :] .* nind[p]
            N[i, :] += @views post[p, :] .* nind[p]
        end
    end
    return N, r
end

@inline function mvn(x)::eltype(x)
    length(x) == 1 ? pdf(Normal(), x)[1] : pdf(MvNormal(length(q), 1.0), x)[1] # nomal density
end

"""
Expected log likelihood function. For DN mode.
"""
function ell(newpar, pars, bgh::DNq, r, N, q)::eltype(newpar)
    lnp::eltype(newpar) = zero(eltype(newpar))
    for k in axes(r, 1)
        for l in axes(bgh.n, 1)
            p = @views irf(newpar, pars, bgh.n[l, :], k-1, q)
            if p < 0.0
                @error "Calculated probability is negative ($p)! See your model setting or data."
            end
            lnp += log(p) * r[k, l] / N[l]
        end
    end
    if typeof(pars) <: guessing
        lnp += log(pdf(Beta(2, 8), newpar[end]))
    end
    if typeof(pars) <: LPE
        lnp += log(pdf(Gamma(2.0, 1.0), exp(newpar[end])))
    end
    return lnp
end

"""
Optimize by Newton method.
"""
function Newton(pars, g, H, j, q; atol = 1e-5, N = 5)
    pars′ = take(pars[j], q) - H(j)\g(j)
    t = 0
    while !all(abs.(take(pars[j], q) .- pars′) .< atol) && t ≤ N
        t += 1
        distribute!(pars′, pars[j], q)
        pars′ = take(pars[j], q) - H(j)\g(j)
    end
    distribute!(pars′, pars[j], q)
end

ext: jl


text:
export AbLoss, MSE, MAD, LOGIT, altloss

abstract type AbLoss end
# Notation
#       y = Observed target
#       z = fit Σ_j ϕ(X, ω_j)
# Methods
#       Loss evaluation
#       Generalized Residual
valnderv(l::AbLoss, y, z) = l(y, z), derv(l, y, z)
altloss(l::AbLoss) =  (_...) -> NaN

################################################################################
type MSELoss <: AbLoss
end

const MSE = MSELoss()

(l::MSELoss)(y, z) = sum(abs2, z-y)/2
derv(::MSELoss, y, z) = z-y
derv2(::MSELoss, y, z) = 1.
function valnderv(::MSELoss, y, z)
  diff = z-y
  sum(abs2, diff)/2, diff
end

################################################################################
type MADLoss <: AbLoss
end

const MAD = MADLoss()

(l::MADLoss)(y, z) =  sum(abs, z-y)
derv(::MADLoss, y, z) = sign.(z-y)
derv2(::MADLoss, y, z) = 0.
function valnderv(::MADLoss, y, z)
    diff = z-y
    sum(abs, diff), sign.(diff)
end

################################################################################
type LogitLoss <: AbLoss
end

const LOGIT = LogitLoss()

log1plusexp(x) = (x<0.) ? log1p(exp(x)) : x+log1p(exp(-x))
sigmoid(x) = (x<0.) ? exp(x)/(1.+exp(x)) : 1./(1.+exp(-x))
(l::LogitLoss)(y, z) =  sum(log1plusexp.(z)) - y'*z
derv(::LogitLoss, y, z) = sigmoid.(z) - y
derv2(::LogitLoss, y, z) = (s=sigmoid.(z); (1-s).*s)
"""
bits(sigmoid(-64)) = "0011101000101001011010011101010001110011001000011110010011001100"
bits(1-(1-sigmoid(-64))) = "0000000000000000000000000000000000000000000000000000000000000000"
"""

misclassification(y, z, threshold=.5) = sum(y .!= (sigmoid.(z) .≥ threshold))
altloss(::LogitLoss) =  misclassification
################################################################################


text:
##########################
## iteration by columns ##
##########################

## column iterator without dates
##------------------------------

## iterate over columns without dates equal to DataFrames:
## (nam, col)
immutable TdColumnIterator
    td::AbstractTimedata
end
import DataFrames.eachcol
eachcol(td::AbstractTimedata) = TdColumnIterator(td)

Base.start(itr::TdColumnIterator) = 1
Base.done(itr::TdColumnIterator, j::Int) = j > size(itr.td, 2)
function Base.next(itr::TdColumnIterator, j::Int)
    return ((names(itr.td)[j], itr.td.vals[j]), j + 1)
end
Base.size(itr::TdColumnIterator) = (size(itr.td, 2), )
Base.length(itr::TdColumnIterator) = size(itr.td, 2)
Base.getindex(itr::TdColumnIterator, j::Any) = itr.td.vals[j]

## column iterator with dates
##---------------------------

## iterate over variable as TimeData object
immutable TdVariableIterator
    td::AbstractTimedata
end
eachvar(td::AbstractTimedata) = TdVariableIterator(td)

Base.start(itr::TdVariableIterator) = 1
Base.done(itr::TdVariableIterator, j::Int) = j > size(itr.td, 2)
function Base.next(itr::TdVariableIterator, j::Int)
    return (itr.td[:, j], j + 1)
end
Base.size(itr::TdVariableIterator) = (size(itr.td, 2), )
Base.length(itr::TdVariableIterator) = size(itr.td, 2)
Base.getindex(itr::TdVariableIterator, j::Any) = itr.td[:, j]

#######################
## iteration by rows ##
#######################

## row without date information
##-----------------------------

## iterate over rows as DataFrames without date information:
immutable TdRowIterator
    td::AbstractTimedata
end
import DataFrames.eachrow
eachrow(td::AbstractTimedata) = TdRowIterator(td)

Base.start(itr::TdRowIterator) = 1
Base.done(itr::TdRowIterator, i::Int) = i > size(itr.td, 1)
function Base.next(itr::TdRowIterator, i::Int)
    return (itr.td.vals[i, :], i + 1)
end
Base.size(itr::TdRowIterator) = (size(itr.td, 1), )
Base.length(itr::TdRowIterator) = size(itr.td, 1)
Base.getindex(itr::TdRowIterator, i::Any) = itr.td.vals[i, :]

## row with date information
##--------------------------

## iterator over rows as TimeData object
immutable TdDateIterator
    td::AbstractTimedata
end
eachdate(td::AbstractTimedata) = TdDateIterator(td)

Base.start(itr::TdDateIterator) = 1
Base.done(itr::TdDateIterator, i::Int) = i > size(itr.td, 1)
function Base.next(itr::TdDateIterator, i::Int)
    return (itr.td[i, :], i + 1)
end
Base.size(itr::TdDateIterator) = (size(itr.td, 1), )
Base.length(itr::TdDateIterator) = size(itr.td, 1)
Base.getindex(itr::TdDateIterator, i::Any) = itr.td[i, :]

#################################
## iteration over each element ##
#################################

## to be implemented
## column-first iterator over individual entries without metadata
immutable TdEntryIterator
    td::AbstractTimedata
end
eachentry(td::AbstractTimedata) = TdEntryIterator(td)

Base.start(itr::TdEntryIterator) = 1
function Base.done(itr::TdEntryIterator, i::Int)
    return i > size(itr.td, 1)*size(itr.td, 2)
end
function Base.next(itr::TdEntryIterator, i::Int)
    rowInd, colInd = ind2sub(size(itr.td), i)
    return (get(itr.td, rowInd, colInd), i + 1)
end
Base.size(itr::TdEntryIterator) = (size(itr.td, 1)*size(itr.td, 2), )
Base.length(itr::TdEntryIterator) = size(itr.td, 1)*size(itr.td, 2)
function Base.getindex(itr::TdEntryIterator, i::Any)
    rowInd, colInd = ind2sub(size(itr.td), i)
    return get(itr.td, rowInd, colInd)
end

## column-first iterator over individual entries with metadata
immutable TdObsIterator
    td::AbstractTimedata
end
eachobs(td::AbstractTimedata) = TdObsIterator(td)

Base.start(itr::TdObsIterator) = 1
function Base.done(itr::TdObsIterator, i::Int)
    return i > size(itr.td, 1)*size(itr.td, 2)
end
function Base.next(itr::TdObsIterator, i::Int)
    rowInd, colInd = ind2sub(size(itr.td), i)
    return (itr.td[rowInd, colInd], i + 1)
end
Base.size(itr::TdObsIterator) = (size(itr.td, 1)*size(itr.td, 2), )
Base.length(itr::TdObsIterator) = size(itr.td, 1)*size(itr.td, 2)
function Base.getindex(itr::TdObsIterator, i::Any)
    rowInd, colInd = ind2sub(size(itr.td), i)
    return itr.td[rowInd, colInd]
end


text:
using TypeDBClient_Speed
using Documenter

DocMeta.setdocmeta!(TypeDBClient_Speed, :DocTestSetup, :(using TypeDBClient_Speed); recursive=true)

makedocs(;
    modules=[TypeDBClient_Speed],
    authors="Frank Urbach",
    repo="https://GitHub.com/FrankUrbach/TypeDBClient_Speed.jl/blob/{commit}{path}#{line}",
    sitename="TypeDBClient_Speed.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://FrankUrbach.github.io/TypeDBClient_Speed.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="GitHub.com/FrankUrbach/TypeDBClient_Speed.jl",
)


text:
module SparseUtilsTests

using SparseArrays
using SparseMatricesCSR
using PartitionedArrays
using LinearAlgebra
using Test

function test_mat(T)

  Tv = eltype(T)
  Ti = indextype(T)

  I = Ti[1,2,5,4,1]
  J = Ti[3,6,1,1,3]
  V = Tv[4,5,3,2,5]
  m = 7
  n = 6

  A = sparse(I,J,V,m,n)
  B = compresscoo(A,I,J,V,m,n)
  @test typeof(B) == typeof(A)
  @test A == B

  B = compresscoo(SparseMatrixCSC{Float64,Int},I,J,V,m,n)
  @test typeof(B) == SparseMatrixCSC{Float64,Int}
  @test A == B

  B = compresscoo(T,I,J,V,m,n)
  @test typeof(B) == T
  @test A == B

  i,j,v = findnz(B)
  for (k,(ki,kj,kv)) in enumerate(nziterator(B))
    @test i[k] == ki
    @test j[k] == kj
    @test v[k] == kv
    @test nzindex(B,ki,kj) == k
  end

  rows = [4,2,3]
  inv_rows = zeros(Int,size(B,1))
  inv_rows[rows] = 1:length(rows)
  cols = [6,2,5,1]
  inv_cols = zeros(Int,size(B,2))
  inv_cols[cols] = 1:length(cols)
  indices = (rows,cols)
  inv_indices = (inv_rows,inv_cols)
  flag = (1,1)
  C = PartitionedArrays.SubSparseMatrix(B,indices,inv_indices,flag)

  x = rand(Tv,size(C,2))
  y = similar(x,Tv,size(C,1))
  mul!(y,C,x)
  @test y ≈ A[rows,cols]*x
  @test C*x ≈ A[rows,cols]*x

end

test_mat(SparseMatrixCSC{Float64,Int})
test_mat(SparseMatrixCSC{Float32,Int32})
test_mat(SparseMatrixCSR{1,Float64,Int})
test_mat(SparseMatrixCSR{1,Float32,Int32})
test_mat(SparseMatrixCSR{0,Float64,Int})
test_mat(SparseMatrixCSR{0,Float32,Int32})

end # module


text:
n = 300
x = randexp(n)
hist(x, nbins=30, fill="lemonchiffon", ecol="navy", norm="pdf")
plot!(x->exp(-x), 0, maximum(x), lw=0.05)


text:
using Compose: x_measure, y_measure

function check_arguments(arg::Vector, len)
    if length(arg)>1
        @assert length(arg) == len
    else
        arg = fill(arg[1], len)
    end
    arg
end
check_arguments(arg, len) = fill(arg, len)

struct HLineGeometry <: Gadfly.GeometryElement
    color::Union{Vector, Color, (Void)}
    size::Union{Vector, Measure, (Void)}
    style::Union{Vector, Symbol, (Void)}
    tag::Symbol

    HLineGeometry(color, size, style, tag) =
            new(color === nothing ? nothing : Gadfly.parse_colorant(color), size, style, tag)
end
HLineGeometry(; color=nothing, size=nothing, style=nothing, tag=empty_tag) =
        HLineGeometry(color, size, style, tag)

const hline = HLineGeometry

element_aesthetics(::HLineGeometry) = [:yintercept]

# Generate a form for the hline geometry
function render(geom::HLineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.hline", aes, :yintercept)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size
    style = geom.style === nothing ? theme.line_style : geom.style

    color = check_arguments(color, length(aes.yintercept))
    size = check_arguments(size, length(aes.yintercept))
    style = check_arguments(style, length(aes.yintercept))

    style = map(Gadfly.get_stroke_vector, style)

    root = compose(context(), svgclass("xfixed"))
    for (idx,y) in enumerate(aes.yintercept)
        compose!(root, (context(),
            Compose.line([(0w, y_measure(y)), (1w, y_measure(y))], geom.tag),
            stroke(color[idx]),
            linewidth(size[idx]),
            strokedash(style[idx])))
    end
    root
end

struct VLineGeometry <: Gadfly.GeometryElement
    color::Union{Vector, Color, (Void)}
    size::Union{Vector, Measure, (Void)}
    style::Union{Vector, Symbol, (Void)}
    tag::Symbol

    VLineGeometry(color, size, style, tag) =
            new(color === nothing ? nothing : Gadfly.parse_colorant(color), size, style, tag)
end
VLineGeometry(; color=nothing, size=nothing, style=nothing, tag=empty_tag) =
        VLineGeometry(color, size, style, tag)

const vline = VLineGeometry

element_aesthetics(::VLineGeometry) = [:xintercept]

# Generate a form for the vline geometry
function render(geom::VLineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.vline", aes, :xintercept)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size
    style = geom.style === nothing ? theme.line_style : geom.style

    color = check_arguments(color, length(aes.xintercept))
    size = check_arguments(size, length(aes.xintercept))
    style = check_arguments(style, length(aes.xintercept))

    style = map(Gadfly.get_stroke_vector, style)

    root = compose(context(), svgclass("yfixed"))
    for (idx,x) in enumerate(aes.xintercept)
        compose!(root, (context(),
                Compose.line([(x_measure(x), 0h), (x_measure(x), 1h)], geom.tag),
                stroke(color[idx]),
                linewidth(size[idx]),
                strokedash(style[idx])))
    end
    root
end

struct ABLineGeometry <: Gadfly.GeometryElement
    color::Union{Vector, Color, (Void)}
    size::Union{Vector, Measure, (Void)}
    style::Union{Vector, Symbol, (Void)}
    tag::Symbol

    ABLineGeometry(color, size, style, tag) =
            new(color === nothing ? nothing : Gadfly.parse_colorant(color), size, style, tag)
end
ABLineGeometry(; color=nothing, size=nothing, style=nothing, tag::Symbol=empty_tag) =
        ABLineGeometry(color, size, style, tag)

abline = ABLineGeometry

element_aesthetics(geom::ABLineGeometry) = [:intercept, :slope]

function render(geom::ABLineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    if aes.intercept == nothing && aes.slope == nothing
        aes.intercept = [0]
        aes.slope = [1]
    elseif aes.intercept == nothing
        aes.intercept = fill(0,length(aes.slope))
    elseif aes.slope == nothing
        aes.slope = fill(1,length(aes.intercept))
    end
    Gadfly.assert_aesthetics_equal_length("Geom.line", aes, element_aesthetics(geom)...)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size
    style = geom.style === nothing ? theme.line_style : geom.style

    color = check_arguments(color, length(aes.intercept))
    size = check_arguments(size, length(aes.intercept))
    style = check_arguments(style, length(aes.intercept))

    style = map(Gadfly.get_stroke_vector, style)

    # it would've been nice to just say x0, x1 = realmin(Float64), realmax(Float64).
    # but SVG() and PDF() have silent overflow errors when plotting way outside the
    # context's bounding box.  so instead, use the extrema of the data, and hope that
    # the line extends to the edges of the graph.

    if typeof(aes.y) <: Array{Function}
        low, high = aes.xmin[1], aes.xmax[1]
    else
        xextrema = extrema(aes.x)
        yextrema = extrema(aes.y)
        low = min(xextrema[1], yextrema[1])
        high = max(xextrema[2], yextrema[2])
    end

    range = high-low
    x0 = low-range
    x1 = high+range

    y0s = [x0 * m + b for (m,b) in zip(aes.slope, aes.intercept)]
    y1s = [x1 * m + b for (m,b) in zip(aes.slope, aes.intercept)]

    root = compose(context(), svgclass("geometry"))
    for (idx,(y0,y1)) in enumerate(zip(y0s,y1s))
        compose!(root, (context(),
                Compose.line([(x_measure(x0),y_measure(y0)),
                              (x_measure(x1),y_measure(y1))], geom.tag),
                stroke(color[idx]),
                linewidth(size[idx]),
                strokedash(style[idx])))
    end
    root
end


text:
# an example on mixed integer programming
#
#   maximize x + 2 y + 5 z
#
#   s.t.  x + y + z <= 10
#         x + 2 y + z <= 15
#
#         x is continuous: 0 <= x <= 5
#         y is integer: 0 <= y <= 10
#         z is binary
#

using Gurobi

env = Gurobi.Env()

model = Gurobi.Model(env, "mip_01", :maximize)

add_cvar!(model, 1., 0., 5.)  # x
add_ivar!(model, 2., 0, 10)   # y
add_bvar!(model, 5.)          # z
update_model!(model)

add_constr!(model, ones(3), '<', 10.)
add_constr!(model, [1., 2., 1.], '<', 15.)

println(model)

optimize(model)

println("sol = $(get_solution(model))")
println("objv = $(get_objval(model))")


text:
using Reactive

function standard_push_test(non_input::Signal)
    m = map(x->2x, non_input)
    pval = number()
    push!(non_input, pval)
    step()

    @test (value(non_input)) == (pval)
    @test (value(m)) == (2pval)
end

@testset "Push to non-input nodes" begin

    a = Signal(number(); name="a")
    b = map(x -> x*x, a; name="b")

    @testset "push to map" begin
        m = map(x->2x, b)

        push!(b, 10.0)
        step()
        @test (value(b)) == (10.0)
        @test (value(m)) == (2value(b))

        push!(a, 2.0)
        step()
        @test (value(b)) == (4.0)
        @test (value(m)) == (2value(b))

        push!(b, 3.0)
        step()
        @test (value(b)) == (3.0)
        @test (value(m)) == (2value(b))
    end

    @testset "push to merge" begin
        ## Merge
        d = Signal(number(); name="d")
        e = merge(b, d, a; name="e")
        m = map(x->2x, e)
        # precedence to d
        @test (value(e)) == (value(d))
        @test (value(m)) == (2value(e))

        standard_push_test(e)
    end

    @testset "push to foldp" begin
        GC.gc()
        x = Signal(number())
        f = foldp(+, 0, x)

        standard_push_test(f)

        # Check for correct continuation
        pval = value(f)
        sval = number()
        push!(x, sval)
        step()

        @test (value(f)) == (pval + sval)
    end

    @testset "push to filter" begin
        g = Signal(0)
        pred = x -> x % 2 != 0
        h = filter(pred, 1, g)
        standard_push_test(h)
    end

    @testset "push to sampleon" begin
        # sampleon
        g = Signal(0)
        nv = number()
        push!(g, nv)
        step()
        i = Signal(true)
        j = sampleon(i, g)
        standard_push_test(j)
    end

    @testset "push to droprepeats" begin
        # droprepeats
        k = Signal(1)
        l = droprepeats(k)

        standard_push_test(l)
    end

    @testset "push to filterwhen" begin
        # filterwhen
        b = Signal(false)
        n = Signal(1)
        dw = filterwhen(b, 0, n)
        standard_push_test(dw)
    end

    @testset "push to previous" begin
        x = Signal(0)
        y = previous(x)
        standard_push_test(y)
    end

    @testset "push to delay" begin
        x = Signal(0)
        y = delay(x)
        standard_push_test(y)
    end

    @testset "bind non-input" begin
        s = Signal(1; name="sig 1")
        m = map(x->2x, s; name="m")
        s2 = Signal(3; name="sig 2")
        push!(m, 10)
        step()
        @test (value(m)) == (10)

        bind!(m, s2) #two-way bind
        @test (value(m)) == (3)
        @test (value(s2)) == (3)

        push!(m, 6)
        step()
        @test (value(m)) == (6)
        @test (value(s2)) == (6)

        push!(s2, 10)
        step()
        @test (value(m)) == (10)
        @test (value(s2)) == (10)
    end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "GR"
version = v"0.52.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/sciapp/gr.git", "90ee73c741271d964f241f098b0dd3ba18eae271"),
]

# Bash recipe for building across all platforms
script = raw"""
if test -f "$prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake"; then
    sed -i 's/_qt5gui_find_extra_libs.*AGL.framework.*//' $prefix/lib/cmake/Qt5Gui/Qt5GuiConfigExtras.cmake
fi

if [[ $target == *"mingw"* ]]; then
    winflags=-DCMAKE_C_FLAGS="-D_WIN32_WINNT=0x0f00"
fi

mkdir build
cd build
cmake $winflags -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_FIND_ROOT_PATH=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release ../gr*
VERBOSE=ON cmake --build . --config Release --target install -- -j${nproc}
install_license $WORKSPACE/srcdir/gr*/LICENSE.md

if [[ $target == *"apple-darwin"* ]]; then
    cd $prefix/lib
    ln -s libGR.so libGR.dylib
    ln -s libGR3.so libGR3.dylib
    ln -s libGRM.so libGRM.dylib
    ln -s libGKS.so libGKS.dylib
    cd ../bin
    ln -s ../Applications/gksqt.app/Contents/MacOS/gksqt ./
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("armv7l", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "windows")
]
platforms = expand_cxxstring_abis(platforms)
push!(platforms, Platform("x86_64", "macos"))

# The products that we will ensure are always built
products = [
    LibraryProduct("libGR", :libGR),
    LibraryProduct("libGR3", :libGR3),
    LibraryProduct("libGRM", :libGRM),
    LibraryProduct("libGKS", :libGKS),
    ExecutableProduct("gksqt", :gksqt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Bzip2_jll"),
    Dependency("Cairo_jll"),
    Dependency("FFMPEG_jll"),
    Dependency("Fontconfig_jll"),
    Dependency("GLFW_jll"),
    Dependency("JpegTurbo_jll"),
    Dependency("libpng_jll"),
    Dependency("Libtiff_jll"),
    Dependency("Pixman_jll"),
    Dependency("Qhull_jll"),
    Dependency("Qt_jll"),
    BuildDependency("Xorg_libX11_jll"),
    BuildDependency("Xorg_xproto_jll"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
# GCC version 7 because of ffmpeg, but building against Qt requires v8 on Windows.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8")

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
using JuMP
using HiGHS
using SCS
using LinearAlgebra

using Graphs
using GenericTensorNetworks
using Random
#g = smallgraph(:petersen)
g = random_diagonal_coupled_graph(10, 10, 0.8)

# MIT course: http://people.csail.mit.edu/moitra/854.html
function maxcut_ilp(g::SimpleGraph; optimizer=HiGHS.Optimizer)
    ws = ones(Int,ne(g))
    edgs = collect(edges(g))
    model = Model(optimizer)
    @variable model xs[1:nv(g)] Bin
    @variable model zs[1:ne(g)] Bin
    @constraint model [i=1:ne(g)] zs[i] <= xs[edgs[i].src] + xs[edgs[i].dst]
    @constraint model [i=1:ne(g)] zs[i] <= 2 - xs[edgs[i].src] - xs[edgs[i].dst]
    @objective model Max ws' * zs
    optimize!(model)
    return objective_value(model), value.(xs)
end

# see also: https://jump.dev/JuMP.jl/stable/tutorials/conic/max_cut_sdp/
function maxcut_relaxed_sdp(g::SimpleGraph; optimizer=SCS.Optimizer)
    L = laplacian_matrix(g)
    model = Model(optimizer)
    # Start with X as the identity matrix to avoid numerical issues.
    @variable model X[i=1:nv(g), j=1:nv(g)] PSD start=(i==j ? 1.0 : 0.0)
    @constraint model diag(X) .== 1
    @objective model Max 0.25 * dot(L, X)
    optimize!(model)
    @assert termination_status(model) == MOI.OPTIMAL
    return objective_value(model), value.(X)
end

function goemans_williamson_postprocess(X)
    num_vertex = size(X, 1)
    # small positive number to avoid zero
    V = cholesky(X + 1e-3 * I).U
    # Generate random vector on unit sphere.
    r = rand(size(V, 1))
    r /= LinearAlgebra.norm(r)
    # Iterate over vertices, and assign each vertex to a side of cut.
    cut = ones(num_vertex)
    for i in 1:num_vertex
        if LinearAlgebra.dot(r, V[:, i]) <= 0
            cut[i] = -1
        end
    end
    return cut
end

function goemans_williamson(g::SimpleGraph)
    loss, X = maxcut_relaxed_sdp(g)
    @info "Got loss = $(loss)"
    return goemans_williamson_postprocess(X)
end

@time solve(MaxCut(g; optimizer=TreeSA(ntrials=3, βs=0.1:0.1:30, niters=10, sc_target=25)), SizeMax())
#function addvar!(model::Model, ::Type{T}, name::String; lower_bound::Union{T,Nothing}=nothing, upper_bound::Union{T,Nothing}=nothing, fixed_value::Union{T,Nothing}=nothing, start) where T
#end
#JuMP.add_variable(model, JuMP.build_variable(error, JuMP.VariableInfo(true, 0, false, NaN, false, NaN, false, NaN, false, true)), "x")
#model[:x] = var"#25###291"


text:
using Revise
using JuMP, EAGO

m = Model(optimizer_with_attributes(EAGO.Optimizer, "verbosity" => 1,
                                                    "output_iterations" => 1000,
                                                    "iteration_limit" => 100000,

                                                    "cp_depth" => -1,
                                                    "cp_repetitions" => -1,
                                                    "cp_forward_reverse_limit" => 2,

                                                    "cut_min_iterations" => 2,
                                                    "cut_max_iterations" => 2,
                                                    "objective_cut_on" => true,
                                                    "subgrad_tighten" => true,

                                                    "obbt_depth" => 4,
                                                    "obbt_repetitions" => 4,
                                                    "obbt_aggressive_on" => true,

                                                    "upper_bounding_depth" => 4,

                                                    "fbbt_lp_depth" => 100000,
                                                    "fbbt_lp_repetitions" => 3))

# ----- Variables ----- #
x_Idx = Any[2, 3, 4, 5]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[2], 1.0e-7)
JuMP.set_upper_bound(x[2], 0.5)
JuMP.set_lower_bound(x[3], 1.0e-7)
JuMP.set_upper_bound(x[3], 0.5)
JuMP.set_lower_bound(x[4], 1.0e-7)
JuMP.set_upper_bound(x[4], 0.5)
JuMP.set_lower_bound(x[5], 1.0e-7)
JuMP.set_upper_bound(x[5], 0.5)

# ----- Constraints ----- #
@constraint(m, e2, x[2]+x[3] == 0.5)
@constraint(m, e3, x[4]+x[5] == 0.5)

# ----- Objective ----- #
@NLobjective(m, Min, (x[2]*(log(x[2]/(x[2]+x[4]))+log(x[2]/(x[2]+0.095173*x[4])))+x[4]*(log(x[4]/(x[2]+x[4]))+log(x[4]/(0.30384*x[2]+x[4])))+(x[2]+2.6738*x[4])*log(x[2]+2.6738*x[4])+(0.374*x[2]+x[4])*log(0.374*x[2]+x[4])+2.6738*x[4]*log(x[4]/(x[2]+2.6738*x[4]))+0.374*x[2]*log(x[2]/(0.374*x[2]+x[4]))+x[3]*(log(x[3]/(x[3]+x[5]))+log(x[3]/(x[3]+0.095173*x[5])))+x[5]*(log(x[5]/(x[3]+x[5]))+log(x[5]/(0.30384*x[3]+x[5])))+(x[3]+2.6738*x[5])*log(x[3]+2.6738*x[5])+(0.374*x[3]+x[5])*log(0.374*x[3]+x[5])+2.6738*x[5]*log(x[5]/(x[3]+2.6738*x[5]))+0.374*x[3]*log(x[3]/(0.374*x[3]+x[5]))-3.6838*x[2]*log(x[2])-1.59549*x[4]*log(x[4])-3.6838*x[3]*log(x[3])-1.59549*x[5]*log(x[5])))

JuMP.optimize!(m)

backend_opt = backend(m).optimizer.model.optimizer
last_relaxed_opt = backend_opt.relaxed_optimizer
run_time = backend_opt._run_time
println("run time: $run_time")


text:
# This file has runs the normal tests and also adds tests that can only be run
# locally on a machine with a sound card. It's mostly to put the library through
# its paces assuming a human is listening.

include("runtests.jl")

# these default values are specific to my machines
if is_windows()
    default_indev = "Microphone Array (Realtek High "
    default_outdev = "Speaker/Headphone (Realtek High"
elseif is_apple()
    default_indev = "Built-in Microph"
    default_outdev = "Built-in Output"
elseif is_linux()
    default_indev = "default"
    default_outdev = "default"
end

@testset ExtendedTestSet "Local Tests" begin
    @testset "Open Default Device" begin
        println("Recording...")
        stream = PortAudioStream(2, 0)
        buf = read(stream, 5s)
        close(stream)
        @test size(buf) == (round(Int, 5 * samplerate(stream)), nchannels(stream.source))
        println("Playing back recording...")
        stream = PortAudioStream(0, 2)
        write(stream, buf)
        println("flushing...")
        flush(stream)
        close(stream)
        println("Testing pass-through")
        stream = PortAudioStream(2, 2)
        write(stream, stream, 5s)
        flush(stream)
        close(stream)
        println("done")
    end
    @testset "Samplerate-converting writing" begin
        stream = PortAudioStream(0, 2)
        write(stream, SinSource(eltype(stream), samplerate(stream)*0.8, [220, 330]), 3s)
        write(stream, SinSource(eltype(stream), samplerate(stream)*1.2, [220, 330]), 3s)
        flush(stream)
        close(stream)
    end
    @testset "Open Device by name" begin
        stream = PortAudioStream(default_indev, default_outdev)
        buf = read(stream, 0.001s)
        @test size(buf) == (round(Int, 0.001 * samplerate(stream)), nchannels(stream.source))
        write(stream, buf)
        io = IOBuffer()
        show(io, stream)
        @test String(take!(io)) == """
        PortAudio.PortAudioStream{Float32}
          Samplerate: 44100.0Hz
          Buffer Size: 4096 frames
          2 channel sink: "$default_outdev"
          2 channel source: "$default_indev\""""
        close(stream)
    end
    @testset "Error on wrong name" begin
        @test_throws ErrorException PortAudioStream("foobarbaz")
    end
    # no way to check that the right data is actually getting read or written here,
    # but at least it's not crashing.
    @testset "Queued Writing" begin
        stream = PortAudioStream(0, 2)
        buf = SampleBuf(rand(eltype(stream), 48000, nchannels(stream.sink))*0.1, samplerate(stream))
        t1 = @async write(stream, buf)
        t2 = @async write(stream, buf)
        @test wait(t1) == 48000
        @test wait(t2) == 48000
        flush(stream)
        close(stream)
    end
    @testset "Queued Reading" begin
        stream = PortAudioStream(2, 0)
        buf = SampleBuf(rand(eltype(stream), 48000, nchannels(stream.source))*0.1, samplerate(stream))
        t1 = @async read!(stream, buf)
        t2 = @async read!(stream, buf)
        @test wait(t1) == 48000
        @test wait(t2) == 48000
        close(stream)
    end
end


text:
#=
   6 digit numbers

   From
   Gunnar Blom, Lars Holst, Dennis Sandell:
   "Problems and Snapshots from the World of Probability"
   Page 19f, Problem 2.5 Problems concerning random numbers

   Given the 6 digits numbers:

   a) Problem 1
      find the probability that at least one of the digits 0..9
      appears exactly twice.
      Answer: 2943/4000 ~ 0.7358

      Distributions of variable k_occ (num:0)
      1.00000 =>    4921  (0.492100)
      0.00000 =>    2656  (0.265600)
      2.00000 =>    2322  (0.232200)
      3.00000 =>     101  (0.010100)

      Distributions of variable k_occ_m (num:0)
      1.00000 =>    7344  (0.734400)
      0.00000 =>    2656  (0.265600)

   b) Problem 2
      find the probability that at least two of the digits 0..9
      appears exactly once.
      Answer: 1179/1250 ~ 0.9432

      Distributions of variable k_occ (num:0)
      4.00000 =>    4450  (0.445000)
      2.00000 =>    2404  (0.240400)
      6.00000 =>    1488  (0.148800)
      3.00000 =>    1044  (0.104400)
      1.00000 =>     503  (0.050300)
      0.00000 =>     111  (0.011100)

      Distributions of variable k_occ_m (num:0)
      1.00000 =>    9386  (0.938600)
      0.00000 =>     614  (0.061400)

=#

using Turing # , StatsPlots, DataFrames
include("jl_utils.jl")

#
# Find the probability that at least m digits appears exactly k times
#
@model function six_digits(n, k, m)
    d ~ filldist(DiscreteUniform(0,9),n)
    digits = tzeros(10)
    for digit in 0:9
        digits[digit+1] ~ Dirac(sum([d[i] == digit ? 1 : 0 for i in 1:n]))
    end
    # Count of exactly k occurrences
    k_occ ~ Dirac(sum([digits[i] == k ? 1 : 0 for i in 1:10]))
    # Is k_occ >= m
    k_occ_m ~ Dirac(k_occ >= m)
end

# 1) Find the probability that at least one of the digits 0..9
#    appears exactly twice.
println("1) Find the probability that at least one of the digits 0..9\nappears exactly twice.")
model = six_digits(6,2,1)
num_chains = 4
num_samples = 10_000

chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(20), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), num_samples)

# display(chns)
show_var_dist_pct(chns,:k_occ)
show_var_dist_pct(chns,:k_occ_m)

# 2) Find the probability that at least two of the digits 0..9
#    appears exactly once.
println("\n2) Find the probability that at least two of the digits 0..9\nappears exactly once.")
model = six_digits(6,1,2)

chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(20), 10_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), num_samples)

# display(chns)

show_var_dist_pct(chns,:k_occ)
show_var_dist_pct(chns,:k_occ_m)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module PupperSim    #  9.716722 seconds (11.72 M allocations: 650.717 MiB, 1.66% gc time)
                    # 10.893923 seconds (13.63 M allocations: 756.021 MiB, 2.05% gc time)
                    # 13.839607 seconds (16.71 M allocations: 918.154 MiB, 1.52% gc time)

# modified from https://github.com/klowrey/MujocoSim.jl/

export loadmodel, pupper, simulate

@time using GLFW                    # 0.418646 seconds (564.20 k allocations: 33.920 MiB)
@time using MuJoCo                  # 2.657531 seconds (10.89 M allocations: 579.276 MiB, 6.15% gc time), was 0.705183 seconds (2.45 M allocations: 162.621 MiB)
      using MuJoCo.MJCore           # 0.000602 seconds (463 allocations: 28.750 KiB)
      using StaticArrays            # 0.000598 seconds (482 allocations: 29.875 KiB)
@time using FixedPointNumbers       # 0.058240 seconds (121.71 k allocations: 7.553 MiB)
@time using ColorTypes              # 0.352808 seconds (366.70 k allocations: 22.285 MiB)

const use_VideoIO = true

@static if use_VideoIO
    @time using VideoIO             # 2.877718 seconds (6.32 M allocations: 345.427 MiB, 5.85% gc time)

    const max_video_duration = 60   # max video duration in seconds
    const video_fps = 30            # frames per second determined by GLFW.GetPrimaryMonitor refresh rate
    const max_video_frames = video_fps * max_video_duration
else
    @time using FFMPEG              # 1.175542 seconds (3.00 M allocations: 160.817 MiB, 3.63% gc time)
end

@time using QuadrupedController     # 1.068143 seconds (3.38 M allocations: 171.183 MiB, 10.19% gc time)
@time using CImGui
      using CImGui.GLFWBackend
      using CImGui.OpenGLBackend
      using CImGui.GLFWBackend.GLFW
      using CImGui.OpenGLBackend.ModernGL

include("Sim.jl")

include("util.jl")
include("record_video.jl")
include("setup.jl")

"""
    loadmodel(modelfile = "model/Pupper.xml", width = 1920, height = 1080)

Loads MuJoCo XML model and starts the simulation
"""
function loadmodel(
        modelfile = normpath(joinpath(dirname(pathof(@__MODULE__)), "../model/Pupper.xml")),
        width = 1920, height = 1080
    )
    m = jlModel(modelfile)
    d = jlData(m)

    s = PupperSim.setup(m, d, width, height)

    s.modelfile = modelfile
    @info("Model file: $modelfile")

    GLFW.SetWindowRefreshCallback(s.window, (w)->render(s,w))    # Called on window resize

    return s
end

"""
    pupper(velocity = 0.4, yaw_rate = 0.0)

Creates a Robot controller with specified initial velocity and yaw_rate
"""
function pupper(velocity = 0.2, yaw_rate = 0.0)
    config = Configuration()
    config.z_clearance = 0.02     # height to pick up each foot during trot

    command = Command([velocity, 0], yaw_rate, crouch_height)
    # command.pitch = 0.1

    # Create the robot (controller and controller state)
    Robot(config, command)
end

include("simstep.jl")
include("render.jl")

# Run the simulation
"""
    simulate([s::mjSim[, robot::Robot]])

Run the simulation loop
"""
function simulate(s::mjSim = loadmodel(), robot::Union{Robot, Nothing} = pupper())
    # Assign the passed in robot to our simulator
    s.robot = robot
    robot !== nothing && toggle_activate(robot)

    # Loop until the user closes the window
    while !GLFW.WindowShouldClose(s.window)
        simstep(s)

        render(s, s.window)
        GLFW.PollEvents()
    end

    GLFW.DestroyWindow(s.window)

    return
end

"""
    simulate(modelpath::String, width = 0, height = 0, robot = nothing)

Run the simulation loop
"""
function simulate(modelpath::String, width = 0, height = 0, robot = nothing)
    simulate(loadmodel(modelpath, width, height), robot)
end

end


text:
using BackgroundSegmenter
using Base.Test

K = Gaussian(0.0, 1.0)
@test neglog_likelihood(K, 0.0) ≈ 0
@test neglog_likelihood(K, 0) ≈ 0

L = Gaussian(0, 2)
@test neglog_likelihood(L, 1) ≈ 0.5

M = Gaussian(0, 4, 4)
@test M.sortkey == 2

bg_input = [0x90 0x10 ; 0x50 0x120]
fg_input = [0x30 0x70; 0x120 0x20]
output = [0x0 0x0; 0x0 0x0]
model = MixtureModel(3, 2, 2)
apply!(model, bg_input)
for i in 1:BackgroundSegmenter.INITIALIZATION_WINDOW
    output .+= apply!(model, bg_input)
end
@test output == [0x0 0x0; 0x0 0x0]
output .= apply!(model, fg_input)
@test output == [0x0 0x1; 0x1 0x0]


text:
import Dates

@testset "Runway.jl" begin
    filename = "data\\test.runway"
    df = DDR2import.Runway.read(filename)
    @test df[1, :AIRPORT] == "BIAR"
    @test df[2, :RWY] == "19"
    @test Dates.year(df[3, :DATEACTIVE]) == 2018
    @test Dates.minute(df[4, :TIMEOPEN]) == 0
    @test Dates.hour(df[5, :TIMECLOSED]) == 23
end


text:
## Signal creation

# With no elements
sig = Signal(1:100)
@test length(sig) == 0
@test isempty(collect(sig))
@test sig.time == 1s:100s
@test samplingfreq(sig) == 1/s
@test samplingrate(sig) == 1s
@test_throws BoundsError sig[0]
@test_throws BoundsError sig[1]

# A matrix of 8 channels
sig = Signal(1:100, reshape(1:100*8, 100, 8))
@test length(sig) == 8
@test sig[1] == [1:100]
@test all([sig[i] == [(1:100)+100*(i-1)] for i=1:length(sig)])
@test length(collect(sig)) == 8
@test_throws BoundsError sig[0]
@test_throws BoundsError sig[9]
sig2 = reshape(sig, 4, 2)
@test collect(sig) == collect(sig2)

sig = Signal(1:100, [1:100 [1.:100.]*2 [1:100]*4])
@test length(sig) == 3
@test sig[1]*4 == sig[2]*2 == sig[3] == [1:100]*4
@test_throws BoundsError sig[0]
@test_throws BoundsError sig[4]

# A Signal with a group of nested Signals (like spikes with snippets)
t1 = (-.2:.025:.8)./1e3
snips = Signal(t1, randn(length(t1), 100))
@test length(snips) == 100
@test length(snips[1]) == length(snips.time) == length(t1)
@test isregular(snips)
@test samplingfreq(snips) == 40_000.0/s
@test samplingrate(snips) == (2.5e-5)s
spikes = Signal(cumsum(rand(100)), snips)
@test length(spikes) == 1
@test spikes[1] == snips
@test !isregular(spikes)

# Ensure we cannot build a Signal with bad time vectors
@test_throws ArgumentError Signal([1,2,3,2])
@test_throws ArgumentError Signal([1.,2.,2.])

## Regularization:
sig = Signal([0:1/100:2pi], [sin, cos])
@test !isregular(sig)
@test isregular(regularize(sig))

# Indexing by vectors returns Signals
sig = Signal(1:100, reshape(1:100*8, 100, 8))
subsig = sig[2:3]
@test length(subsig) == 2
@test subsig[1] == sig[2]
@test subsig[2] == sig[3]


text:
"""
    write2disk(filename::AbstractString, wv::WordVectors [;kind=:binary])

Writes embeddings to disk.

# Arguments
  * `filename::AbstractString` the embeddings file name
  * `wv::WordVectors` the embeddings

# Keyword arguments
  * `kind::Symbol` specifies whether the embeddings file is textual (`:text`)
or binary (`:binary`); default `:binary`
"""
function write2disk(fid::IO, wv::WordVectors{S,T,H};
                    kind::Symbol=:binary
                   ) where {S<:AbstractString, T<:Real, H<:Integer}
    if kind == :binary
        _write2disk_binary(fid, wv)
    elseif kind == :text
        _write2disk_text(fid, wv)
    else
        throw(ErrorException("Supported values for the kind keyword"*
                             "argument are :text and :binary."))
    end

end

function write2disk(file::AbstractString, wv::WordVectors{S,T,H};
                    kind::Symbol=:binary
                   ) where {S<:AbstractString, T<:Real, H<:Integer}
    open(file, "w") do fid
        write2disk(fid, wv, kind=kind)
    end
end

function _write2disk_binary(fid::IO, wv::WordVectors{S,T,H}
                           ) where {S<:AbstractString, T<:Real, H<:Integer}
    vector_size, vocab_size = size(wv.vectors)
    println(fid, "$vocab_size $vector_size")
    for i in 1:vocab_size
        write(fid, wv.vocab[i])
        write(fid, ' ')
        write(fid, Float32.(wv.vectors[:,i]))
        write(fid, '\n')
    end
end

function _write2disk_text(fid::IO, wv::WordVectors{S,T,H}
                         ) where {S<:AbstractString, T<:Real, H<:Integer}
    vector_size, vocab_size = size(wv.vectors)
    println(fid, "$vocab_size $vector_size")
    for i in 1:vocab_size
        nstr = join(map(string, wv.vectors[:,i]), " ")
        println(fid, "$(wv.vocab[i]) $(nstr)")
    end
end

"""
    write2disk(filename::AbstractString, wv::CompressedWordVectors [;kind=:binary])

Writes compressed embeddings to disk.

# Arguments
  * `filename::AbstractString` the embeddings file name
  * `wv::CompressedWordVectors` the embeddings

# Keyword arguments
  * `kind::Symbol` specifies whether the embeddings file is textual (`:text`)
or binary (`:binary`); default `:binary`
"""
function write2disk(fid::IO, cwv::CompressedWordVectors{Q,U,D,T,S,H};
                    kind::Symbol=:binary) where {Q,U,D,T,S,H}
    if kind == :binary
        _write2disk_binary(fid, cwv)
    elseif kind == :text
        _write2disk_text(fid, cwv)
    else
        throw(ErrorException("Supported values for the kind keyword"*
                             "argument are :text and :binary."))
    end
end

function write2disk(file::AbstractString,
                    cwv::CompressedWordVectors{Q,U,D,T,S,H};
                    kind::Symbol=:binary) where {Q,U,D,T,S,H}
    open(file, "w") do fid
        write2disk(fid, cwv, kind=kind)
    end
end

function _write2disk_binary(fid::IO,
                            cwv::CompressedWordVectors{Q,U,D,T,S,H}
                           ) where {Q,U,D,T,S,H}
    # Initialize all variables needed to write
    nrows, vocab_size = cwv.vectors.quantizer.dims
    d = size(cwv.vectors.quantizer.codebooks[1].vectors, 1)
    m = length(cwv.vectors.quantizer.codebooks)
    k = cwv.vectors.quantizer.k
    quanttype_str = string(Q)
    disttype_str = string(D)
    quanteltype_str = string(U)
    origeltype_str = string(T)

    # Start writing text information
    println(fid, "$nrows $vocab_size")
    println(fid, "$d $k $m")
    println(fid, "$quanttype_str")
    println(fid, "$disttype_str")
    println(fid, "$quanteltype_str")
    println(fid, "$origeltype_str")

    # Write vocabulary and compressed vectors
    for i in 1:vocab_size
        write(fid, cwv.vocab[i])
        write(fid, ' ')
        write(fid, cwv.vectors.data[:,i])
    end

    # Write codebooks
    for i in 1:m
        write(fid, cwv.vectors.quantizer.codebooks[i].codes)
        for j in 1:d
            write(fid, cwv.vectors.quantizer.codebooks[i].vectors[j,:])
        end
    end
    for i in 1:nrows  # rotation matrix
        write(fid, cwv.vectors.quantizer.rot[:,i])
    end
end

function _write2disk_text(fid::IO,
                          cwv::CompressedWordVectors{Q,U,D,T,S,H}
                         ) where {Q,U,D,T,S,H}
    # Initialize all variables needed to write
    nrows, vocab_size = cwv.vectors.quantizer.dims
    d = size(cwv.vectors.quantizer.codebooks[1].vectors, 1)
    m = length(cwv.vectors.quantizer.codebooks)
    k = cwv.vectors.quantizer.k
    quanttype_str = string(Q)
    disttype_str = string(D)
    quanteltype_str = string(U)
    origeltype_str = string(T)

    # Start writing text information
    println(fid, "$nrows $vocab_size")
    println(fid, "$d $k $m")
    println(fid, "$quanttype_str")
    println(fid, "$disttype_str")
    println(fid, "$quanteltype_str")
    println(fid, "$origeltype_str")

    # Write vocabulary and compressed vectors
    for i in 1:vocab_size
        nstr = join(map(string, cwv.vectors.data[:,i]), " ")
        println(fid, "$(cwv.vocab[i]) $nstr")
    end

    # Write codebooks
    for i in 1:m
        nstr = join(map(string, cwv.vectors.quantizer.codebooks[i].codes), " ")
        println(fid, "$nstr")
        for j in 1:d
            nstr2 = join(map(string, cwv.vectors.quantizer.codebooks[i].vectors[j,:]), " ")
            println(fid, "$nstr2")
        end
    end
    for i in 1:nrows  # rotation matrix
        nstr = join(map(string, cwv.vectors.quantizer.rot[:,i]), " ")
        println(fid, "$nstr")
    end
end

max_stars_repo_name: zgornel/EmbeddingsAnalysis

max_issues_repo_name: zgornel/EmbeddingsAnalysis

max_forks_repo_name: zgornel/EmbeddingsAnalysis


text:
# Test CachedElement
let
    a = CachedElement(1)
    @test get_element(a) == 1
    set_up_to_date!(a,false)
    @test !is_up_to_date(a)
    update_element!(a,2)
    @test is_up_to_date(a)
    @test get_element(a) == 2
    b = CachedElement(1)
    update_element!(b,4)
    @test time_stamp(b) > time_stamp(a)
    update_element!(a,2)
    @test time_stamp(b) < time_stamp(a)
end


text:
include("diffusion_mnist.jl")
using Images
using ProgressMeter
using DifferentialEquations
using Plots

"""
Helper function yielding the diffusion coefficient from a SDE.
"""
diffusion_coeff(t, sigma=convert(eltype(t), 25.0f0)) = sigma .^ t

"""
Helper function that produces images from a batch of images.
"""
function convert_to_image(x, y_size)
    Gray.(permutedims(vcat(reshape.(chunk(x |> cpu, y_size), 28, :)...), (2, 1)))
end

"""
Helper to make an animation from a batch of images.
"""
function convert_to_animation(x)
    frames = size(x)[end]
    batches = size(x)[end-1]
    animation = @animate for i = 1:frames+frames÷4
        if i <= frames
            heatmap(
                convert_to_image(x[:, :, :, :, i], batches),
                title="Iteration: $i out of $frames"
            )
        else
            heatmap(
                convert_to_image(x[:, :, :, :, end], batches),
                title="Iteration: $frames out of $frames"
            )
        end
    end
    return animation
end

"""
Helper function that generates inputs to a sampler.
"""
function setup_sampler(device, num_images=5, num_steps=500, ϵ=1.0f-3)
    t = ones(Float32, num_images) |> device
    init_x = (
        randn(Float32, (28, 28, 1, num_images)) .*
        expand_dims(marginal_prob_std(t), 3)
    ) |> device
    time_steps = LinRange(1.0f0, ϵ, num_steps)
    Δt = time_steps[1] - time_steps[2]
    return time_steps, Δt, init_x
end

"""
Sample from a diffusion model using the Euler-Maruyama method.

# References
https://yang-song.github.io/blog/2021/score/#how-to-solve-the-reverse-sde
"""
function Euler_Maruyama_sampler(model, init_x, time_steps, Δt)
    x = mean_x = init_x
    @showprogress "Euler-Maruyama Sampling" for time_step in time_steps
        batch_time_step = fill!(similar(init_x, size(init_x)[end]), 1) .* time_step
        g = diffusion_coeff(batch_time_step)
        mean_x = x .+ expand_dims(g, 3) .^ 2 .* model(x, batch_time_step) .* Δt
        x = mean_x .+ sqrt(Δt) .* expand_dims(g, 3) .* randn(Float32, size(x))
    end
    return mean_x
end

"""
Sample from a diffusion model using the Predictor-Corrector method.

# References
https://yang-song.github.io/blog/2021/score/#how-to-solve-the-reverse-sde
"""
function predictor_corrector_sampler(model, init_x, time_steps, Δt, snr=0.16f0)
    x = mean_x = init_x
    @showprogress "Predictor Corrector Sampling" for time_step in time_steps
        batch_time_step = fill!(similar(init_x, size(init_x)[end]), 1) .* time_step
        # Corrector step (Langevin MCMC)
        grad = model(x, batch_time_step)
        num_pixels = prod(size(grad)[1:end-1])
        grad_batch_vector = reshape(grad, (size(grad)[end], num_pixels))
        grad_norm = mean(sqrt, sum(abs2, grad_batch_vector, dims=2))
        noise_norm = Float32(sqrt(num_pixels))
        langevin_step_size = 2 * (snr * noise_norm / grad_norm)^2
        x += (
            langevin_step_size .* grad .+
            sqrt(2 * langevin_step_size) .* randn(Float32, size(x))
        )
        # Predictor step (Euler-Maruyama)
        g = diffusion_coeff(batch_time_step)
        mean_x = x .+ expand_dims((g .^ 2), 3) .* model(x, batch_time_step) .* Δt
        x = mean_x + sqrt.(expand_dims((g .^ 2), 3) .* Δt) .* randn(Float32, size(x))
    end
    return mean_x
end

"""
Helper to create a SDEProblem with DifferentialEquations.jl

# Notes
The reverse-time SDE is given by:
𝘥x = -σ²ᵗ 𝘚₀(𝙭, 𝘵)𝘥𝘵 + σᵗ𝘥𝘸
⟹ `f(u, p, t)` = -σ²ᵗ 𝘚₀(𝙭, 𝘵)
⟹ `g(u, p, t` = σᵗ
"""
function DifferentialEquations_problem(model, init_x, time_steps, Δt)
    function f(u, p, t)
        batch_time_step = fill!(similar(u, size(u)[end]), 1) .* t
        return (
            -expand_dims(diffusion_coeff(batch_time_step), 3) .^ 2 .*
            model(u, batch_time_step)
        )
    end

    function g(u, p, t)
        batch_time_step = fill!(similar(u), 1) .* t
        diffusion_coeff(batch_time_step)
    end
    tspan = (time_steps[begin], time_steps[end])
    SDEProblem(f, g, init_x, tspan), ODEProblem(f, init_x, tspan)
end

function plot_result(unet, args)
    args = Args(; args...)
    args.seed > 0 && Random.seed!(args.seed)
    device = args.cuda && CUDA.has_cuda() ? gpu : cpu
    unet = unet |> device
    time_steps, Δt, init_x = setup_sampler(device)

    # Euler-Maruyama
    euler_maruyama = Euler_Maruyama_sampler(unet, init_x, time_steps, Δt)
    sampled_noise = convert_to_image(init_x, size(init_x)[end])
    save(joinpath(args.save_path, "sampled_noise.jpeg"), sampled_noise)
    em_images = convert_to_image(euler_maruyama, size(euler_maruyama)[end])
    save(joinpath(args.save_path, "em_images.jpeg"), em_images)

    # Predictor Corrector
    pc = predictor_corrector_sampler(unet, init_x, time_steps, Δt)
    pc_images = convert_to_image(pc, size(pc)[end])
    save(joinpath(args.save_path, "pc_images.jpeg"), pc_images)

    # Setup an SDEProblem and ODEProblem to input to `solve()`.
    # Use dt=Δt to make the sample paths comparable to calculating "by hand".
    sde_problem, ode_problem = DifferentialEquations_problem(unet, init_x, time_steps, Δt)

    @info "Euler-Maruyama Sampling w/ DifferentialEquations.jl"
    diff_eq_em = solve(sde_problem, EM(), dt=Δt)
    diff_eq_em_end = diff_eq_em[:, :, :, :, end]
    diff_eq_em_images = convert_to_image(diff_eq_em_end, size(diff_eq_em_end)[end])
    save(joinpath(args.save_path, "diff_eq_em_images.jpeg"), diff_eq_em_images)
    diff_eq_em_animation = convert_to_animation(diff_eq_em)
    gif(diff_eq_em_animation, joinpath(args.save_path, "diff_eq_em.gif"), fps=50)
    em_plot = plot(diff_eq_em, title="Euler-Maruyama", legend=false, ylabel="x", la=0.25)
    plot!(time_steps, diffusion_coeff(time_steps), xflip=true, ls=:dash, lc=:red)
    plot!(time_steps, -diffusion_coeff(time_steps), xflip=true, ls=:dash, lc=:red)
    savefig(em_plot, joinpath(args.save_path, "diff_eq_em_plot.png"))

    @info "Probability Flow ODE Sampling w/ DifferentialEquations.jl"
    diff_eq_ode = solve(ode_problem, dt=Δt, adaptive=false)
    diff_eq_ode_end = diff_eq_ode[:, :, :, :, end]
    diff_eq_ode_images = convert_to_image(diff_eq_ode_end, size(diff_eq_ode_end)[end])
    save(joinpath(args.save_path, "diff_eq_ode_images.jpeg"), diff_eq_ode_images)
    diff_eq_ode_animation = convert_to_animation(diff_eq_ode)
    gif(diff_eq_ode_animation, joinpath(args.save_path, "diff_eq_ode.gif"), fps=50)
    ode_plot = plot(diff_eq_ode, title="Probability Flow ODE", legend=false, ylabel="x", la=0.25)
    plot!(time_steps, diffusion_coeff(time_steps), xflip=true, ls=:dash, lc=:red)
    plot!(time_steps, -diffusion_coeff(time_steps), xflip=true, ls=:dash, lc=:red)
    savefig(ode_plot, joinpath(args.save_path, "diff_eq_ode_plot.png"))
end

if abspath(PROGRAM_FILE) == @__FILE__
    ############################################################################
    # Issue loading function closures with BSON:
    # https://github.com/JuliaIO/BSON.jl/issues/69
    #
    BSON.@load "output/model.bson" unet args
    #
    # BSON.@load does not work if defined inside plot_result(⋅) because
    # it contains a function closure, GaussFourierProject(⋅), containing W.
    ###########################################################################
    plot_result(unet, args)
end


text:
# for profiling run with julia --track-allocation=user

using Profile
using GeometricIntegrators
using GeometricIntegrators.TestProblems.KuboOscillatorProblem
using GeometricIntegrators.TestProblems.KuboOscillatorProblem: Δt

const ntime = 10000

set_config(:nls_stol_break, 1E3)

function profile_kubo(tableau, sde, prof=false)
    int = Integrator(sde, tableau, Δt)
    sol = Solution(sde, Δt, ntime)

    integrate!(int, sol)
    set_initial_conditions!(sol, sde)

    @time integrate!(int, sol)

    if prof
        set_initial_conditions!(sol, sde)

        Profile.clear()
        Profile.clear_malloc_data()

        @profile integrate!(int, sol)
    end
end

# SERK
profile_kubo(getTableauBurrageE1(), kubo_oscillator_sde_1())
profile_kubo(getTableauBurrageE1(), kubo_oscillator_sde_2())
#
# SIRK
profile_kubo(getTableauStochasticGLRK(1), kubo_oscillator_sde_1())
profile_kubo(getTableauStochasticGLRK(1), kubo_oscillator_sde_2())
#
# WERK
profile_kubo(getTableauRosslerRS1(), kubo_oscillator_sde_1())
profile_kubo(getTableauRosslerRS1(), kubo_oscillator_sde_2())

# WIRK
profile_kubo(getTableauSRKw1(), kubo_oscillator_sde_1())
profile_kubo(getTableauSRKw1(), kubo_oscillator_sde_2())

# SIPRK
profile_kubo(getTableauStochasticStormerVerlet(), kubo_oscillator_psde_1())
profile_kubo(getTableauStochasticStormerVerlet(), kubo_oscillator_psde_2())

# SISPRK
profile_kubo(getTableauModifiedStochasticStormerVerlet(), kubo_oscillator_spsde_1())
profile_kubo(getTableauModifiedStochasticStormerVerlet(), kubo_oscillator_spsde_2())

# SERK
# profile_kubo(getTableauBurrageE1(), kubo_oscillator_sde_1(), true)
# profile_kubo(getTableauBurrageE1(), kubo_oscillator_sde_2(), true)

# SIRK
# profile_kubo(getTableauStochasticGLRK(1), kubo_oscillator_sde_1(), true)
# profile_kubo(getTableauStochasticGLRK(1), kubo_oscillator_sde_2(), true)

# WERK
# profile_kubo(getTableauRosslerRS1(), kubo_oscillator_sde_1(), true)
# profile_kubo(getTableauRosslerRS1(), kubo_oscillator_sde_2(), true)

# WIRK
# profile_kubo(getTableauSRKw1(), kubo_oscillator_sde_1(), true)
# profile_kubo(getTableauSRKw1(), kubo_oscillator_sde_2(), true)

# SIPRK
# profile_kubo(getTableauStochasticStormerVerlet(), kubo_oscillator_psde_1(), true)
# profile_kubo(getTableauStochasticStormerVerlet(), kubo_oscillator_psde_2(), true)

# SISPRK
# profile_kubo(getTableauModifiedStochasticStormerVerlet(), kubo_oscillator_spsde_1(), true)
# profile_kubo(getTableauModifiedStochasticStormerVerlet(), kubo_oscillator_spsde_2(), true)


text:
using RunLengthArrays
using Test

empty_array = RunLengthArray{Int,String}()
@test isempty(empty_array)

##################################################################################

original = Int8[3, 3, 3, 7, 7, 7, 7, 7, 7, 4]
compressed = RunLengthArray{Int,Int8}(original)

@test runs(compressed) == Int[3, 6, 1]
@test values(compressed) == Int8[3, 7, 4]

##################################################################################
# Check that looping works

total = Int8(0)
for elem in compressed
    global total
    total += elem
end
@test total == sum(original)

##################################################################################
# Check that the original array and the compressed array behave identically

@test length(compressed) == length(original)
@test sum(compressed) == sum(original)

for i in eachindex(original)
    @test compressed[i] == original[i]
end
@test compressed[end] == 4

@test collect(compressed) == original

##################################################################################
# Modify the two arrays and check that they are still identical

push!(compressed, Int8(4))
push!(original, 4)
@test length(compressed) == length(original)
@test sum(compressed) == sum(original)

push!(compressed, (4, Int8(6)))
append!(original, Int8[6, 6, 6, 6])
@test length(compressed) == length(original)
@test sum(compressed) == sum(original)

for i in eachindex(original)
    @test compressed[i] == original[i]
end

append!(compressed, Int8[3, 4, 5])
append!(original, Int8[3, 4, 5])

for i in eachindex(original)
    @test compressed[i] == original[i]
end

@test minimum(compressed) == minimum(original)
@test maximum(compressed) == maximum(original)
@test extrema(compressed) == extrema(original)

##################################################################################
# Finally, check that sorting works

@test collect(sort(compressed)) == sort(original)

sort!(compressed)
@test collect(compressed) == sort(original)


text:
#solves the EulerLagrane equations setup in DiscreteShape and DiscreteShapeBound.
function solve_EulerLagrange(ρ, I, iterations, p, q, R, λ, E, Z)
    #initilize common variables
    M, N = size(E)
    δp = zeros(Complex{Float64},axes(E))
    δq = zeros(Complex{Float64},axes(E))
    pq = zeros(Complex{Float64}, size(p))
    dRδp = zeros(Complex{Float64}, size(R))
    dRδq = zeros(Complex{Float64}, size(R))
    C = zeros(Complex{Float64}, size(R))
    w = centered(0.25*[0.0 1.0 0.0;1.0 0.0 1.0;0.0 1.0 0.0])
    wx, wy = setup_transform_values(M, N)
    for i = 1:iterations
        #calculate second derivatives of p and q
        δp = imfilter(p, reflect(w), "replicate")
        δq = imfilter(q, reflect(w), "replicate")

        for i in CartesianIndices(R)
            #calculate reflectance map
            R[i] = (ρ * (-I[1] * p[i] - I[2] * q[i] + I[3])) / sqrt(1 + p[i]^2
                + q[i]^2)
            #calculate partial derivatives of p and q in respect to p and q
            pq[i] = (1 + p[i]^2 + q[i]^2)
            dRδp[i] = ((-ρ * I[1]) / sqrt(pq[i])) + ((-I[1] * ρ) * p[i] - I[2]
                * ρ * q[i] + I[3] * ρ) * (-1 * p[i] * (pq[i]^(-3/2)))

            dRδq[i] = ((-ρ * I[2]) / sqrt(pq[i])) + (-I[1] * ρ * p[i] - I[2]
                * ρ * q[i] + I[3] * ρ) * (-1 * q[i] * (pq[i]^(-3/2)))
            #calculate new surface normals
            p[i] = δp[i] + (1 / (4 * λ)) * (E[i] - R[i]) * dRδp[i]
            q[i] = δq[i] + (1 / (4 * λ)) * (E[i] - R[i]) * dRδq[i]
        end

        #calculate Transform of surface normals
        fft!(p)
        fft!(q)

        for i in CartesianIndices(C)
            #compute C and prepare to update p and q
            C[i] = (-1im * wx[i] * p[i] + 1im * wy[i] * q[i]) / (2*π*(wx[i]^2 + wy[i]^2 + eps()))
            p[i] = 1im * wx[i] * C[i]
            q[i] = 1im * wy[i] * C[i]
        end

        #update p and q
        ifft!(p)
        ifft!(q)

        #recover Z by taking inverse transform of C
        ifft!(C)

        #take abs of C to get Z
        for i in CartesianIndices(Z)
            Z[i] = abs(C[i])
        end
    end
    return Z, p, q
end

function solve_EulerLagrange2(ρ, I, iterations, p, q, R, λ, E, Z)
    bound = (2:(first(size(R))-1),(2:last(size(R))-1))
    for i = 1:iterations
        pᵏ⁺¹ = copy(p)
        qᵏ⁺¹ = copy(q)
        for i in CartesianIndices(bound)
            #calculate reflectance map
            R[i] = (ρ * (-I[1] * p[i] - I[2] * q[i] + I[3])) / sqrt(1 + p[i]^2 + q[i]^2)
            pq = (1 + p[i]^2 + q[i]^2)
            dRδp = ((-ρ * I[1]) / sqrt(pq)) + ((-I[1] * ρ) * p[i] - I[2] * ρ * q[i] + I[3] * ρ) * (-1 * p[i] * (pq^(-3/2)))
            dRδq = ((-ρ * I[2]) / sqrt(pq)) + (-I[1] * ρ * p[i] - I[2] * ρ * q[i] + I[3] * ρ) * (-1 * q[i] * (pq^(-3/2)))
            pbar = (p[i[1]+1, i[2]] + p[i[1]-1, i[2]]) / 2
            qbar = (q[i[1], i[2]+1] + q[i[1], i[2]-1]) / 2
            pt = (p[i[1]+1, i[2]+1] + p[i[1]-1, i[2]-1] - p[i[1]+1, i[2]-1] - p[i[1]-1, i[2]+1]) / 4
            qt = (q[i[1]+1, i[2]+1] + q[i[1]-1, i[2]-1] - q[i[1]+1, i[2]-1] - q[i[1]-1, i[2]+1]) / 4
            pᵏ⁺¹[i] = pbar - qt/2 + (1 / (2 * λ)) * (E[i] - R[i]) * dRδq
            qᵏ⁺¹[i] = qbar - pt/2 + (1 / (2 * λ)) * (E[i] - R[i]) * dRδp
        end
        p = pᵏ⁺¹
        q = qᵏ⁺¹
    end
    return p, q
end

# WIP only works for odd dims
function setup_transform_values(N, M)
    wx = [j for i in 1:M, j in -0.5:(1/(N-1)):0.5]
    wy = [j for j in -0.5:(1/(M-1)):0.5, i in 1:N]
    wx = ifftshift(wx)
    wy = ifftshift(wy)
    return wx, wy
end

function setup_transform_values_pentland(M, N)
    wx = zeros(M, N)
    wy = zeros(M, N)
    for i in CartesianIndices(wx)
        wx[i] = (2 * π * i[2]) / M
        wy[i] = (2 * π * (N - i[1] + 1)) / N
    end
    return wx, wy
end

# Sets up xy grid as array
function setup_xy(dim)
    x = zeros(dim , dim)
    y = zeros(dim , dim)
    xyrange = ceil(Int, -dim/2):floor(Int, (dim/2))
    for i in CartesianIndices(x)
        x[i] = xyrange[i[2]]
        y[i] = -xyrange[i[1]]
    end
    return x,y
end


text:
using Documenter, Lerche

makedocs(sitename="Lerche documentation",
	  format = Documenter.HTML(
				   prettyurls = get(ENV,"CI",nothing) == "true"
				   ),
         pages = [
             "Overview and Guide" => "index.md",
	     "Grammar reference" => "grammar.md",
             "API" => "api.md"
             ],
         #doctest = :fix
	  )

deploydocs(
    repo = "github.com/jamesrhester/Lerche.jl.git",
)


text:
using Test
import GeneralizedSVD:GenSVD, gsvd
using LinearAlgebra
using DataFrames
using CSV

"""
    test(m, p, n)

This function tests the time performance
and backward stability of Generalized Singular Value Decomposition (GSVD)

First, we generate a m by n random matrix A and p by n random matrix B,
Create a copy of A and B, respectively (original A and B will be overwritten)

Then, we call gsvd function, return U, V, Q, C, S, R, k, l.
Calculate and print time elapsed for computing GSVD.

Compute:

      res_a = ||U' * A_ * Q - C * R||/(max(m,n)* ||A_|| * eps)

      res_b = ||V' * B_ * Q - S * R||/(max(p,n)* ||B_|| * eps)

      orthog_u = ||U' * U - I||/(m * eps)

      orthog_v = ||V' * V - I||/(p * eps)

      orthog_q = ||Q' * Q - I||/(n * eps)

# Arguments
- m: num of rows in A
- p: num of rows in B
- n: num of cols in A, B

no return
"""
function computemetric(F::GenSVD, A, B)
    e = eps(Float64)
    m, n = size(A)
    p = size(B)[1]
    res_a = opnorm(F.U'*A*F.Q - F.C*F.R, 1)/(max(m,n)*opnorm(A, 1)*e)
    res_b = opnorm(F.V'*B*F.Q - F.S*F.R, 1)/(max(p,n)*opnorm(B, 1)*e)
    orthog_cs = opnorm(I - F.C'*F.C - F.S'*F.S, 1)/((F.k+F.l)*e)
    orthog_u = opnorm(I - F.U'*F.U, 1)/(m*e)
    orthog_v = opnorm(I - F.V'*F.V, 1)/(p*e)
    orthog_q = opnorm(I - F.Q'*F.Q, 1)/(n*e)
    return res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q
end

function computemetric(F::GeneralizedSVD, A, B)
    e = eps(Float64)
    m, n = size(A)
    p = size(B)[1]
    res_a = opnorm(F.U'*A*F.Q - F.D1*F.R0, 1)/(max(m,n)*opnorm(A, 1)*e)
    res_b = opnorm(F.V'*B*F.Q - F.D2*F.R0, 1)/(max(p,n)*opnorm(B, 1)*e)
    orthog_cs = opnorm(I - F.D1'*F.D1 - F.D2'*F.D2, 1)/((F.k+F.l)*e)
    orthog_u = opnorm(I - F.U'*F.U, 1)/(m*e)
    orthog_v = opnorm(I - F.V'*F.V, 1)/(p*e)
    orthog_q = opnorm(I - F.Q'*F.Q, 1)/(n*e)
    return res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q
end

function test(m, p, n)
    A = randn(m, n)
    B = randn(p, n)
    F = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F, A, B)
    return F.k, F.l, res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q
end

function test()
    df = DataFrame(m = Int64[], p = Int64[], n = Int64[], k = Int64[], l = Int64[], res_a = Float64[],
    res_b = Float64[], orthog_cs = Float64[], orthog_u = Float64[],
    orthog_v = Float64[], orthog_q = Float64[])

    # case 1: m >= n & p >= n
    println("=======================Case 1============================")
    for i = 1:20
        a, b, c, d, e, f, g, h = test(60, 50, 40)
        push!(df, [60 50 40 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(300, 250, 200)
        push!(df, [300 250 200 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(900, 750, 600)
        push!(df, [900 750 600 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(1500, 1250, 1000)
        push!(df, [1500 1250 1000 a b c d e f g h])
        println("---------------------------------------------------------------")
    end

    # case 2: m >= n > p
    println("=======================Case 2================================")
    for i = 1:20
        a, b, c, d, e, f, g, h = test(60, 40, 50)
        push!(df, [60 40 50 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(300, 200, 250)
        push!(df, [300 200 250 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(900, 600, 750)
        push!(df, [900 600 750 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(1500, 1000, 1250)
        push!(df, [1500 1000 1250 a b c d e f g h])
        println("---------------------------------------------------------------")
    end

    # case 3: p >= n > m
    println("=======================Case 3==============================")
    for i = 1:20
        a, b, c, d, e, f, g, h = test(40, 60, 50)
        push!(df, [40 60 50 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(200, 300, 250)
        push!(df, [200 300 250 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(600, 900, 750)
        push!(df, [600 900 750 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(1000, 1500, 1250)
        push!(df, [1000 1500 1250 a b c d e f g h])
        println("---------------------------------------------------------------")
    end

    # case 4: n > m & n > p
    println("=====================Case 4============================")
    for i = 1:20
        a, b, c, d, e, f, g, h = test(20, 30, 60)
        push!(df, [20 30 60 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(200, 300, 600)
        push!(df, [200 300 600 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(400, 600, 1200)
        push!(df, [400 600 1200 a b c d e f g h])
        println("---------------------------------------------------------------")
    end
    for i = 1:20
        a, b, c, d, e, f, g, h = test(1000, 1500, 3000)
        push!(df, [1000 1500 3000 a b c d e f g h])
        println("---------------------------------------------------------------")
    end

    CSV.write("gsvd_rand_matrix_test.csv", df, append=true)
    return df
end

# Example 1 from Ch.4 of Jenny's Thesis
function test1()
    A  = Float64[1 2 1 0; 2 3 1 1; 3 4 1 2;4 5 1 3;5 6 1 4]
    B = Float64[6 7 1 5;7 1 -6 13;-4 8 9 -2]
    @time F1 = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F1, A, B)
    println("--------Residual of gsvd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    println([A; B] ≈ [F1.U*F1.C; F1.V*F1.S]*F1.R*F1.Q')

    F2 = svd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F2, A, B)
    println("--------Residual of svd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    return F1, F2
end

# Example 2 from Ch.4 of Jenny's Thesis
# It turns out this example is wrong

# This thus is used as Example 2 in working notes
function test2()
    A = [1. 2 1 0;2 3 1 1;3 4 1 2]
    B = [4. 5 1 3;5 6 1 4;6 7 1 5;7 1 -6 13]
    F1 = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F1, A, B)
    println("--------Residual of gsvd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    println([A; B] ≈ [F1.U*F1.C; F1.V*F1.S]*F1.R*F1.Q')

    F2 = svd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F2, A, B)
    println("--------Residual of svd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    return F1, F2
end

function test2_wrapper()
    A = [1. 2 1 0;2 3 1 1;3 4 1 2];
    B = [4. 5 1 3;5 6 1 4;6 7 1 5;7 1 -6 13];
    return U, V, Q, D1, D2, k, l, R = gsvdwrapper(A, B)
end

# Example 1 in working notes
function test3()
    A = [1. 2 3 0; 5 4 2 1; 0 3 5 2; 2 1 3 3; 2 0 5 3];
    B = [1. 0 3 -1; -2 5 0 1; 4 2 -1 2];
    F1 = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F1, A, B)
    println("--------Residual of gsvd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    println([A; B] ≈ [F1.U*F1.C; F1.V*F1.S]*F1.R*F1.Q')

    F2 = svd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F2, A, B)
    println("--------Residual of svd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    return F1, F2
end

function test3_wrapper()
    A = [1. 2 3 0; 5 4 2 1; 0 3 5 2; 2 1 3 3; 2 0 5 3];
    B = [1. 0 3 -1; -2 5 0 1; 4 2 -1 2];
    return U, V, Q, D1, D2, k, l, R = gsvdwrapper(A, B)
end

# Example 4 in working notes
function test4()
    A = [1. 4 2 3 0;3 4 0 -2 1;4 7 5 6 3];
    B = [1. 4 2 3 0;2 5 3 4 1;3 6 4 5 2;0 1 -1 3 1];
    F1 = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F1, A, B)
    println("--------Residual of gsvd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    println([A; B] ≈ [F1.U*F1.C; F1.V*F1.S]*F1.R*F1.Q')

    F2 = svd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F2, A, B)
    println("--------Residual of svd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    return F1, F2
end

# Early exit after preprocessing
function test5()
    # rank(A) = 3, rank(B) = 1, rank([A;B]) = 4
    A = [1. 2 3 4;3 5 2 1;4 6 2 1];
    B = [0 0 3. 0;0 0 0 0;0 0 2. 0; 0 0 0 0];
    F1 = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F1, A, B)
    println("--------Residual of gsvd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    println([A; B] ≈ [F1.U*F1.C; F1.V*F1.S]*F1.R*F1.Q')

    F2 = svd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F2, A, B)
    println("--------Residual of svd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    return F1, F2
end

# Early exit in CSD
function test6()
    # rank(A) = 4, rank(B) = 1, rank([A;B]) = 4
    A = [1.0 2 3 4;2 5 2 1;4 2 3 1;2 4 5 1];
    B = [1. 1 0 0;0 0 0 0; 0 0 0 0; 0 0 0 0];
    F1 = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F1, A, B)
    println("--------Residual of gsvd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    println([A; B] ≈ [F1.U*F1.C; F1.V*F1.S]*F1.R*F1.Q')

    F2 = svd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F2, A, B)
    println("--------Residual of svd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    return F1, F2
end

# Example in the document GSVD.jl VS. gsvd.m
function test7()
    U = [1/2 0 sqrt(3)/2 0;
        0 -sqrt(3)/2 0 1/2;
        sqrt(3)/2 0 -1/2 0;
        0 1/2 0 sqrt(3)/2]
    V = [1/2 0 sqrt(3)/2 0;
        0 -sqrt(3)/2 0 1/2;
        sqrt(3)/2 0 -1/2 0;
        0 1/2 0 sqrt(3)/2]
    Q = [1/2 0 sqrt(3)/2 0;
        0 -sqrt(3)/2 0 1/2;
        sqrt(3)/2 0 -1/2 0;
        0 1/2 0 sqrt(3)/2]
    R = [0 1 0 0;0 0 1 0;0 0 0 1]
    C = [sqrt(3)/2 0 0;0 1/2 0;0 0 0;0 0 0]
    S = [1/2 0 0;0 sqrt(3)/2 0;0 0 1;0 0 0]
    A = U*C*R*Q'
    B = V*S*R*Q'
    F1 = gsvd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F1, A, B)
    println("--------Residual of gsvd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    println([A; B] ≈ [F1.U*F1.C; F1.V*F1.S]*F1.R*F1.Q')

    F2 = svd(A, B)
    res_a, res_b, orthog_cs, orthog_u, orthog_v, orthog_q = computemetric(F2, A, B)
    println("--------Residual of svd---------")
    println("res_a: ", res_a)
    println("res_b: ", res_b)
    println("orthog_cs: ", orthog_cs)
    println("orthog_u: ", orthog_u)
    println("orthog_v: ", orthog_v)
    println("orthog_q: ", orthog_q)
    return F1, F2
end

# function profiling_function(m, p, n)
#     maxI = 50
#     t_pre = 0.0
#     t_qr = 0.0
#     t_csd = 0.0
#     t_all = 0.0
#     for i = 1:maxI
#         println("---$i-th run---")
#         A = randn(m, n)
#         B = randn(p, n)
#         ans = @timed gsvd(A, B)
#         t_pre = t_pre + ans[1][1]
#         t_qr = t_qr + ans[1][2]
#         t_csd = t_csd + ans[1][3]
#         t_all = t_all + ans[2]
#     end
#     t_pre = t_pre/maxI
#     t_qr = t_qr/maxI
#     t_csd = t_csd/maxI
#     t_all = t_all/maxI
#     p_pre = 100*t_pre/t_all
#     p_qr = 100*t_qr/t_all
#     p_csd = 100*t_csd/t_all
#     return m, p, n, t_pre, p_pre, t_qr, p_qr, t_csd, p_csd, t_all
# end
#
# function test_profiling()
#     df = DataFrame(m = Int64[], p = Int64[], n = Int64[],
#     t_pre = Float64[], p_pre = Float64[],
#     t_qr = Float64[], p_qr = Float64[],
#     t_csd = Float64[], p_csd = Float64[],
#     t_all = Float64[])
#
#     # case 1: m >= n & p >= n
#     println("=======================Case 1============================")
#     t1a = profiling_function(1500, 1200, 1000)
#     push!(df, [t1a[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t1b = profiling_function(500, 500, 500)
#     push!(df, [t1b[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t1c = profiling_function(650, 310, 230)
#     push!(df, [t1c[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t1d = profiling_function(430, 610, 210)
#     push!(df, [t1d[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#
#     # case 2: m >= n > p
#     println("=======================Case 2================================")
#     t2a = profiling_function(1500, 1000, 1200)
#     push!(df, [t2a[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t2b = profiling_function(720, 220, 540)
#     push!(df, [t2b[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t2c = profiling_function(440, 180, 440)
#     push!(df, [t2c[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t2d = profiling_function(370, 290, 350)
#     push!(df, [t2d[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#
#     # case 3: p >= n > m
#     println("=======================Case 3==============================")
#     # Catch LAPACK Exception: DBDSDC did not converge
#     t3a = profiling_function(1000, 1500, 1200)
#     push!(df, [t3a[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t3b = profiling_function(250, 300, 300)
#     push!(df, [t3b[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t3c = profiling_function(360, 660, 600)
#     push!(df, [t3c[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t3d = profiling_function(130, 520, 480)
#     push!(df, [t3d[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#
#     # case 4: n > m & n > p
#     println("=====================Case 4============================")
#     # Catch LAPACK Exception: DBDSDC did not converge
#     # t4a = profiling_function(1000, 1200, 1500)
#     t4a = profiling_function(1300, 1200, 1500)
#     push!(df, [t4a[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     t4b = profiling_function(260, 600, 770)
#     push!(df, [t4b[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     # Might early exit
#     # t4c = profiling_function(370, 250, 700)
#     t4c = profiling_function(370, 500, 700)
#     push!(df, [t4c[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#     # Might early exit
#     # t4d = profiling_function(120, 120, 400)
#     t4d = profiling_function(220, 220, 400)
#     push!(df, [t4d[i] for i in 1:10])
#     println("---------------------------------------------------------------")
#
#     CSV.write("LAPACK-JuliaX.X-profiling.csv", df, append=true)
#     return df
# end
#
# function wrapper_profiling_function(m, p, n)
#     maxI = 20
#     t_dggsvp3 = 0.0
#     t_dtgsja = 0.0
#     t_all = 0.0
#     for i = 1:maxI
#         println("---$i-th run---")
#         A = randn(m, n)
#         B = randn(p, n)
#         ans = @timed gsvdwrapper(A, B)
#         t_dggsvp3 = t_dggsvp3 + ans[1][1]
#         t_dtgsja = t_dtgsja + ans[1][2]
#         t_all = t_all + ans[2]
#     end
#     t_dggsvp3 = t_dggsvp3/maxI
#     t_dtgsja = t_dtgsja/maxI
#     t_all = t_all/maxI
#     p_dggsvp3 = 100*t_dggsvp3/t_all
#     p_dtgsja = 100*t_dtgsja/t_all
#     return m, p, n, t_dggsvp3, p_dggsvp3, t_dtgsja, p_dtgsja, t_all
# end
#
# function test_wrapper_profiling()
#     df = DataFrame(m = Int64[], p = Int64[], n = Int64[],
#     t_dggsvp3 = Float64[], p_dggsvp3 = Float64[],
#     t_dtgsja = Float64[], p_dtgsja = Float64[],
#     t_all = Float64[])
#
#     # case 1: m >= n & p >= n
#     println("=======================Case 1============================")
#     t1a = wrapper_profiling_function(1500, 1200, 1000)
#     push!(df, [t1a[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t1b = wrapper_profiling_function(500, 500, 500)
#     push!(df, [t1b[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t1c = wrapper_profiling_function(650, 310, 230)
#     push!(df, [t1c[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t1d = wrapper_profiling_function(430, 610, 210)
#     push!(df, [t1d[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#
#     # case 2: m >= n > p
#     println("=======================Case 2================================")
#     t2a = wrapper_profiling_function(1500, 1000, 1200)
#     push!(df, [t2a[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t2b = wrapper_profiling_function(720, 220, 540)
#     push!(df, [t2b[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t2c = wrapper_profiling_function(440, 180, 440)
#     push!(df, [t2c[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t2d = wrapper_profiling_function(370, 290, 350)
#     push!(df, [t2d[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#
#     # case 3: p >= n > m
#     println("=======================Case 3==============================")
#     t3a = wrapper_profiling_function(1000, 1500, 1200)
#     push!(df, [t3a[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t3b = wrapper_profiling_function(250, 300, 300)
#     push!(df, [t3b[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t3c = wrapper_profiling_function(360, 660, 600)
#     push!(df, [t3c[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t3d = wrapper_profiling_function(130, 520, 480)
#     push!(df, [t3d[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#
#     # case 4: n > m & n > p
#     println("=====================Case 4============================")
#     t4a = wrapper_profiling_function(1300, 1200, 1500)
#     push!(df, [t4a[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t4b = wrapper_profiling_function(260, 600, 770)
#     push!(df, [t4b[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t4c = wrapper_profiling_function(370, 500, 700)
#     push!(df, [t4c[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#     t4d = wrapper_profiling_function(220, 220, 400)
#     push!(df, [t4d[i] for i in 1:8])
#     println("---------------------------------------------------------------")
#
#     CSV.write("Julia1.3-profiling.csv", df, append=true)
#     return df
# end


text:
## Predictive Information Criteria

#function log_sum_exp(x)
#    xmax = maximum(x)
#    xmax + log(sum(exp.(x .- xmax)))
#end

AIC(x::T, data::Data) where T <: AbstractFit = AIC(merge(x.chains...), data)

## Akaike Information Criterion (AIC)
function AIC(x::T, data::Data) where T <: Chain
    map = MAP(x, data)
    N, q = size(x.theta[1].b)
    ## number of parameters without random effects and latent class indicators
    k = length(vectorize(map)) - N*q - N
    -2*(loglik(map, data) - k)
end

function mean_param(x::Fit{T}, data::Data) where T <: AbstractParam
    par = T(data, x.theta[1].G)
    for prop in fieldnames(typeof(par))
        y = mean([getproperty(w, prop) for w in x.theta])

        ## check for Int to avoid InexactError
        if eltype(getproperty(par, prop)) <: Int
            setproperty!(par, prop, Int.(round.(y, digits=0)))
        else
            setproperty!(par, prop, y)
        end
    end
    par
end

BIC(x::T, data::Data) where T <: AbstractFit = BIC(merge(x.chains...), data)

## Bayesian Information Criterion (BIC)
function BIC(x::T, data::Data) where T <: Chain
    post_mean = mean_param(x, data)

    N, q = size(x.theta[1].b)
    ## number of parameters without random effects and latent class indicators
    k = length(vectorize(post_mean)) - N*q - N

    -2*loglik(post_mean, data) + k*log(N)
end

DIC(x::T, data::Data) where T <: AbstractFit = DIC(merge(x.chain...), data)

## Deviance Information Criterion (DIC)
function DIC(x::T, data::Data) where T <: Chain
    post_mean = mean_param(x, data)
    lppd = loglik(post_mean, data)
    pDIC = 2*(lppd - mean(map(par -> loglik(par, data), x.theta)))
    -2*(lppd - pDIC)
end

WAIC(x::T, data::Data) where T <: AbstractFit = WAIC(merge(x.chain...), data)

## Watanabe Akaike Information Criterion (WAIC)
function WAIC(x::T, data::Data) where T <: Chain
    N = length(data)
    pd = zeros(Float64, N)
    lpd = zeros(Float64, N)
    lpd_sq = zeros(Float64, N)

    M = n_samples(x)
    lp = get_loglik_p(x) ##dropdims(loglik_array(x, data), dims=(2))

    lppd = -N*log(M)
    for c in eachcol(lp) ##i in 1:N
        lppd += logsumexp(c) ##lp[:,i])
    end

    pWAIC1 = 2*(lppd - sum(mean(lp, dims=1)))
    pWAIC2 = sum(var(lp, dims=1))

    -2*(lppd - pWAIC1), -2*(lppd - pWAIC2)
end

function loo(x::T, data::Data; cores=1) where T <: AbstractFit
    ll = loglik_array(x, data)

    R"library(loo)"
    @rput ll
    @rput cores
    R"reff <- relative_eff(ll)"
    R"res <- loo(ll, r_eff=reff, cores=cores)"
    R"print(res)"
    R"w <- waic(ll)"
    R"print(w)"
end