text:
"""
    LSAModel{S<:AbstractString, T<:AbstractFloat, A<:AbstractMatrix{T}, H<:Integer}

LSA (latent semantic analysis) model. It constructs from a document term matrix (dtm)
a model that can be used to embed documents in a latent semantic space pertaining to
the data. The model requires that the document term matrix be a
`DocumentTermMatrix{T<:AbstractFloat}` because the elements of the matrices resulted
from the SVD operation are floating point numbers and these have to match or be
convertible to type `T`.

# Fields
  * `vocab::Vector{S}` a vector with all the words in the corpus
  * `vocab_hash::OrderedDict{S,H}` a word to index in word embeddings matrix mapping
  * `Σinv::A` diagonal of the inverse singular value matrix
  * `Uᵀ::A` transpose of the word embedding matrix
  * `stats::Symbol` the statistical measure to use for word importances in documents. Available values are: `:count` (term count), `:tf` (term frequency), `:tfidf` (default, term frequency-inverse document frequency) and `:bm25` (Okapi BM25)
  * `idf::Vector{T}` inverse document frequencies for the words in the vocabulary
  * `nwords::T` averge number of words in a document
  * `ngram_complexity::Int` ngram complexity
  * `κ::Int` the `κ` parameter of the BM25 statistic
  * `β::Float64` the `β` parameter of the BM25 statistic
  * `tol::T` minimum size of the vector components (default `T(1e-15)`)

# SVD matrices `U`, `Σinv` and `V`:
  If `X` is a `m`×`n` document-term-matrix with `n` documents and `m` words so that
`X[i,j]` represents a statistical indicator of the importance of term `i` in document `j`
then:
  * `U, Σ, V = svd(X)`
  * `Σinv = diag(inv(Σ))`
  * `Uᵀ = U'`
  * `X ≈ U * Σ * V'`
  The matrix `V` of document embeddings is not actually stored in the model.

# Examples
```
julia> using StringAnalysis

       doc1 = StringDocument("This is a text about an apple. There are many texts about apples.")
       doc2 = StringDocument("Pears and apples are good but not exotic. An apple a day keeps the doctor away.")
       doc3 = StringDocument("Fruits are good for you.")
       doc4 = StringDocument("This phrase has nothing to do with the others...")
       doc5 = StringDocument("Simple text, little info inside")

       crps = Corpus(AbstractDocument[doc1, doc2, doc3, doc4, doc5])
       prepare!(crps, strip_punctuation)
       update_lexicon!(crps)
       dtm = DocumentTermMatrix{Float32}(crps, collect(keys(crps.lexicon)))

       ### Build LSA Model ###
       lsa_model = LSAModel(dtm, k=3, stats=:tf)

       query = StringDocument("Apples and an exotic fruit.")
       idxs, corrs = cosine(lsa_model, crps, query)

       println("Query: \"\$(query.text)\"")
       for (idx, corr) in zip(idxs, corrs)
           println("\$corr -> \"\$(crps[idx].text)\"")
       end
Query: "Apples and an exotic fruit."
0.9746108 -> "Pears and apples are good but not exotic  An apple a day keeps the doctor away "
0.870703 -> "This is a text about an apple  There are many texts about apples "
0.7122063 -> "Fruits are good for you "
0.22725986 -> "This phrase has nothing to do with the others "
0.076901935 -> "Simple text  little info inside "
```

# References:
  * [The LSA wiki page](https://en.wikipedia.org/wiki/Latent_semantic_analysis)
  * [Deerwester et al. 1990](http://lsa.colorado.edu/papers/JASIS.lsi.90.pdf)

"""
struct LSAModel{S<:AbstractString, T<:AbstractFloat, A<:AbstractMatrix{T}, H<:Integer}
    vocab::Vector{S}        # vocabulary
    vocab_hash::OrderedDict{S,H}   # term to column index in U
    Σinv::Vector{T}         # diagonal of inverse of Σ
    Uᵀ::A                   # word vectors (transpose of U)
    stats::Symbol           # term/document importance
    idf::Vector{T}          # inverse document frequencies
    nwords::T               # average words/document in corpus
    ngram_complexity::Int   # ngram complexity
    κ::Int                  # κ parameter for Okapi BM25 (used if stats==:bm25)
    β::Float64              # β parameter for Okapi BM25 (used if stats==:bm25)
    tol::T                  # Minimum size of vector elements
end

function LSAModel(dtm::DocumentTermMatrix{T}; kwargs...) where T<:Integer
    throw(ErrorException(
        """A LSA model requires a that the document term matrix
        be a DocumentTermMatrix{<:AbstractFloat}!"""))
end

function LSAModel(dtm::DocumentTermMatrix{T};
                  k::Int=size(dtm.dtm, 2),
                  stats::Symbol=:tfidf,
                  ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
                  tol::T=T(1e-15),
                  κ::Int=BM25_KAPPA,
                  β::Float64=BM25_BETA
                 ) where T<:AbstractFloat
    p, n = size(dtm.dtm)
    k = clamp(k, 1, min(k, n))
    zeroval = zero(T)
    minval = T(tol)
    # Checks
    length(dtm.terms) != p &&
        throw(DimensionMismatch("Dimensions inside dtm are inconsistent."))
    if !(stats in [:count, :tf, :tfidf, :bm25])
        @warn "stats has to be either :tf, :tfidf or :bm25; defaulting to :tfidf"
        stats = :tfidf
    end
    # Calculate inverse document frequency, mean document size
    documents_containing_term = vec(sum(dtm.dtm .> 0, dims=2)) .+ one(T)
    idf = log.(n ./ documents_containing_term) .+ one(T)
    nwords = mean(sum(dtm.dtm, dims=1))
    # Get X
    if stats == :count
        X = dtm.dtm
    elseif stats == :tf
        X = tf(dtm.dtm)
    elseif stats == :tfidf
        X = tf_idf(dtm.dtm)
    elseif stats == :bm25
        X = bm_25(dtm.dtm, κ=κ, β=β)
    end
    # Decompose document-word statistic
    local U₀::Matrix{T}, Σ::Vector{T}
    try
        U₀, Σ, _ = tsvd(X, k)
    catch
        # Try regular svd as a
        # last resort
        U₀, Σ, _ = svd(Matrix(X))
        U₀ = U₀[:, 1:k]
        Σ = Σ[1:k]
    end
    # Build model components
    Σinv = 1 ./ Σ
    Σinv[abs.(Σinv) .< minval] .= zeroval
    U₀[abs.(U₀).< minval] .= zeroval
    U = Matrix{T}(U₀')
    # Return the model
    return LSAModel(dtm.terms, dtm.row_indices,
                    Σinv, U, stats, idf, nwords,
                    ngram_complexity, κ, β, minval)
end

function Base.show(io::IO, lm::LSAModel{S,T,A,H}) where {S,T,A,H}
    len_vecs, num_terms = size(lm.Uᵀ)
    print(io, "LSA Model ($(lm.stats)), $(num_terms) terms, " *
          "dimensionality $(len_vecs), $(T) vectors")
end

"""
    lsa(X [;k=<num documents>, stats=:tfidf, κ=2, β=0.75, tol=1e-15])

Constructs a LSA model. The input `X` can be a `Corpus` or a `DocumentTermMatrix`.
Use `?LSAModel` for more details. Vector components smaller than `tol` will be
zeroed out.
"""
function lsa(dtm::DocumentTermMatrix{T};
             k::Int=size(dtm.dtm, 2),
             stats::Symbol=:tfidf,
             ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
             tol::T=T(1e-15),
             κ::Int=BM25_KAPPA,
             β::Float64=BM25_BETA) where T<:AbstractFloat
    LSAModel(dtm, k=k, stats=stats, ngram_complexity=ngram_complexity, κ=κ, β=β, tol=tol)
end

function lsa(crps::Corpus,
             ::Type{T} = DEFAULT_FLOAT_TYPE;
             k::Int=length(crps),
             stats::Symbol=:tfidf,
             ngram_complexity::Int=DEFAULT_NGRAM_COMPLEXITY,
             tol::T=T(1e-15),
             κ::Int=BM25_KAPPA,
             β::Float64=BM25_BETA) where T<:AbstractFloat
    lex = create_lexicon(crps, ngram_complexity)
    lsa(DocumentTermMatrix{T}(crps, lex, ngram_complexity=ngram_complexity),
        k=k, stats=stats, ngram_complexity=ngram_complexity,
        κ=κ, β=β, tol=tol)
end

"""
    vocabulary(lm)

Return the vocabulary as a vector of words of the LSA model `lm`.
"""
vocabulary(lm::LSAModel) = lm.vocab

"""
    in_vocabulary(lm, word)

Return `true` if `word` is part of the vocabulary of the LSA model `lm` and
`false` otherwise.
"""
in_vocabulary(lm::LSAModel, word::AbstractString) = word in lm.vocab

"""
    size(lm)

Return a tuple containin input and output dimensionalities of the LSA model `lm`.
"""
size(lm::LSAModel) = size(lm.Uᵀ,2), length(lm.Σinv)

"""
    index(lm, word)

Return the index of `word` from the LSA model `lm`.
"""
index(lm::LSAModel, word) = lm.vocab_hash[word]

"""
    get_vector(lm, word)

Returns the vector representation of `word` from the LSA model `lm`.
"""
function get_vector(lm::LSAModel{S,T,A,H}, word) where {S,T,A,H}
    default = zeros(T, length(lm.Σinv))
    idx = get(lm.vocab_hash, word, 0)
    if idx == 0
        return default
    else
        return lm.Uᵀ[:, idx]
    end
end

"""
    embed_document(lm, doc)

Return the vector representation of `doc`, obtained using the LSA model `lm`.
`doc` can be an `AbstractDocument`, `Corpus` or DTV or DTM.
"""
embed_document(lm::LSAModel{S,T,A,H}, doc::AbstractDocument) where {S,T,A,H} =
    # Hijack vocabulary hash to use as lexicon (only the keys needed)
    embed_document(lm, dtv(doc, lm.vocab_hash, T,
                           ngram_complexity=lm.ngram_complexity,
                           lex_is_row_indices=true))

embed_document(lm::LSAModel{S,T,A,H}, doc::AbstractString) where {S,T,A,H} =
    embed_document(lm, NGramDocument{S}(doc, DocumentMetadata(), lm.ngram_complexity))

embed_document(lm::LSAModel{S,T,A,H}, doc::Vector{S2}) where {S,T,A,H,S2<:AbstractString} =
    embed_document(lm, TokenDocument{S}(doc))

# Actual embedding function: takes as input the LSA model `lm` and a document
# term vector `dtv`. Returns the representation of `dtv` in the embedding space.
function embed_document(lm::LSAModel{S,T,A,H}, dtv::AbstractVector{T}) where {S,T,A,H}
    words_in_document = sum(dtv)
    # Calculate document vector
    if lm.stats == :count
        v = dtv
    elseif lm.stats == :tf
        v = sqrt.(dtv ./ max(words_in_document, one(T)))
    elseif lm.stats == :tfidf
        v = sqrt.(dtv ./ max(words_in_document, one(T))) .* lm.idf
    elseif lm.stats == :bm25
        k = T(lm.κ)
        b = T(lm.β)
        tf = sqrt.(dtv ./ max(words_in_document, one(T)))
        v = lm.idf .* (k + 1) .* tf ./
                      (k * (one(T) - b + b * words_in_document/lm.nwords) .+ tf)
    end
    # Embed
    d̂ = (lm.Σinv .* lm.Uᵀ) * v      # embed
    d̂ = d̂ ./ (norm(d̂,2) .+ eps(T))  # normalize
    d̂[abs.(d̂) .< lm.tol] .= zero(T) # zero small elements
    return d̂
end

function embed_document(lm::LSAModel{S,T,A,H}, dtm::DocumentTermMatrix{T}) where {S,T,A,H}
    n = size(dtm.dtm,1)
    k = size(lm.Uᵀ, 1)
    if lm.stats == :count
        X = dtm.dtm
    elseif lm.stats == :tf
        X = tf(dtm)
    elseif lm.stats == :tfidf
        X = tf_idf(dtm)
    elseif lm.stats == :bm25
        X = bm_25(dtm, κ=lm.κ, β=lm.β)
    end
    V = (lm.Σinv .* lm.Uᵀ) * X
    V ./= (sqrt.(sum(V.^2, dims=1)) .+ eps(T))
    V[abs.(V) .< lm.tol] .= zero(T)
    return V
end

function embed_document(lm::LSAModel{S,T,A,H}, crps::Corpus) where {S,T,A,H}
    lex = create_lexicon(crps, lm.ngram_complexity)
    embed_document(lm, DocumentTermMatrix{T}(crps, lex, ngram_complexity=lm.ngram_complexity))
end

"""
    embed_word(lm, word)

Return the vector representation of `word` using the LSA model `lm`.
"""
function embed_word(lm::LSAModel, word)
    throw(ErrorException(
        """Word embedding is not supported as it would require storing
        all documents in the model in order to determine the counts
        of `word` across the corpus."""))
end

"""
    cosine(model, docs, doc, n=10)

Return the positions of the `n` closest neighboring documents to `doc`
found in `docs`. `docs` can be a corpus or document term matrix.
The vector representations of `docs` and `doc` are obtained with the
`model` which can be either a `LSAModel` or `RPModel`.
"""
function cosine(model, docs, doc, n=10)
    metrics = embed_document(model, docs)' * embed_document(model, doc)
    n = min(n, length(metrics))
    topn_positions = sortperm(metrics[:], rev = true)[1:n]
    topn_metrics = metrics[topn_positions]
    return topn_positions, topn_metrics
end

"""
    similarity(model, doc1, doc2)

Return the cosine similarity value between two documents `doc1` and `doc2`
whose vector representations have been obtained using the `model`,
which can be either a `LSAModel` or `RPModel`.
"""
function similarity(model, doc1, doc2)
    return embed_document(model, doc1)' * embed_document(model, doc2)
end

"""
    save(lm, filename)

Saves an LSA model `lm` to disc in file `filename`.
"""
function save_lsa_model(lm::LSAModel{S,T,A,H}, filename::AbstractString) where {S,T,A,H}
    k, nwords = size(lm.Uᵀ)
    open(filename, "w") do fid
        println(fid, "LSA Model saved at $(Dates.now())")
        println(fid, "$nwords $k")  # number of documents, words, k
        println(fid, lm.stats)
        writedlm(fid, lm.idf', " ")
        println(fid, lm.nwords)
        println(fid, lm.ngram_complexity)
        println(fid, lm.κ)
        println(fid, lm.β)
        println(fid, lm.tol)
        # Σinv diagonal
        writedlm(fid, lm.Σinv', " ")
        # Word embeddings
        idxs = [lm.vocab_hash[word] for word in lm.vocab]
        writedlm(fid, [lm.vocab lm.Uᵀ[:, idxs]'], " ")
    end
end

"""
    load_lsa_model(filename, eltype; [sparse=false])

Loads an LSA model from `filename` into an LSA model object. The embeddings matrix
element type is specified by `eltype` (default `DEFAULT_FLOAT_TYPE`) while the keyword argument
`sparse` specifies whether the matrix should be sparse or not.
"""
function load_lsa_model(filename::AbstractString, ::Type{T}=DEFAULT_FLOAT_TYPE;
                        sparse::Bool=false) where T<: AbstractFloat
    # Matrix type for LSA model
    A = ifelse(sparse, SparseMatrixCSC{T, Int}, Matrix{T})
    # Define parsed variables local to outer scope of do statement
    local vocab, vocab_hash, Σinv, Uᵀ, stats, idf, nwords, ngram_complexity, κ, β, tol
    open(filename, "r") do fid
        readline(fid)  # first line, header
        line = readline(fid)
        vocab_size, k = map(x -> parse(Int, x), split(line, ' '))
        # Preallocate
        vocab = Vector{String}(undef, vocab_size)
        vocab_hash = OrderedDict{String, Int}()
        if sparse
            Uᵀ = SparseMatrixCSC{T, Int}(UniformScaling(0), k, vocab_size)
            Σinv = spzeros(T, k)
        else
            Uᵀ = Matrix{T}(undef, k, vocab_size)
            Σinv = zeros(T, k)
        end
        # Start parsing the rest of the file
        stats = Symbol(strip(readline(fid)))
        idf = map(x->parse(T, x), split(readline(fid), ' '))
        nwords = parse(T, readline(fid))
        ngram_complexity = parse(Int, readline(fid))
        κ = parse(Int, readline(fid))
        β = parse(Float64, readline(fid))
        tol = parse(Float64, readline(fid))
        Σinv = map(x->parse(T, x), split(readline(fid), ' '))
        for i in 1:vocab_size
            line = strip(readline(fid))
            parts = split(line, ' ')
            word = parts[1]
            vector = map(x-> parse(T, x), parts[2:end])
            vocab[i] = word
            push!(vocab_hash, word=>i)
            Uᵀ[:, i] = vector
        end
    end
    return LSAModel{String,T,A,Int}(vocab, vocab_hash, Σinv, Uᵀ, stats, idf,
                                    nwords, ngram_complexity, κ, β, tol)
end


text:
#Lambert-W

## Lambert-W package on github does not work with julia 0.7
## Native implementations
function lambertw_lb(z)
	@assert -exp(-1) <= z <= 0 "Lambertw only defined on [-1/e, 0] : $z"
	x = -exp(1) * z
	-1 - sqrt(2 * log(1/x)) - log(1/x)
end

function lambertw_ub(z)
	@assert -exp(-1) <= z <= 0 "Lambertw only defined on [-1/e, 0] : $z"
	x = -exp(1) * z
	-1 - sqrt(2 * log(1/x)) - 2 * log(1/x) / 3
end

#the -1 branch
function lambertw(x)
	@assert -exp(-1) <= x <= 0 "Lambertw -1 only defined on [-1/e, 0] : $x"
	#use the Chaterizighou bound to get a guesstimate
	lb = lambertw_lb(x)
	f(z) = z * exp(z) - x
	find_zero(f, (lb, -1.))
end

#the 0 branch
function lambertw0(x)
	@assert x >= -exp(-1) "Lambertw0 only defined on [-1/e, Infty] : $x"

	#Really loose upper bound
	if x <= exp(1)
		ub = 1.
	else
		ub = log(x)
	end

	f(z) = z * exp(z) - x
	find_zero(f, (-1., ub))
end

ext: jl


text:
using CUDAdrv, CuArrays

include("../utils.jl")
include("../read/read_datasets.jl")
include("../initializations.jl")
include("../encodings/cuda/cudaUtilsModule.jl")

function get_new_B(
  B::Matrix{Int16},
  m::Int,
  n::Int)

  newB = Matrix{Int16}(m, n); # codes
  @inbounds @simd for j = 1:m*n
    newB[j] = B[j];
  end

  return newB
end

"Encodes a database with ILS in cuda"
function encode_icm_cuda_single(
  RX::Matrix{Float32},         # in. The data to encode
  B::Matrix{Int16},            # in. Initial list of codes
  C::Vector{Matrix{Float32}},  # in. Codebooks
  ilsiters::Vector{Int64},     # in. ILS iterations to record Bs and obj function. Its max is the total number of iteration we will run.
  icmiter::Integer,            # in. Number of ICM iterations
  npert::Integer,              # in. Number of entries to perturb
  randord::Bool,               # in. Whether to randomize the order in which nodes are visited in ILS
  V::Bool=false)               # in. Whether to print progress

  d, n = size( RX )
  m    = length( C )
  _, h = size( C[1] )

  # Number of results to keep track of
  nr   = length( ilsiters )

  # Make space for the outputs
  Bs   = Vector{Matrix{Int16}}(nr)
  objs = Vector{Float32}(nr)

  # === Compute binary terms (products between all codebook pairs) ===
  binaries, cbi = get_binaries( C )
  _, ncbi       = size( cbi )

  # Create a transposed copy of the binaries for cache-friendliness
  binaries_t = similar( binaries )
  for j = 1:ncbi
    binaries_t[j] = binaries[j]'
  end

  # Create an index from codebook pairs to indices
  cbpair2binaryidx   = zeros(Int32, m, m)
  for j = 1:ncbi
    cbpair2binaryidx[ cbi[1,j], cbi[2,j] ] = j
  end

  dev = CuDevice(0)
  ctx = CuContext(dev)

  # Initialize the cuda module, and choose the GPU
  gpuid = 0
  CudaUtilsModule.init( gpuid, "src/encodings/cuda/cudautils.ptx" )
  # CUDArt.device( devlist[gpuid] )

  # === Create a state for random number generation ===
  if V; @printf("Creating %d random states... ", n); tic(); end
  d_state = CUDAdrv.Mem.alloc( n*64 )
  CudaUtilsModule.setup_kernel( cld(n, 1024), 1024, Cint(n), d_state )

  CUDAdrv.synchronize(ctx)
  if V; @printf("done in %.2f seconds\n", toq()); end

  # Measure time for encoding
  if V; tic(); end

  # Copy X and C to the GPU
  d_RX = CuArrays.CuArray(RX)
  d_C  = CuArrays.CuArray(cat(2, C... ))

  # === Get unaries in the gpu ===
  d_prevcost = CuArrays.CuArray{Cfloat}(n)
  d_newcost  = CuArrays.CuArray{Cfloat}(n)

  # CUBLAS.cublasCreate_v2( CUBLAS.cublashandle )
  d_unaries   = Vector{CuArrays.CuArray{Float32}}(m)
  d_codebooks = Vector{CuArrays.CuArray{Float32}}(m)
  for j = 1:m
    d_codebooks[j] = CuArrays.CuArray(C[j])
    # -2 * C' * X
    d_unaries[j] = CuArrays.BLAS.gemm('T', 'N', -2.0f0, d_codebooks[j], d_RX)
    # d_unaries[j] = -2.0f0 * d_codebooks[j]' * d_RX <-- thus runs out of memory real fast

    # Add self-codebook interactions || C_{i,i} ||^2
    CudaUtilsModule.vec_add( n, (1,h), d_unaries[j].buf, CuArrays.CuArray(diag( C[j]' * C[j] )).buf, Cint(n), Cint(h) )
  end

  # === Get binaries to the GPU ===
  d_binaries  = Vector{CuArrays.CuArray{Float32}}( ncbi )
  d_binariest = Vector{CuArrays.CuArray{Float32}}( ncbi )
  for j = 1:ncbi
    d_binaries[j]  = CuArrays.CuArray( binaries[j] )
    d_binariest[j] = CuArrays.CuArray( binaries_t[j] )
  end

  # Allocate space for temporary results
  bbs = Vector{Matrix{Cfloat}}(m-1)

  # Initialize the previous cost
  prevcost = veccost( RX, B, C )
  IDX = 1:n;

  # For codebook i, we have to condition on these codebooks
  to_look      = 1:m
  to_condition = zeros(Int32, m-1, m)
  for j = 1:m
    tmp = collect(1:m)
    splice!( tmp, j )
    to_condition[:,j] = tmp
  end

  # Loop for the number of requested ILS iterations
  for i = 1:maximum( ilsiters )

    # @show i, ilsiters
    # @time begin

    to_look_r      = to_look
    to_condition_r = to_condition

    # Compute the cost of the previous assignments
    # CudaUtilsModule.veccost(n, (1, d), d_RX, d_C, CuArray( convert(Matrix{Cuchar}, (B')-1) ), d_prevcost, Cint(m), Cint(n))
    CudaUtilsModule.veccost2(n, (1, d), d_RX.buf, d_C.buf, CuArrays.CuArray( convert(Matrix{Cuchar}, (B')-1) ).buf, d_prevcost.buf, Cint(d), Cint(m), Cint(n))
    CUDAdrv.synchronize(ctx)
    prevcost = Array( d_prevcost )

    newB = copy(B)

    # Randomize the visit order in ICM
    if randord
      to_look_r      = randperm( m )
      to_condition_r = to_condition[:, to_look_r]
    end

    d_newB = CuArrays.CuArray( convert(Matrix{Cuchar}, newB-1 ) )

    # Perturn npert entries in each code
    CudaUtilsModule.perturb( n, (1,m), d_state, d_newB.buf, Cint(n), Cint(m), Cint(npert) )

    newB = Array( d_newB )
    Bt   = newB'
    d_Bs = CuArrays.CuArray( Bt )

    CUDAdrv.synchronize(ctx)

    # Run the number of ICM iterations requested
    for j = 1:icmiter

      kidx = 1;
      # Loop through each MRF node
      for k = to_look_r

        lidx = 1;
        for l = to_condition_r[:, kidx]
          # Determine the pairwise tables that we'll use in this conditioning
          if k < l
            binariidx = cbpair2binaryidx[ k, l ]
            bbs[lidx] = binaries[ binariidx ]
          else
            binariidx = cbpair2binaryidx[ l, k ]
            bbs[lidx] = binaries_t[ binariidx ]
          end
          lidx = lidx+1;
        end

        # Transfer pairwise tables to the GPU
        d_bbs = CuArrays.CuArray( convert(Matrix{Cfloat}, cat(2,bbs...)) )
        # Sum binaries (condition) and minimize
        CudaUtilsModule.condition_icm3(
          n, (1, h),
          d_unaries[k].buf, d_bbs.buf, d_Bs.buf, Cint(k-1), Cint(m), Cint(n) )
        CUDAdrv.synchronize(ctx)

        kidx = kidx + 1;
      end # for k = to_look_r
    end # for j = 1:icmiter

    newB = Array( d_Bs )
    newB = convert(Matrix{Int16}, newB')
    newB .+= 1

    # Keep only the codes that improved
    CudaUtilsModule.veccost2(n, (1, d), d_RX.buf, d_C.buf, d_Bs.buf, d_newcost.buf, Cint(d), Cint(m), Cint(n))
    CUDAdrv.synchronize(ctx)

    newcost = Array( d_newcost )

    areequal = newcost .== prevcost
    if V; @printf(" ILS iteration %d/%d done. ", i, maximum(ilsiters)); end
    if V; @printf("%5.2f%% new codes are equal. ", 100*sum(areequal)/n ); end

    arebetter = newcost .< prevcost
    if V; @printf("%5.2f%% new codes are better.\n", 100*sum(arebetter)/n ); end

    newB[:, .~arebetter] = B[:, .~arebetter]
    B = copy( newB )
    # end # @time

    # Check if this # of iterations was requested
    if i in ilsiters

      ithidx = find( i .== ilsiters )[1]

      # Compute and save the objective
      obj = qerror( RX, B, C )
      # @show obj
      objs[ ithidx ] = obj
      # @show size(B)
      Bs[ ithidx ] = B

    end # end if i in ilsiters
  end # end for i=1:max(ilsiters)

  # CUBLAS.cublasDestroy_v2( CUBLAS.cublashandle )
  CudaUtilsModule.finit()

  destroy!(ctx)
  # end # do devlist

  if V; @printf(" Encoding done in %.2f seconds\n", toq()); end

  return Bs, objs
end

# OLD API
# function encode_icm_cuda(
#   RX::Matrix{Float32},
#   B::Matrix{Int16},
#   C::Vector{Matrix{Float32}},
#   cbnorms::Vector{Float32},
#   ilsiters::Vector{Int64},
#   icmiter::Integer,
#   npert::Integer,
#   randord::Bool,
#   qdbnorms::Bool
#   )
#
#   return Bs, objs
# end

"Encodes a database with ILS in cuda"
function encode_icm_cuda(
  RX::Matrix{Float32},         # in. The data to encode
  B::Matrix{Int16},            # in. Initial list of codes
  C::Vector{Matrix{Float32}},  # in. Codebooks
  ilsiters::Vector{Int64},     # in. ILS iterations to record Bs and obj function. Its max is the total number of iteration we will run.
  icmiter::Integer,            # in. Number of ICM iterations
  npert::Integer,              # in. Number of entries to perturb
  randord::Bool,               # in. Whether to randomize the order in which nodes are visited in ILS
  nsplits::Integer=2,          # in. Number of splits of the data (for limited memory GPUs)
  V::Bool=false)

  # TODO check that splits >= 1
  if nsplits == 1
    Bs, objs =  encode_icm_cuda_single(RX, B, C, ilsiters, icmiter, npert, randord, V)
    return Bs, objs
  end

  # Split the data
  d, n = size( RX )
  splits = splitarray(1:n, nsplits)
  nr = length(ilsiters)

  Bs   = Vector{Matrix{Int16}}(nr)
  objs = zeros(Float32,nr)

  # Create storage space for the codes
  for i = 1:nr; Bs[i] = Matrix{Int16}(size(B)); end

  # Run encoding in the GPU for each split
  for i = 1:nsplits
    aaBs, _ = encode_icm_cuda_single(RX[:,splits[i]], B[:,splits[i]], C, ilsiters, icmiter, npert, randord, V)
    for j = 1:nr
      # Save the codes
      Bs[j][:,splits[i]] = aaBs[j]
    end
  end

  # Compute the cost again
  for j = 1:nr
    objs[j] = qerror(RX, Bs[j], C)
  end

  return Bs, objs
end

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
function moma(flux_soln_array::Array{Float64,1}, data_dictionary::Dict{String,Any})

    # what is the size of the array?
    number_of_fluxes = length(flux_soln_array)

    # Setup the Gurobi solver parameters -
    env = Gurobi.Env()

    # set presolve to 0
    setparam!(env, "Presolve", 0)

    # Setup objective function -
    HM = Matrix{Float64}(I,number_of_fluxes,number_of_fluxes)
    fV = -1*flux_soln_array

    # get the bounds -
    flux_bounds_array = data_dictionary["flux_bounds_array"]
    lba = flux_bounds_array[:,1]
    uba = flux_bounds_array[:,2]

    # Get constraints -
    stoichiometric_matrix = data_dictionary["stoichiometric_matrix"]
    bV = data_dictionary["species_bounds_array"][:,1]

    # setup the model -
    model = gurobi_model(env; name = "QP_01", sense = :minimize, H = HM, f = fV, Aeq = stoichiometric_matrix, beq = bV, lb = lba, ub = uba)

    # run optimization
    Gurobi.optimize(model)

    # show results
    sol = get_solution(model)

    # return -
    return sol
end

function objective_function_sweep(sweep_flux_index_array::Array{Int64,1},data_dictionary::Dict{String,Any}, number_of_divisions::Int64)

    # how many fluxes do we have?
    number_of_fluxes = data_dictionary["number_of_reactions"]
    number_of_fluxes_to_sample = length(sweep_flux_index_array)

    # make a copy of the flux_bounds_array -
    fba_original = copy(data_dictionary["flux_bounds_array"])

    # initialize -
    flux_ensemble = zeros(number_of_fluxes,1)

    # do max -
    for flux_index in sweep_flux_index_array

        # initialize -
        c_vector = zeros(number_of_fluxes)
        c_vector[flux_index] = -1.0

        # update the data d -
        data_dictionary["objective_coefficient_array"] = c_vector

        # evaluate lp max and min -
        (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag,status_flag) = calculate_optimal_flux_distribution(data_dictionary)

        # add this solution to the ensemble -
        flux_ensemble = [flux_ensemble calculated_flux_array]

        # what is the original max solution?
        flux_solution_max = calculated_flux_array

        # sample the between the bounds for this flux -
        theta = 0.0
        fba_new = fba_original
        for sample_index = 1:number_of_divisions

            theta = (number_of_divisions - sample_index)*(1/number_of_divisions)
            upper_bound = theta*flux_solution_max[flux_index]
            fba_new[flux_index,2] = upper_bound
            data_dictionary["flux_bounds_array"] = fba_new

            # call the lp again -
            (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag,status_flag) = calculate_optimal_flux_distribution(data_dictionary)

            # cache -
            flux_ensemble = [flux_ensemble calculated_flux_array]
        end
    end

    return flux_ensemble[:,2:end]
end

function calculate_optimal_flux_distribution(data_dictionary)

    # set some constants -
    TIME_RESTART_LIM = 60

    # Get the stoichiometric_matrix from data_dictionary -
    tmp_stoichiometric_matrix = data_dictionary["stoichiometric_matrix"]
    additional_constraint_array = data_dictionary["additional_constraint_array"];
    stoichiometric_matrix = [tmp_stoichiometric_matrix ; additional_constraint_array];
    (number_of_species,number_of_fluxes) = size(stoichiometric_matrix);

    # # Setup the GLPK problem -
    lp_problem = GLPK.Prob();
    GLPK.set_prob_name(lp_problem, "sample");
    GLPK.set_obj_name(lp_problem, "objective")

    # Set solver parameters
    solver_parameters = GLPK.SimplexParam();
    GLPK.init_smcp(solver_parameters);
    solver_parameters.msg_lev = GLPK.MSG_OFF;

    # Are we doing min -or- max?
    min_flag = data_dictionary["is_minimum_flag"];
    if min_flag == true
    	GLPK.set_obj_dir(lp_problem, GLPK.MIN);
    else
    	GLPK.set_obj_dir(lp_problem, GLPK.MAX);
    end

    # Set the number of constraints and fluxes -
    GLPK.add_rows(lp_problem, number_of_species);
    GLPK.add_cols(lp_problem, number_of_fluxes);

    # Setup flux bounds, and objective function -
    default_bounds_array = data_dictionary["flux_bounds_array"]
    (number_of_fluxes,number_of_bounds) = size(default_bounds_array)
    for flux_index = 1:number_of_fluxes

    	flux_lower_bound = default_bounds_array[flux_index,1]
    	flux_upper_bound = default_bounds_array[flux_index,2]

    	# Check bounds type ... default is DB -
    	if (flux_upper_bound == flux_lower_bound)
    		flux_constraint_type = GLPK.FX
    	else
    		flux_constraint_type = GLPK.DB
    	end

        # check -
        if (flux_upper_bound < flux_lower_bound)
            println("Bound violation at $(flux_index)")
        end

    	# flux symbol? (later use name - for now, fake it)
    	flux_symbol = "R_"*string(flux_index)

    	# Set the bounds in GLPK -
    	GLPK.set_col_name(lp_problem, flux_index, flux_symbol);
    	GLPK.set_col_bnds(lp_problem, flux_index, flux_constraint_type, flux_lower_bound, flux_upper_bound);
    end

    # Setup objective function -
    objective_coefficient_array = data_dictionary["objective_coefficient_array"]
    for (flux_index,obj_coeff) in enumerate(objective_coefficient_array)

    	# Set the objective function value in GLPK -
    	GLPK.set_obj_coef(lp_problem, flux_index, obj_coeff);
    end

    # Setup problem constraints for the metabolites -
    species_bounds_array = data_dictionary["species_bounds_array"]
    for species_index = 1:number_of_species

    	species_lower_bound = species_bounds_array[species_index,1]
    	species_upper_bound = species_bounds_array[species_index,2]

    	# defualt
    	species_constraint_type = GLPK.FX
    	if (species_lower_bound != species_upper_bound)
    		species_constraint_type = GLPK.DB
    	end

    	# set the symbol -
    	species_symbol = "x_"*string(species_index)

    	# Set the species bounds in GLPK -
    	GLPK.set_row_name(lp_problem, species_index, species_symbol);
    	GLPK.set_row_bnds(lp_problem, species_index, species_constraint_type, species_lower_bound, species_upper_bound);

    end

    # Setup the stoichiometric array -
    counter = 1;
    row_index_array = zeros(Int,number_of_species*number_of_fluxes);
    col_index_array = zeros(Int,number_of_species*number_of_fluxes);
    species_index_vector = collect(1:number_of_species);
    flux_index_vector = collect(1:number_of_fluxes);
    flat_stoichiometric_array = zeros(Float64,number_of_species*number_of_fluxes);
    for species_index in species_index_vector
    	for flux_index in flux_index_vector
    		row_index_array[counter] = species_index;
    		col_index_array[counter] = flux_index;
    		flat_stoichiometric_array[counter] = stoichiometric_matrix[species_index,flux_index];
    		counter+=1;
    	end
    end
    GLPK.load_matrix(lp_problem, number_of_species*number_of_fluxes, row_index_array, col_index_array, flat_stoichiometric_array);

    # Call the solver -
    exit_flag = GLPK.simplex(lp_problem, solver_parameters);

    # Get the objective function value -
    objective_value = GLPK.get_obj_val(lp_problem);

    # Get the calculated flux values from GLPK -
    calculated_flux_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	calculated_flux_array[flux_index] = GLPK.get_col_prim(lp_problem, flux_index);
    end

    # Get the dual values -
    dual_value_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	dual_value_array[flux_index] = GLPK.get_col_dual(lp_problem, flux_index);
    end

    # is this solution optimal?
    status_flag = GLPK.get_status(lp_problem)

    # Calculate the uptake array -
    uptake_array = stoichiometric_matrix*calculated_flux_array;

    # Formulate the return tuple -
    return (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag,status_flag);
end

function flux_variabilty_analysis(data_dictionary, additional_flux_constraints_array)

    # --- Phase 0: GET SOME STUFF ---------------------------------------------------#
    # Get the stoichiometric_matrix from data_dictionary -
    tmp_stoichiometric_matrix = data_dictionary["stoichiometric_matrix"]
    stoichiometric_matrix = Float64[]

    # do we have additional constraints?
    if haskey(data_dictionary,"additional_constraint_array") == true
        additional_constraint_array = data_dictionary["additional_constraint_array"]
        stoichiometric_matrix = [tmp_stoichiometric_matrix ; additional_constraint_array]
    else
        stoichiometric_matrix = tmp_stoichiometric_matrix
    end

    # get the size -
    (number_of_species,number_of_fluxes) = size(stoichiometric_matrix);

    # initialize the output array -
    calculated_flux_array = zeros(number_of_fluxes,3)

    # set some constants -
    TIME_RESTART_LIM = 60
    EXIT_FLAG_SUCCESS = 0
    TOL = 1e-6
    # ------------------------------------------------------------------------------- #

    # --- Phase 1: SOLVE THE ORIGINAL PROBLEM --------------------------------------- #

    # Setup the GLPK problem -
    lp_problem = GLPK.Prob();
    GLPK.set_prob_name(lp_problem, "sample");
    GLPK.set_obj_name(lp_problem, "objective")

    # Are we doing min -or- max?
    min_flag = data_dictionary["is_minimum_flag"];
    if min_flag == true
    	GLPK.set_obj_dir(lp_problem, GLPK.MIN);
    else
    	GLPK.set_obj_dir(lp_problem, GLPK.MAX);
    end

    # Set the number of constraints and fluxes -
    GLPK.add_rows(lp_problem, number_of_species);
    GLPK.add_cols(lp_problem, number_of_fluxes);

    # Set solver parameters
    solver_parameters = GLPK.SimplexParam();
    solver_parameters.msg_lev = GLPK.MSG_ERR;
    solver_parameters.presolve = GLPK.ON;

    # Setup flux bounds, and objective function -
    default_bounds_array = data_dictionary["flux_bounds_array"]
    (number_of_fluxes,number_of_bounds) = size(default_bounds_array)
    for flux_index = 1:number_of_fluxes

    	flux_lower_bound = default_bounds_array[flux_index,1]
    	flux_upper_bound = abs(default_bounds_array[flux_index,2])

    	# Check bounds type ... default is DB -
    	if (flux_upper_bound == flux_lower_bound)
    		flux_constraint_type = GLPK.FX
    	else
    		flux_constraint_type = GLPK.DB
    	end

    	# flux symbol? (later use name - for now, fake it)
    	flux_symbol = "R_"*string(flux_index)

    	# Set the bounds in GLPK -
    	GLPK.set_col_name(lp_problem, flux_index, flux_symbol);
    	GLPK.set_col_bnds(lp_problem, flux_index, flux_constraint_type, flux_lower_bound, flux_upper_bound);
    end

    # Setup objective function -
    objective_coefficient_array = data_dictionary["objective_coefficient_array"]
    for (flux_index,obj_coeff) in enumerate(objective_coefficient_array)

    	# Set the objective function value in GLPK -
    	GLPK.set_obj_coef(lp_problem, flux_index, obj_coeff);
    end

    # Setup problem constraints for the metabolites -
    species_bounds_array = data_dictionary["species_bounds_array"]
    for species_index = 1:number_of_species

    	species_lower_bound = species_bounds_array[species_index,1]
    	species_upper_bound = species_bounds_array[species_index,2]

    	# defualt
    	species_constraint_type = GLPK.FX
    	if (species_lower_bound != species_upper_bound)
    		species_constraint_type = GLPK.DB
    	end

    	# set the symbol -
    	species_symbol = "x_"*string(species_index)

    	# Set the species bounds in GLPK -
    	GLPK.set_row_name(lp_problem, species_index, species_symbol);
    	GLPK.set_row_bnds(lp_problem, species_index, species_constraint_type, species_lower_bound, species_upper_bound);

    end

    # Setup the stoichiometric array -
    counter = 1;
    row_index_array = zeros(Int,number_of_species*number_of_fluxes);
    col_index_array = zeros(Int,number_of_species*number_of_fluxes);
    species_index_vector = collect(1:number_of_species);
    flux_index_vector = collect(1:number_of_fluxes);
    flat_stoichiometric_array = zeros(Float64,number_of_species*number_of_fluxes);
    for species_index in species_index_vector
    	for flux_index in flux_index_vector
    		row_index_array[counter] = species_index;
    		col_index_array[counter] = flux_index;
    		flat_stoichiometric_array[counter] = stoichiometric_matrix[species_index,flux_index];
    		counter+=1;
    	end
    end
    GLPK.load_matrix(lp_problem, number_of_species*number_of_fluxes, row_index_array, col_index_array, flat_stoichiometric_array);

    # solve --
    # Call the solver -
    exit_flag = GLPK.simplex(lp_problem, solver_parameters);

    # Get the objective function value -
    objective_value = GLPK.get_obj_val(lp_problem);

    # Get the calculated flux values for the optimal distrubution from GLPK -
    for flux_index in flux_index_vector
    	calculated_flux_array[flux_index,1] = GLPK.get_col_prim(lp_problem, flux_index);
    end

    # Get the dual values -
    dual_value_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	dual_value_array[flux_index] = GLPK.get_col_dual(lp_problem, flux_index);
    end
    # ------------------------------------------------------------------------------- #

    # --- Phase 2: SOLVE THE FVA PROBLEM -------------------------------------------- #
    target_value = 0.0
    if (GLPK.get_obj_dir(lp_problem) == GLPK.MAX)
        target_value = floor(objective_value/TOL)*(TOL);
    else
        target_value = ceil(objective_value/TOL)*(TOL);
    end

    # do we have any additional flux constraints?
    number_of_extra_flux_constraints = length(additional_flux_constraints_array)

    # Add an extra row for the objective -
    m = GLPK.add_rows(lp_problem,number_of_extra_flux_constraints+1)
    if (GLPK.get_obj_dir(lp_problem) == GLPK.MAX)
        GLPK.set_row_bnds(lp_problem, m, GLPK.LO,target_value,0.0);
    else
        GLPK.set_row_bnds(lp_problem, m, GLPK.UP,0.0,target_value);
    end

    # populate new row, add the object value -
    index_array = Int[]
    value_array = Float64[]
    for flux_index = 1:number_of_fluxes
        push!(index_array,flux_index)
        push!(value_array,GLPK.get_obj_coef(lp_problem, flux_index))
    end
    GLPK.set_mat_row(lp_problem,m,number_of_fluxes,index_array,value_array)

    # Add additional rows for the extra flux constraints -
    for flux_constraint_index = 1:number_of_extra_flux_constraints

        # what rate do I want to constrain?
        target_flux_index = additional_flux_constraints_array[flux_constraint_index]

        # What is my current row index?
        current_row_index = m + flux_constraint_index

        # What is my target value?
        target_value = calculated_flux_array[target_flux_index]

        # add the row bound -
        if (GLPK.get_obj_dir(lp_problem) == GLPK.MAX)
            GLPK.set_row_bnds(lp_problem, current_row_index, GLPK.LO,target_value,0.0);
        else
            GLPK.set_row_bnds(lp_problem, current_row_index, GLPK.UP,0.0,target_value);
        end

        # Add the extra row to the st array -
        index_array = Int[]
        value_array = Float64[]
        for flux_index = 1:number_of_fluxes
            push!(index_array,flux_index)

            if (flux_index == target_flux_index)
                push!(value_array,1.0)
            else
                push!(value_array,0.0)
            end
        end
        GLPK.set_mat_row(lp_problem,current_row_index,number_of_fluxes,index_array,value_array)
    end

    # zero out the objective coefficients -
    for flux_index = 1:number_of_fluxes
        GLPK.set_obj_coef(lp_problem, flux_index, 0.0);
    end

    solver_parameters.presolve = GLPK.OFF;
    solver_parameters.msg_lev = GLPK.MSG_OFF;
    solver_parameters.tm_lim = 1000*TIME_RESTART_LIM;

    # main FVA loop -
    for round_index = 1:2

        # are we min are maximizing the flux?
        if round_index == 1
            GLPK.set_obj_dir(lp_problem, GLPK.MIN);
        else
            GLPK.set_obj_dir(lp_problem, GLPK.MAX);
        end

        for flux_index = 1:number_of_fluxes

            # run the calc -
            GLPK.set_obj_coef(lp_problem, flux_index, 1.0);
            exit_flag = GLPK.simplex(lp_problem, solver_parameters);

            # reset the objective coefficient -
            GLPK.set_obj_coef(lp_problem,flux_index,0.0)

            if (exit_flag != EXIT_FLAG_SUCCESS)

                solver_parameters.tm_lim = 10000*TIME_RESTART_LIM;
                solver_parameters.presolve = GLPK.ON
                GLPK.adv_basis(lp_problem, 0);
                exit_flag = GLPK.simplex(lp_problem, solver_parameters);
                if (exit_flag != EXIT_FLAG_SUCCESS)
                    # ok, something muy malo is happening. throw an error ...
                    error("ERROR: Unable to solve modified problem. Exiting ... ")
                end
            end

            # reset the problem options -
            solver_parameters.presolve = GLPK.OFF;
            solver_parameters.tm_lim = 1000*TIME_RESTART_LIM;

            if (GLPK.get_obj_dir(lp_problem) == GLPK.MIN)
                calculated_flux_array[flux_index,2] = GLPK.get_obj_val(lp_problem)
            else
                calculated_flux_array[flux_index,3] = GLPK.get_obj_val(lp_problem)
            end
        end
    end
    # ------------------------------------------------------------------------------- #

    # return the array -
    return (calculated_flux_array,dual_value_array)
end

function calculate_optimal_flux_distribution2(data_dictionary) #modified, doesn't
    #include the additional flux array

    # set some constants -
    TIME_RESTART_LIM = 60

    # Get the stoichiometric_matrix from data_dictionary -
    stoichiometric_matrix = data_dictionary["stoichiometric_matrix"]
    #additional_constraint_array = data_dictionary["additional_constraint_array"];
    #stoichiometric_matrix = [tmp_stoichiometric_matrix ; additional_constraint_array];
    (number_of_species,number_of_fluxes) = size(stoichiometric_matrix);

    # # Setup the GLPK problem -
    lp_problem = GLPK.Prob();
    GLPK.set_prob_name(lp_problem, "sample");
    GLPK.set_obj_name(lp_problem, "objective")

    # Set solver parameters
    solver_parameters = GLPK.SimplexParam();
    GLPK.init_smcp(solver_parameters);
    solver_parameters.msg_lev = GLPK.MSG_OFF;

    # Are we doing min -or- max?
    min_flag = data_dictionary["is_minimum_flag"];
    if min_flag == true
    	GLPK.set_obj_dir(lp_problem, GLPK.MIN);
    else
    	GLPK.set_obj_dir(lp_problem, GLPK.MAX);
    end

    # Set the number of constraints and fluxes -
    GLPK.add_rows(lp_problem, number_of_species);
    GLPK.add_cols(lp_problem, number_of_fluxes);

    # Setup flux bounds, and objective function -
    default_bounds_array = data_dictionary["flux_bounds_array"]
    (number_of_fluxes,number_of_bounds) = size(default_bounds_array)
    for flux_index = 1:number_of_fluxes

    	flux_lower_bound = default_bounds_array[flux_index,1]
    	flux_upper_bound = default_bounds_array[flux_index,2]

    	# Check bounds type ... default is DB -
    	if (flux_upper_bound == flux_lower_bound)
    		flux_constraint_type = GLPK.FX
    	else
    		flux_constraint_type = GLPK.DB
    	end

        # check -
        if (flux_upper_bound < flux_lower_bound)
            println("Bound violation at $(flux_index)")
        end

    	# flux symbol? (later use name - for now, fake it)
    	flux_symbol = "R_"*string(flux_index)

    	# Set the bounds in GLPK -
    	GLPK.set_col_name(lp_problem, flux_index, flux_symbol);
    	GLPK.set_col_bnds(lp_problem, flux_index, flux_constraint_type, flux_lower_bound, flux_upper_bound);
    end

    # Setup objective function -
    objective_coefficient_array = data_dictionary["objective_coefficient_array"]
    for (flux_index,obj_coeff) in enumerate(objective_coefficient_array)

    	# Set the objective function value in GLPK -
    	GLPK.set_obj_coef(lp_problem, flux_index, obj_coeff);
    end

    # Setup problem constraints for the metabolites -
    species_bounds_array = data_dictionary["species_bounds_array"]
    for species_index = 1:number_of_species

    	species_lower_bound = species_bounds_array[species_index,1]
    	species_upper_bound = species_bounds_array[species_index,2]

    	# defualt
    	species_constraint_type = GLPK.FX
    	if (species_lower_bound != species_upper_bound)
    		species_constraint_type = GLPK.DB
    	end

    	# set the symbol -
    	species_symbol = "x_"*string(species_index)

    	# Set the species bounds in GLPK -
    	GLPK.set_row_name(lp_problem, species_index, species_symbol);
    	GLPK.set_row_bnds(lp_problem, species_index, species_constraint_type, species_lower_bound, species_upper_bound);

    end

    # Setup the stoichiometric array -
    counter = 1;
    row_index_array = zeros(Int,number_of_species*number_of_fluxes);
    col_index_array = zeros(Int,number_of_species*number_of_fluxes);
    species_index_vector = collect(1:number_of_species);
    flux_index_vector = collect(1:number_of_fluxes);
    flat_stoichiometric_array = zeros(Float64,number_of_species*number_of_fluxes);
    for species_index in species_index_vector
    	for flux_index in flux_index_vector
    		row_index_array[counter] = species_index;
    		col_index_array[counter] = flux_index;
    		flat_stoichiometric_array[counter] = stoichiometric_matrix[species_index,flux_index];
    		counter+=1;
    	end
    end
    GLPK.load_matrix(lp_problem, number_of_species*number_of_fluxes, row_index_array, col_index_array, flat_stoichiometric_array);

    # Call the solver -
    exit_flag = GLPK.simplex(lp_problem, solver_parameters);

    # Get the objective function value -
    objective_value = GLPK.get_obj_val(lp_problem);

    # Get the calculated flux values from GLPK -
    calculated_flux_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	calculated_flux_array[flux_index] = GLPK.get_col_prim(lp_problem, flux_index);
    end

    # Get the dual values -
    dual_value_array = zeros(Float64,number_of_fluxes);
    for flux_index in flux_index_vector
    	dual_value_array[flux_index] = GLPK.get_col_dual(lp_problem, flux_index);
    end

    # is this solution optimal?
    status_flag = GLPK.get_status(lp_problem)

    # Calculate the uptake array -
    uptake_array = stoichiometric_matrix*calculated_flux_array;

    # Formulate the return tuple -
    return (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag,status_flag);
end


text:
__precompile__(false)

module TextureManager

using GraphicsManager
const GPU =  GraphicsManager
GPU.addListener(@__MODULE__)

using Images
using ModernGL

mutable struct Texture
  id::Symbol
  refID::GLuint
  created::Bool

  function Texture(id::Symbol)
    this = new(id,0,true)
    set(id,this)
    this
  end
end

LIST = Dict{Symbol, Texture}()

get(id::Symbol) = haskey(LIST, id) ? LIST[id] : nothing
set(id::Symbol, tex::Texture) = LIST[id] = tex

function cleanUp()
  global LIST = typeof(LIST)()
end

init(this::Texture) = this.refID = GPU.create(:TEXTURE, this.id)

function load(id::Symbol)
  this = get(id)
  if this == nothing this = Texture(id) else this.created = false end
  init(this)
  this
end

getRefIDs() = [this.refID for (_,this) in LIST]

function bindTexture(id::Symbol, active::Union{Nothing,Integer}=nothing)
  if active != nothing glActiveTexture(GL_TEXTURE0 + UInt32(active)) end
  glBindTexture(GL_TEXTURE_2D, get(id).refID)
end

export bindTexture

"""
uploads a texture by given file path
"""
function createTextureGray(id::Symbol, path::String)
  this = load(id)
  if !this.created return this.refID end

  img = Images.load(path)
  (imgwidth, imgheight) = size(img)
  gray = Array{Float32}(Gray.(img)) # do not vec

  textureID = GPU.create(:TEXTURE, id)
  glActiveTexture(GL_TEXTURE0)
  glBindTexture(GL_TEXTURE_2D, this.refID)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST) #GL_LINEAR
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST) #GL_LINEAR

  glTexImage2D(GL_TEXTURE_2D, 0, GL_LUMINANCE, imgwidth, imgheight, 0, GL_LUMINANCE, GL_FLOAT, gray)
  #glGenerateMipmap(GL_TEXTURE_2D)
  #glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR)
  #glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR_MIPMAP_LINEAR)
  #glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_LOD_BIAS, -0.8f0)
  glBindImageTexture(0, this.refID, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F)
  glBindTexture(GL_TEXTURE_2D, 0)
  glCheckError("texture")

  this.refID
end

createTextureGray(path::String) = createTextureGray(Symbol(path*" (gray)"), path)

export createTextureGray

"""
uploads a texture by given file path
"""
function createTexture(id::Symbol, path::String)
  this = load(id)
  if !this.created return this.refID end

  img = Images.load(path)
  (imgwidth, imgheight) = size(img)
  imga = reinterpret.(channelview(img)) # do not vec

  glActiveTexture(GL_TEXTURE0)
  glBindTexture(GL_TEXTURE_2D, this.refID)
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, imgwidth, imgheight, 0, GL_RGBA, GL_UNSIGNED_BYTE, imga)
  glGenerateMipmap(GL_TEXTURE_2D)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR) #GL_LINEAR
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_LOD_BIAS, -0.4f0)
  largest_supported_anisotropy = GLfloat[0]
  GL_MAX_TEXTURE_MAX_ANISOTROPY = 0x84FF
  GL_TEXTURE_MAX_ANISOTROPY = 0x84FE
  glGetFloatv(GL_MAX_TEXTURE_MAX_ANISOTROPY, largest_supported_anisotropy)
  glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_MAX_ANISOTROPY, largest_supported_anisotropy[])
  #glGenerateMipmap(GL_TEXTURE_2D)
  glBindTexture(GL_TEXTURE_2D, 0)
  glCheckError("texture")

  this.refID
end

"""
uploads a texture by given file path
"""
createTexture(path::String) = createTexture(Symbol(path), path)

"""
TODO
"""
function createTexture(id::Symbol, sz::Tuple{Integer,Integer})
  this = load(id) #Symbol(id, string(sz[1])*"x"*string(sz[2]))
  if !this.created return this.refID end

  width, height = sz

  glActiveTexture(GL_TEXTURE0)
  glBindTexture(GL_TEXTURE_2D,  this.refID)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEA
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR) #GL_NEAREST
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR) #GL_NEAREST
  glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, width, height, 0, GL_RGBA, GL_FLOAT, C_NULL) #GL_RGBA8, ..., GL_RGBA, GL_UNSIGNED_BYTE  GL_FLOAT
  glBindImageTexture(0,  this.refID, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F)
  glGenerateMipmap(GL_TEXTURE_2D)
  #glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR_MIPMAP_LINEAR)
  #glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR_MIPMAP_LINEAR)
  #glTexParameterf(GL_TEXTURE_2D, GL_TEXTURE_LOD_BIAS, -0.8f0)
  glBindTexture(GL_TEXTURE_2D, 0)
  glCheckError("texture")

  this.refID
end

"""
TODO
"""
function createTextureStorage(id::Symbol, sz::Tuple{Integer,Integer}; level=1)
  this = load(id) #Symbol(string(sz[1])*"x"*string(sz[2])*"x"*string(level))
  if !this.created return this.refID end

  width, height = sz

  glActiveTexture(GL_TEXTURE0)
  glBindTexture(GL_TEXTURE_2D, this.refID)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEAT
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, GL_CLAMP_TO_EDGE) #GL_CLAMP_TO_EDGE,GL_REPEA
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST) #GL_NEAREST
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST_MIPMAP_NEAREST) #GL_NEAREST, GL_NEAREST_MIPMAP_NEAREST
  #glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_FUNC, GL_LEQUAL)
  #glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_COMPARE_MODE, GL_NONE)
  # Useful for debugging purposes so depth shows up as graytone and not just red.
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_SWIZZLE_R, GL_RED)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_SWIZZLE_G, GL_RED)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_SWIZZLE_B, GL_RED)
  glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_SWIZZLE_A, GL_ONE)
  #glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F, width, height, 0, GL_RGBA, GL_FLOAT, C_NULL)
  #glTexImage2D(GL_TEXTURE_2D, 0, GL_DEPTH_COMPONENT24, width, height, 0, GL_DEPTH_COMPONENT, GL_FLOAT, C_NULL)
  #glTexStorage2D(GL_TEXTURE_2D, 1, GL_DEPTH_COMPONENT24, width, height)
  glTexStorage2D(GL_TEXTURE_2D, level, GL_DEPTH24_STENCIL8, width, height)
  glBindTexture(GL_TEXTURE_2D, 0)
  #glBindImageTexture(0, this.refID, 0, GL_FALSE, 0, GL_WRITE_ONLY, GL_RGBA32F)
  glCheckError("texture")

  this.refID
end

"""
TODO
"""
function createTextureMultiSample(id::Symbol, sz::Tuple{Integer,Integer})
  this = load(id) #Symbol(string(sz[1])*"x"*string(sz[2])*"x"*string(level))
  if !this.created return this.refID end

  width, height = sz

  glBindTexture(GL_TEXTURE_2D_MULTISAMPLE, this.refID)
  glTexImage2DMultisample(GL_TEXTURE_2D_MULTISAMPLE, 4, GL_RGB, width, height, GL_TRUE)
  glBindTexture(GL_TEXTURE_2D_MULTISAMPLE, 0)

  this.refID
end

export createTextureMultiSample

export createTexture
export createTextureStorage
export uploadTexture

end #TextureManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module EMLstoic
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...reactors.EMLtank_basic
	import EMLtypes.length
	include("stoic/stoic_vap.jl")
	include("stoic/stoic_liq.jl")
	include("stoic/stoic_extent_vap.jl")
	include("stoic/stoic_extent_liq.jl")
	include("stoic/stoic_conv_vap.jl")
	include("stoic/stoic_conv_liq.jl")
end


text:
using UnPack
using LinearAlgebra

export  Moire,
        lattice_img

"""
    Moire(Gmax,nmax,lattice₁,lattice₂)
    Gmax    : the maximum norm of folding vector
    nmax    : the range to search folding vectors
    lattice : constructors of lattice
"""
struct Moire
    fold :: AbstractArray{Int}
    G₁ :: AbstractArray{Real}               # reciprocal lattice of layer 1
    G₂ :: AbstractArray{Real}               # reciprocal lattice of layer 1
    numfold :: Int                          # fold k space of layer 2 for numfold times
    M :: Int                          # fold k space of layer 2 for numfold times
    #matrix :: AbstractArray{ComplexF64}     # matrix for moire hamiltonian

    function Moire(Gmax,nmax,lattice₁,lattice₂)
        function make_fold(Gmax,nmax,b)
            G = Array{Float64,2}(undef,(2nmax+1)^2,2)
            F = Array{Int,2}(undef,(2nmax+1)^2,2)
            index = 0
            @inbounds for g₁ in -nmax:nmax
                @inbounds for g₂ in -nmax:nmax
                    g = [g₁ g₂]*b
                    if norm(g) ≤ Gmax
                        index += 1
                        F[index,:] = [g₁ g₂]
                        G[index,:] = g
                    end
                end
            end
            return G[1:index,:] , F[1:index,:]
        end
        G₁ , fold=make_fold(Gmax,nmax,lattice₁.b)
        new(fold,
            G₁,
            similar(G₁),
            size(G₁,1),
            lattice₁.M+lattice₂.M*size(G₁,1)
        )
    end
end

"""
    lattice_img(Moire,lattice₁)
    lattice₂    : lattice constructor of layer1
"""
function lattice_img(Moire,lattice₁)
    @unpack a, τ, b, Ham!, M= lattice₁
    @unpack fold, G₁, G₂, numfold , M = Moire

    p1 = plot(title = "Lattice") ; p2 = plot(title="Reciprocal lattice");
    lat = fold*a' ; G = fold*b
    scatter!(p1,lat[:,1] , lat[:,2]) ; scatter!(p2,G[:,1],G[:,2])
    plot(p1,p2,aspect_ratio=:equal,msw = 0,ms=2,mc = :blue)
end


text:
"""
    glcm = glcm(img, distance, angle, mat_size=16)
    glcm = glcm(img, distances, angle, mat_size=16)
    glcm = glcm(img, distance, angles, mat_size=16)
    glcm = glcm(img, distances, angles, mat_size=16)

Calculates the GLCM (Gray Level Co-occurence Matrix) of an image. The `distances` and `angles` arguments may be
a single integer or a vector of integers if multiple GLCMs need to be calculated. The `mat_size` argument is used
to define the granularity of the GLCM.
"""
function glcm(img::AbstractArray{T, 2}, distance::Integer, angle::Real, mat_size::Integer = 16) where T<:Real
    max_img = maximum(img)
    img_rescaled = map(i -> max(1, Int(ceil(i * mat_size / max_img))), img)
    _glcm(img_rescaled, distance, angle, mat_size)
end

function glcm(img::AbstractArray{T, 2}, distances::Array{Int, 1}, angles::Array{A, 1}, mat_size::Integer = 16) where {T<:Real, A<:Real}
    max_img = maximum(img)
    img_rescaled = map(i -> max(1, Int(ceil(i * mat_size / max_img))), img)
    glcm_matrices = [_glcm(img_rescaled, d, a, mat_size) for d in distances, a in angles]
    glcm_matrices
end

function glcm(img::AbstractArray{T, 2}, distance::Integer, angle::Real, mat_size::Integer = 16) where T<:Colorant
    img_rescaled = map(i -> max(1, Int(ceil((reinterpret(gray(convert(Gray{N0f8}, i)))) * mat_size / 256))), img)
    _glcm(img_rescaled, distance, angle, mat_size)
end

function glcm(img::AbstractArray{T, 2}, distances::Array{Int, 1}, angles::Array{A, 1}, mat_size::Integer = 16) where {T<:Colorant, A<:Real}
    img_rescaled = map(i -> max(1, Int(ceil((reinterpret(gray(convert(Gray{N0f8}, i)))) * mat_size / 256))), img)
    glcm_matrices = [_glcm(img_rescaled, d, a, mat_size) for d in distances, a in angles]
    glcm_matrices
end

glcm(img::AbstractArray{T, 2}, distances::Int, angles::Array{A, 1}, mat_size::Integer = 16) where {T<:Union{Colorant, Real}, A<:Real} = glcm(img, [distances], angles, mat_size)
glcm(img::AbstractArray{T, 2}, distances::Array{Int, 1}, angles::Real, mat_size::Integer = 16) where {T<:Union{Colorant, Real}} = glcm(img, distances, [angles], mat_size)

function _glcm(img::AbstractArray{T, 2}, distance::Integer, angle::Number, mat_size::Integer) where T
    co_oc_matrix = zeros(Int, mat_size, mat_size)
    sin_angle = sin(angle)
    cos_angle = cos(angle)
    for j = 1:size(img, 1)
        j_off = j + Int(round(sin_angle * distance))
        checkbounds(Bool, img, j_off, 1) || continue
        for i = 1:size(img, 2)
            int_one = img[j, i]
            i_off = i + Int(round(cos_angle * distance))

            if checkbounds(Bool, img, j_off, i_off)
                int_two = img[j_off, i_off]
                co_oc_matrix[int_one, int_two] += 1
            end
        end
    end
    co_oc_matrix
end

"""
    glcm = glcm_symmetric(img, distance, angle, mat_size=16)
    glcm = glcm_symmetric(img, distances, angle, mat_size=16)
    glcm = glcm_symmetric(img, distance, angles, mat_size=16)
    glcm = glcm_symmetric(img, distances, angles, mat_size=16)

Symmetric version of the [`glcm`](@ref) function.
"""
function glcm_symmetric(img::AbstractArray, distance::Integer, angle::Real, mat_size=16)
    co_oc_matrix = glcm(img, distance, angle, mat_size)
    co_oc_matrix_trans = co_oc_matrix'
    co_oc_matrix_symm = co_oc_matrix + co_oc_matrix_trans
    co_oc_matrix_symm
end

function glcm_symmetric(img::AbstractArray, distances, angles, mat_size=16)
    co_oc_matrices = glcm(img, distances, angles, mat_size)
    co_oc_matrices_sym = map(gmat -> gmat + gmat', co_oc_matrices)
    co_oc_matrices_sym
end

"""
    glcm = glcm_norm(img, distance, angle, mat_size)
    glcm = glcm_norm(img, distances, angle, mat_size)
    glcm = glcm_norm(img, distance, angles, mat_size)
    glcm = glcm_norm(img, distances, angles, mat_size)

Normalised version of the [`glcm`](@ref) function.
"""
function glcm_norm(img::AbstractArray, distance::Integer, angle::Real, mat_size=16)
    co_oc_matrix = glcm(img, distance, angle, mat_size)
    co_oc_matrix_norm = co_oc_matrix / Float64(sum(co_oc_matrix))
    co_oc_matrix_norm
end

function glcm_norm(img::AbstractArray, distances, angles, mat_size=16)
    co_oc_matrices = glcm(img, distances, angles, mat_size)
    co_oc_matrices_norm = map(gmat -> gmat /= Float64(sum(gmat)), co_oc_matrices)
    co_oc_matrices_norm
end

"""
Multiple properties of the obtained GLCM can be calculated by using the `glcm_prop` function which calculates the
property for the entire matrix. If grid dimensions are provided, the matrix is divided into a grid and the property
is calculated for each cell resulting in a height x width property matrix.
```julia
prop = glcm_prop(glcm, property)
prop = glcm_prop(glcm, height, width, property)
```
Various properties can be calculated like `mean`, `variance`, `correlation`, `contrast`, `IDM` (Inverse Difference Moment),
 `ASM` (Angular Second Moment), `entropy`, `max_prob` (Max Probability), `energy` and `dissimilarity`.
"""
function glcm_prop(gmat::Array{T, 2}, window_height::Integer, window_width::Integer, property::Function) where T<:Real
    k_h = Int(floor(window_height / 2))
    k_w = Int(floor(window_width / 2))
    glcm_size = size(gmat)
    R = CartesianIndices(glcm_size)
    prop_mat = zeros(Float64, glcm_size)
    for I in R
        prop_mat[I] = property(gmat[max(1, I[1] - k_h) : min(glcm_size[1], I[1] + k_h), max(1, I[2] - k_w) : min(glcm_size[2], I[2] + k_w)])
    end
    prop_mat
end

glcm_prop(gmat::Array{T, 2}, window_size::Integer, property::Function) where {T<:Real} = glcm_prop(gmat, window_size, window_size, property)

function glcm_prop(gmat::Array{T, 2}, property::Function) where T<:Real
    property(gmat)
end

function contrast(glcm_window::Array{T, 2}) where T<:Real
    sum([(id[1] - id[2]) ^ 2 * glcm_window[id] for id in CartesianIndices(size(glcm_window))])
end

function dissimilarity(glcm_window::Array{T, 2}) where T<:Real
    sum([glcm_window[id] * abs(id[1] - id[2]) for id in CartesianIndices(size(glcm_window))])
end

function glcm_entropy(glcm_window::Array{T, 2}) where T<:Real
    -sum(map(i -> i * log(i), glcm_window))
end

function ASM(glcm_window::Array{T, 2}) where T<:Real
    sum(map(i -> i ^ 2, glcm_window))
end

function IDM(glcm_window::Array{T, 2}) where T<:Real
    sum([glcm_window[id] / (1 + (id[1] - id[2]) ^ 2) for id in CartesianIndices(size(glcm_window))])
end

function glcm_mean_ref(glcm_window::Array{T, 2}) where T<:Real
    sumref = sum(glcm_window, dims=2)
    meanref = sum([id * sumref[id] for id = 1:size(glcm_window)[1]])
    meanref
end

function glcm_mean_neighbour(glcm_window::Array{T, 2}) where T<:Real
    sumneighbour = sum(glcm_window, dims=1)
    meanneighbour = sum([id * sumneighbour[id] for id = 1:size(glcm_window)[2]])
    meanneighbour
end

function glcm_var_ref(glcm_window::Array{T, 2}) where T<:Real
    mean_ref = glcm_mean_ref(glcm_window)
    sumref = sum(glcm_window, dims=2)
    var_ref = sum([((id - mean_ref) ^ 2) * sumref[id] for id = 1:size(glcm_window)[1]])
    var_ref ^ 0.5
end

function glcm_var_neighbour(glcm_window::Array{T, 2}) where T<:Real
    mean_neighbour = glcm_mean_neighbour(glcm_window)
    sumneighbour = sum(glcm_window, dims=1)
    var_neighbour = sum([((id - mean_neighbour) ^ 2) * sumneighbour[id] for id = 1:size(glcm_window)[2]])
    var_neighbour ^ 0.5
end

function correlation(glcm_window::Array{T, 2}) where T<:Real
    mean_ref = glcm_mean_ref(glcm_window)
    var_ref = glcm_var_ref(glcm_window)
    mean_neighbour = glcm_mean_neighbour(glcm_window)
    var_neighbour = glcm_var_neighbour(glcm_window)
    if var_ref == 0 || var_neighbour == 0
        return 1
    end
    sum([glcm_window[i, j] * (i - mean_ref) * (j - mean_neighbour) for i = 1:size(glcm_window)[1], j = 1:size(glcm_window)[2]]) / (var_ref * var_neighbour)
end

function max_prob(glcm_window::Array{T, 2}) where T<:Real
    maxfinite(glcm_window)
end

function energy(glcm_window::Array{T, 2}) where T<:Real
    ASM(glcm_window) ^ 0.5
end


text:
using Chemostat_Human1
using Test

@testset "Chemostat_Human1.jl" begin
    # Write your tests here.
end


text:
function del_rxn!(metnet::MetNet, rxn::IDER_TYPE)
    rxni = rxnindex(metnet, rxn)
    metnet.rxns[rxni] = EMPTY_SPOT
    z = zero(eltype(metnet.S))
    metnet.S[:, rxni] .= z
    metnet.lb[rxni] = z
    metnet.ub[rxni] = z
    metnet.c[rxni] = z

    return metnet
end

function del_rxn!(metnet::MetNet, rxns::Vector)
    for rxn in rxns
        del_rxn!(metnet, rxn)
    end
    return metnet
end


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

if !isdefined(@__MODULE__, :DefsEnterpriseId)
    const DefsEnterpriseId = String
else
    @warn("Skipping redefinition of DefsEnterpriseId to String")
end


text:
@testset "Ambiguities" begin
    if VERSION.prerelease == () && !Sys.iswindows() && VERSION < v"1.8.0"
        mbs = Test.detect_ambiguities(ManifoldsBase)
        # Interims solution until we follow what was proposed in
        # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2
        fmbs = filter(x -> !any(has_type_in_signature.(x, Identity)), mbs)
        FMBS_LIMIT = 15
        println("Number of ManifoldsBase.jl ambiguities: $(length(fmbs))")
        @test length(fmbs) <= FMBS_LIMIT
        if length(fmbs) > FMBS_LIMIT
            for amb in fmbs
                println(amb)
                println()
            end
        end
        ms = Test.detect_ambiguities(Manifolds)
        # Interims solution until we follow what was proposed in
        # https://discourse.julialang.org/t/avoid-ambiguities-with-individual-number-element-identity/62465/2
        fms = filter(x -> !any(has_type_in_signature.(x, Identity)), ms)
        FMS_LIMIT = 20
        println("Number of Manifolds.jl ambiguities: $(length(fms))")
        if length(fms) > FMS_LIMIT
            for amb in fms
                println(amb)
                println()
            end
        end
        @test length(fms) <= FMS_LIMIT
        # this test takes way too long to perform regularly
        # @test length(our_base_ambiguities()) <= 4
    else
        @info "Skipping Ambiguity tests for pre-release versions"
    end
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

# module: MyGame

import FlatBuffers

mutable struct InParentNamespace
end
FlatBuffers.@ALIGN(InParentNamespace, 1)

InParentNamespace(buf::AbstractVector{UInt8}) = FlatBuffers.read(InParentNamespace, buf)
InParentNamespace(io::IO) = FlatBuffers.deserialize(io, InParentNamespace)


text:
using Test
@isdefined(TestHelpers) || include("TestHelpers.jl")

TestHelpers.@timed_testset "integration/" begin
    include("integration/sequential.jl")
end


text:
using sem, Feather, ModelingToolkit, Statistics, LinearAlgebra,
    Optim, SparseArrays, Test, Zygote, LineSearches, ForwardDiff

## Observed Data
three_path_dat = Feather.read("test/comparisons/three_path_dat.feather")
three_path_par = Feather.read("test/comparisons/three_path_par.feather")

semobserved = SemObsCommon(data = Matrix(three_path_dat))

loss = Loss([SemML(semobserved)])

diff_for = SemForwardDiff(LBFGS(), Optim.Options())

## Model definition
@variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

fun_A = eval(ModelingToolkit.build_function(A, x)[2])
fun_S = eval(ModelingToolkit.build_function(S, x)[2])

S = sparse(S)
#F
F = sparse(F)
#A
A = sparse(A)

start_val = vcat(
    vec(var(Matrix(three_path_dat), dims = 1))./2,
    fill(0.05, 3),
    fill(0.0, 6),
    fill(1.0, 8),
    fill(0, 3)
    )

imply = ImplySymbolic(A, S, F, x, start_val)

imply_alloc = sem.ImplySymbolicAlloc(A, S, F, x, start_val)

model_for = Sem(semobserved, imply_alloc, loss, diff_for)

ForwardDiff.gradient(model_for, start_val)[1]

D = model_for.observed.obs_cov

size(start_val)

start_val

function mygrad(y, D)

    S =[y[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
        0     y[2]  0     0     0     0     0     0     0     0     0     0     0     0
        0     0     y[3]  0     0     0     0     0     0     0     0     0     0     0
        0     0     0     y[4]  0     0     0     y[15] 0     0     0     0     0     0
        0     0     0     0     y[5]  0     y[16] 0     y[17] 0     0     0     0     0
        0     0     0     0     0     y[6]  0     0     0     y[18] 0     0     0     0
        0     0     0     0     y[16] 0     y[7]  0     0     0     y[19] 0     0     0
        0     0     0     y[15] 0     0     0     y[8]  0     0     0     0     0     0
        0     0     0     0     y[17] 0     0     0     y[9]  0     y[20] 0     0     0
        0     0     0     0     0     y[18] 0     0     0     y[10] 0     0     0     0
        0     0     0     0     0     0     y[19] 0     y[20] 0     y[11] 0     0     0
        0     0     0     0     0     0     0     0     0     0     0     y[12] 0     0
        0     0     0     0     0     0     0     0     0     0     0     0     y[13] 0
        0     0     0     0     0     0     0     0     0     0     0     0     0     y[14]]

    F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
        0 1 0 0 0 0 0 0 0 0 0 0 0 0
        0 0 1 0 0 0 0 0 0 0 0 0 0 0
        0 0 0 1 0 0 0 0 0 0 0 0 0 0
        0 0 0 0 1 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 1 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 1 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 1 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 1 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 1 0 0 0 0
        0 0 0 0 0 0 0 0 0 0 1 0 0 0]

    A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
        0  0  0  0  0  0  0  0  0  0  0     y[21] 0     0
        0  0  0  0  0  0  0  0  0  0  0     y[22] 0     0
        0  0  0  0  0  0  0  0  0  0  0     0     1     0
        0  0  0  0  0  0  0  0  0  0  0     0     y[23] 0
        0  0  0  0  0  0  0  0  0  0  0     0     y[24] 0
        0  0  0  0  0  0  0  0  0  0  0     0     y[25] 0
        0  0  0  0  0  0  0  0  0  0  0     0     0     1
        0  0  0  0  0  0  0  0  0  0  0     0     0     y[26]
        0  0  0  0  0  0  0  0  0  0  0     0     0     y[27]
        0  0  0  0  0  0  0  0  0  0  0     0     0     y[28]
        0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     y[29] 0     0
        0  0  0  0  0  0  0  0  0  0  0     y[30] y[31] 0]

    B = inv(I-A)
    E = B*S*B'
    Σ_inv = inv(F*E*F')

    parnum = size(y, 1)
    d = zeros(parnum)
    gradvec = zeros(parnum)
    for i = 1:parnum
        d .= 0.0
        d[i] = 1.0

        S_der =[d[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
                0     d[2]  0     0     0     0     0     0     0     0     0     0     0     0
                0     0     d[3]  0     0     0     0     0     0     0     0     0     0     0
                0     0     0     d[4]  0     0     0     d[15] 0     0     0     0     0     0
                0     0     0     0     d[5]  0     d[16] 0     d[17] 0     0     0     0     0
                0     0     0     0     0     d[6]  0     0     0     d[18] 0     0     0     0
                0     0     0     0     d[16] 0     d[7]  0     0     0     d[19] 0     0     0
                0     0     0     d[15] 0     0     0     d[8]  0     0     0     0     0     0
                0     0     0     0     d[17] 0     0     0     d[9]  0     d[20] 0     0     0
                0     0     0     0     0     d[18] 0     0     0     d[10] 0     0     0     0
                0     0     0     0     0     0     d[19] 0     d[20] 0     d[11] 0     0     0
                0     0     0     0     0     0     0     0     0     0     0     d[12] 0     0
                0     0     0     0     0     0     0     0     0     0     0     0     d[13] 0
                0     0     0     0     0     0     0     0     0     0     0     0     0     d[14]]

        A_der =[0  0  0  0  0  0  0  0  0  0  0     0     0     0
            0  0  0  0  0  0  0  0  0  0  0     d[21] 0     0
            0  0  0  0  0  0  0  0  0  0  0     d[22] 0     0
            0  0  0  0  0  0  0  0  0  0  0     0     0     0
            0  0  0  0  0  0  0  0  0  0  0     0     d[23] 0
            0  0  0  0  0  0  0  0  0  0  0     0     d[24] 0
            0  0  0  0  0  0  0  0  0  0  0     0     d[25] 0
            0  0  0  0  0  0  0  0  0  0  0     0     0     0
            0  0  0  0  0  0  0  0  0  0  0     0     0     d[26]
            0  0  0  0  0  0  0  0  0  0  0     0     0     d[27]
            0  0  0  0  0  0  0  0  0  0  0     0     0     d[28]
            0  0  0  0  0  0  0  0  0  0  0     0     0     0
            0  0  0  0  0  0  0  0  0  0  0     d[29] 0     0
            0  0  0  0  0  0  0  0  0  0  0     d[30] d[31] 0]

        term = F*B*A_der*E*F'

        Σ_der =  F*B*S_der*B'F' + term + term'

        gradvec[i] = tr(Σ_inv*Σ_der) + tr((-Σ_inv)*Σ_der*Σ_inv*D)
    end

    return gradvec
end

@btime mygrad(start_val, D)

all(mygrad(start_val, D) .≈ ForwardDiff.gradient(model_for, start_val))

@btime model_fin(start_val)

In = [1, 4, 3, 5]; J = [4, 7, 18, 9]; V = [1, 2, -5, 3];

S = sparse(In,J,V)

S[1,4] = 0.0

parnum = size(start_val, 1)
d = zeros(parnum)
gradvec = zeros(parnum)

S_ind_vec = Vector{Tuple{Array{Int64,1},Array{Int64,1},Array{Float64,1}}}()
A_ind_vec = Vector{Tuple{Array{Int64,1},Array{Int64,1},Array{Float64,1}}}()

for i = 1:parnum
    d .= 0.0
    d[i] = 1.0

    S_der =[d[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
            0     d[2]  0     0     0     0     0     0     0     0     0     0     0     0
            0     0     d[3]  0     0     0     0     0     0     0     0     0     0     0
            0     0     0     d[4]  0     0     0     d[15] 0     0     0     0     0     0
            0     0     0     0     d[5]  0     d[16] 0     d[17] 0     0     0     0     0
            0     0     0     0     0     d[6]  0     0     0     d[18] 0     0     0     0
            0     0     0     0     d[16] 0     d[7]  0     0     0     d[19] 0     0     0
            0     0     0     d[15] 0     0     0     d[8]  0     0     0     0     0     0
            0     0     0     0     d[17] 0     0     0     d[9]  0     d[20] 0     0     0
            0     0     0     0     0     d[18] 0     0     0     d[10] 0     0     0     0
            0     0     0     0     0     0     d[19] 0     d[20] 0     d[11] 0     0     0
            0     0     0     0     0     0     0     0     0     0     0     d[12] 0     0
            0     0     0     0     0     0     0     0     0     0     0     0     d[13] 0
            0     0     0     0     0     0     0     0     0     0     0     0     0     d[14]]

    A_der =[0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     d[21] 0     0
        0  0  0  0  0  0  0  0  0  0  0     d[22] 0     0
        0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     0     d[23] 0
        0  0  0  0  0  0  0  0  0  0  0     0     d[24] 0
        0  0  0  0  0  0  0  0  0  0  0     0     d[25] 0
        0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     0     0     d[26]
        0  0  0  0  0  0  0  0  0  0  0     0     0     d[27]
        0  0  0  0  0  0  0  0  0  0  0     0     0     d[28]
        0  0  0  0  0  0  0  0  0  0  0     0     0     0
        0  0  0  0  0  0  0  0  0  0  0     d[29] 0     0
        0  0  0  0  0  0  0  0  0  0  0     d[30] d[31] 0]

    S_der = sparse(S_der)
    A_der = sparse(A_der)

    S_ind = findnz(S_der)
    A_ind = findnz(A_der)

    push!(S_ind_vec, S_ind)
    push!(A_ind_vec, A_ind)
end

test = sparse(S_ind_vec[11]..., matsize...)

matsize = size(A)

function mygrad2(y, D, A, S, F, S_ind_vec, A_ind_vec)

    fun_S(S, y)

    fun_A(A, y)

    B = inv(I-A)
    E = B*S*B'
    Σ_inv = inv(F*E*F')

    parnum = size(y, 1)
    d = zeros(parnum)
    gradvec = zeros(parnum)
    for i = 1:parnum
        S_der = sparse(S_ind_vec[i]..., matsize...)
        A_der = sparse(A_ind_vec[i]..., matsize...)

        term = F*B*A_der*E*F'

        Σ_der =  F*B*S_der*B'F' + term + term'

        gradvec[i] = tr(Σ_inv*Σ_der) + tr((-Σ_inv)*Σ_der*Σ_inv*D)
    end

    return gradvec
end

A_pre = zeros(14,14)
S_pre = zeros(14,14)

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

all(mygrad2(start_val, D, A_pre, S_pre, F, S_ind_vec, A_ind_vec) .≈
    ForwardDiff.gradient(model_for, start_val))

testfun2 = eval(ModelingToolkit.build_function(A, x)[2])

In = [1, 4, 3, 5]; J = [4, 7, 18, 9]; V = [1, 2, -5, 3];
mat = sparse(In,J,V)

testfun2(mat, start_val)

testfun2()

## optim test

function fg!(F,G,x)
  if G != nothing
    G .= nest1()
  end
  if F != nothing
    return res
  end
end

struct nest1{grad}
    G::grad
end

struct nest2{grad}
    G::grad
end

mynest1 = nest1([0.0])

mynest2 = nest2([0.0])

function (nest1::nest1)(F, G, x)
    res = x[1]^2
    if G != nothing
      ForwardDiff.gradient!(nest1.G, x -> x[1]^2, x)
    end
    if F != nothing
      return res
    end
end

function (nest2::nest2)(F, G, x)
    res = x[1]^2
    if G != nothing
      ForwardDiff.gradient!(nest2.G, x -> x[1]^2, x)
    end
    if F != nothing
      return res
    end
end

function nestedfg!(F, G, x)
    if G != nothing
        mynest1(F, G, x)
        mynest2(F, G, x)
        G .= mynest1.G + mynest2.G
    end
    if F != nothing
        return mynest1(F, G, x) + mynest2(F, G, x)
    end
end

Optim.optimize(Optim.only_fg!(nestedfg!), [2.0], Optim.LBFGS())

ForwardDiff.gradient(model_for, fill(0, 31))

function testfun(mat)
    let A, B = mat, mat
    end
end

mat = rand(10,10)

@btime testfun(mat)

## test structs

sem.SemAnalyticDiff(LBFGS(), Optim.Options(), A, S, F, x, start_val)


text:
function [D,J,X]=JacobiTri(ksi,F,Grid)
dX1dksi1=(Grid.Nodes(F.N(2)).P(1)-Grid.Nodes(F.N(1)).P(1));
dX1dksi2=(Grid.Nodes(F.N(3)).P(1)-Grid.Nodes(F.N(1)).P(1));
dX2dksi1=(Grid.Nodes(F.N(2)).P(2)-Grid.Nodes(F.N(1)).P(2));
dX2dksi2=(Grid.Nodes(F.N(3)).P(2)-Grid.Nodes(F.N(1)).P(2));
dX3dksi1=(Grid.Nodes(F.N(2)).P(3)-Grid.Nodes(F.N(1)).P(3));
dX3dksi2=(Grid.Nodes(F.N(3)).P(3)-Grid.Nodes(F.N(1)).P(3));
J=zeros(3,2);

J(1,1)=dX1dksi1;
J(2,1)=dX2dksi1;
J(3,1)=dX3dksi1;
J(1,2)=dX1dksi2;
J(2,2)=dX2dksi2;
J(3,2)=dX3dksi2;

D=norm(cross(J(:,1),J(:,2)),2);
D=J(1,1)*J(2,2)-J(2,1)*J(1,2);
if nargout > 2
  X=zeros(1,3);
  X(1)=Grid.Nodes(F.N(1)).P(1)...
    +(Grid.Nodes(F.N(2)).P(1)-Grid.Nodes(F.N(1)).P(1))*ksi(1)...
    +(Grid.Nodes(F.N(3)).P(1)-Grid.Nodes(F.N(1)).P(1))*ksi(2);

  X(2)=Grid.Nodes(F.N(1)).P(2)...
    +(Grid.Nodes(F.N(2)).P(2)-Grid.Nodes(F.N(1)).P(2))*ksi(1)...
    +(Grid.Nodes(F.N(3)).P(2)-Grid.Nodes(F.N(1)).P(2))*ksi(2);
  X(3)=Grid.Nodes(F.N(1)).P(3)...
    +(Grid.Nodes(F.N(2)).P(3)-Grid.Nodes(F.N(1)).P(3))*ksi(1)...
    +(Grid.Nodes(F.N(3)).P(3)-Grid.Nodes(F.N(1)).P(3))*ksi(2);
end
end


text:
## connect to NLopt as backend

function (model::Sem{A, B, C, D} where {A, B, C, D <: SemFiniteDiff})(par::Vector, grad::Vector)
    if length(grad) > 0
        FiniteDiff.finite_difference_gradient!(grad, model, par)
    end
    return model(par)
end

function grad_nlopt(model, par, grad)
    if length(grad) > 0
        model(par, grad)
    end
    return model(par)
end

function sem_fit_nlopt(
        model::Sem{O, I, L, D};
        ftol_rel = 1e-10,
        xtol_rel = 1.5e-8,
        lower = nothing,
        upper = nothing,
        local_algo = nothing,
        maxeval = 200) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemDiff}

    opt = NLopt.Opt(model.diff.algorithm, length(model.imply.start_val))
    #cache = FiniteDiff.GradientCache(start)
    opt.min_objective = (x,y) -> grad_nlopt(model, x, y)
    opt.ftol_rel = ftol_rel
    opt.xtol_rel = xtol_rel
    opt.maxeval = maxeval
    !isnothing(lower) ? opt.lower_bounds = lower : nothing
    !isnothing(upper) ? opt.upper_bounds = upper : nothing
    !isnothing(local_algo) ? opt.local_optimizer = NLopt.Opt(local_algo, length(model.imply.start_val)) : nothing
    result = NLopt.optimize(opt, model.imply.start_val)

    return result
end


text:
using CSV
using LsqFit
using Dates

function calc(dates::Array{Date}, data::Array{Float64}, n::Int64)

    function j_m_exp(t,p)
        J = Array{Float64}(undef, length(t), length(p))
        J[:,1] = exp.(p[2] .* t)       #df/dp[1]
        J[:,2] = t .* p[1] .* J[:,1]   #df/dp[2]
        J
    end

    function j_m_pow2(t,p)
        J = Array{Float64}(undef, length(t), length(p))
        J[:,1] = 2.0 .^ (t ./ p[2])       #df/dp[1]

        c = -log(2.0)*p[1]/(p[2] * p[2])
        J[:,2] =  t .* c .* J[:,1]        #df/dp[2]
        J
    end

    datesWindow = dates[end-n+1 : end]
    date1 = dates[end-n+1]

    x = map(d -> (d-date1).value, datesWindow)
    y = data[end-n+1 : end] .- data[end-n+1];

    m_exp(t, p) = p[1] * exp.(p[2] * t)
    p0_exp = [x[1], 0.5]
    fit_exp = curve_fit(m_exp, j_m_exp, x, y, p0_exp)

    m_pow2(t, p) = p[1] * 2.0 .^ (t/p[2])
    p0_pow2 = [x[1], 0.5]
    fit_pow2 = curve_fit(m_pow2, j_m_pow2, x, y, p0_pow2)

    confidence_intervals = confidence_interval(fit_pow2, 0.1)

    p = fit_pow2.param

    # y_fitted = m_pow2(x,p)
    # diff = 100.0*(abs.(y.-y_fitted))./y
    # println(diff)

    print(p)
    print("\t")
    println(confidence_intervals)
    println("---")

end

df = CSV.File("CoronaData2.csv")

dates = df.Date

dataDE = convert(Array{Float64}, df.Germany)
dataCH = convert(Array{Float64}, df.Switzerland)

print("Germany: ")
# calc(dates[1:end], dataDE, 15)
calc(dates[1:end], dataDE, 10)
# calc(dates[1:end], dataDE, 5)

print("Switzerland: ")
# calc(dates[1:end], dataCH, 15)
calc(dates[1:end], dataCH, 10)
# calc(dates[1:end], dataCH, 5)


text:
function op!(Op::ITensor, ::OpName"Id", ::SiteType"Generic", s::Index)
  for n in 1:dim(s)
    Op[n, n] = 1.0
  end
end

op(::OpName"I", t::SiteType"Generic", s::Index) = op("Id", s)

function op!(Op::ITensor, ::OpName"F", st::SiteType"Generic", s::Index)
  return op!(Op, OpName("Id"), st, s)
end


text:
function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    solution_info = model.inner.solution_info
    if solution_info === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    solver_status = solution_info.solver_status
    model_status = solution_info.model_status

    if solver_status == NORMAL_COMPLETION
        if model_status == OPTIMAL
            return MOI.OPTIMAL
        elseif model_status == INFEASIBLE
            return MOI.INFEASIBLE
        elseif model_status == UNBOUNDED
            return MOI.DUAL_INFEASIBLE
        elseif model_status == INTERMEDIATE_FEASIBLE
            return LOCALLY_SOLVED
        elseif model_status == UNKNOWN
            return MOI.OTHER_ERROR
        end
    elseif solver_status == INSUFFICIENT_MEMORY_FOR_NODES
        return MOI.MEMORY_LIMIT
    elseif solver_status == ITERATION_LIMIT
        return MOI.ITERATION_LIMIT
    elseif solver_status == TIME_LIMIT
        return MOI.TIME_LIMIT
    elseif solver_status == NUMERICAL_SENSITIVITY
        return MOI.NUMERICAL_ERROR
    elseif solver_status == INSUFFICIENT_MEMORY_FOR_SETUP
        return MOI.MEMORY_LIMIT
    elseif solver_status == RESERVED
        return MOI.OTHER_ERROR
    elseif solver_status == TERMINATED_BY_BARON
        return MOI.OTHER_ERROR
    elseif solver_status == SYNTAX_ERROR
        return MOI.INVALID_MODEL
    elseif solver_status == LICENSING_ERROR
        return MOI.OTHER_ERROR
    elseif solver_status == USER_HEURISTIC_TERMINATION
        return MOI.OTHER_LIMIT
    end

    error("Unrecognized Baron status: $solver_status, $model_status")
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return (model.inner.solution_info.feasible_point === nothing) ? 0 : 1
end

function MOI.get(model::Optimizer, ::MOI.DualStatus)
    MOI.NO_SOLUTION
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if attr.result_index != 1
        return MOI.NO_SOLUTION
    end
    solution_info = model.inner.solution_info
    if solution_info === nothing || solution_info.feasible_point === nothing
        return MOI.NO_SOLUTION
    else
        return solution_info.model_status == UNBOUNDED ? MOI.INFEASIBILITY_CERTIFICATE : MOI.FEASIBLE_POINT
    end
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    # if model.inner.solution_info.model_status == UNBOUNDED
    #     # BARON will set the same large number
    #     # for both abj and variables in case of unbounded
    #     # this would make it retur MOI consistent values
    #     # but getters are no define and NL one would be messy
    #     return MOIU.get_fallback(model, attr)
    # else
        model.inner.solution_info.objective_value
    # end
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(model, attr)
    # MOI.throw_if_not_valid(model, vi) # TODO implement is_valid
    solution_info = model.inner.solution_info
    if solution_info === nothing || solution_info.feasible_point === nothing
        error("VariablePrimal not available.")
    end
    check_variable_indices(model, vi)
    return solution_info.feasible_point[vi.value]
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveBound)
    solution_info = model.inner.solution_info
    return solution_info.dual_bound
end

function MOI.get(model::Optimizer, ::MOI.SolveTimeSec)
    solution_info = model.inner.solution_info
    return solution_info.wall_time
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    info = model.inner.solution_info
    return "solver: $(info.solver_status), model: $(info.model_status)"
end

# TODO: desirable?
function MOI.get(model::MOIU.CachingOptimizer{BARON.Optimizer}, attr::MOI.ConstraintPrimal, ci::MOI.ConstraintIndex)
    return MOIU.get_fallback(model, attr, ci)
end

function MOI.supports(::Optimizer, ::MOI.SolverVersion)
    false
end


text:
using NestedEnvironments
using DifferentialEquations
using Transducers

using Random
using Plots
using Test

function dynamics(env::NestedEnvironments.InputAffineQuadraticCostEnv)
    return function (x, p, t)
        u = command(env, x)
        NestedEnvironments.ẋ(env, x, t, u)
    end
end
command(env, x) = NestedEnvironments.u_optimal(env, x)
NestedEnvironments.initial_condition(env::NestedEnvironments.InputAffineQuadraticCostEnv) = 2*(rand(2) .- 0.5)
# register envs
__env = NestedEnvironments.InputAffineQuadraticCostEnv()
__x0 = NestedEnvironments.initial_condition(__env)
@reg_env __env __x0

function test()
    env = NestedEnvironments.InputAffineQuadraticCostEnv()
    x0 = NestedEnvironments.initial_condition(env)
    @show x0
    @show _x0 = @raw(env, x0)
    @test _x0 == @raw x0
    @test x0 == @readable(env, _x0)
    t0 = 0.0
    tf = 10.0
    Δt = 0.01
    tspan = (t0, tf)
    ts = t0:Δt:tf
    prob = ODEProblem(env, dynamics(env), x0, tspan)
    @time sol = solve(prob, Tsit5(), saveat=ts)
    xs = sol.u |> Map(u -> @readable u) |> collect
    @show xs
end


text:
export hmm_generate

import ..HMMModel

using mlpack._Internal.io

import mlpack_jll
const hmm_generateLibrary = mlpack_jll.libmlpack_julia_hmm_generate

# Call the C binding of the mlpack hmm_generate binding.
function hmm_generate_mlpackMain()
  success = ccall((:hmm_generate, hmm_generateLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module hmm_generate_internal
  import ..hmm_generateLibrary

import ...HMMModel

# Get the value of a model pointer parameter of type HMMModel.
function IOGetParamHMMModel(paramName::String)::HMMModel
  HMMModel(ccall((:IO_GetParamHMMModelPtr, hmm_generateLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type HMMModel.
function IOSetParamHMMModel(paramName::String, model::HMMModel)
  ccall((:IO_SetParamHMMModelPtr, hmm_generateLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeHMMModel(stream::IO, model::HMMModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeHMMModelPtr, hmm_generateLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeHMMModel(stream::IO)::HMMModel
  buffer = read(stream)
  HMMModel(ccall((:DeserializeHMMModelPtr, hmm_generateLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    hmm_generate(length, model; [seed, start_state, verbose])

This utility takes an already-trained HMM, specified as the `model` parameter,
and generates a random observation sequence and hidden state sequence based on
its parameters. The observation sequence may be saved with the `output` output
parameter, and the internal state  sequence may be saved with the `state` output
parameter.

The state to start the sequence in may be specified with the `start_state`
parameter.

For example, to generate a sequence of length 150 from the HMM `hmm` and save
the observation sequence to `observations` and the hidden state sequence to
`states`, the following command may be used:

```julia
julia> observations, states = hmm_generate(150, hmm)
```

# Arguments

 - `length::Int`: Length of sequence to generate.
 - `model::HMMModel`: Trained HMM to generate sequences with.
 - `seed::Int`: Random seed.  If 0, 'std::time(NULL)' is used.  Default
      value `0`.

 - `start_state::Int`: Starting state of sequence.  Default value `0`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Float64, 2}`: Matrix to save observation sequence to.
 - `state::Array{Int, 2}`: Matrix to save hidden state sequence to.

"""
function hmm_generate(length::Int,
                      model::HMMModel;
                      seed::Union{Int, Missing} = missing,
                      start_state::Union{Int, Missing} = missing,
                      verbose::Union{Bool, Missing} = missing,
                      points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, hmm_generateLibrary), Nothing, ());

  IORestoreSettings("Hidden Markov Model (HMM) Sequence Generator")

  # Process each input argument before calling mlpackMain().
  IOSetParam("length", length)
  hmm_generate_internal.IOSetParamHMMModel("model", convert(HMMModel, model))
  if !ismissing(seed)
    IOSetParam("seed", convert(Int, seed))
  end
  if !ismissing(start_state)
    IOSetParam("start_state", convert(Int, start_state))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  IOSetPassed("state")
  # Call the program.
  hmm_generate_mlpackMain()

  return IOGetParamMat("output", points_are_rows),
         IOGetParamUMat("state", points_are_rows)
end