text:
module QRDowndate
# copyright (c) 2017 Ralph Smith
# SPDX-License-Identifier: MIT-Expat
# License-Filename: LICENSE.md

export qrdowndate, qrdowndate!
using LinearAlgebra

# classical Gramm-Schmidt
function gs_step(U,::Type{Val{:CGS}},v)
    g = U' * v
    r = v - U * g
    β = norm(r,2)
    rmul!(r,1/β)
    r, g, β
end

# modified Gramm-Schmidt
function gs_step(U,::Type{Val{:MGS}},v)
    n = size(U,2)
    r = copy(v)
    g = zeros(eltype(U),n)
    for k=1:n
        # paper has a temp vector q; we just use r here
        ζ = dot(view(U,:,k),r)
        r = r - ζ * view(U,:,k)
        g[k] = ζ
    end
    β = norm(r,2)
    rmul!(r,1/β)
    r, g, β
end

# backwards modified Gramm-Schmidt
function gs_step(U,::Type{Val{:BMGS}},v)
    n = size(U,2)
    r = copy(v)
    g = zeros(eltype(U),n)
    for k=n:-1:1
        ζ = dot(view(U,:,k),r)
        r = r - ζ * view(U,:,k)
        g[k] = ζ
    end
    β = norm(r,2)
    rmul!(r,1/β)
    r, g, β
end

"""
Given `U[m,n]`, find unit vector `u0[m]` and vector `d[n+1]` such that
`U' * u0 = 0` and `hcat(u0,U) * d = e1`.

Algorithm from [^Barlow2005]

[Barlow2005]: Barlow et al., BIT 45, 259 (2005)
"""
function getorth(U::AbstractMatrix{T},M1,M2) where T
    m,n = size(U)
    @assert m >= n
    e1 = zeros(T,m); e1[1] = 1
    v1,g1,β1 = gs_step(U,M1,e1)
    βmin = sqrt(1/2)
    if β1 == 0
        β2 = α = zero(T)
        f = copy(g1)
        d = vcat(zeros(eltype(g1),1),g1)
        branch = 10
    else
        v2,g2,β2 = gs_step(U,M2,v1)
        if β1 >= βmin
            f = copy(g1)
            α = β1
            # β2 = zero(T) # early guess since paper fails to specify
            u0 = copy(v1)
            d = vcat([α],g1)
            branch = 20
        else
            f = g1 + β1 * g2
            α = β1 * β2
            if β2 >= βmin
                u0 = v2
                d = vcat([α],f)
            end
            branch = 30
        end
    end
    if β2 < βmin
        # FIXME: inelegant
        cval = Inf
        k = 0
        for j=1:m
            ek = zeros(m); ek[j] = 1
            tval = norm(U' * ek)
            if tval < cval
                k = j
                cval = tval
            end
        end

        ek = zeros(m); ek[k] = 1
        w1,k1,γ1 = gs_step(U,M1,ek)
        u0,k2,γ2 = gs_step(U,M2,w1)
        if α > 0
            d1 = vcat(zeros(T,1),f)
            Uhat = hcat(u0, U)
            # confusingly this is named $v_2$ in the paper,
            # but is NOT the same as the above v2.
            v2b = -U * f; v2b[1] += one(T)
            rmul!(v2,1/α)
            w2,d2,γ3 = gs_step(Uhat,M1,v2b)
            d = d1 + α * d2
            branch += 1
        else
            branch += 2
        end
    end

    return u0, d, branch
end

"""
    qrdowndate(Q,R) -> Qbar, Rbar

downdate the QR decomposition of a matrix (by the top row).

Given a QR decomposition `X=Q*R`, computes `Qbar, Rbar` such that
`Xbar = Qbar * Rbar` where `Xbar == X[2:end,:]`.
"""
function qrdowndate(U,R)
    Unew,Rnew = copy(U),copy(R)
    qrdowndate!(Unew,Rnew)
    Unew[2:end,:],Rnew
end

"""
    qrdowndate!(Q,R)

downdate the QR decomposition of a matrix.

Given a QR decomposition `X=Q*R`, computes `Qbar, Rbar` such that
`Xbar = Qbar * Rbar` where `Xbar == X[2:end,:]`, overwriting
the bottom of `Q` and `R`.
"""
function qrdowndate!(U::AbstractMatrix{T},R::AbstractMatrix{T}) where T
    m,n = size(U)
    (m > n) || throw(ArgumentError("U must be tall"))
    # check size of R?
    u0, d, b = getorth(U,Val{:MGS},Val{:BMGS})

    ρ = d[1]
    r = d[2:end]
    u = zeros(T,1,n)
    # find V st
    #     V is a product of Givens rotations
    #     V' * d = norm(d)*e1
    #     V' * vcat(zeros(1,n),R) = vcat(z',Rbar) for Up.Tri. Rhat
    # Ubar = hcat(u0,U) * V
    for k=n:-1:1
        G,ρ = givens(ρ,r[k],1,2)
        tmp = vcat(Transpose(u[k:n]),Transpose(R[k,k:n]))
        lmul!(G,tmp)
        u[k:n] = tmp[1,:]
        R[k,k:n] = tmp[2,:]
        tmp2 = hcat(u0,U[:,k])
        rmul!(tmp2,adjoint(G))
        u0 = tmp2[:,1]
        U[:,k] = tmp2[:,2]
    end
    # unceremoniously return indicator of branch used in orthogonalizer
    b
end

end


text:
###
### LongSequence specific specializations of src/biosequence/indexing.jl
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# assumes `i` is positive and `bitsof(A)` is a power of 2

@inline function bitindex(seq::LongSubSeq, i::Integer)
    return bitindex(BitsPerSymbol(seq), encoded_data_eltype(seq), i + first(seq.part) - 1)
end

# More efficient due to copyto!
function Base.getindex(seq::LongSequence, part::UnitRange{<:Integer})
	@boundscheck checkbounds(seq, part)
	newseq = typeof(seq)(length(part))
	return copyto!(newseq, 1, seq, first(part), length(part))
end

# More efficient due to copyto!
function Base.setindex!(seq::SeqOrView{A},
                        other::SeqOrView{A},
                        locs::UnitRange{<:Integer}) where {A <: Alphabet}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    return copyto!(seq, locs.start, other, 1, length(locs))
end

@inline function encoded_setindex!(seq::SeqOrView{A},
				   bin::Unsigned, i::Integer) where {A <: Alphabet}
	return encoded_setindex!(seq, bin, bitindex(seq, i))
end

@inline function encoded_setindex!(s::SeqOrView, v::Unsigned, i::BitIndex)
    vi, off = i
    data = encoded_data(s)
    bits = @inbounds data[vi]
	v_ = v % encoded_data_eltype(s)
    @inbounds data[vi] = (v_ << off) | (bits & ~(bindata_mask(s) << off))
    return s
end

function Base.setindex!(seq::SeqOrView, x, locs::AbstractVector{<:Integer})
	@boundscheck checkbounds(seq, locs)
	unsafe_setindex!(seq, x, locs)
end

function Base.setindex!(seq::SeqOrView, x::BioSequence, locs::AbstractVector{<:Integer})
	@boundscheck checkbounds(seq, locs)
	unsafe_setindex!(seq, x, locs)
end

@inline function unsafe_setindex!(seq::SeqOrView, x, locs::AbstractVector{<:Integer})
	bin = encode(Alphabet(seq), convert(eltype(seq), x))
    for i in locs
        encoded_setindex!(seq, bin, bitindex(seq, i))
    end
    return seq
end

# Only to avoid ambiguity
function Base.setindex!(seq::SeqOrView, x::SeqOrView, locs::AbstractVector{Bool})
	@boundscheck checkbounds(seq, locs)
	checkdimension(x, locs)
	j = 0
	@inbounds for i in eachindex(locs)
		if locs[i]
			j += 1
			seq[i] = x[j]
		end
	end
	return seq
end

# To save encoding.
function unsafe_setindex!(seq::SeqOrView, x, locs::AbstractVector{Bool})
    bin = encode(Alphabet(seq), convert(eltype(seq), x))
    i = j = 0
    while true
        i = findnext(locs, i + 1)
        if i === nothing
            break
        end
        encoded_setindex!(seq, bin, i)
    end
    return seq
end

# To avoid ambiguity errors
function Base.setindex!(seq::SeqOrView{A},
                        other::BioSequence{A},
                        locs::AbstractVector{<:Integer}) where {A <: Alphabet}
    @boundscheck checkbounds(seq, locs)
    checkdimension(other, locs)
    unsafe_setindex!(seq, other, locs)
end

function unsafe_setindex!(seq::SeqOrView{A},
                        other::BioSequence,
                        locs::AbstractVector{<:Integer}) where {A <: Alphabet}
	@inbounds for (i, n) in zip(locs, other)
		seq[i] = n
	end
	return seq
end

ext: jl


text:
# Helper function for the PCA demo

mutable struct pPCAParams
    W::Matrix{Float64}
    m::Vector{Float64}
    v::Float64
end

function readDataSet(filename)
    f = open(filename);
    data = nothing
    for line in eachline(f)
        row = map((s) -> parse(Int, s), split(replace(line, "  " => " "), " "))
        data = (data == nothing) ? row' : vcat(data, row')
    end
    close(f)

    return data
end

function pPCA(X::Matrix, M::Int64)
    # Implement probabilistic PCA using the EM algorithm
    # X collects the data points on its columns
    # X may contain missing values, indicated by NaN

    # Probabilistic model:
    # x[n] = W*z[n] + m + ε, where W is a matrix, m is a vector and ε is a noise term
    # z[n] ~ N(0,I), where I is the identity matrix
    # ε ~ N(0, v*I), where v is a scalar
    N = size(X, 2); D = size(X, 1)
    missing_values = isnan.(X)
    has_missing_values = sum(missing_values) > 0
    # Find row means of X, ignoring missing values
    if has_missing_values
        X[findall(missing_values)] .= 0
        m = ( sum(X, dims=2) ./ (N*ones(D) - sum(missing_values,dims=2)) )[:,1] # [:,1] makes sure this is a vector
        X[findall(missing_values)] .= NaN
        observed_d = [findall(.!missing_values[:,n]) for n=1:N] # observed dimensions for every data vector
        observed_n = [findall(.!missing_values[d,:]) for d=1:D] # observed value indexes per dimension
        x_observed = [X[observed_d[n],n] for n=1:N]  # data vectors without the missing data
    else
        m = (sum(X, dims=2) ./ N)[:,1] # [:,1] makes sure this is a vector
        _X = X .- m # shift data to be zero-mean
    end
    # Initialize variables
    θ = pPCAParams(Matrix{Float64}(undef,D,M), m, 1.0) # θ holds the parameters
    Mz = randn(M,N)              # Mz[:,n] = MEAN[z[n]]
    Σz = [zeros(M,M) for n=1:N]  # Σz[n] = COV[z[n]]

    function eStep!(Mz, Σz)
        # Update Mz and Σz
        W = θ.W; v = θ.v

        if has_missing_values
            for n=1:N
                _W = W[observed_d[n],:]
                _x = x_observed[n] - θ.m[observed_d[n]]
                M_inv = inv(_W'*_W + v*Matrix{Float64}(I, M, M))
                Σz[n] = v * M_inv
                Mz[:,n] = M_inv * _W' * _x
            end
        else

            M_inv = inv(W'*W + v*Matrix{Float64}(I, M, M))
            Mz[:] = M_inv * W' * _X
            Σz[1] = v * M_inv # COV[z[n]] is equal for all n. We only update the one for n=1
        end
    end

    function mStep!(θ)
        # Update W, m, v
        if has_missing_values
            m = zero(θ.m)
            for n=1:N
                m[observed_d[n]] += x_observed[n] .- (θ.W[observed_d[n],:] * Mz[:,n])[:,1]
            end
            m = m ./ (N .- sum(missing_values, dims=2)[:,1])
            θ.m = m

            for d=1:D
                A = inv(Mz[:,observed_n[d]]*Mz[:,observed_n[d]]' + sum(Σz[observed_n[d]]))
                b = zeros(M)
                for n in observed_n[d]
                    b += Mz[:,n]*(X[d,n] .- m[d])
                end
                θ.W[d,:] = ( A * b )'
            end
            W = θ.W

            θ.v = 0.0
            for n=1:N
                for d in observed_d[n]
                    θ.v +=  (X[d,n] .- (W[d,:]'*Mz[:,n])[1] - m[d])^2 + (W[d,:]' * Σz[n] * W[d,:])[1]
                end
            end
            θ.v = θ.v / (N*D-sum(missing_values))

        else
            # no need to update m
            W = _X * Mz' * inv(Mz*Mz' + N*Σz[1])
            θ.W = W
            θ.v = sum((_X - W*Mz).^2)/(N*D) + tr(W*Σz[1]*W')/D
        end
    end

    for i=1:30 # Here should be a decent convergence check, but hey, this also works...
        mStep!(θ)
        eStep!(Mz, Σz)
    end

    return (θ, Mz) # Mz is actually redundant, but we return it for convenient plotting
end

max_stars_repo_name: Yikeru/BMLIP

max_issues_repo_name: Yikeru/BMLIP

max_forks_repo_name: Yikeru/BMLIP


text:
using ChainRulesTestUtils, DifferenceEquations, Distributions, LinearAlgebra, Test, Zygote
using DelimitedFiles
using DiffEqBase
using FiniteDiff: finite_difference_gradient

function joint_likelihood_1(A, B, C, u0, noise, observables, D; kwargs...)
    problem = LinearStateSpaceProblem(A, B, u0, (0, size(observables, 2)); C, observables_noise = D,
                                      noise, observables, kwargs...)
    return solve(problem, DirectIteration()).logpdf
end

# CRTU has problems with generating random MvNormal, so just testing diagonals
function kalman_likelihood(A, B, C, u0, observables, D; kwargs...)
    problem = LinearStateSpaceProblem(A, B, u0, (0, size(observables, 2)); C, observables_noise = D,
                                      u0_prior = MvNormal(u0, diagm(ones(length(u0)))),
                                      noise = nothing, observables, kwargs...)
    return solve(problem).logpdf
end

# Matrices from RBC
A_rbc = [0.9568351489231076 6.209371005755285;
         3.0153731819288737e-18 0.20000000000000007]
B_rbc = reshape([0.0; -0.01], 2, 1) # make sure B is a matrix
C_rbc = [0.09579643002426148 0.6746869652592109; 1.0 0.0]
D_rbc = MvNormal(Diagonal(abs2.([0.1, 0.1])))
u0_rbc = zeros(2)

observables_rbc = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/RBC_observables.csv"),
                          ',')' |> collect
noise_rbc = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/RBC_noise.csv"), ',')' |>
            collect
# Data and Noise
T = 5
observables_rbc = observables_rbc[:, 1:T]
noise_rbc = noise_rbc[:, 1:T]

@testset "basic inference" begin
    prob = LinearStateSpaceProblem(A_rbc, B_rbc, u0_rbc, (0, size(observables_rbc, 2)); C = C_rbc,
                                   observables_noise = D_rbc, noise = noise_rbc,
                                   observables = observables_rbc)
    @inferred LinearStateSpaceProblem(A_rbc, B_rbc, u0_rbc, (0, size(observables_rbc, 2));
                                      C = C_rbc, observables_noise = D_rbc, noise = noise_rbc,
                                      observables = observables_rbc)

    sol = solve(prob)
    @inferred solve(prob)

    DiffEqBase.get_concrete_problem(prob, false)
    @inferred DiffEqBase.get_concrete_problem(prob, false)

    joint_likelihood_1(A_rbc, B_rbc, C_rbc, u0_rbc, noise_rbc, observables_rbc, D_rbc)
    @inferred joint_likelihood_1(A_rbc, B_rbc, C_rbc, u0_rbc, noise_rbc, observables_rbc, D_rbc)
end

@testset "basic kalman inference" begin
    prob = LinearStateSpaceProblem(A_rbc, B_rbc, u0_rbc, (0, size(observables_rbc, 2)); C = C_rbc,
                                   observables_noise = D_rbc, observables = observables_rbc,
                                   u0_prior = MvNormal(u0_rbc, diagm(ones(length(u0_rbc)))))
    @inferred LinearStateSpaceProblem(A_rbc, B_rbc, u0_rbc, (0, size(observables_rbc, 2));
                                      C = C_rbc, observables_noise = D_rbc,
                                      observables = observables_rbc,
                                      u0_prior = MvNormal(u0_rbc, diagm(ones(length(u0_rbc)))))

    sol = solve(prob)
    @inferred solve(prob)

    prob_concrete = DiffEqBase.get_concrete_problem(prob, false)
    @inferred DiffEqBase.get_concrete_problem(prob, false)

    kalman_likelihood(A_rbc, B_rbc, C_rbc, u0_rbc, observables_rbc, D_rbc)
    @inferred kalman_likelihood(A_rbc, B_rbc, C_rbc, u0_rbc, observables_rbc, D_rbc)
end

gradient((args...) -> joint_likelihood_1(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc,
         u0_rbc, noise_rbc)

@testset "linear rbc joint likelihood" begin
    @test joint_likelihood_1(A_rbc, B_rbc, C_rbc, u0_rbc, noise_rbc, observables_rbc, D_rbc) ≈
          -690.9407412360038
    @inferred joint_likelihood_1(A_rbc, B_rbc, C_rbc, u0_rbc, noise_rbc, observables_rbc, D_rbc) #
    gradient((args...) -> joint_likelihood_1(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc,
             u0_rbc, noise_rbc)

    test_rrule(Zygote.ZygoteRuleConfig(),
               (args...) -> joint_likelihood_1(args..., observables_rbc, D_rbc), A_rbc, B_rbc,
               C_rbc, u0_rbc, noise_rbc; rrule_f = rrule_via_ad, check_inferred = false)
end

gradient((args...) -> kalman_likelihood(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc,
         u0_rbc)

@testset "linear rbc kalman likelihood" begin
    @test kalman_likelihood(A_rbc, B_rbc, C_rbc, u0_rbc, observables_rbc, D_rbc) ≈
          -607.3698273765538
    @inferred kalman_likelihood(A_rbc, B_rbc, C_rbc, u0_rbc, observables_rbc, D_rbc) # would this catch inference problems in the solve?
    test_rrule(Zygote.ZygoteRuleConfig(),
               (args...) -> kalman_likelihood(args..., observables_rbc, D_rbc), A_rbc, B_rbc, C_rbc,
               u0_rbc; rrule_f = rrule_via_ad, check_inferred = false)
end

# Load FVGQ data for checks
A_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_A.csv"), ',')
B_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_B.csv"), ',')
C_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_C.csv"), ',')
D_FVGQ = MvNormal(Diagonal(abs2.(ones(6) * 1e-3)))

observables_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations),
                                    "test/data/FVGQ20_observables.csv"), ',')' |> collect

noise_FVGQ = readdlm(joinpath(pkgdir(DifferenceEquations), "test/data/FVGQ20_noise.csv"), ',')' |>
             collect
u0_FVGQ = zeros(size(A_FVGQ, 1))

@testset "linear FVGQ joint likelihood" begin
    @test joint_likelihood_1(A_FVGQ, B_FVGQ, C_FVGQ, u0_FVGQ, noise_FVGQ, observables_FVGQ,
                             D_FVGQ) ≈ -1.4648817357717388e9
    @inferred joint_likelihood_1(A_FVGQ, B_FVGQ, C_FVGQ, u0_FVGQ, noise_FVGQ, observables_FVGQ,
                                 D_FVGQ)
    test_rrule(Zygote.ZygoteRuleConfig(),
               (args...) -> joint_likelihood_1(args..., observables_FVGQ, D_FVGQ), A_FVGQ, B_FVGQ,
               C_FVGQ, u0_FVGQ, noise_FVGQ; rrule_f = rrule_via_ad, check_inferred = false)
end

@testset "linear FVGQ Kalman" begin
    # Note: set rtol to be higher than the default case because of huge gradient numbers
    @test kalman_likelihood(A_FVGQ, B_FVGQ, C_FVGQ, u0_FVGQ, observables_FVGQ, D_FVGQ) ≈
          -108.52706300389917
    gradient((args...) -> kalman_likelihood(args..., observables_FVGQ, D_FVGQ), A_FVGQ, B_FVGQ,
             C_FVGQ, u0_FVGQ)

    # TODO: this is not turned on because the numbers explode.  Need better unit test data to be interior
    # test_rrule(Zygote.ZygoteRuleConfig(), (args...) -> kalman_likelihood(args..., observables, D),
    #            A, B, C, u0; rrule_f = rrule_via_ad, check_inferred = false, rtol = 1e-5)
end


text:
using TestFunctionRunner
TestFunctionRunner.@run


text:
# https://storopoli.io/Bayesian-Julia/pages/4_Turing/

using Turing
using Statistics: mean, std
using Random:seed!
seed!(123)

@model linreg(X, y; predictors=size(X, 2)) = begin
    #priors
    α ~ Normal(mean(y), 2.5 * std(y))
    β ~ filldist(TDist(3), predictors)
    σ ~ Exponential(1)

    #likelihood
    y ~ MvNormal(α .+ X * β, σ)
end;

using DataFrames, CSV, HTTP

url = "https://raw.githubusercontent.com/storopoli/Bayesian-Julia/master/datasets/kidiq.csv"
kidiq = CSV.read(HTTP.get(url).body, DataFrame)
describe(kidiq)

X = Matrix(select(kidiq, Not(:kid_score)))
y = kidiq[:, :kid_score]
model = linreg(X, y);

chain = sample(model, NUTS(), MCMCThreads(), 2_000, 4)
# typeof(chain)
summarystats(chain)
quantile(chain)


text:
################################################################################################
# core. The Core Functionality
################################################################################################

# 1.  Basic Structures
#     Point
#     Size
#     Rect
#     RotatedRect
#     TermCriteria
#     TermCriteria::TermCriteria
#     Vec
#     Scalar
#     Range
#     Mat::Mat
#     Mat::~Mat
#     Mat::operator =
#     Mat::row
#     Mat::col
#     Mat::rowRange
#     Mat::colRange
#     Mat::diag
#     Mat::clone
#     Mat::copyTo
#     Mat::convertTo
#     Mat::assignTo
#     Mat::setTo
#     Mat::reshape
#     Mat::t
#     Mat::inv
#     Mat::mul
#     Mat::cross
#     Mat::dot
#     Mat::zeros
#     Mat::ones
#     Mat::eye
#     Mat::create
#     Mat::addref
#     Mat::release
#     Mat::resize
#     Mat::reserve
#     Mat::push_back
#     Mat::pop_back
#     Mat::locateROI
#     Mat::adjustROI
#     Mat::operator()
#     Mat::total
#     Mat::isContinuous
#     Mat::elemSize
#     Mat::elemSize1
#     Mat::type
#     Mat::depth
#     Mat::channels
#     Mat::step1
#     Mat::size
#     Mat::empty
#     Mat::ptr
#     Mat::at
#     Mat::begin
#     Mat::end
#     Mat::forEach
#     Mat_
#     InputArray
#     OutputArray
#     MatIterator

# 2. Operations on Arrays
#     abs
#     absdiff
#     add
#     addWeighted
#     bitwise_and
#     bitwise_not
#     bitwise_or
#     bitwise_xor
#     calcCovarMatrix
#     cartToPolar
#     checkRange
#     compare
#     completeSymm
#     convertScaleAbs
#     countNonZero
#     cvarrToMat
#     dct
#     dft
#     divide
#     determinant
#     eigen
#     exp
#     extractImageCOI
#     insertImageCOI
#     flip
#     gemm
#     getOptimalDFTSize
#     idct
#     idft
#     inRange
#     invert
#     log
#     LUT
#     magnitude
#     Mahalanobis
#     max
#     mean
#     meanStdDev
#     merge
#     min
#     minMaxIdx
#     minMaxLoc
#     mixChannels
#     mulSpectrums
#     multiply
#     mulTransposed
#     norm
#     normalize
#     PCA
#     PCA::PCA
#     PCA::operator ()
#     PCA::project
#     PCA::backProject
#     perspectiveTransform
#     phase
#     polarToCart
#     pow
#     RNG
#     RNG::RNG
#     RNG::next
#     RNG::operator T
#     RNG::operator ()
#     RNG::uniform
#     RNG::gaussian
#     RNG::fill
#     randu
#     randn
#     randShuffle
#     reduce
#     repeat
#     scaleAdd
#     setIdentity
#     solve
#     solveCubic
#     solvePoly
#     sort
#     sortIdx
#     split
#     sqrt
#     subtract
#     SVD
#     SVD::SVD
#     SVD::operator ()
#     SVD::compute
#     SVD::solveZ
#     SVD::backSubst
#     sum
#     theRNG
#     trace
#     transform
#     transpose
#     borderInterpolate
#     copyMakeBorder

# 3.  Clustering
#     kmeans
#     partition

# 4.  Utility and System Functions and Macros
#     alignPtr
#     alignSize
#     allocate
#     deallocate
#     fastAtan2
#     cubeRoot
#     Ceil
#     Floor
#     Round
#     IsInf
#     IsNaN
#     CV_Assert
#     error
#     Exception
#     fastMalloc
#     fastFree
#     format
#     getBuildInformation
#     checkHardwareSupport
#     getNumberOfCPUs
#     getNumThreads
#     getThreadNum
#     getTickCount
#     getTickFrequency
#     getCPUTickCount
#     saturate_cast
#     setNumThreads
#     setUseOptimized
#     useOptimized

# import base functions to explicitly extend
import Base.copy,
       Base.convert,
       Base.split

#-------------------------------------------------------------------------------------------------------------------#
# 1. Basic structures

cxx"""
// UInt8 => Cuchar
cv::Vec2b vec2b(uchar a, uchar b) { return cv::Vec2b(a,b); }
cv::Vec3b vec3b(uchar a, uchar b, uchar c) { return cv::Vec3b(a,b,c); }
cv::Vec4b vec4b(uchar a, uchar b, uchar c, uchar d) { return cv::Vec4b(a,b,c,d); }

// Int32 => Cint
cv::Vec2i vec2i(int a, int b) { return cv::Vec2i(a,b); }
cv::Vec3i vec3i(int a, int b, int c) { return cv::Vec3i(a,b,c); }
cv::Vec4i vec4i(int a, int b, int c, int d) { return cv::Vec4i(a,b,c,d); }

// Float32
cv::Vec2f vec2f(float a, float b) { return cv::Vec2f(a,b); }
cv::Vec3f vec3f(float a, float b, float c) { return cv::Vec3f(a,b,c); }
cv::Vec4f vec4f(float a, float b, float c, float d) { return cv::Vec4f(a,b,c,d); }
cv::Vec6f vec6f(float a, float b, float c, float d, float e, float f) { return cv::Vec6f(a,b,c,d,e,f); }

// Float64 => Cdouble
cv::Vec2d vec2d(double a, double b) { return cv::Vec2d(a,b); }
cv::Vec3d vec3d(double a, double b, double c) { return cv::Vec3d(a,b,c); }
cv::Vec4d vec4d(double a, double b, double c, double d) { return cv::Vec4d(a,b,c,d); }
cv::Vec6d vec6d(double a, double b, double c, double d, double e, double f) { return cv::Vec6d(a,b,c,d,e,f); }
"""

vec2b(a,b) = @cxx vec2b(a, b)
vec3b(a,b,c) = @cxx vec3b(a, b, c)
vec4b(a,b,c,d) = @cxx vec4b(a, b, c, d)

vec2i(a,b) = @cxx vec2i(cint(a), cint(b))
vec3i(a,b,c) = @cxx vec3i(cint(a), cint(b), cint(c))
vec4i(a,b,c,d) = @cxx vec4i(cint(a), cint(b), cint(c), cint(d))

vec2f(a,b) = @cxx vec2f(float32(a), float32(b))
vec3f(a,b,c) = @cxx vec3f(float32(a), float32(b), float32(c))
vec4f(a,b,c,d) = @cxx vec4f(float32(a), float32(b), float32(c), float32(d))
vec6f(a,b,c,d,e,f) = @cxx vec6f(float32(a),float32(b), float32(c), float32(d), float32(e), float32(f))

vec2d(a,b) = @cxx vec2d(float(a), float(b))
vec3d(a,b,c) = @cxx vec3d(float(a), float(b), float(c))
vec4d(a,b,c,d) = @cxx vec4d(float(a), float(b), float(c), float(d))
vec6d(a,b,c,d,e,f) = @cxx vec6d(float(a), float(b), float(c), float(d), float(e), float(f))

# Point
cvPoint(x, y) = @cxx cv::Point(x,y)
cvPoint2f(x, y) = @cxx cv::Point2f(float32(x),float32(y))
cvPoint2d(x, y) = @cxx cv::Point2d(float(x),float(y))
cvPoint3i(x, y, z) = @cxx cv::Point3f(cint(x),cint(y),cint(z))
cvPoint3f(x, y, z) = @cxx cv::Point3f(float32(x),float32(y),float32(z))
cvPoint3d(x, y, z) = @cxx cv::Point3d(float64(x),float64(y),float64(z))

# Size
cvSize(width, height) = @cxx cv::Size(width,height)
cvSize2f(width, height) = @cxx cv::Size2f(width,height)

# Scalar
cvScalar(blue::Int, green::Int, red::Int) = @cxx cv::Scalar(blue,green,red)

# Range
cvRange(start::Int, tend::Int) = @cxx cv::Range(start, tend)
cvRange_all(range) = @cxx range->all()

# Rectangle
cvRect(x::Int, y::Int, width::Int, height::Int) = @cxx cv::Rect(x, y, width, height)

# Rotated Rectangle
cvRotatedRect() =  @cxx cv::RotatedRect()
# const center = cvPoint2f(x,y), const size = cv::Size2f(width,height)
cvRotatedRect(center, size, angle::Float32) = @cxx cv::RotatedRect(center, size, angle)
# const point1 = cvPoint2f(x,y), const point2 = cvPoint2f(x,y), const point3 = cvPoint2f(x,y)
cvRotatedRect(point1, point2, point3) = @cxx cv::RotatedRect(point1, point2, point3)
# arrpts(x,y) = [cvPoint2f(x,y)]
cvRotatedRectPoints(arrpts) = @cxx rRect->points(arrpts)
# rRect = cvRotatedRect(pts, size, angle)
cvRotatedRectBoundingRect(rRect) = @cxx rRect->boundingRect()

#cv::String class
cvString() = @cxx cv::String()
cvString(str) = @cxx cv::String(str)
cvString(str, pos, len) = @cxx cv::String(str, csize_t(pos), csize_t(len))
cvString(s::String) = @cxx cv::String(pointer(s))
cvString(s::String, n) = @cxx cv::String(s::String, csize_t(n))
cvString(first::String, last::String) = @cxx cv::String(pointer(first), pointer(last))

#TermCriteria
# type     – TermCriteria::COUNT, TermCriteria::EPS or TermCriteria::COUNT + TermCriteria::EPS
# maxCount – The maximum number of iterations or elements to compute
# epsilon  – The desired accuracy or change in parameters

TermCriteria()      = @cxx cv::TermCriteria::TermCriteria()
TermCriteria(rtype::Int, maxCount::Int, epsilon::Float64) =
     @cxx cv::TermCriteria::TermCriteria(rtype, maxCount, epsilon)
TermCriteriaCount() = @cxx cv::TermCriteria::COUNT  # ::CppEnum  (not int)!
TermCriteriaEPS()   = @cxx cv::TermCriteria::EPS      # ::CppEnum  (not int)!

# to get the value of TermCriteriaCount and TermCriteriaEPS, e.g.,
# TermCriteriaCount().val or TermCriteriaEPS().val
# e.g., criteria=TermCriteria(int(TermCriteriaCount().val+TermCriteriaEPS().val), 30, 0.01)

###################################################
# Mat: The core image array structure in OpenCV
###################################################
# Mat array constructors
# Parameters:
# ndims    –>   dimension
# rws      –>   Number of rows in a 2D array
# cols     –>   Number of columns in a 2D array
# roi      –>   Region of interest
# size     –>   2D array size: Size(cols, rows)
# sizes    –>   Integer array
# type     –>   # select using CV_MatType in OpenCV_hpp.jl
                # Array type CV_8UC(n), ..., CV_64FC(n)
                # CV_8U       8-bit unsigned (0:255)
                # CV_8S       8-bit signed (-128:127)
                # CV_16U     16-bit unsigned
                # CV_16S     16-bit signed
                # CV_32S     32-bit signed
                # CV_32F     32-bit float
                # CV_64F     64-bit float
# s        –>   l value to initialize each matrix element with.
#              To set all the matrix elements to the particular value after the construction,
#              use the assignment operator Mat::operator=(const Scalar& value).
# data     –>   Pointer to the user data.
# step     –>   Number of bytes each matrix row occupies, if missing set to AUTO_STEP.
#              cols*elemSize(). See Mat::elemSize().
#              The value should include the padding bytes at the end of each row.
# steps    –>   Array of ndims-1 steps in case of a multi-dimensional array
#              (the last step is always set to the element size).
#              If not specified, the matrix is assumed to be continuous
# m        –>  Array that (as a whole or partly) is assigned to the constructed matrix.
#              No data is copied by these constructors.
#              Instead, the header pointing to m data or its sub-array is constructed and associated with it.
#              The reference counter, if any, is incremented.
# img      –>  Pointer to the old-style IplImage image structure.
# vec      –>  STL vector whose elements form the matrix.
#              1-column matrix =>  n rows = n vector elem
#              matrix Type => matches type of vector elem
#              arbitrary types OK, as long as declared DataType.
#              elems => primitive numbers or uni-type numerical tuples of numbers.
#              ***Mixed-type structures are not supported.***
#              "Mat(vec)" declared explicitly.
#              if copyData=true (see next), no new elements will be added to the vector
# copyData  –> Flag for STL vector, true = copied, false =  shared the newly constructed matrix.
#              (old-style CvMat or IplImage)
# rowRange  –> Range of the m rows to take.    Range::all()=> take all the rows.
# colRange  –> Range of the m columns to take. Range::all()=> take all the columns.
# ranges    –> Array of selected ranges of m along each dimensionality.

# Mat::Mat()
Mat() = @cxx cv::Mat()

# Mat::Mat(int rows, int cols, int type)
Mat(rows::Int, cols::Int, matType::CV_MatType) = @cxx cv::Mat(rows, cols, matType)

# Mat::Mat(Size size, int type)
#cxx""" cv::Mat Mat(int rows, int cols, int matType){ cv::Mat img(cv::Size(rows,cols), matType); return(img); }"""
Mat(size, matType::CV_MatType) = @cxx cv::Mat(size, matType)

# Mat::Mat(int rows, int cols, int type, const Scalar& s)
#cxx""" cv::Mat Mat(int rows, int cols, int matType){ cv::Mat img(cv::Size(rows,cols), matType); return(img); }"""
Mat(rows::Int, cols::Int, matType::CV_MatType, s) = @cxx cv::Mat(rows, cols, matType, s)

# Mat::Mat(Size size, int type, const Scalar& s)
Mat(size, matType::CV_MatType, s) = @cxx cv::Mat(size, matType, s)

# Mat::Mat(const Mat& m)
Mat(m) = @cxx cv::Mat(m)

# Mat::Mat(int ndims, const int* sizes, int type)
# const psizes(sizes) = pointer([sizes::Int])
Mat(ndims::Int, psizes::Ptr{Int}, matType::CV_MatType) = @cxx cv::Mat(ndims, psizes, matType)

# Mat::Mat(int ndims, const int* sizes, int type, const Scalar& s)
Mat(ndims::Int, psizes::Ptr{Int}, matType::CV_MatType, s) = @cxx cv::Mat(ndims, psizes, matType, s)

# Mat::Mat(const Mat& m, const Rect& roi)
Mat(img, roi) = @cxx cv::Mat(img, roi)

# Mat::Mat(const Mat& m, const Range* ranges)
# const ranges = pointer(range)
Mat(img, ranges) = @cxx cv::Mat(img, ranges)

# Mat class methods
# addition
cxx""" cv::Mat add(cv::Mat img1, cv::Mat img2) { return(img1 + img2); } """
imadd(img1, img2) = @cxx add(img1, img2)

# substract
cxx""" cv::Mat substract(cv::Mat img1, cv::Mat img2) { return(img1 - img2); } """
imsubstract(img1, img2) = @cxx substract(img1, img2)

# Mat::row
cxx""" cv::Mat row(cv::Mat img, int y) { return(img.row(y)); } """
row(img, x::Int) = @cxx row(img, x)

# Mat::col
cxx""" cv::Mat col(cv::Mat img, int x) { return(img.col(x)); } """
col(img, y::Int) = @cxx col(img, y)

# Mat::rowRange
cxx""" cv::Mat rowRange(cv::Mat img, int startrow, int endrow) { return(img.rowRange(startrow, endrow)); } """
rowRange(img, startrow::Int, endrow::Int) = @cxx rowRange(img, startrow, endrow)

cxx""" cv::Mat rowRange(cv::Mat img, const cv::Range& r) { return(img.rowRange(r)); } """
rowRange(img, range) = @cxx rowRange(img, range)

# Mat::colRange
# const range = cvRange(start::Int, tend::Int)
cxx""" cv::Mat colRange(cv::Mat img, int startcol, int endcol) { return(img.colRange(startcol, endcol)); } """
colRange(img, startcol::Int, endcol::Int) = @cxx colRange(img, startcol, endcol)

cxx""" cv::Mat colRange(cv::Mat img, const cv::Range& r) { return(img.colRange(r)); } """
colRange(img, range) = @cxx colRange(img, range)

# Mat::diag
# d=0 is the main diagonal
# d>0 is a diagonal from the lower half
# d=1 below the main one
# d<0 is a diagonal from the upper half

cxx""" cv::Mat diag(cv::Mat img, int d=0) { return(img.diag(d)); } """
diag(img, d::Int) = @cxx diag(img, d)

cxx""" cv::Mat diag(cv::Mat img, const cv::Mat& m) { return(img.diag(m)); } """
diag(img, m) = @cxx diag(img, m)

# Mat::clone()
cxx""" cv::Mat clone(cv::Mat img) { return(img.clone()); } """
clone(img) = @cxx clone(img)

# Mat::copyTo
cxx""" void copy(cv::Mat out, cv::Mat img) { img.copyTo(out); } """
cxx""" void copyTomask(cv::Mat img, cv::Mat mask, cv::Mat out) { img.copyTo(out, mask); } """
copy(out, img) = @cxx copy(out, img)
copyTomask(img, mask, out) = @cxx copyTomask(img, mask, out)

cxx""" cv::Mat imageROI(cv::Mat img, cv::Rect roi) { return img(roi); } """
ROImage(img, ROI) = @cxx imageROI(img, ROI)

# Mat::convertTo
# rtype - depth, rtype < 0 output = input type
# alpha – optional scale factor.
# beta – optional delta added to the scaled values.
convert(img1, img2, rtype::Int, alpha=1, beta=0) = @cxx img1->convertTo(img2, rtype, alpha, beta)

# Mat::assignTo
cxx""" void  assignTo(cv::Mat img, cv::Mat& m, int type) { img.assignTo(m, type); } """
assignTo(img, m, rtype=-1) = @cxx assignTo(img, m, rtype)

# Mat::setTo
# value(cvScalar), mask (same size as img)
set(img, value) = @cxx img->setTo(value, mask)

# Mat::reshape
# rows = 0 (no change)
cxx""" cv::Mat reshape(cv::Mat img, int cn, int rows) { return(img.reshape(cn, rows)); } """
reshape(img, ch::Int, rows=0) = @cxx reshape(img, ch, rows)

# Mat::t() (transpose)
cxx""" cv::Mat transpose(cv::Mat img, double lambda) { return(img.t()*lambda); } """
transpose(img, lambda) = @cxx transpose(img, lambda)

# Mat::inv (invert)
cxx""" cv::Mat inv(cv::Mat img, int method) { return(img.inv(method)); } """
inv(img, method=DECOMP_LU) = @cxx inv(img, method)

# Mat::mul (mutiply)   # element-wise only
cxx""" cv::Mat mul(cv::Mat img, double scale) { return(img.mul(scale)); } """
mul(img, scale=1) = @cxx mul(img, scale)

# Mat::cross (cross-product of 2 Vec<float,3>)
cxx""" cv::Mat cross(cv::Mat img, cv::Mat m) { return(img.cross(m)); } """
cross(img, m) = @cxx cross(img, m)

# Mat::dot (Computes a dot-product of two equally sized matrices)
cxx"""  double dot(cv::Mat img, cv::Mat m) { return(img.dot(m)); } """
dot(img, m) = @cxx dot(img, m)

# Empty array (NULL pointer)
noArray() = @cxx cv::noArray()

# Mat::zeros()
# ndims –> Array dimensionality.
# rows  –> Number of rows.
# cols  –> Number of columns.
# size  –> Alternative to the matrix size specification Size(cols, rows) .
# sz    –> Array of integers specifying the array shape.
# type  –> Created matrix type.

cxx""" cv::Mat zerosM(int rows, int cols, int matType) {cv::Mat A = cv::Mat::zeros(rows, cols, matType); return (A); }"""
cxx""" cv::Mat zerosM(cv::Size size, int matType) {cv::Mat A = cv::Mat::zeros(size, matType); return (A); }"""
cxx""" cv::Mat zerosM(int ndims, int* sz, int matType) {cv::Mat A = cv::Mat::zeros(ndims, sz, matType); return (A); }"""
zerosM(rows::Int, cols::Int, matType::CV_MatType) = @cxx zerosM(rows, cols, matType)
zerosM(size, matType::CV_MatType) = @cxx zerosM(size, matType)
zerosM(ndims::Int, sz::Ptr{Int32}, matType::CV_MatType) = @cxx zerosM(ndims, sz, matType)

# Mat::ones()
ones(rows::Int, cols::Int, matType::CV_MatType) = @cxx cv::Mat::ones(rows, cols, matType)
ones(size, matType::CV_MatType) = @cxx cv::Mat::ones(size, matType)
ones(ndims::Int, sz::Ptr{Int32}, matType::CV_MatType) = @cxx cv::Mat::ones(ndims, sz, matType)

# Mat::eye()
eye(rows::Int, cols::Int, matType::CV_MatType) = @cxx cv::Mat::eye(rows, cols, matType)
eye(size, matType::CV_MatType) = @cxx cv::Mat::eye(size, matType)

# Mat::addref
addref(img) = @cxx img->addref()

# Mat::release
destroy(img) = @cxx img->release()

# Mat::resize
resizeMat(img, sz) = @cxx img->resize(csize_t(sz))         # sz – new n rows (UInt64), s = cvScalar()
resizeMat(img, sz, s) = @cxx img->resize(csize_t(sz), s)

# Mat::reserve
# Reserves space for the certain number of rows
reserve(img,sz) = @cxx img->reserve(csize_t(sz))

# Mat::push_back
# Adds elements to the bottom of the matrix:
# const Mat& m, const T& elem
push(img, m) = @cxx img->push_back(m)
push(img, elem) = @cxx img->push_back(elem)

# Mat::pop_back
pop(img, nelems::Csize_t) = @cxx img->pop_back(nelems)

# Mat::locateROI
# void Mat::locateROI(Size& wholeSize, Point& ofs)
locateROI(img, wholeSize, ofs) = @cxx img->locateROI(wholeSize, ofs)

# Mat::adjustROI
# Mat& Mat::adjustROI(int dtop, int dbottom, int dleft, int dright)
# dtop – Shift of the top submatrix boundary upwards.
# dbottom – Shift of the bottom submatrix boundary downwards.
# dleft – Shift of the left submatrix boundary to the left.
# dright – Shift of the right submatrix boundary to the right.
adjustROI(img, dtop::Int, dbottom::Int, dleft::Int, dright::Int) = @cxx img->adjustROI(dtop, dbottom, dleft, dright)

total(img) = @cxx img->total()                               # returns size_t: number of array elements
dims(img) = @cxx img->dims                                   # ndims
rows(img) = @cxx img->rows                                   # rows
cols(img) = @cxx img->cols                                   # cols
isContinuous(img) = @cxx img->isContinuous()                 # bool
elemSize(img) = @cxx img->elemSize()                         # step (size_t), e.g., CV_16SC3 = 6 (3*sizeof(short))
elemSize1(img) = @cxx img->elemSize1()                       # size of element e.g., CV_16SC3 = 2 (sizeof(short))
depth(img) = @cxx img->depth()
step(img) = elemSize(img)*cols(img)                          # assumes no padding
step1(img, i) = @cxx img->step1(i)                           # normalized step, default i = 1
flags(img) = @cxx img->flags                                 # array dimensionality >= 2
data(img) = @cxx img->data                                   # pointer to the user data
refcount(img) = @cxx img->refcount                           # reference count, pointer is NULL (if user-allocated data)
empty(img) = @cxx img->empty()                               # bool
ptr(img, row) = @cxx img->ptr(row)                           # return uchar* or typed pointer for matrix row

cxx""" cv::Size sizeofimage(cv::Mat img) {  return(img.size()); } """
sizeofimage(img) = @cxx sizeofimage(img)                                   # returns Size(cols, rows), if matrix > 2d, size = (-1,-1)

cxx""" int channels(cv::Mat img) { return(img.channels()); } """
channels(img) = @cxx channels(img)                           # number of matrix channels

findlabel(value::Int) = CV_MAT_TYPE[value]                   # e.g., CV_8UC1 type CV_MatType in OpenCV_hpp.jl
cxx""" int cvtype(cv::Mat img) { return(img.type()); } """
cvtypeval(img) = @cxx cvtype(img)
cvtypelabel(img) = findlabel(int(cvtypeval(img)))

# SparseMat: multi-dimensional sparse numerical arrays
# Parameters:
# m – Source matrix for copy constructor
#    If m is dense matrix (ocv:class:Mat) -> sparse representation
# dims – Array dimensionality
# _sizes – Sparce matrix size on all dementions, i.e.,  const int*
# _type – Sparse matrix data type, i.e., const SparseMat&

SparseMat() = @cxx cv::SparseMat()
SparseMat(dims::Int, _sizes::Ptr{Int32}, _type::Int) = @cxx cv::SparseMat::SparseMat(dims, _sizes, _type)
SparseMat(sparseM) = @cxx cv::SparseMat::SparseMat(sparseM)
SparseMat(M) = @cxx cv::SparseMat::SparseMat(M)

# clone
cxx""" cv::SparseMat cloneSparse(cv::SparseMat img) { return(img.clone()); } """
cloneSparse(SparseM) = @cxx cloneSparse(SparseM)

# SparseMat::copyTo
cxx""" void copyToSparse(cv::SparseMat img, cv::SparseMat copy) { img.copyTo(copy); } """
cxx""" void copyToSparse(cv::SparseMat img, cv::Mat copy) { img.copyTo(copy); } """
copySparse(SparseM, copy) = @cxx copySparse(SparseM, copy)
copySparse(M, copy) = @cxx copySparse(M, copy)

# SparseMat::convertTo
convertSparse(sparseM1, sparseM2, rtype::Int, alpha=1, beta=0) = @cxx img1->convertToSparse(sparseM1, sparseM2, rtype, alpha, beta)
SparseToMat(sparseM, M, rtype::Int, alpha=1, beta=0) = @cxx img1->SparseToMat(sparseM, M, rtype, alpha, beta)

# SparseMat::create
createSparse(sparseM) = @cxx sparseM->create(dims, _sizes, _type)

# SparseMat::clear
clearSparse(sparseM) = @cxx sparseM->clear()
addrefSparse(sparseM) = @cxx sparseM->addref()
destroySparse(sparseM) = @cxx sparseM->release()              # ~SparseMat()
elemSizeSparse(sparseM) = @cxx sparseM->elemSize()
elemSize1Sparse(sparseM) = @cxx sparseM->elemSize1()

cxx""" int typeSparse(cv::SparseMat img) { return(img.type()); } """
typeSparse(sparseM) = @cxx typeSparse(sparseM)
typeSparselabel(sparseM) = findlabel(int(typeSparse(sparseM)))
depthSparse(sparseM) = @cxx sparseM->depth()

cxx""" int channelsSparse(cv::SparseMat img) { return(img.channels()); } """
channelsSparse(sparseM) = @cxx channelsSparse(sparseM)         # number of matrix channels

#-------------------------------------------------------------------------------------------------------------------#
# Command Line Parser (incomplete)
# The CommandLineParser class is designed for command line arguments parsing

# CommandLineParser::CommandLineParser(int argc, const char* const argv[], const String& keys)
# Parameters:
# argc –
# argv –
# keys –
CommandLineParser(argc::Int, argv::Ptr{Ptr{UInt8}}, keys::String) =
       @cxxnew cv::CommandLineParser(argc, argv, pointer(keys))

# bool CommandLineParser::has(const String& name)
has(parser, name::String) = @cxx parser->has(stdstring(name))

#-------------------------------------------------------------------------------------------------------------------#
# 2. Operations on Arrays

abs(img) = @cxx cv::abs(img)
absdiff(src1, src2, dst) = @cxx absdiff(src1, src2, dst)

add(src1, src2, dst, mask = noArray(), dtype=-1) = @cxx cv::add(src1, src2, dst, mask, dtype)
# mask – optional operation mask - 8-bit single channel array, default mask = noArray()
# dtype – optional depth of the output array

addWeighted(src1, alpha, src2, beta::Float64, gamma::Float64, dst, dtype=-1) =
    @cxx cv::addWeighted(src1, alpha, src2, beta, gamma, dst, dtype)
# beta – weight of the second array elements
# gamma – scalar added to each sum

bitwise_and(src1, src2, dst,  mask = noArray()) = @cxx cv::bitwise_and(src1, src2, dst, mask)
bitwise_not(src, dst, mask = noArray()) = @cxx cv::bitwise_not(src, dst, mask)
bitwise_or(src1, src2, dst, mask = noArray()) = @cxx cv::bitwise_or(src1, src2, dst,mask)
bitwise_xor(src1, src2, dst, mask = noArray()) = @cxx cv::bitwise_xor(src1, src2, dst,mask)
# mask – optional operation mask - 8-bit single channel array, default mask = noArray()

calcCovarMatrix(sample, nsamples::Int, covar, mean, flags::Int, ctype=CV_64F) =
    @cxx cv::calcCovarMatrix(samples, nsamples, covar, mean, flags, ctype)
# samples – samples stored either as separate matrices or as rows/columns of a single matrix - pointer(mat)
# covar   – output covariance matrix of the type ctype and square size
# mean    – input or output (depending on the flags) array
# vects   – a set of vectors
# flags:
#    CV_COVAR_SCRAMBLED
#    CV_COVAR_USE_AVG
#    CV_COVAR_SCALE
#    CV_COVAR_ROWS
#    CV_COVAR_COLS

cartToPolar(x, y, magnitude, angle, angleInDegrees=false) =
    @cxx cv::cartToPolar(x, y, magnitude, angle, angleInDegrees=false)
polarToCart(magnitude, angle, x, y, angleInDegrees=false) =
    @cxx cv::polarToCart(magnitude, angle, x, y, angleInDegrees=false)
# x              – array of x-coordinates; single-precision or double-precision floating-point array
# y              – array of y-coordinates, same size and same type as x
# magnitude      – output array of magnitudes of the same size and type as x
# angle          – output array of angles that has the same size and type as x;
#                   in radians (from 0 to 2*Pi) or in degrees (0 to 360 degrees)
# angleInDegrees – bool flag = false (radians), flag = true (degrees)

checkRange(a, point, quiet=true, minVal=-DBL_MAX, maxVal=DBL_MAX) =
    @cxx cv::checkRange(a, quiet, point, minVal, maxVal)
# a       – input array
# quiet   – bool flag to throw error if array elements are out of range (true/false)
# pos     – optional output parameter, pos of first outlier is stored
#           when not NULL,it must be a pointer to array of src.dims elements
# minVal  – inclusive lower boundary of valid values range
# maxVal  – exclusive upper boundary of valid values range

compare(src1, src2, dst, cmpop::Int) = @cxx cv::compare(rc2, dst, cmpop)
# cmpop     –  flag, that specifies correspondence between the arrays:
# CMP_EQ if src1 = src2
# CMP_GT if src1 > src2
# CMP_GE src1 >= src2
# CMP_LT src1 < src2
# CMP_LE src1 =< src2
# CMP_NE src1 != src2
# e.g., Mat dst2 = src1 < 8;

completeSymm(mtx, lowerToUpper=false) = @cxx cv::completeSymm(mtx, lowerToUpper)
# mtx          – input-output floating-point square matrix
# lowerToUpper – operation flag; if true, the lower half is copied to the upper half
#                Otherwise, the upper half is copied to the lower half

convertScaleAbs(src, dst, alpha=1.0, beta=float(0)) = @cxx cv::convertScaleAbs(src, dst, alpha, beta)
# src   – input array
# dst   – output array
# alpha – optional scale factor
# beta  – optional delta added to the scaled values

countNonZero(src) = @cxx cv::countNonZero(src)

dct(src, dst, flags=0) = @cxx cv::dct(src, dst, flags)
# forward or inverse discrete Cosine transform of 1D or 2D array
# flags  - transformation flags
# DCT_INVERSE
# DCT_ROWS

dft(src, dst, flags=0, nonzeroRows=0) = @cxx cv::dft(src, dst, flags, nonzeroRows)
# forward or inverse Discrete Fourier transform of a 1D or 2D floating-point array

divide(src1, src2, dst, scale=1.0, dtype=-1) = @cxx cv::divide(src1, src2, dst, scale, dtype)
# per-element division of two arrays or a scalar by an array

determinant(mtx) = @cxx cv::determinant(mtx)
# determinant of a square floating-point matrix

eigen(src, eigenvalues, eigenvectors=noArray()) = @cxx cv::eigen(src, eigenvalues,eigenvectors)
# src           – input matrix that must have CV_32FC1 or CV_64FC1
# eigenvectors  – default =noArray()
# eigenvalues   – output vector [eigenvalues, type = src], stored in descending order
# lowindex      – optional index of largest eigenvalue/-vector to calculate
# highindex     – optional index of smallest eigenvalue/-vector to calculate

exp(src, dst) = @cxx cv::exp(src, dst)

flip(src, dst, flipCode) = @cxx cv::flip(src, dst, flipCode)
# flipCode – a flag to specify how to flip the array
#   0 -> x-axis
#  >0 -> y-axis
#  <0 -> x,y axes

gemm(src1, src2, alpha, src3, beta, dst, flags=0) = @cxx cv::gemm(src1, src2, alpha, src3, beta, dst, flags)
# src1      – first multiplied input matrix that could be real(CV_32FC1, CV_64FC1) or complex(CV_32FC2, CV_64FC2)
# src2      – second multiplied input matrix of the same type as src1
# alpha     – weight of the matrix product
# src3      – third optional delta matrix added to the matrix product;
# beta      – weight of src3
# dst       – output matrix
# flags     –
# GEMM_1_T transposes src1
# GEMM_2_T transposes src2
# GEMM_3_T transposes src3

getOptimalDFTSize(vecsize::Int) = @cxx cv::getOptimalDFTSize(vecsize)
# optimal DFT size for a given vector size

idct(src, dst, flags=0) = @cxx cv::idct(src, dst, flags)
# inverse Discrete Cosine Transform of a 1D or 2D array
# flags – operation flags

idft(src, dst, flags=0, nonzeroRows=0) = @cxx cv::idft(src, dst, flags, nonzeroRows)
# flags – operation flags
# nonzeroRows – number of dst rows to process

inRange(src, lowerb, upperb, dst) = @cxx cv::inRange(src, lowerb, upperb, dst)
# src    – first input array
# lowerb – inclusive lower boundary array or a scalar
# upperb – inclusive upper boundary array or a scalar
# dst    – output array of the same size as src and CV_8U type

# same as inv for cv::Mat
invert(src, dst, flags=DECOMP_LU) = @cxx cv::invert(src, dst, flags)

# natural logarithm of every array element
log(src, dst) = @cxx cv::log(src, dst)

# look-up table transform of an array
LUT(src, lut, dst) = @cxx cv::LUT(src, lut, dst)
# src – input array of 8-bit elements
# dst – output array of the same size
# lut – look-up table of 256 elements
# multi-channel input array, single channel (same lut table for all channel)
# or the same number of channels as in the input array.

# magnitude of 2D vectors
magnitude(x, y, magnitude) = @cxx cv::magnitude(x, y, magnitude)
# x         – floating-point array of x-coordinates of the vectors
# y         – floating-point array of y-coordinates of the vectors; it must have the same size as x
# magnitude – output array of the same size and type as x

# Mahalanobis distance between two vectors
Mahalanobis(v1, v2, icovar) = @cxx cv::Mahalanobis(v1, v2, icovar)
# vec1   – first 1D input vector
# vec2   – second 1D input vector
# icovar – inverse covariance matrix

max(a, b) = @cxx cv::max(a,b)
max(a, scalar) = @cxx cv::max(a,scalar)

mean(src,mask=noArray()) = @cxx cv::mean(src, mask)
# mask  – optional operation mask, default mask=noArray()

meanStdDev(src, mean, stddev, mask=noArray()) = @cxx cv::meanStdDev(src, mean, stddev, mask)
# mask  – optional operation mask, default mask=noArray()

merge(mv, count::UInt64, dst) = @cxx cv::merge(mv, count, dst)
# mv  => const pointer(cv::Mat) or const Mat* mv
# dst – output array of the same size

min(a, b) = @cxx cv::min(a,b)
min(a, scalar) = @cxx cv::min(a,scalar)

# global minimum and maximum in an array
minMaxIdx(src, minVal::Ptr{Float64}, maxVal=pointer([float(0)]), iminIdx=pointer([0]), maxIdx=pointer([0]), mask=noArray()) =
    @cxx cv::minMaxIdx(src, minVal, maxVal, iminIdx, maxIdx, mask)
# maxVal::Ptr{Float64}
# minIdx::Ptr{Int}
# maxIdx::Ptr{Int}
# src – input single-channel array.
# minVal – pointer to the returned minimum value; NULL is used if not required.
# maxVal – pointer to the returned maximum value; NULL is used if not required.
# minIdx – minIdx is not NULL, it must have at least 2 elements (as well as maxIdx),

# global minimum and maximum in an array
minMaxLoc(src, minVal::Ptr{Float64}, maxVal=pointer([float(0)]), minLoc=pointer([0]), maxLoc=pointer([0]), mask=noArray()) =
    @cxx cv::minMaxLoc(src, minVal, maxVal, minLoc, maxLoc, mask)
# minLoc = convert(C_NULL, pointer(cvPoint))
# maxLoc = convert(C_NULL, pointer(cvPoint))
# mask  – optional operation mask, default mask=noArray()
# minVal – pointer to the returned minimum value; NULL is used if not required.
# maxVal – pointer to the returned maximum value; NULL is used if not required.
# minLoc – pointer to the returned minimum location (in 2D case); NULL is used if not required.
# maxLoc – pointer to the returned maximum location (in 2D case); NULL is used if not required.
# mask – optional mask used to select a sub-array.

mixChannels(src, dst, fromTo::Ptr{Int}, npairs::UInt64) = @cxx cv::mixChannels(src, dst, fromTo, npairs)
# const fromTo::Ptr{Int}
# fromTo = const std::vector<int>&

# per-element multiplication of two Fourier spectrums
mulSpectrums(a, b, c, flags=DFT_ROWS, conjB=false) = @cxx cv::mulSpectrums(a, b, c, flags, conjB)
# flags – currently, the only supported flag is DFT_ROWS

# per-element scaled product of two arrays
multiply(src1, src2, dst, scale=1.0, dtype=-1) = @cxx cv::multiply(src1, src2, dst, scale, dtype)

# product of a matrix and its transposition
mulTransposed(src, dst, aTa::Bool, delta=noArray(), scale=1.0, dtype=-1) = @cxx cv::mulTransposed(src, dst, aTa,delta, scale, dtype)
# InputArray delta=noArray()

# absolute array norm, an absolute difference norm, or a relative difference norm
norm(src1, normType=NORM_L2,mask=noArray()) = @cxx cv::norm(src1, normType=NORM_L2, mask)
norm(src1, src2, normType=NORM_L2,mask=noArray()) = @cxx cv::norm(src1, src2, normType=NORM_L2, mask)
# InputArray mask=noArray()

# norm or value range of an array
normalize(src, dst, alpha=1.0, beta=0, norm_type=NORM_L2, dtype=-1, mask=noArray()) = @cxx cv::normalize(src, dst, alpha, beta, norm_type, dtype, mask)
# InputArray mask=noArray()

# PCA: Principal Component Analysis class
# Constructors
pca() = cv::PCA::PCA()
pca(data, mean=noArray(), flags=CV_PCA_DATA_AS_ROW, maxComponents=0) =  @cxx cv::PCA::PCA(data, mean, flags, maxComponents)
pca(data, mean=noArray(), flags=CV_PCA_DATA_AS_ROW, retainedVariance=2.0) = @cxx cv::PCA::PCA(data, mean, flags, retainedVariance)
# Run PCA
pca_operator() = @cxx cv::PCA::PCA->operator()
# data             – input samples stored as matrix rows or matrix columns.
# mean             – optional mean value; if the matrix is empty (noArray()), the mean is computed from the data.
# flags            – operation flags; currently the parameter is only used to specify the data layout:
#                     CV_PCA_DATA_AS_ROW stored as matrix rows
#                     CV_PCA_DATA_AS_COL stored as matrix columns
# maxComponents    – maximum number of components that PCA should retain; by default, all the components are retained.
# retainedVariance – Percentage of variance that PCA should retain.
#                    Using this parameter will let the PCA decided how many components to retain
#                    but it will always keep at least 2.

# Projects vector(s) to the principal component subspace
pca_project(vec) =  @cxx cv::PCA::project(vec)                  # returns cv::Mat
pca_project(vec, result) =  @cxx cv::PCA::project(vec, result)  # void
# vec    – input vector(s), if CV_PCA_DATA_AS_ROW then vec.cols==data.cols
# result – output vectors, if CV_PCA_DATA_AS_COL then result.cols==vec.cols

# Reconstructs vectors from their PC projections
pca_backProject(vec) = @cxx cv::PCA::backProject(vec)
pca_backProject(vec, result) = @cxx cv::PCA::backProject(vec, result)

# Performs the perspective matrix transformation of vectors
perspectiveTransform(src, dst, m) = @cxx cv::perspectiveTransform(src, dst, m)
# src   – input two-channel or three-channel floating-point array
# dst   – output array of the same size and type as src
# m     – 3x3 or 4x4 floating-point transformation matrix

# rotation angle of 2D vectors
phase(x, y, angle, angleInDegrees=false) = @cxx cv::phase(x, y, angle, angleInDegrees)

# Raises every array element to a power
pow(src, power::Float64, dst) = @cxx cv::pow(src, power, dst)

# Random number generator
rng() = @cxx cv::RNG::RNG()
rng(state::UInt64) = @cxx cv::RNG::RNG(state)
rng_next() = @cxx cv::RNG::RNG->next()
rng_uniform(a::Int, b::Int) = @cxx cv::RNG::RNG->uniform(a, b)
rng_uniform(a::Float32, b::Float32) = @cxx cv::RNG::RNG->uniform(a, b)
rng_uniform(a::Float64, b::Float64) = @cxx cv::RNG::RNG->uniform(a, b)
rng_gaussian(sigma::Float64) = @cxx cv::RNG::RNG->gaussian(sigma)
rng_fill(mat, distType, a, b, saturateRange=false) = @cxx cv::RNG::RNG->fill(mat, distType, a, b, saturateRange)
# distType = RNG::UNIFORM or RNG::NORMAL

# array with normally distributed random numbers
randn(dst, mean, stddev) = @cxx cv::randn(dst, mean, stddev)

# Shuffles the array elements randomly
theRNG() = @cxx cv::theRNG()
randShuffle(dst, iterFactor=1.0, rng = theRNG()) = cv::randShuffle(dst, iterFactor, rng)
# RNG* rng=0

# Reduces a matrix to a vector
reduce(src, dst, dim::Int, rtype::Int, dtype=-1) = @cxx cv::reduce(src, dst, dim, rtype, dtype=-1)

# output array with repeated copies of the input array
repeat(src, ny::Int, nx::Int, dst) = @cxx cv::repeat(src, ny, nx, dst)
repeat(src, ny::Int, nx::Int) = @cxx cv::repeat(src, ny, nx)

# sum of a scaled array and another array
scaleAdd(src1, alpha::Float64, src2, dst) = @cxx cv::scaleAdd(src1, alpha, src2, dst)

# Initializes a scaled identity matrix
setIdentity(mtx, s=cvScalar(1, 0, 0)) =  @cxx cv::setIdentity(mtx, s)

# Solves one or more linear systems or least-squares problems
solve(src1, src2, dst, flags=DECOMP_LU) = @cxx cv::solve(src1, src2, dst, flags)
# DECOMP_LU
# DECOMP_CHOLESKY
# DECOMP_EIG
# DECOMP_SVD
# DECOMP_QR
# DECOMP_NORMAL

# real roots of a cubic equation
solveCubic(coeffs, roots) = @cxx cv::solveCubic(coeffs, roots)
#coeffs – equation coefficients, an array of 3 or 4 elements
#roots – output array of real roots that has 1 or 3 elements

# real or complex roots of a polynomial equation
solvePoly(coeffs, roots, maxIters=300) = @cxx cv::solvePoly(coeffs, roots, maxIters)

# Sorts each row or each column of a matrix
# CV_SORT_EVERY_ROW
# CV_SORT_EVERY_COLUMN
# CV_SORT_ASCENDING
# CV_SORT_DESCENDING
sort(src, dst, flags::Int) = @cxx cv::sort(src, dst, flags)
sortIdx(src, dst, flags::Int) = @cxx cv::sortIdx(src, dst, flags)

# Divides a multi-channel array into several single-channel arrays
split(src, mvbegin) = @cxx cv::split(src, mvbegin)    #mvbegin = Mat*
split(m, mv) = @cxx cv::split(src, mvbegin)

# square root of array elements
sqrt(src, dst) = @cxx cv::sqrt(src, dst)
# src - input floating-point array
# dst – output array

subtract(src1, src2, dst, mask=noArray(), dtype=-1) = @cxx cv::subtract(src1, src2, dst, mask, dtype)
# InputArray mask=noArray()

# Singular Value Decomposition of a floating-point matrix
svd() = cv::SVD::SVD()
svd(src, flags=0) = @cxx cv::SVD::SVD(src, flags)
# src – decomposed matrix
# flags -
# SVD::MODIFY_A  # modify the decomposed matrix
# SVD::NO_UV     # use only singular values
# SVD::FULL_UV   # u and vt are full-size square orthogonal matrices

# Perform SVD of a matrix
svd_compute(src, w, u, vt, flags=0) = @cxx cv::SVD::compute(src, w, u, vt, flags)
svd_compute(src, w, flags=0) = @cxx cv::SVD::compute(src, w, flags)
# src – decomposed matrix
# w   – calculated singular values
# u   – calculated left singular vectors
# V   – calculated right singular vectors
# vt  – transposed matrix of right singular values

# Solves an under-determined singular linear system
svd_solveZ(src, dst) = @cxx cv::SVD::solveZ(src, dst)

# sum of array elements
sum(src) = @cxx cv::sum(src)

# trace of a matrix
trace(mtx) = @cxx cv::trace(mtx)

# matrix transformation of every array element
transform(src, dst, m) = @cxx cv::transform(src, dst, m)

# Transpose a matrix
transposeMat(src, dst) = @cxx cv::transpose(src, dst)

# Computes the source location of an extrapolated pixel
borderInterpolate(p::Int, len::Int, borderType::Int) = @cxx cv::borderInterpolate(p, len, borderType)
# Forms a border around an image
copyMakeBorder(src, dst, top::Int, bottom::Int, left::Int, right::Int, borderType::Int, value=cvScalar()) =
    @cxx cv::copyMakeBorder(src, dst, top, bottom, left, right, borderType, value)
# borderType ->
# BORDER_CONSTANT
# BORDER_TRANSPARENT
# BORDER_ISOLATED

#-------------------------------------------------------------------------------------------------------------------#
# 3. Clustering

# Finds centers of clusters and groups input samples around the clusters
kmeans(data, K::Int, bestLabels, criteria, attempts::Int, flags::Int,centers=noArray()) =
    @cxx cv::kmeans(data, K, bestLabels, criteria, attempts, flags, centers)
# OutputArray centers=noArray()
# TermCriteria
# partition

#-------------------------------------------------------------------------------------------------------------------#
# 4. Utility and System Functions and Macros

# angle of a 2D vector in degrees
fastAtan2(y::Float64, x::Float64) = @cxx cv::fastAtan2(y, x)

# cube root of an argument
cubeRoot(val::Float64) =  @cxx cv::cubeRoot(val)

# fastMalloc
fastMalloc(bufSize::UInt64) = @cxx cv::fastMalloc(bufSize)  # returns void*

# fastFree
fastFree(ptr::Ptr{Nothing}) = @cxx cv::fastFree(ptr)
# ptr – Pointer to the allocated buffer

# format strings
cvString(text::String) = @cxx cv::format(pointer(text)) # returns cv::String

# convert cv::String to Julia String
julString(cvstr) = bytestring(@cxx cvstr->c_str())

# getBuildInformation
getBuildInformation() = @cxx cv::getBuildInformation()

# heckHardwareSupport(int feature)
checkHardwareSupport(feature::Int) = @cxx cv::checkHardwareSupport(feature)
# CV_CPU_MMX
# CV_CPU_SSE
# CV_CPU_SSE2
# CV_CPU_SSE3
# CV_CPU_SSSE3
# CV_CPU_SSE4_1
# CV_CPU_SSE4_2
# CV_CPU_POPCNT
# CV_CPU_AVX

getNumberOfCPUs() = @cxx cv::getNumberOfCPUs()
getNumThreads() = @cxx cv::getNumThreads()
setNumThreads(nthreads::Int) = @cxx cv::setNumThreads(nthreads)

# Returns the number of ticks after a certain event
getTickCount() = @cxx cv::getTickCount()
getTickFrequency() = @cxx cv::getTickFrequency()
getCPUTickCount() = @cxx cv::getCPUTickCount()

# Enables or disables the optimized code
cvUseOptimized(on_off::Int) = @cxx cv::cvUseOptimized(on_off)
useOptimized() = @cxx cv::useOptimized()  #true/false status


text:
# ----------------------------------------------------------------------
# compute Wigner U-functions for one neighbor
# ------------------------------------------------------------------------- */
function compute_ui(
    ielem::Int,
    snap::SNAPParams,
    runtime_arrays :: RuntimeArrays)

    if ~snap.chem_flag
        ielem = 1
    end

    zero_uarraytot(ielem, snap, runtime_arrays)

    num_of_interactions = length(runtime_arrays.indij)
    for ind = 1:num_of_interactions
        jj = runtime_arrays.indij[ind][2]
        rij = runtime_arrays.rij[ind]
        rcutij = runtime_arrays.rcutij[ind]
        wj = runtime_arrays.wj[ind]
        elementj = runtime_arrays.element[ind]

        x = rij[1]
        y = rij[2]
        z = rij[3]
        rsq = x*x + y*y + z*z
        r   = sqrt(rsq)
        theta0 = (r - snap.rmin0) * snap.rfac0 * pi / (rcutij - snap.rmin0)
        z0 = r / tan(theta0)
        compute_uarray(x, y, z, z0, r, jj,
                snap,
                runtime_arrays)

        if (snap.chem_flag)
            add_uarraytot(r, rcutij, wj,
                    jj, elementj, snap, runtime_arrays)
        else
            add_uarraytot(r, rcutij, wj, jj, 1, snap, runtime_arrays)
        end
    end
end

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

function compute_uarray(x::AbstractFloat, y::AbstractFloat, z::AbstractFloat,
    z0::AbstractFloat, r::AbstractFloat, jj::Int,
    snap::SNAPParams, runtime_arrays::RuntimeArrays)

    # compute Cayley-Klein parameters for unit quaternion
    r0inv = 1.0 / sqrt(r * r + z0 * z0);
    a_r = r0inv * z0;
    a_i = -r0inv * z;
    b_r = r0inv * y;
    b_i = -r0inv * x;
    # // VMK Section 4.8.2
    ulist_r = runtime_arrays.ulist_r_ij[jj, :];
    ulist_i = runtime_arrays.ulist_i_ij[jj, :];
    ulist_r[1] = 1.0;
    ulist_i[1] = 0.0;

    for j = 1:(snap.twojmax)
        jju = snap.prebuilt_arrays.idxu_block[j+1];
        jjup = snap.prebuilt_arrays.idxu_block[j];

        # fill in left side of matrix layer from previous layer
        mb = 0
        while 2*mb <= j
            ulist_r[jju+1] = 0.0;
            ulist_i[jju+1] = 0.0;

            for ma = 0:(j-1)
                rootpq = snap.prebuilt_arrays.rootpqarray[j - ma+1, j - mb+1];
                ulist_r[jju+1] +=
                    rootpq *
                    (a_r * ulist_r[jjup+1] +
                    a_i * ulist_i[jjup+1]);
                ulist_i[jju+1] +=
                    rootpq *
                    (a_r * ulist_i[jjup+1] -
                    a_i * ulist_r[jjup+1]);
                rootpq = snap.prebuilt_arrays.rootpqarray[ma + 1+1, j - mb+1];
                ulist_r[jju+1+1] =
                        -rootpq *
                        (b_r * ulist_r[jjup+1] +
                        b_i * ulist_i[jjup+1]);
                ulist_i[jju+1+1] =
                        -rootpq *
                        (b_r * ulist_i[jjup+1] -
                        b_i * ulist_r[jjup+1]);
                jju+=1;
                jjup+=1;
            end
            jju+=1;
            mb +=1;
        end

        #  // copy left side to right side with inversion symmetry VMK 4.4(2)
        #  // u[ma-j][mb-j] = (-1)^(ma-mb)*Conj([u[ma][mb])
        jju = snap.prebuilt_arrays.idxu_block[j+1];
        jjup = jju+(j+1)*(j+1) - 1 ;
        mbpar = 1;
        mb = 0;
        while 2*mb <= j
            mapar = mbpar;
            for ma = 0:j
                if (mapar == 1)
                    ulist_r[jjup+1] = ulist_r[jju+1];
                    ulist_i[jjup+1] = -ulist_i[jju+1];
                else
                    ulist_r[jjup+1] = -ulist_r[jju+1];
                    ulist_i[jjup+1] = ulist_i[jju+1];
                end
                mapar = -mapar;
                jju+=1;
                jjup-=1;
            end
            mbpar = -mbpar;
            mb += 1;
        end
    end
    runtime_arrays.ulist_r_ij[jj, :] = ulist_r
    runtime_arrays.ulist_i_ij[jj, :] = ulist_i
end

function zero_uarraytot(ielem::Int, snap::SNAPParams, runtime_arrays::RuntimeArrays)
    num_elements = snap.prebuilt_arrays.num_elements
    for jelem = 0:(num_elements-1)
        for j = 0:snap.twojmax
            jju = snap.prebuilt_arrays.idxu_block[j+1];
            for mb = 0:j
                for ma = 0:j
                    runtime_arrays.ulisttot_r[jelem*snap.prebuilt_arrays.idxu_max+jju+1] = 0.0;
                    runtime_arrays.ulisttot_i[jelem*snap.prebuilt_arrays.idxu_max+jju+1] = 0.0;

                    # // utot(j,ma,ma) = wself, sometimes
                    if ((jelem+1) == ielem) || (snap.wselfall_flag)
                        if (ma==mb)
                            runtime_arrays.ulisttot_r[jelem*snap.prebuilt_arrays.idxu_max+jju+1] = snap.prebuilt_arrays.wself; # ///// double check this
                        end
                    end
                    jju+=1;
                end
            end
        end
    end
end

function add_uarraytot(r::AbstractFloat, rcut::AbstractFloat, wj::AbstractFloat,
        jj::Int, jelem :: Int,
        snap::SNAPParams,
        runtime_arrays::RuntimeArrays)
    sfac = compute_sfac(snap, r, rcut);
    sfac *= wj;
    ulist_r = runtime_arrays.ulist_r_ij[jj, :];
    ulist_i = runtime_arrays.ulist_i_ij[jj, :];
    for j = 0:snap.twojmax
        jju = snap.prebuilt_arrays.idxu_block[j+1];
        for mb = 0:j
            for ma = 0:j
                runtime_arrays.ulisttot_r[(jelem-1)*snap.prebuilt_arrays.idxu_max+jju+1] +=
                        sfac * ulist_r[jju+1];
                runtime_arrays.ulisttot_i[(jelem-1)*snap.prebuilt_arrays.idxu_max+jju+1] +=
                        sfac * ulist_i[jju+1];
                jju+=1
            end
        end
    end
end

# /* ----------------------------------------------------------------------
# Compute derivatives of Wigner U-functions for one neighbor
# see comments in compute_uarray()
# ------------------------------------------------------------------------- */

function compute_duarray(x::AbstractFloat, y::AbstractFloat, z::AbstractFloat,
     z0::AbstractFloat, r::AbstractFloat, dz0dr::AbstractFloat,
     wj::AbstractFloat, rcut::AbstractFloat, jj::Int,
     snap::SNAPParams, runtime_arrays::RuntimeArrays)

    rinv = 1.0/r;
    ux   = x * rinv;
    uy   = y * rinv;
    uz   = z * rinv;

    r0inv = 1.0 / sqrt(r * r + z0 * z0);
    a_r   = z0 * r0inv;
    a_i   = -z * r0inv;
    b_r   = y * r0inv;
    b_i   = -x * r0inv;

    dr0invdr = -r0inv*r0inv*r0inv * (r + z0 * dz0dr);

    dr0inv = [dr0invdr * ux, dr0invdr * uy, dr0invdr * uz];
    dz0    = [dz0dr * ux, dz0dr * uy, dz0dr * uz];

    da_r   = [dz0[k] * r0inv + z0 * dr0inv[k] for k = 1:3]
    da_i   = [-z * dr0inv[k] for k = 1:3]

    da_i[3] += -r0inv;

    db_r   = [ y * dr0inv[k] for k = 1:3]
    db_i   = [ -x * dr0inv[k] for k = 1:3]

    db_i[1] += -r0inv;
    db_r[2] += r0inv;

    ulist_r = runtime_arrays.ulist_r_ij[jj, :];
    ulist_i = runtime_arrays.ulist_i_ij[jj, :];

    runtime_arrays.dulist_r[1, :] = [0.0, 0.0, 0.0]
    runtime_arrays.dulist_i[1, :] = [0.0, 0.0, 0.0]

    for j = 1:snap.twojmax
        jju = snap.prebuilt_arrays.idxu_block[j+1];
        jjup = snap.prebuilt_arrays.idxu_block[j];
        mb = 0
        while 2*mb <= j

            runtime_arrays.dulist_r[jju+1, :] = [0.0, 0.0, 0.0]
            runtime_arrays.dulist_i[jju+1, :] = [0.0, 0.0, 0.0]

            for ma = 0:(j-1)
                rootpq = snap.prebuilt_arrays.rootpqarray[j - ma+1, j - mb+1];
                for k = 1:3
                    runtime_arrays.dulist_r[jju+1, k] +=
                        rootpq * (da_r[k] * ulist_r[jjup+1] +
                        da_i[k] * ulist_i[jjup+1] +
                        a_r * runtime_arrays.dulist_r[jjup+1, k] +
                        a_i * runtime_arrays.dulist_i[jjup+1, k]);
                    runtime_arrays.dulist_i[jju+1, k] +=
                        rootpq * (da_r[k] * ulist_i[jjup+1] -
                        da_i[k] * ulist_r[jjup+1] +
                        a_r * runtime_arrays.dulist_i[jjup+1, k] -
                        a_i * runtime_arrays.dulist_r[jjup+1, k]);
                end

                rootpq = snap.prebuilt_arrays.rootpqarray[ma + 1 + 1,j - mb + 1];
                for k = 1:3
                    runtime_arrays.dulist_r[jju+1+1, k] =
                        -rootpq * (db_r[k] * ulist_r[jjup+1] +
                        db_i[k] * ulist_i[jjup+1] +
                        b_r * runtime_arrays.dulist_r[jjup+1, k] +
                        b_i * runtime_arrays.dulist_i[jjup+1, k]);
                    runtime_arrays.dulist_i[jju+1+1, k] =
                        -rootpq * (db_r[k] * ulist_i[jjup+1] -
                        db_i[k] * ulist_r[jjup+1] +
                        b_r * runtime_arrays.dulist_i[jjup+1, k] -
                        b_i * runtime_arrays.dulist_r[jjup+1, k]);
                end

                jju+=1;
                jjup+=1;
            end
            jju+=1
            mb +=1;
        end

        #  // copy left side to right side with inversion symmetry VMK 4.4(2)
        #  // u[ma-j][mb-j] = (-1)^(ma-mb)*Conj([u[ma][mb])

        jju = snap.prebuilt_arrays.idxu_block[j+1];
        jjup = jju+(j+1)*(j+1)-1;
        mbpar = 1;
        mb = 0
        while 2*mb <= j
            mapar = mbpar;
            for ma = 0:j
                if (mapar == 1)
                    for k = 1:3
                        runtime_arrays.dulist_r[jjup+1, k] = runtime_arrays.dulist_r[jju+1, k];
                        runtime_arrays.dulist_i[jjup+1, k] = -runtime_arrays.dulist_i[jju+1, k];
                    end
                else
                    for k = 1:3
                        runtime_arrays.dulist_r[jjup+1, k] = -runtime_arrays.dulist_r[jju+1, k];
                        runtime_arrays.dulist_i[jjup+1, k] = runtime_arrays.dulist_i[jju+1, k];
                    end
                end
                mapar = -mapar;
                jju+=1;
                jjup-=1;
            end
            mbpar = -mbpar;
            mb +=1;
        end
    end

    sfac = compute_sfac(snap, r, rcut);
    dsfac = compute_dsfac(snap, r, rcut);
    sfac *= wj;
    dsfac *= wj;
    for j = 0:snap.twojmax
        jju = snap.prebuilt_arrays.idxu_block[j+1];
        mb = 0;
        while 2*mb <= j
            for ma = 0:j
                runtime_arrays.dulist_r[jju+1,1] = dsfac * ulist_r[jju+1] * ux +
                                sfac * runtime_arrays.dulist_r[jju+1, 1];
                runtime_arrays.dulist_i[jju+1,1] = dsfac * ulist_i[jju+1] * ux +
                                sfac * runtime_arrays.dulist_i[jju+1, 1];
                runtime_arrays.dulist_r[jju+1,2] = dsfac * ulist_r[jju+1] * uy +
                                sfac * runtime_arrays.dulist_r[jju+1, 2];
                runtime_arrays.dulist_i[jju+1,2] = dsfac * ulist_i[jju+1] * uy +
                                sfac * runtime_arrays.dulist_i[jju+1, 2];
                runtime_arrays.dulist_r[jju+1,3] = dsfac * ulist_r[jju+1] * uz +
                                sfac * runtime_arrays.dulist_r[jju+1, 3];
                runtime_arrays.dulist_i[jju+1,3] = dsfac * ulist_i[jju+1] * uz +
                                sfac * runtime_arrays.dulist_i[jju+1, 3];
                jju+=1;
            end
            mb +=1;
        end
    end
end

function compute_duidrj(rij::SVector{3, <:AbstractFloat}, wj::AbstractFloat,
        rcut::AbstractFloat, jj::Int,
        snap::SNAPParams, runtime_arrays::RuntimeArrays)
    x = rij[1];
    y = rij[2];
    z = rij[3];
    rsq = x * x + y * y + z * z;
    r = sqrt(rsq);
    rscale0 = snap.rfac0 * pi / (rcut - snap.rmin0);
    theta0 = (r - snap.rmin0) * rscale0;
    cs = cos(theta0);
    sn = sin(theta0);
    z0 = r * cs / sn;
    dz0dr = z0 / r - (r*rscale0) * (rsq + z0 * z0) / rsq;

    compute_duarray(x, y, z,
        z0, r, dz0dr,
        wj, rcut, jj,
        snap, runtime_arrays)
end

function compute_sfac(snap::SNAPParams, r::AbstractFloat, rcut::AbstractFloat)
    if ~snap.switch_flag
        sfac = 1.0
    else
        if r <= snap.rmin0
            sfac = 1.0
        elseif r > rcut
            sfac =  0.0
        else
            rcutfac = pi / (rcut - snap.rmin0)
            sfac = 0.5 * ( cos( (r - snap.rmin0) * rcutfac ) + 1.0 )
        end
    end
    return sfac
end

function compute_dsfac(snap::SNAPParams, r::AbstractFloat, rcut::AbstractFloat)
    if ~snap.switch_flag
        dsfac = 0.0
    else
        if r <= snap.rmin0
            dsfac = 0.0
        elseif rcut > rcut
            dsfac = 0.0
        else
            rcutfac = pi / (rcut - snap.rmin0)
            dsfac = -0.5 * ( sin( (r - snap.rmin0) * rcutfac ) * rcutfac )
        end
    end
    return dsfac
end


text:
@snn_kw struct HHParameter{FT=Float32}
    Cm::FT = 1uF * cm^(-2) * 20000um^2
    gl::FT = 5e-5siemens * cm^(-2) * 20000um^2
    El::FT = -65mV
    Ek::FT = -90mV
    En::FT = 50mV
    gn::FT = 100msiemens * cm^(-2) * 20000um^2
    gk::FT = 30msiemens * cm^(-2) * 20000um^2
    Vt::FT = -63mV
    τe::FT = 5ms
    τi::FT = 10ms
    Ee::FT = 0mV
    Ei::FT = -80mV
end

@snn_kw mutable struct HH{VFT=Vector{Float32},VBT=Vector{Bool}}
    param::HHParameter = HHParameter()
    N::Int32 = 100
    v::VFT = param.El .+ 5(randn(N) .- 1)
    m::VFT = zeros(N)
    n::VFT = zeros(N)
    h::VFT = ones(N)
    ge::VFT  = (1.5randn(N) .+ 4) .* 10nS
    gi::VFT  = (12randn(N) .+ 20) .* 10nS
    fire::VBT = zeros(Bool, N)
    I::VFT = zeros(N)
    records::Dict = Dict()
end

"""
[Hodgkin–Huxley Neuron](https://en.wikipedia.org/wiki/Hodgkin%E2%80%93Huxley_model)
"""
HH

function integrate!(p::HH, param::HHParameter, dt::Float32)
    @unpack N, v, m, n, h, ge, gi, fire, I = p
    @unpack Cm, gl, El, Ek, En, gn, gk, Vt, τe, τi, Ee, Ei = param
    @inbounds for i = 1:N
        m[i] += dt * (0.32f0 * (13f0 - v[i] + Vt) / (exp((13f0 - v[i] + Vt) / 4f0) - 1f0) * (1f0 - m[i]) -
        0.28f0 * (v[i] - Vt - 40f0) / (exp((v[i] - Vt - 40f0) / 5f0) - 1f0) * m[i])
        n[i] += dt * (0.032f0 * (15f0 - v[i] + Vt) / (exp((15f0 - v[i] + Vt) / 5f0) - 1f0) * (1f0 - n[i]) -
        0.5f0 * exp((10f0 - v[i] + Vt) / 40f0) * n[i])
        h[i] += dt * (0.128f0 * exp((17f0 - v[i] + Vt) / 18f0) * (1f0 - h[i]) -
        4f0 / (1f0 + exp((40f0 - v[i] + Vt) / 5f0)) * h[i])
        v[i] += dt / Cm * ( I[i] + gl * (El - v[i]) + ge[i] * (Ee - v[i]) + gi[i] * (Ei - v[i]) +
        gn * m[i]^3 * h[i] * (En - v[i]) + gk * n[i]^4 * (Ek - v[i]) )
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > -20f0
    end
end


text:
module BoundaryCellFields

using Gridap
using Gridap.Helpers
using Gridap.CellFields: NonIterableCellFieldLike

export restrict
import Base: length
import Gridap: evaluate
import Gridap: gradient

function restrict(cf::IndexCellFieldLike,desc::BoundaryDescriptor)
  BoundaryCellFieldLike(cf,desc)
end

function restrict(cf::IndexCellFieldLike{D},trian::BoundaryTriangulation{Z}) where {D,Z}
  @assert D == Z + 1
  restrict(cf,trian.descriptor)
end

struct BoundaryCellFieldLike{Z,T,N,D} <: NonIterableCellFieldLike{Z,T,N}
  cellfield::CellFieldLike{D,T,N}
  descriptor::BoundaryDescriptor
end

function BoundaryCellFieldLike(
  cellfield::IndexCellFieldLike{D,T,N}, descriptor::BoundaryDescriptor) where {D,T,N}
  Z = D - 1
  BoundaryCellFieldLike{Z,T,N,D}(cellfield,descriptor)
end

function gradient(cf::BoundaryCellFieldLike)
  cfg = gradient(cf.cellfield)
  BoundaryCellFieldLike(cfg,cf.descriptor)
end

length(bcf::BoundaryCellFieldLike) = length(bcf.descriptor.facet_to_cell)

function evaluate(bcf::BoundaryCellFieldLike{Z}, q::CellPoints{Z}) where Z

  facet_qpoint_to_qcoord = q
  facet_to_lfacet = bcf.descriptor.facet_to_lfacet
  cell_to_polytope = bcf.descriptor.cell_to_polytope
  facet_to_cell = bcf.descriptor.facet_to_cell

  r = coordinates_in_ref_cells(
    facet_qpoint_to_qcoord,
    facet_to_lfacet,
    cell_to_polytope,
    facet_to_cell)

  cf = reindex(bcf.cellfield,facet_to_cell)

  evaluate(cf,r)

end

function coordinates_in_ref_cells(
  facet_qpoint_to_qcoord,
  facet_to_lfacet,
  cell_to_polytope,
  facet_to_cell)

  _coordinates_in_ref_cells(
    facet_qpoint_to_qcoord,
    facet_to_lfacet,
    cell_to_polytope,
    facet_to_cell)

end

function  _coordinates_in_ref_cells(
  facet_qpoint_to_qcoord,
  facet_to_lfacet,
  cell_to_polytope,
  facet_to_cell)
  @notimplemented
end

function  _coordinates_in_ref_cells(
  facet_qpoint_to_qcoord::ConstantCellValue,
  facet_to_lfacet,
  cell_to_polytope::ConstantCellValue,
  facet_to_cell)

  polytope = cell_to_polytope.value
  D = length(polytope.extrusion)
  grid = Grid(polytope,D-1)
  trian = Triangulation(grid)
  phi = CellGeomap(trian)
  nlfacets = length(phi)

  _qpoint_to_qcoord = facet_qpoint_to_qcoord.value
  qpoint_to_qcoord = ConstantCellValue(_qpoint_to_qcoord,nlfacets)
  _lfacet_qpoint_to_rcoord = evaluate(phi,qpoint_to_qcoord)

  lfacet_qpoint_to_rcoord = collect(_lfacet_qpoint_to_rcoord)
  facet_qpoint_to_rcoord = CompressedCellArray(lfacet_qpoint_to_rcoord,facet_to_lfacet)

  facet_qpoint_to_rcoord
end

end # module


text:
using BioinformaticsBISC195
using Test

@testset "BioinformaticsBISC195" begin

@testset "Using Strings" begin

    @testset "normalizeDNA" begin
        @test normalizeDNA("aatgn") == "AATGN"
        @test normalizeDNA("C") isa String
        @test_throws Exception normalizeDNA("ZCA")
        @test_throws Exception normalizeDNA(42)
        @test_throws Exception normalizeDNA('C')
    end # normalizeDNA

    @testset "composition" begin
        seq = rand(['A','T','G','C','N'], 20) |> join
        bc = composition(seq)
        @test bc isa Dict

        @test bc['A'] == count(x-> x == 'A', seq)
        @test bc['C'] == count(x-> x == 'C', seq)
        @test bc['G'] == count(x-> x == 'G', seq)
        @test bc['T'] == count(x-> x == 'T', seq)
        @test bc['N'] == count(x-> x == 'N', seq)

        bc = composition(lowercase(seq))

        @test bc['A'] == count(x-> x == 'A', seq)
        @test bc['C'] == count(x-> x == 'C', seq)
        @test bc['G'] == count(x-> x == 'G', seq)
        @test bc['T'] == count(x-> x == 'T', seq)
        @test bc['N'] == count(x-> x == 'N', seq)
    end # composition

    @testset "gc_content" begin
        @test gc_content("ANTG") == 1/3
        @test gc_content("cccggg") * 100 == 100.0
        @test gc_content("ATta") == 0.0
        @test_throws Exception gc_content("ATtq")
    end # gc_content

    @testset "complement" begin
        @test complement("ATTAN") == "TAATN"
        @test complement("gcta") == "CGAT"
        @test complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception complement("AxC")
    end # complement

    @testset "reverse_complement" begin
        @test reverse_complement("ATTAN") == "NTAAT"
        @test reverse_complement("gcta") == "TAGC"
        @test reverse_complement("nnnnnnn") == "NNNNNNN"
        @test_throws Exception reverse_complement("AxC")
    end # reverse_complement

    @testset "parse_fasta" begin
        testpath = normpath(joinpath(@__DIR__, "..", "data"))
        genomes = joinpath(testpath, "cov2_genomes.fasta")
        ex1_path = joinpath(testpath, "ex1.fasta")
        ex2_path = joinpath(testpath, "ex2.fasta")

        ex1 = parse_fasta(ex1_path)
        @test ex1 isa Tuple
        @test all(x-> x isa String, ex1[1])
        @test all(x-> x isa String, ex1[2])

        @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
        @test ex1[2] == ["AATTATAGC", "CGCCCCCCAGTCGGATT"]

        @test_throws Exception parse_fasta(ex2_path)

        cov2 = parse_fasta(genomes)
        @test length(cov2[1]) == 8
        @test length(cov2[2]) == 8
    end #parse_fasta

    @testset "getKmers" begin
        @test getKmers("ggg", 3) == Set(["GGG"])
        @test getKmers("ATATATATA", 4) == Set(["TATA", "ATAT"])
        @test getKmers("GTAGAGCTGT", 6) == Set(["GTAGAG", "TAGAGC", "AGAGCT", "GAGCTG", "AGCTGT"])
        @test getKmers("GTAGAGCTGT", 8) == Set(["GTAGAGCT", "TAGAGCTG", "AGAGCTGT"])
        @test_throws Exception getKmers("A", 2)
        @test_throws Exception getKmers("CXG", 2)
    end # getKmers

    @testset "getKmerDist" begin
        kmer1, kmer2, kmer3 = Set(["GGG", "CCAT"]), Set(["ATAT", "GGG", "ATCG"]), Set(["TAG"])
        @test getKmerDist(kmer1, kmer2) == 0.75
        @test getKmerDist(kmer1, kmer3) == 1
        @test getKmerDist(kmer2, kmer2) == 0
    end # getKmerDist

    @testset "getKmerCount" begin
        @test getKmerCount("ggg", 3) == Dict("GGG" => 1)
        @test getKmerCount("ATATATATA", 4) == Dict("TATA" => 3, "ATAT" => 3)
        @test getKmerCount("GTAGAGCTGT", 6) == Dict("GTAGAG" => 1, "TAGAGC" => 1, "AGAGCT" => 1, "GAGCTG" => 1, "AGCTGT" => 1)
        @test getKmerCount("GTAGAGCTGT", 8) == Dict("GTAGAGCT" => 1, "TAGAGCTG" => 1, "AGAGCTGT" => 1)
        @test_throws Exception getKmerCount("A", 2)
    end # getKmerCount

    @testset "getHeaderAttrib" begin
        @test getHeaderAttrib("NC_045512.2 |Severe acute respiratory syndrome-related coronavirus||China|2019-12", "|", [4,5]) == ["China", "2019-12"]
        @test getHeaderAttrib("MZ413975.1 |Severe acute respiratory syndrome-related coronavirus|oronasopharynx|Bangladesh: Sylhet|2020-12-31", "|", [4,5]) == ["Bangladesh: Sylhet", "2020-12-31"]
        @test getHeaderAttrib("MZ042984.1 |Severe acute respiratory syndrome-related coronavirus|oronasopharynx|Egypt|2021-04-09", "|", [1,3,5]) == ["MZ042984.1", "oronasopharynx", "2021-04-09"]
        @test_throws Exception getHeaderAttrib("A|B|C", 6)
    end # getHeaderAttrib

    @testset "filterByEmptyAttrib" begin
        testHeader = ["T|T|T||T", "T|T||T|T", "T||T||T"]
        testSeq = ["1", "2", "3"]
        noDateHeader = ["MZ562753.1 |Severe acute respiratory syndrome-related coronavirus||India: Madhya Pradesh|2021-02-08", "AY278554.2 |Severe acute respiratory syndrome-related coronavirus||China: Hong Kong, Prince of Wales Hospital|"]
        noDateSeq = ["AATG", "GTAC"]

        @test filterByEmptyAttrib(testHeader, testSeq, "|", 1) == (["T|T|T||T", "T|T||T|T", "T||T||T"], ["1", "2", "3"])
        @test filterByEmptyAttrib(testHeader, testSeq, "|", 2) == (["T|T|T||T", "T|T||T|T"], ["1", "2"])
        @test filterByEmptyAttrib(testHeader, testSeq, "|", 3) == (["T|T|T||T", "T||T||T"], ["1", "3"])
        @test filterByEmptyAttrib(testHeader, testSeq, "|", 4) == (["T|T||T|T"], ["2"])
        @test filterByEmptyAttrib(noDateHeader, noDateSeq, "|", 5) == (["MZ562753.1 |Severe acute respiratory syndrome-related coronavirus||India: Madhya Pradesh|2021-02-08"], ["AATG"])
        @test_throws Exception filterByEmptyAttrib(noDateHeader, testSeq, "|", 3)
    end # filterByEmptyAttrib

end # strings

# @testset "Using BioSequences" begin

#     @testset "normalizeDNA" begin
#         @test normalizeDNA("aatgn") == dna"AATGN"
#         @test_throws Exception normalizeDNA("ZCA")
#         @test_throws Exception normalizeDNA(42)
#         c = normalizeDNA('C')
#         @test c == dna"c"
#         @test c isa LongSequence
#     end #  normalizeDNA

#     @testset "gc_content" begin
#         @test gc_content(dna"ANTG") == 0.25
#         @test gc_content(dna"cccggg") * 100 == 100.0
#         @test gc_content(dna"ATta") == 0.0
#     end #  composition

#     @testset "composition" begin
#         seq = rand(['A','T','G','C','N'], 20) |> join |> LongDNASeq
#         bc = composition(seq)

#         @test bc[DNA_A] == count(==(DNA_A), collect(seq))
#         @test bc[DNA_C] == count(==(DNA_C), collect(seq))
#         @test bc[DNA_G] == count(==(DNA_G), collect(seq))
#         @test bc[DNA_T] == count(==(DNA_T), collect(seq))
#         @test bc[DNA_N] == count(==(DNA_N), collect(seq))
#     end #  gc_content

#     @testset "complement" begin
#         @test complement(dna"ATTAN") == dna"TAATN"
#         @test complement(dna"gcta") == dna"CGAT"
#         @test complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  complement

#     @testset "reverse_complement" begin
#         @test reverse_complement(dna"ATTAN") == dna"NTAAT"
#         @test reverse_complement(dna"gcta") == dna"TAGC"
#         @test reverse_complement(dna"nnnnnnn") == dna"NNNNNNN"
#     end #  reverse_complement

#     @testset "parse_fasta" begin
#         testpath = normpath(joinpath(@__DIR__, "..", "data"))
#         genomes = joinpath(testpath, "cov2_genomes.fasta")
#         ex1_path = joinpath(testpath, "ex1.fasta")
#         ex2_path = joinpath(testpath, "ex2.fasta")

#         ex1 = parse_fasta(ex1_path)
#         @test ex1 isa Tuple
#         @test all(x-> x isa String, ex1[1])
#         @test all(x-> x isa LongSequence, ex1[2])

#         @test ex1[1] == ["ex1.1 | easy", "ex1.2 | multiline"]
#         @test ex1[2] == [dna"AATTATAGC", dna"CGCCCCCCAGTCGGATT"]

#         @test_throws Exception parse_fasta(ex2_path)

#         cov2 = parse_fasta(genomes)
#         @test length(cov2[1]) == 8
#         @test length(cov2[2]) == 8
#     end # parse_fasta

# end # BioSequences

end # BioinformaticsBISC195


text:
_process_sp(t,tN,data,encParams,sp,display_evaluation=identity)= begin
  # z: encoder activation
  z,bucket= encode_powerDay(data.power_hourly_kw[t], data.hour[t], data.is_weekend[t];
                       encParams...)
  # a: SP activation
  a= step!(sp,z)
  t%(tN÷10)==0 && display_evaluation(t,sp,a,sp.params.szₛₚ)
  return z,a,bucket
end
_process_tm(t,tN, tm,a)= step!(tm,a)
identity(a...)= a
read_gympower(dataset_filename= "test/test_data/gym_power_benchmark.csv")= begin
  data= CSV.read(dataset_filename, DataFrame)#; truestrings=["1"], falsestrings=["0"], types=(String,Float16,Int,String))
  ((power_hourly_kw= data.power_hourly_kw,
    hour= data.hour,
    is_weekend= data.is_weekend
   ), size(data,1))
end
"""
Find the top-k% most similar SDRs to the last one from the entire history.
"""
top_similar_sdr(sdrHist,k)= begin
  k= size(sdrHist,2)==1 ? 0 : ceil(Int,k/100*size(sdrHist,2))
  # Shame: I can't find ANY high-level construct to even approach the efficiency of this!
  overlap= zeros(Int32,size(sdrHist,2)-1)
  for t in 1:length(overlap) for i in 1:size(sdrHist,1)
      overlap[t]+= sdrHist[i,end] & sdrHist[i,t] end end
  topIdx= @> overlap partialsortperm(1:k,rev=true)
  (sdrHist[:,topIdx], topIdx)
end

# Calculate prediction "mean absolute scaled error"
# (https://en.wikipedia.org/wiki/Mean_absolute_scaled_error)
function mase(data,prediction,predict_timesteps)
  aligned_data= data[1+predict_timesteps:end]
  aligned_pred= prediction[1:end-predict_timesteps]
  mean(abs.(aligned_data .- aligned_pred)) / mean(abs.(diff(aligned_data)))
end

plot_ts_similarEncSp(t,ts,encOnly,spOnly,encANDspHistory)=
    if t>1 _plot_ts_similarEncSP(t,ts,encOnly,spOnly,encANDspHistory)
    end

function _plot_ts_similarEncSP(t,ts,encOnly,spOnly,encANDspHistory)
  graph_ts()= begin
    graph_ts= plot(ts, label="power", ylims=(0.9minimum(ts),1.05maximum(ts)),
              xlims=(1,1.3length(ts)))
    scatter!(graph_ts, encOnly, ts[encOnly], label="encoding only")
    scatter!(graph_ts, spOnly, ts[spOnly], label="SP only")
    scatter!(graph_ts, encANDspHistory[t][1], ts[encANDspHistory[t][1]],
              label="overlapping", legendfontsize=6)
    vline!(graph_ts,[t], label="")
    title!(graph_ts, "Spatial Pooler mapping property evaluation", titlefont=font(10))
    return graph_ts
  end
  graph_ovp(overlapFraction)= begin
    graph_ovp= plot(1:t,overlapFraction*100, label="", linestyle=:dash,
        ylabel="overlap%",ylims=(-1,102),xlims=(1,1.3length(ts)))
    title!(graph_ovp, "Percentage of overlapping encoder and SP SDRs", titlefont=font(10))
    return graph_ovp
  end
  overlapFraction= map(x->length(x.encANDsp)/x.Nenc, encANDspHistory[1:t])
  overlapFraction[isnan.(overlapFraction)].= 1.0
  plot(graph_ts(),graph_ovp(overlapFraction), layout=(2,1))|> display
end


text:
mutable struct Timer
    idx::Int
    last5::Vector{DateTime}

    Timer() = new(1, zeros(DateTime, 5))
end

function _ratelimit!(t::Timer)
    n = now()
    e = t.last5[t.idx]
    rl = n - e
    if rl < Millisecond(1000)
        @warn "Pausing for rate-limiting"
        sleep(Millisecond(1000) - rl)
        n = now()
    end

    t.last5[t.idx] = n
    t.idx = 1 + (t.idx % 5)
end

const _ratelimiter = Timer()

@testset "Rate Limiting" begin
    tab = AirTable("Table 1", AirBase("appphImnhJO8AXmmo"))
    # Make > 5 requests in a second, and be sure no error is thrown
    for i in 1:7
        q = Airtable.query(tab)
        if i == 7
            @test q isa Vector{AirRecord}
        end
    end

end


text:
using XPORTA: IEQ

export AbstractGame, Game, BellGame

"""
An `AbstractGame` is the abstract type that is parent to type representing a cost
function for strategies.
A `Game` is a linear inequality dual to `Strategy` matrices, it is represented by
a matrix ``G`` containing scalar coefficients and a bound ``\\beta``.
The bound of the linear inequality is typically a maximum score attainable by Bell
scenario using a certain set of resources.
A `Game` ``G`` is played by a `Strategy` ``S`` to achieve a score computed as,

```math
\\beta \\geq \\langle G, S\\rangle = \\sum_{x,y} G_{y,x}S(y|x).
```

The game is "won" if the strategy scores greater than ``\\beta``, that is, the above
inequality is violated.
Since ``\\beta`` is the maximum score for a set of resources, the game is won only
if the tested strategy used a set of resources of greater operational value than
than considered when computing the bound ``\\beta``.

Conveniently, `Strategy` matrices have normalized columns and non-negative elements.
This means that any game inequality can be converted into a form where game matrix
``G`` has positive elements and ``\\beta`` designates a positive upper bound.
"""
abstract type AbstractGame{T} <: AbstractMatrix{T} end
Base.size(G::AbstractGame) = size(G.game)
Base.getindex(G::AbstractGame, I::Vararg{Int64,2}) = getindex(G.game, I...)
Base.setindex!(G::AbstractGame, v, I::Vararg{Int64,2}) = (G.game[I...] = v)

"""
    Game(game::Matrix{T}, β::Real) <: AbstractGame{T}

A `Game` is represented by a `Matrix` of coefficients and a scalar bound `β`.

Type parameter `T` can be either `Int64` or `Float64`.
"""
struct Game{T} <: AbstractGame{T}
    game :: Matrix{T}
    β :: Real
    scenario :: Scenario
    Game(game::Matrix{T}, β::Real) where T = new{T}(game, β)
end

"""
    BellGame(game::Matrix{Int64}, β::Int64, scenario::Scenario)

A `BellGame` represents a Bell inequality or tight facet of the local polytope.
Since the vertices of the local polytope are deterministic strategies with 0,1 elements,
the linear inequalities describing facets of the local polytope have rational coefficients.
Therefore, if a inequality tightly bounds the local polytope, it can be represented by
a game with integer coefficients.
"""
struct BellGame <: AbstractGame{Int64}
    game :: Matrix{Int64}
    β :: Int64
    scenario :: Scenario
    BellGame(game::Matrix{Int64}, β::Int64) = new(game, β)
end


text:
using OrdinaryDiffEq
using Trixi

# define new structs inside a module to allow re-evaluating the file
module TrixiExtensionExample

using Trixi
using DiffEqCallbacks: DiscreteCallback, u_modified!

# This is an example implementation for a simple stage callback (i.e., a callable
# that is executed after each Runge-Kutta *stage*), which records some values
# each time it is called. Its sole purpose here is to showcase how to implement
# a stage callback for Trixi.
struct ExampleStageCallback
  times::Vector{Float64}
  min_values::Vector{Float64}
  max_values::Vector{Float64}

  # You can optionally define an inner constructor like the one below to set up
  # some required stuff. You can also create outer constructors (not demonstrated
  # here) for further customization options.
  function ExampleStageCallback()
    new(Float64[], Float64[], Float64[])
  end
end

# This method is called when the `ExampleStageCallback` is used as `stage_limiter!`
# which gets called after every RK stage. There is no specific initialization
# method for such `stage_limiter!`s in OrdinaryDiffEq.jl.
function (example_stage_callback::ExampleStageCallback)(u_ode::AbstractVector, _, semi, t)

  min_val, max_val = extrema(u_ode)
  push!(example_stage_callback.times, t)
  push!(example_stage_callback.min_values, min_val)
  push!(example_stage_callback.max_values, max_val)

  return nothing
end

# This is an example implementation for a simple step callback (i.e., a callable
# that is potentially executed after each Runge-Kutta *step*), which records
# some values each time it is called. Its sole purpose here is to showcase
# how to implement a step callback for Trixi.
struct ExampleStepCallback
  message::String
  times::Vector{Float64}
  min_values::Vector{Float64}
  max_values::Vector{Float64}

  # You can optionally define an inner constructor like the one below to set up
  # some required stuff. You can also create outer constructors (not demonstrated
  # here) for further customization options.
  function ExampleStepCallback(message::String)
    new(message, Float64[], Float64[], Float64[])
  end
end

# This method is called when the `ExampleStepCallback` is used as callback
# which gets called after RK steps.
function (example_callback::ExampleStepCallback)(integrator)
  u_ode = integrator.u
  t = integrator.t
  # You can also access semi = integrator.p

  min_val, max_val = extrema(u_ode)
  push!(example_callback.times, t)
  push!(example_callback.min_values, min_val)
  push!(example_callback.max_values, max_val)

  # avoid re-evaluating possible FSAL stages
  u_modified!(integrator, false)
  return nothing
end

# This method is used to wrap an `ExampleStepCallback` inside a `DiscreteCallback`
# which gets called after every RK step. You can pass an additional initialization
# method and a separate condition specifying whether the callback shall be called.
function ExampleStepCallback(; message::String)
  # Call the `ExampleStepCallback` after every RK step.
  condition = (u_ode, t, integrator) -> true

  # You can optionally pass an initialization method. There, you can access the
  # `ExampleStepCallback` as `cb.affect!`.
  initialize = (cb, u_ode, t, integrator) -> println(cb.affect!.message)

  example_callback = ExampleStepCallback(message)

  DiscreteCallback(condition, example_callback,
                   save_positions=(false,false),
                   initialize=initialize)
end

end # module TrixiExtensionExample

import .TrixiExtensionExample

###############################################################################
# semidiscretization of the linear advection equation

advectionvelocity = (1.0, 1.0)
# advectionvelocity = (0.2, -0.3)
equations = LinearScalarAdvectionEquation2D(advectionvelocity)

initial_condition = initial_condition_convergence_test

surface_flux = flux_lax_friedrichs
solver = DGSEM(3, surface_flux)

coordinates_min = (-1, -1)
coordinates_max = ( 1,  1)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                n_cells_max=30_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan);

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_integrals=(entropy, energy_total))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2cons)

example_callback = TrixiExtensionExample.ExampleStepCallback(message="안녕하세요?")

stepsize_callback = StepsizeCallback(cfl=1.6)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution,
                        example_callback,
                        stepsize_callback)

# In OrdinaryDiffEq.jl, the `step_limiter!` is called after every Runge-Kutta step
# but before possible RHS evaluations of the new value occur. Hence, it's possible
# to modify the new solution value there without impacting the performance of FSAL
# methods.
# The `stage_limiter!` is called additionally after computing a Runge-Kutta stage
# value but before evaluating the corresponding stage derivatives.
# Hence, if a limiter should be called before each RHS evaluation, it needs to be
# set both as `stage_limiter!` and as `step_limiter!`.
example_stage_callback! = TrixiExtensionExample.ExampleStageCallback()
stage_limiter! = example_stage_callback!
step_limiter!  = example_stage_callback!

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, step_limiter!, williamson_condition=false),
                  dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                  save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

# Check whether we recorded the same values.
# Remember that CarpenterKennedy2N54 has five stages per step.
@assert example_stage_callback!.times[5:5:end] ≈ example_callback.affect!.times
@assert example_stage_callback!.min_values[5:5:end] ≈ example_callback.affect!.min_values
@assert example_stage_callback!.max_values[5:5:end] ≈ example_callback.affect!.max_values


text:
module WingBiology

    using LinearAlgebra
    using ForwardDiff
    using SparseArrays
    using SparseDiffTools
    using LinearMaps
    using Arpack

    using PyCall

    include("FuselageAero.jl")
    include("WingAero.jl")
    include("Engine.jl")

    include("StructuralAnalysis.jl")

    include("DataTypes.jl")
    include("Plotting.jl")

    include("StateSpace.jl")
    include("Newton.jl")

    _linv_reg(A::AbstractMatrix, λ::Real = 1e-3) = let D = A' * A
        inv(
            D .+ Matrix{Float64}(I * λ, size(D)...)
        ) * A'
    end

    include("Arnoldi.jl")
    include("AssumedModes.jl")

end # module

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology