text:
export @ein_str, @ein, IndexSize
"""
    ein"ij,jk -> ik"(A,B)

String macro interface which understands `numpy.einsum`'s notation.
Translates strings into `EinCode`-structs that can be called to evaluate
an `einsum`.
To control evaluation order, use parentheses - instead of an `EinCode`,
a `NestedEinsumStable` is returned which evaluates the expression
according to parens.
The valid character ranges for index-labels are `a-z` and `α-ω`.

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b, c = rand(10,10), rand(10,10), rand(10,1);

julia> ein"ij,jk,kl -> il"(a,b,c) ≈ ein"(ij,jk),kl -> il"(a,b,c) ≈ a * b * c
true
```
"""
macro ein_str(s::AbstractString)
    s = replace(s, " " => "")
    m = match(r"([\(\)a-z,α-ω]*)->([a-zα-ω]*)", s)
    m == nothing && throw(ArgumentError("invalid einsum specification $s"))
    sixs, siy = m.captures
    if '(' in sixs
        return parse_nested(sixs, collect(siy))
    else
        iy  = Tuple(siy)
        ixs = Tuple(Tuple(ix) for ix in split(sixs,','))
        return EinCode(ixs, iy)
    end
end

function (code::EinCode)(xs...; size_info=nothing) where {T, N}
    size_dict = get_size_dict(getixs(code), xs)
    if !(size_info isa Nothing)
        size_dict += size_info
    end
    einsum(code, xs, size_dict)
end

einsum(code::EinCode{ixs, iy}, xs) where {ixs, iy} = einsum(code, xs, nothing)

function einsum(code::EinCode{ixs, iy}, xs, ::Nothing) where {ixs, iy}
    einsum(code, xs, get_size_dict(ixs, xs))
end

@doc raw"
    get_size_dict(ixs, xs)

return the `IndexSize` struct that is used to get the size of an index-label
in the einsum-specification with input-indices `ixs` and tensors `xs` after
consistency within `ixs` and between `ixs` and `xs` has been verified.
"
function get_size_dict(ixs::NTuple{N, NTuple{M, T} where M}, xs::NTuple{X, AbstractArray}) where {N,T,X}
    # check size of input tuples
    N != X && throw(ArgumentError("$X tensors labelled by $N indices"))

    # check tensor orders
    foreach(ixs, xs) do ix, x
        length(ix) == ndims(x) || throw(
            ArgumentError("indices $ix invalid for tensor with ndims = $(ndims(x))"))
    end
    sd = getindexsize(ixs, xs)
    check_dimensions(sd)
    return sd
end

@doc raw"
    IndexSize{N,T}(k::NTuple{N,T},v::NTuple{N,Int})

struct to hold the size of indices specified by their labels.
Note that while a dict would work, for the small sizes we usually
have, a tuple of keys and values is much faster to construct
and competitive for lookup.
"
struct IndexSize{N,T}
    indices::NTuple{N,T}
    sizes::NTuple{N,Int}
end

IndexSize(sizes::Pair...) = IndexSize(first.(sizes), last.(sizes))
IndexSize(indices::Tuple{}, sizes::Tuple{}) = IndexSize{0,Any}(indices, sizes)

Base.:+(x::IndexSize, y::IndexSize) = IndexSize((x.indices..., y.indices...), (x.sizes..., y.sizes...))

function getindexsize(ixs, xs)
    indices = TupleTools.flatten(ixs)
    sizes = TupleTools.flatten(map(size,xs))
    IndexSize(indices,sizes)
end

Base.getindex(inds::IndexSize{N,T},i::T) where {N,T} = inds.sizes[findfirst(==(i), inds.indices)]

@doc raw"
    check_dimensions(inds::IndexSize)

check whether all non-unique indexlabels point to the same
dimension - otherwise throw an error.
"
function check_dimensions(inds::IndexSize)
    for (c,i) in enumerate(inds.indices)
        j = findnext(==(i), inds.indices, c+1)
        j != nothing && inds.sizes[c] != inds.sizes[j] && return throw(
            DimensionMismatch("index $i has incosistent sizes"))
    end
    return true
end

using MacroTools
"""
    @ein A[i,k] := B[i,j] * C[j,k]     # A = B * C

Macro interface similar to that of other packages.

You may use numbers in place of letters for dummy indices, as in `@tensor`,
and need not name the output array. Thus `A = @ein [1,2] := B[1,ξ] * C[ξ,2]`
is equivalent to the above. This can also be written `A = ein"ij,jk -> ik"(B,C)`
using the numpy-style string macro.

# example

```jldoctest; setup = :(using OMEinsum)
julia> a, b = rand(2,2), rand(2,2);

julia> @ein c[i,k] := a[i,j] * b[j,k];

julia> c ≈ a * b
true

julia> c ≈ ein"ij,jk -> ik"(a,b)
true
```
"""
macro ein(exs...)
    _ein_macro(exs...)
end

primefix!(ind) = map!(i -> @capture(i, (j_)') ? Symbol(j, '′') : i, ind, ind)

function _ein_macro(ex; einsum=:einsum)
    @capture(ex, (left_ := right_)) || throw(ArgumentError("expected A[] := B[]... "))
    @capture(left, Z_[leftind__] | [leftind__] ) || throw(
        ArgumentError("can't understand LHS, expected A[i,j] etc."))
    if Z===nothing
        @gensym Z
    end
    primefix!(leftind)

    rightind, rightpairs = [], []
    @capture(right, *(factors__)) || (factors = Any[right])
    for fact in factors
        @capture(fact, A_[Aind__]) || return _nested_ein_macro(ex)
        primefix!(Aind)
        append!(rightind, Aind)
        push!(rightpairs, (A, Aind) )
    end
    unique!(rightind)
    isempty(setdiff(leftind, rightind)) || throw(
        ArgumentError("some indices appear only on the left"))

    lefttuple = Tuple(indexin(leftind, rightind))
    righttuples = [ Tuple(indexin(ind, rightind)) for (A, ind) in rightpairs ]
    rightnames = [ esc(A) for (A, ind) in rightpairs ]

    return :( $(esc(Z)) = $einsum( EinCode{($(righttuples...),), $lefttuple}(), ($(rightnames...),)) )
end


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 24c00466-e86f-410c-bb60-6dcb2399b63a
# Load the required Julia packages -
using LinearAlgebra

# ╔═╡ ff56afbd-630f-42c3-9171-ccd808d5589f
md"""

### Background
High Fructose Corn Syrup (HFCS) is a common sweetener in many food products. As a sweetener, HFCS is often compared to granulated table sugar (sucrose is made up of one molecule of glucose and one molecule of fructose), but HFCS is easier to handle and cheaper than granulated table sugar. HFCS comes in two varieties, HFCS-42 and HFCS-55. HFCS-42 and HFCS-55 refer to dry weight fructose compositions of 42% and 55% respectively, the rest being glucose. HFCS-42 is mainly used for processed foods and breakfast cereals, whereas HFCS-55 is used mostly for production of soft drinks.

### Manufacturing Process
Corn is milled to extract corn starch and an "acid-enzyme" process is used, in which the corn-starch solution is acidified to begin breaking up the existing carbohydrates. High-temperature enzymes are added to further metabolize the starch and convert the resulting sugars to fructose.
The first enzyme, alpha-amylase, breaks the long starch chains down into shorter sugar chains (oligosaccharides).
Glucoamylase, a second enzyme, converts the oligosaccharides to glucose (a common six carbon sugar for biotechnology applications). The resulting solution is filtered to remove protein, then using activated carbon, and then demineralized using ion-exchange resins.

The purified solution is then run over immobilized xylose isomerase, which converts the sugars to ~50–52% glucose with some unconverted oligosaccharides and 42% fructose (HFCS-42), and again demineralized and again purified using activated carbon. Some is processed into HFCS-90 by liquid chromatography, and then mixed with HFCS-42 to form HFCS-55. The enzymes used in the process are made by microbial fermentation.
"""

# ╔═╡ c7c294a2-43bf-4ae2-a6cc-f917eb6857a3
md"""

### Model and Assumptions
Let's model the immobilized xylose isomerase step of the HCFS-42 process.

__Assumptions__:

* The _system_ is a well-mixed reactor holding the immobilized xylose isomerase with a single input and out stream
* The _system_ operates at steady-state
* The _system_ temperature and pressure are held constant
"""

# ╔═╡ 58aff590-027b-4e63-8375-7ea304d1292d
md"""

#### Total mass balance
The number of streams $\mathcal{S}$ = 2 and we are operating at steady state. Thus, the general total mass balance equation:

$$\sum_{s=1}^{\mathcal{S}}v_{s}\dot{m}_{s} = \frac{dm}{dt}$$

becomes:

$$\dot{m}_{1} - \dot{m}_{2} = 0$$
"""

# ╔═╡ 9ad11963-8b06-42dd-9a92-0310893c03fc
md"""

#### Species mass balance equations
The _system_ involves a single reaction (with rate $r_{1}$), and three chemical species, glucose (1), fructose (2) and unreacted oligosaccharides (3). Glucose is converted to fructose by the enzyme (catalyst) xylose isomerase. Thus, the total number of species $\mathcal{M}$ = 3 and the number of streams $\mathcal{S}$ = 2. Therefore, our weight fraction composition matrix $\mathbf{W}$ is a 3 $\times$ 2 array. The $\dot{\mathbf{m}}$ vector is a 2$\times$ 1 column vector, while the generation terms are a 3 $\times$ 1 column vector. Let the rate of generation of species $i$ be given as:

$$\dot{m}_{i,gen} = \sigma_{i1}r_{1}$$

###### Stoichiometric coefficients
The $\sigma_{ij}$ in the generation term is called a _stoichiometric coefficient_. Stoichiometric coefficients relate chemical species with reaction rates.

* if $\sigma_{ij}>0$: species _i_ is __produced__ by reaction _j_
* if $\sigma_{ij}=0$: species _i_ is not connected with reaction _j_
* if $\sigma_{ij}<0$: species _i_ is __consumed__ by reaction _j_

Putting all these together gives the species material balances in matrix-vector form:

$\begin{pmatrix}
w_{11} & -w_{12} \\
w_{21} & -w_{22} \\
w_{31} & -w_{32} \\
\end{pmatrix}
\begin{pmatrix}
\dot{m}_{1} \\
\dot{m}_{2} \\
\end{pmatrix} = -
\begin{pmatrix}
\sigma_{11}r_{1}\\
\sigma_{21}r_{1}\\
\sigma_{31}r_{1}\\
\end{pmatrix}$

or in compact notion:

$$\mathbf{W}\dot{\mathbf{m}}=\mathbf{b}$$

"""

# ╔═╡ f6dcdae8-e1a1-44da-b085-4859df008dec
# Let's start doing some Julia! (the fun part ...)

# Define mass flow rate vector (kg/hr) -
mdot = [100,100] # kg/hr

# ╔═╡ eb822e57-423a-4662-b402-052a616180d5
# define the composition array -
W = [0.92 -0.50 ; 0.0 -0.42 ; 0.08 -0.08]

# ╔═╡ 2ce2941b-047b-4c19-b757-6a3b3af83d21
# Compute the b-vector (reaction rate vector)
bV = W*mdot

# ╔═╡ 783ce2ee-398e-4d2e-9389-acfaa47b96df
sum(bV)

# ╔═╡ 9ad93173-2896-4aec-bb77-2aa3b51c8b50
md"""

#### Hmmm. Can we increase the flow rate to as fast as we want?
Unfortunately, no. The mass flow in the HFCS-42 process above will be governed by the rate at which xylose isomerase ($E_{1}$, units: kg) can process glucose. The fastest the enzyme can go will be given by:

$$r_{1} = k_{cat}E_{1}$$

where $k_{cat}$ is a _rate constant_, in this case called the turnover rate for xylose isomerase. These constants have been tabulated for many common enzymes in databases such as [BRENDA](https://www.brenda-enzymes.org). Looking at the various records, we find that $k_{cat}$ = 16.2$s^{-1}$.
"""

# ╔═╡ 4444d802-7d3c-4978-b9f3-5dd421c337b6
# How fast can the flow rate be for 10g of enzyme?
# convert to kg
E1 = 10*(1/1000) # convert to kg

# ╔═╡ 146cd344-c295-4d14-aa87-5f0ad4655640
# convert kcat to hr^-1 -
kcat = (16.2)*(3600) # convert to hr^-1

# ╔═╡ 9326629c-679a-4f8b-a06b-4eb39a752aa3
# what is the max rate we can expect for xylose isomerase -
max_rate = kcat*E1 # units: kg/hr

# ╔═╡ bb74b0b6-b209-425c-ad8f-e2cbc16d7670
md"""

#### Solve for flow rate given the reaction rate
If we knew the reaction rate, and the composition matrix were square, we could compute the flow rate by computing the inverse of the composition matrix:

$$\dot{\mathbf{m}} = \mathbf{W}^{-1}\mathbf{b}$$

where $\mathbf{W}^{-1}$ denotes the _inverse_ of the composition matrix $\mathbf{W}$ such that $\mathbf{W}\mathbf{W}^{-1}=\mathbf{W}^{-1}\mathbf{W} = \mathbf{I}$ (the identity matrix).

However, our compostion matrix is _not square_. We'll have to compute the [pseudo inverse](https://en.wikipedia.org/wiki/Moore–Penrose_inverse).
Let's put a pin in how we actually compute $\mathbf{W}^{-1}$ for the moment (and let Julia do it for us using the [pinv](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.pinv) command).

"""

# ╔═╡ 1a010c01-b515-4fd4-9f47-d7595c3a99f2
# compute W inverse -
W_inv = pinv(W)

# ╔═╡ d144a329-ea6c-4b69-80ca-aaf38d793491
W_inv*W

# ╔═╡ c0da6422-1318-49a3-8405-4a62f4323f55
# what is my bV -
bV_max = [max_rate ; -1*max_rate; 0.0]

# ╔═╡ eb90e0aa-df15-44b3-a11d-29742667fc9f
# compute the maximum flow rate -
mdot_max = W_inv*bV_max

# ╔═╡ a0523812-c4f7-461d-8ef0-1533bb86da3a
md"""
 #### Is there a diffent (less complicated) way to do this?
"""

# ╔═╡ 7920ee9f-cab4-4d20-b0a2-61eda1ada306
W_new = W[1:2,:]

# ╔═╡ 2dc5f057-b0d7-469c-a4be-b1ffb9190efe
W_inv_new = inv(W_new)

# ╔═╡ 4910fdf8-306e-4937-afe1-3a6bb099f4f0
W_inv_new*W_new

# ╔═╡ 0b6c30a4-53ea-4dbb-a145-81e37f511c2d
bV_max_new = bV_max[1:2]

# ╔═╡ e316ae95-549c-4fd5-80d6-2b653db355f6
m_dot_new = W_inv_new*bV_max_new

# ╔═╡ dc80bb02-0fb5-11ec-0866-0709291e7ff7
html"""<style>
main {
    max-width: 1200px;
  	width: 80%;
	margin: auto;
	font-family: "Roboto, monospace";
}

a {
	color: blue;
	text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
"""

# ╔═╡ Cell order:
# ╠═24c00466-e86f-410c-bb60-6dcb2399b63a
# ╟─ff56afbd-630f-42c3-9171-ccd808d5589f
# ╟─c7c294a2-43bf-4ae2-a6cc-f917eb6857a3
# ╟─58aff590-027b-4e63-8375-7ea304d1292d
# ╟─9ad11963-8b06-42dd-9a92-0310893c03fc
# ╠═f6dcdae8-e1a1-44da-b085-4859df008dec
# ╠═eb822e57-423a-4662-b402-052a616180d5
# ╠═2ce2941b-047b-4c19-b757-6a3b3af83d21
# ╠═783ce2ee-398e-4d2e-9389-acfaa47b96df
# ╟─9ad93173-2896-4aec-bb77-2aa3b51c8b50
# ╠═4444d802-7d3c-4978-b9f3-5dd421c337b6
# ╠═146cd344-c295-4d14-aa87-5f0ad4655640
# ╠═9326629c-679a-4f8b-a06b-4eb39a752aa3
# ╟─bb74b0b6-b209-425c-ad8f-e2cbc16d7670
# ╠═1a010c01-b515-4fd4-9f47-d7595c3a99f2
# ╠═d144a329-ea6c-4b69-80ca-aaf38d793491
# ╠═c0da6422-1318-49a3-8405-4a62f4323f55
# ╠═eb90e0aa-df15-44b3-a11d-29742667fc9f
# ╟─a0523812-c4f7-461d-8ef0-1533bb86da3a
# ╠═7920ee9f-cab4-4d20-b0a2-61eda1ada306
# ╠═2dc5f057-b0d7-469c-a4be-b1ffb9190efe
# ╠═4910fdf8-306e-4937-afe1-3a6bb099f4f0
# ╠═0b6c30a4-53ea-4dbb-a145-81e37f511c2d
# ╠═e316ae95-549c-4fd5-80d6-2b653db355f6
# ╟─dc80bb02-0fb5-11ec-0866-0709291e7ff7
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
__precompile__(true)

module ToeplitzMatrices

import Compat.view

import StatsBase
include("iterativeLinearSolvers.jl")

import Base: *, \, full, getindex, print_matrix, size, tril, triu, inv, A_mul_B!, Ac_mul_B,
    A_ldiv_B!, convert
import Base.LinAlg: BlasReal, DimensionMismatch

export Toeplitz, SymmetricToeplitz, Circulant, TriangularToeplitz, Hankel,
       chan, strang

# Abstract
abstract AbstractToeplitz{T<:Number} <: AbstractMatrix{T}

size(A::AbstractToeplitz) = (size(A, 1), size(A, 2))
getindex(A::AbstractToeplitz, i::Integer) = A[mod(i, size(A,1)), div(i, size(A,1)) + 1]

convert(::Type{Matrix}, S::AbstractToeplitz) = full(S)
convert{T}(::Type{AbstractMatrix{T}}, S::AbstractToeplitz) = convert(AbstractToeplitz{T}, S)
convert{T}(::Type{AbstractArray{T}}, S::AbstractToeplitz) = convert(AbstractToeplitz{T}, S)

# Convert an abstract Toeplitz matrix to a full matrix
function full{T}(A::AbstractToeplitz{T})
    m, n = size(A)
    Af = Array(T, m, n)
    for j = 1:n
        for i = 1:m
            Af[i,j] = A[i,j]
        end
    end
    return Af
end

# Fast application of a general Toeplitz matrix to a column vector via FFT
function A_mul_B!{T}(α::T, A::AbstractToeplitz{T}, x::StridedVector{T}, β::T,
      y::StridedVector{T})
    m = size(A,1)
    n = size(A,2)
    N = length(A.vcvr_dft)
    if m != length(y)
        throw(DimensionMismatch(""))
    end
    if n != length(x)
        throw(DimensionMismatch(""))
    end
    if N < 512
        y[:] *= β
        for j = 1:n
            tmp = α * x[j]
            for i = 1:m
                y[i] += tmp*A[i,j]
            end
        end
        return y
    end
    for i = 1:n
        A.tmp[i] = x[i]
    end
    for i = n+1:N
        A.tmp[i] = 0
    end
    A_mul_B!(A.tmp, A.dft, A.tmp)
    for i = 1:N
        A.tmp[i] *= A.vcvr_dft[i]
    end
    A.dft \ A.tmp
    for i = 1:m
        y[i] *= β
        y[i] += α * (T <: Real ? real(A.tmp[i]) : A.tmp[i])
    end
    return y
end

# Application of a general Toeplitz matrix to a general matrix
function A_mul_B!{T}(α::T, A::AbstractToeplitz{T}, B::StridedMatrix{T}, β::T,
    C::StridedMatrix{T})
    l = size(B, 2)
    if size(C, 2) != l
        throw(DimensionMismatch("input and output matrices must have same number of columns"))
    end
    for j = 1:l
        A_mul_B!(α, A, view(B, :, j), β, view(C, :, j))
    end
    return C
end

# Translate three to five argument A_mul_B!
A_mul_B!(y, A::AbstractToeplitz, x) = A_mul_B!(one(eltype(A)), A, x, zero(eltype(A)), x)

# * operator
function (*)(A::AbstractToeplitz, B::StridedVecOrMat)
    T = promote_type(eltype(A), eltype(B))
    A_mul_B!(one(T), convert(AbstractMatrix{T}, A), convert(AbstractArray{T}, B), zero(T),
        size(B,2) == 1 ? zeros(T, size(A, 1)) : zeros(T, size(A, 1), size(B, 2)))
end

# Left division of a general matrix B by a general Toeplitz matrix A, i.e. the solution x of Ax=B.
function A_ldiv_B!(A::AbstractToeplitz, B::StridedMatrix)
    if size(A, 1) != size(A, 2)
        error("Division: Rectangular case is not supported.")
    end
    for j = 1:size(B, 2)
        A_ldiv_B!(A, view(B, :, j))
    end
    return B
end

function (\)(A::AbstractToeplitz, b::AbstractVector)
    T = promote_type(eltype(A), eltype(b))
    if T != eltype(A)
        throw(ArgumentError("promotion of Toeplitz matrices not handled yet"))
    end
    bb = similar(b, T)
    copy!(bb, b)
    A_ldiv_B!(A, bb)
end

# General Toeplitz matrix
type Toeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
    vc::Vector{T}
    vr::Vector{T}
    vcvr_dft::Vector{S}
    tmp::Vector{S}
    dft::Base.DFT.Plan{S}
end

# Ctor
function Toeplitz(vc::Vector, vr::Vector)
    m, n = length(vc), length(vr)
    if vc[1] != vr[1]
        error("First element of the vectors must be the same")
    end

    T = promote_type(eltype(vc), eltype(vr), Float32)
    vcp, vrp = Vector{T}(vc), Vector{T}(vr)

    tmp = Array(promote_type(T, Complex{Float32}), m + n - 1)
    copy!(tmp, vcp)
    for i = 1:n - 1
        tmp[i + m] = vrp[n - i + 1]
    end
    dft = plan_fft!(tmp)
    return Toeplitz(vcp, vrp, dft*tmp, similar(tmp), dft)
end

convert{T}(::Type{AbstractToeplitz{T}},A::Toeplitz) = convert(Toeplitz{T},A)
convert{T}(::Type{Toeplitz{T}},A::Toeplitz) = Toeplitz(convert(Vector{T},A.vc),
                                                       convert(Vector{T},A.vr))

# Size of a general Toeplitz matrix
function size(A::Toeplitz, dim::Int)
    if dim == 1
        return length(A.vc)
    elseif dim == 2
        return length(A.vr)
    elseif dim > 2
        return 1
    else
        error("arraysize: dimension out of range")
    end
end

# Retrieve an entry
function getindex(A::Toeplitz, i::Integer, j::Integer)
    m = size(A,1)
    n = size(A,2)
    if i > m || j > n
        error("BoundsError()")
    end

    if i >= j
        return A.vc[i - j + 1]
    else
        return A.vr[1 - i + j]
    end
end

# Form a lower triangular Toeplitz matrix by annihilating all entries above the k-th diaganal
function tril(A::Toeplitz, k = 0)
    if k > 0
        error("Second argument cannot be positive")
    end
    Al = TriangularToeplitz(copy(A.vc), 'L', length(A.vr))
    if k < 0
      for i in -1:-1:k
          Al.ve[-i] = zero(eltype(A))
      end
    end
    return Al
end

# Form a lower triangular Toeplitz matrix by annihilating all entries below the k-th diaganal
function triu(A::Toeplitz, k = 0)
    if k < 0
        error("Second argument cannot be negative")
    end
    Al = TriangularToeplitz(copy(A.vr), 'U', length(A.vc))
    if k > 0
      for i in 1:k
          Al.ve[i] = zero(eltype(A))
      end
    end
    return Al
end

A_ldiv_B!(A::Toeplitz, b::StridedVector) =
    copy!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, strang(A), 1000, 100eps())[1])

# Symmetric
type SymmetricToeplitz{T<:BlasReal} <: AbstractToeplitz{T}
    vc::Vector{T}
    vcvr_dft::Vector{Complex{T}}
    tmp::Vector{Complex{T}}
    dft::Base.DFT.Plan
end
function SymmetricToeplitz{T<:BlasReal}(vc::Vector{T})
    tmp = convert(Array{Complex{T}}, [vc; zero(T); reverse(vc[2:end])])
    dft = plan_fft!(tmp)
    return SymmetricToeplitz(vc, dft*tmp, similar(tmp), dft)
end

convert{T}(::Type{AbstractToeplitz{T}},A::SymmetricToeplitz) = convert(SymmetricToeplitz{T},A)
convert{T}(::Type{SymmetricToeplitz{T}},A::SymmetricToeplitz) = SymmetricToeplitz(convert(Vector{T},A.vc))

function size(A::SymmetricToeplitz, dim::Int)
    if 1 <= dim <= 2
        return length(A.vc)
    else
        error("arraysize: dimension out of range")
    end
end

getindex(A::SymmetricToeplitz, i::Integer, j::Integer) = A.vc[abs(i - j) + 1]

A_ldiv_B!(A::SymmetricToeplitz, b::StridedVector) =
    copy!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1])

# Circulant
type Circulant{T<:Number,S<:Number} <: AbstractToeplitz{T}
    vc::Vector{T}
    vcvr_dft::Vector{S}
    tmp::Vector{S}
    dft::Base.DFT.Plan
end

function Circulant{T<:Real}(vc::Vector{T})
    TT = promote_type(eltype(vc), Float32)
    tmp = zeros(promote_type(TT, Complex{Float32}), length(vc))
    return Circulant(real(vc), fft(vc), tmp, plan_fft!(tmp))
end
function Circulant(vc::Vector)
    T = promote_type(eltype(vc), Float32)
    tmp = zeros(promote_type(T, Complex{Float32}), length(vc))
    return Circulant(vc, fft(vc), tmp, plan_fft!(tmp))
end

convert{T}(::Type{AbstractToeplitz{T}},A::Circulant) = convert(Circulant{T},A)
convert{T}(::Type{Circulant{T}},A::Circulant) = Circulant(convert(Vector{T},A.vc))

function size(C::Circulant, dim::Integer)
    if 1 <= dim <= 2
        return length(C.vc)
    else
        error("arraysize: dimension out of range")
    end
end

function getindex(C::Circulant, i::Integer, j::Integer)
    n = size(C, 1)
    if i > n || j > n
        error("BoundsError()")
    end
    return C.vc[mod(i - j, length(C.vc)) + 1]
end

function Ac_mul_B{T<:Real}(A::Circulant{T}, B::Circulant{T})
    tmp = similar(A.vcvr_dft)
    for i = 1:length(tmp)
        tmp[i] = conj(A.vcvr_dft[i]) * B.vcvr_dft[i]
    end
    return Circulant(real(A.vc)(A.dft \ tmp), tmp, A.tmp, A.dft)
end
function Ac_mul_B(A::Circulant, B::Circulant)
    T = promote_type(eltype(A), eltype(B))
    tmp = similar(A.vcvr_dft, T)
    for i = 1:length(tmp)
        tmp[i] = conj(A.vcvr_dft[i]) * B.vcvr_dft[i]
    end
    tmp2 = A.dft \ tmp
    return Circulant(Vector{T}(tmp2), tmp, Vector{T}(A.tmp), eltype(A) == T ? A.dft : plan_fft!(tmp2))
end

function A_ldiv_B!{T}(C::Circulant{T}, b::AbstractVector{T})
    n = length(b)
    size(C, 1) == n || throw(DimensionMismatch(""))
    for i = 1:n
        C.tmp[i] = b[i]
    end
    C.dft * C.tmp
    for i = 1:n
        C.tmp[i] /= C.vcvr_dft[i]
    end
    C.dft \ C.tmp
    for i = 1:n
        b[i] = (T <: Real ? real(C.tmp[i]) : C.tmp[i])
    end
    return b
end

function inv{T<:Real}(C::Circulant{T})
    vdft = 1 ./ C.vcvr_dft
    return Circulant(real(C.dft \ vdft), copy(vdft), similar(vdft), C.dft)
end
function inv(C::Circulant)
    vdft = 1 ./ C.vcvr_dft
    return Circulant(C.dft \ vdft, copy(vdft), similar(vdft), C.dft)
end

function strang{T}(A::AbstractMatrix{T})
    n = size(A, 1)
    v = Array(T, n)
    n2 = div(n, 2)
    for i = 1:n
        if i <= n2 + 1
            v[i] = A[i,1]
        else
            v[i] = A[1, n - i + 2]
        end
    end
    return Circulant(v)
end
function chan{T}(A::AbstractMatrix{T})
    n = size(A, 1)
    v = Array(T, n)
    for i = 1:n
        v[i] = ((n - i + 1) * A[i, 1] + (i - 1) * A[1, min(n - i + 2, n)]) / n
    end
    return Circulant(v)
end

# Triangular
type TriangularToeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
    ve::Vector{T}
    uplo::Char
    vcvr_dft::Vector{S}
    tmp::Vector{S}
    dft::Base.DFT.Plan
end

function TriangularToeplitz(ve::Vector, uplo::Symbol)
    n = length(ve)

    T = promote_type(eltype(ve), Float32)
    vep = Vector{T}(ve)

    tmp = zeros(promote_type(T, Complex{Float32}), 2n - 1)
    if uplo == :L
        copy!(tmp, vep)
    else
        tmp[1] = vep[1]
        for i = 1:n - 1
            tmp[n + i] = ve[n - i + 1]
        end
    end
    dft = plan_fft!(tmp)
    return TriangularToeplitz(vep, string(uplo)[1], dft * tmp, similar(tmp), dft)
end

function convert(::Type{Toeplitz}, A::TriangularToeplitz)
    if A.uplo == 'L'
        Toeplitz(A.ve, [A.ve[1]; zeros(length(A.ve) - 1)])
    else
        @assert A.uplo == 'U'
        Toeplitz([A.ve[1]; zeros(length(A.ve) - 1)], A.ve)
    end
end

convert{T}(::Type{AbstractToeplitz{T}},A::TriangularToeplitz) = convert(TriangularToeplitz{T},A)
convert{T}(::Type{TriangularToeplitz{T}},A::TriangularToeplitz) =
    TriangularToeplitz(convert(Vector{T},A.ve),A.uplo=='U'?(:U):(:L))

function size(A::TriangularToeplitz, dim::Int)
    if dim == 1 || dim == 2
        return length(A.ve)
    elseif dim > 2
        return 1
    else
        error("arraysize: dimension out of range")
    end
end

function getindex{T}(A::TriangularToeplitz{T}, i::Integer, j::Integer)
    if A.uplo == 'L'
        return i >= j ? A.ve[i - j + 1] : zero(T)
    else
        return i <= j ? A.ve[j - i + 1] : zero(T)
    end
end

function (*)(A::TriangularToeplitz, B::TriangularToeplitz)
    n = size(A, 1)
    if n != size(B, 1)
        throw(DimensionMismatch(""))
    end
    if A.uplo == B.uplo
        return TriangularToeplitz(conv(A.ve, B.ve)[1:n], A.uplo)
    end
    return Triangular(full(A), A.uplo) * Triangular(full(B), B.uplo)
end

Ac_mul_B(A::TriangularToeplitz, b::AbstractVector) =
    TriangularToeplitz(A.ve, A.uplo == 'U' ? :L : :U) * b

# NB! only valid for lower triangular
function smallinv{T}(A::TriangularToeplitz{T})
    n = size(A, 1)
    b = zeros(T, n)
    b[1] = 1 ./ A.ve[1]
    for k = 2:n
        tmp = zero(T)
        for i = 1:k-1
            tmp += A.uplo == 'L' ? A.ve[k - i + 1]*b[i] : A.ve[i + 1] * b[k - i]
        end
        b[k] = -tmp/A.ve[1]
    end
    return TriangularToeplitz(b, symbol(A.uplo))
end

function inv{T}(A::TriangularToeplitz{T})
    n = size(A, 1)
    if n <= 64
        return smallinv(A)
    end
    np2 = nextpow2(n)
    if n != np2
        return TriangularToeplitz(inv(TriangularToeplitz([A.ve, zeros(T, np2 - n)],
            symbol(A.uplo))).ve[1:n], symbol(A.uplo))
    end
    nd2 = div(n, 2)
    a1 = inv(TriangularToeplitz(A.ve[1:nd2], symbol(A.uplo))).ve
    return TriangularToeplitz([a1, -(TriangularToeplitz(a1, symbol(A.uplo)) *
        (Toeplitz(A.ve[nd2 + 1:end], A.ve[nd2 + 1:-1:2]) * a1))], symbol(A.uplo))
end

# A_ldiv_B!(A::TriangularToeplitz,b::StridedVector) = inv(A)*b
A_ldiv_B!(A::TriangularToeplitz, b::StridedVector) =
    copy!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, chan(A), 1000, 100eps())[1])

# extend levinson
StatsBase.levinson!(x::StridedVector, A::SymmetricToeplitz, b::StridedVector) =
    StatsBase.levinson!(A.vc, b, x)
function StatsBase.levinson!(C::StridedMatrix, A::SymmetricToeplitz, B::StridedMatrix)
    n = size(B, 2)
    if n != size(C, 2)
        throw(DimensionMismatch("input and output matrices must have same number of columns"))
    end
    for j = 1:n
        StatsBase.levinson!(view(C, :, j), A, view(B, :, j))
    end
    C
end
StatsBase.levinson(A::AbstractToeplitz, B::StridedVecOrMat) =
    StatsBase.levinson!(zeros(size(B)), A, copy(B))

# BlockTriangular
# type BlockTriangularToeplitz{T<:BlasReal} <: AbstractMatrix{T}
#     Mc::Array{T,3}
#     uplo::Char
#     Mc_dft::Array{Complex{T},3}
#     tmp::Vector{Complex{T}}
#     dft::Base.DFT.Plan
# end
# function BlockTriangularToeplitz{T<:BlasReal}(Mc::Array{T,3}, uplo::Symbol)
#     n, p, _ = size(Mc)
#     tmp = zeros(Complex{T}, 2n)
#     dft = plan_fft!(tmp)
#     Mc_dft = Array(Complex{T}, 2n, p, p)
#     for j = 1:p
#         for i = 1:p
#             Mc_dft[1,i,j] = complex(Mc[1,i,j])
#             for t = 2:n
#                 Mc_dft[t,i,j] = uplo == :L ? complex(Mc[t,i,j]) : zero(Complex{T})
#             end
#             Mc_dft[n+1,i,j] = zero(Complex{T})
#             for t = n+2:2n
#                 Mc_dft[t,i,j] = uplo == :L ? zero(Complex{T}) : complex(Mc[2n-t+2,i,j])
#             end
#             dft(view(Mc_dft, 1:2n, 1:p, 1:p))
#         end
#     end
#     return BlockTriangularToeplitz(Mc, string(uplo)[1], Mc_dft, tmp, dft, idft)
# end

#= Hankel Matrix
 A Hankel matrix is a matrix that is constant across the anti-diagonals:

  [a_0 a_1 a_2 a_3 a_4
   a_1 a_2 a_3 a_4 a_5
   a_2 a_3 a_4 a_5 a_6]

 This is precisely a Toeplitz matrix with the columns reversed:
                             [0 0 0 0 1
  [a_4 a_3 a_2 a_1 a_0        0 0 0 1 0
   a_5 a_4 a_3 a_2 a_1   *    0 0 1 0 0
   a_6 a_5 a_4 a_3 a_2]       0 1 0 0 0
                              1 0 0 0 0]
 We represent the Hankel matrix by wrapping the corresponding Toeplitz matrix.=#

# Hankel Matrix
type Hankel{T<:Number} <: AbstractMatrix{T}
    T::Toeplitz{T}
end

# Ctor: vc is the leftmost column and vr is the bottom row.
function Hankel(vc,vr)
    if vc[end] != vr[1]
        error("First element of rows must equal last element of columns")
    end
    n = length(vr)
    p = [vc; vr[2:end]]
    Hankel(Toeplitz(p[n:end],p[n:-1:1]))
end

convert(::Type{Array},A::Hankel) = convert(Matrix,A)
convert(::Type{Matrix},A::Hankel) = full(A)

convert{T}(::Type{AbstractMatrix{T}},A::Hankel) = convert(Hankel{T},A)
convert{T}(::Type{Hankel{T}},A::Hankel) = Hankel(convert(Toeplitz{T},A.T))

# Size
size(H::Hankel,k...) = size(H.T,k...)

# Full version of a Hankel matrix
function full{T}(A::Hankel{T})
    m, n = size(A)
    Af = Array(T, m, n)
    for j = 1:n
        for i = 1:m
            Af[i,j] = A[i,j]
        end
    end
    return Af
end

# Retrieve an entry by two indices
getindex(A::Hankel, i::Integer, j::Integer) = A.T[i,end-j+1]

# Retrieve an entry by one index
getindex(H::Hankel, i::Integer) = H[mod(i, size(H,1)), div(i, size(H,1)) + 1]

# Fast application of a general Hankel matrix to a general vector
*(A::Hankel,b::AbstractVector) = A.T * reverse(b)

# Fast application of a general Hankel matrix to a general matrix
*(A::Hankel,B::AbstractMatrix) = A.T * flipdim(B, 1)
## BigFloat support

(*){T<:BigFloat}(A::Toeplitz{T}, b::Vector) = irfft(
    rfft([
        A.vc;
        reverse(A.vr[2:end])]
    ) .* rfft([
        b;
        zeros(length(b) - 1)
    ]),
    2 * length(b) - 1
)[1:length(b)]

end #module


text:
using Vega
using Compat
using Patchwork

Patchwork.load_js_runtime()

stock = ["chairs", "tables", "desks", "rugs", "lamps"];

quantity = [15, 10, 10, 5, 20];

piechart(x = stock, y = quantity)

using PyPlot
precipitation = [0,0,0,0,0,0,0,0,0,0,0.12,0.01,0,0,0,0.37,0,0,0,0,0.01,0,0,0,0.01,0.01,0,0.17,0.01,0.11,0.31]
date = collect(1:31)
fig = figure(1, figsize=(4, 4))
plot(date, precipitation, ".")
title("Boston Precipitation")
xlabel("May 2013")
ylabel("Precipitation")

addprocs(1)
srand(111)
r = remotecall(rand, 2, 3, 4)
s = @spawnat 2 1 .+ fetch(r)
fetch(s)

max_stars_repo_path: Section 4/Vega.jl

max_issues_repo_path: Section 4/Vega.jl

max_forks_repo_path: Section 4/Vega.jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

_all_in_ui(X::AbstractArray) = _all_lteq_impl(0, X, 1)

"""
    y = fromui(x::Real, lo::Real, hi::Real)
    y = fromui(x::Real, lo_hi::ClosedInterval{<:Real})

*BAT-internal, not part of stable public API.*

Linear bijective transformation from the unit inverval (i.e. `x ∈ 0..1`) to
`y ∈ lo..hi`.

Use `inv(fromui)` to get the the inverse transformation.

Use `@inbounds` to disable range checking on the input value.
"""
function fromui end

@inline unsafe_fromui(x::Real, lo::Real, hi::Real) = muladd(x, (hi - lo), lo)

@inline function fromui(x::Real, lo::Real, hi::Real)
    @boundscheck x in 0..1 || throw(ArgumentError("Input value not in 0..1"))
    unsafe_fromui(x, lo, hi)
end

Base.@propagate_inbounds fromui(x::Real, lo_hi::ClosedInterval{<:Real}) =
    fromui(x, minimum(lo_hi), maximum(lo_hi))

@inline unsafe_inv_fromui(x::Real, lo::Real, hi::Real) = (x - lo) / (hi - lo)

@inline function inv_fromui(x::Real, lo::Real, hi::Real)
    @boundscheck x in lo..hi || throw(ArgumentError("Input value not in lo..hi"))
    unsafe_inv_fromui(x, lo, hi)
end

Base.@propagate_inbounds inv_fromui(x::Real, lo_hi::ClosedInterval{<:Real}) =
    inv_fromui(x, minimum(lo_hi), maximum(lo_hi))

@inline Base.inv(::typeof(fromui)) = inv_fromui
@inline Base.inv(::typeof(inv_fromui)) = fromui


text:
__precompile__()

module Tools
using LinearAlgebra
using Distributed

#################### utilities
function discount_sum(x::Vector{DT}, gamma::DT, terminal::DT) where DT<:AbstractFloat
   y = similar(x)
   runsum = terminal
   for t=length(x):-1:1
      runsum = x[t] + gamma*runsum
      y[t] = runsum
   end
   return y
end

function discount_sum!(y::Matrix{DT}, x::Matrix{DT},
                       t::Int, gamma::DT, terminal::DT) where DT<:AbstractFloat
   runsum = terminal
   for i=size(x, 1):-1:1
      runsum = x[i,t] + gamma*runsum
      y[i,t] = runsum
   end
end

function discount_sum!(y::Matrix{T}, x::Matrix{DT},
                       t::Int, gamma::DT, terminal::DT) where {T<:AbstractFloat,
                                                             DT<:AbstractFloat}
   runsum = terminal
   for i=size(x, 1):-1:1
      runsum = x[i,t] + gamma*runsum
      y[i,t] = runsum
   end
end

# reward: T x K
function compute_returns!(returns::Vector{DT},
                          rewards::Matrix{Float64},
                          gamma::Float64) where DT<:AbstractFloat
   T, K = size(rewards)
   mat_ret = reshape(returns, T, K)
   for t=1:K
      #mat_ret[:,t] = discount_sum(rewards[:,t], gamma, 0.0)
      discount_sum!(mat_ret, rewards, t, gamma, 0.0)
   end
   #return reshape(returns, T*K)
end

function compute_advantages(returns::Vector{Float64},
                            baseline::Vector{Float64},
                            obs::Matrix{Float64},
                            gamma::Float64)
   return returns - baseline
end

function computeGAEadvantages!(adv::Vector{DT},
                               baseline::Vector,
                               rewards::Matrix{Float64},
                               obs::Matrix{Float64},
                               gae::Float64,
                               gamma::Float64) where DT<:AbstractFloat
   T, K = size(rewards)
   tdsum = Array{DT}(undef, size(rewards))
   mat_adv = reshape(adv, T, K)

   nthread = min(K, Threads.nthreads()) # cant have more threads than data
   Threads.@threads for tid=1:nthread # WHY IS MT SO MUCH SLOWER??
      #nthread = 1
      #for tid=1:nthread # WHY IS MT SO MUCH SLOWER??
      thread_range = Distributed.splitrange(K, nthread)[tid]

      for t=thread_range
         base = view(baseline, (t-1)*T+1:t*T)
         for k=1:T-1
            tdsum[k,t] = rewards[k,t] + gamma * base[k+1] - base[k]
         end
         tdsum[T,t] = rewards[T,t] - base[T]
         #mat_adv[:,t] = discount_sum(tdsum[:,t], DT(gae*gamma), DT(0.0))
         discount_sum!(mat_adv, tdsum, t, DT(gae*gamma), DT(0.0))
      end
   end
end

# solver algorithms
function cpu_cg_solve(fim::Matrix{T}, # ((n+1)*m + m) X T*K
                      vpg::Vector{T}, # ((n+1)*m + m)
                      cg_iters::Integer=10,
                      reg::Float64=1e-4,
                      tol::Float64=1e-10) where T<:AbstractFloat
   # Initialize cg variables
   r = copy(vpg)
   p = copy(vpg)
   x = zeros(T, size(vpg))  # I want x to be same shape as vpg but full of zeros
   rdr = dot(vpg, vpg)  # outputs a scalar
   z = fim*p
   #z += p*reg

   iters = 1
   for i=1:cg_iters
      v = rdr/dot(p, z)      # scalar

      x .+= v.*p
      r .-= v.*z

      rdr_new = dot(r, r)    # scalar
      ratio = rdr_new/rdr
      rdr = rdr_new

      iters = i
      if rdr < tol
         break   # this should break the for loop
      end

      p .*= ratio
      p .+= r

      #println(norm(z))
      mul!(z, fim, p)
   end
   println("    Used $iters cojugate-gradient iterations.");
   return x
end

function hvp_cg_solve(hvpfim,
                      vpg::Vector{DT}, # ((n+1)*m + m)
                      cg_iters::Integer=10,
                      reg::Float64=1e-4,
                      tol::Float64=1e-10) where DT<:AbstractFloat
   # Initialize cg variables
   r = copy(vpg)
   p = copy(vpg)
   x = zeros(DT, size(vpg))
   rdr = dot(vpg, vpg)  # outputs a scalar
   z = zeros(DT, size(vpg))
   hvpfim(z, p)

   iters = 1
   for i=1:cg_iters
      v = rdr/dot(p, z)      # scalar

      x += v*p
      r -= v*z
      #@. x += v*p
      #@. r -= v*z

      rdr_new = dot(r, r)    # scalar
      ratio = rdr_new/rdr
      rdr = rdr_new

      iters = i
      if rdr < tol
         break   # this should break the for loop
      end

      @. p = r + ratio*p
      #p[:] .*= ratio
      #p[:] += r

      #println(norm(z))
      hvpfim(z, p)
   end
   println("    Used $iters HVP-cg iterations.");
   return x
end

function par_cg_solver(fim::Matrix{T}, # ((n+1)*m + m) X T*K
                       vpg::Vector{DT}, # ((n+1)*m + m)
                       toworkers, #::Vector{RemoteChannel},
                       tomaster, #::RemoteChannel,
                       cg_iters::Integer=10,
                       reg::DT=1e-4) where {T<:AbstractFloat,
                                            DT<:AbstractFloat}
   put!(tomaster[myid()], vpg) # send vpg
   vpg[:] = take!(toworkers[myid()])

   z = fim*vpg
   put!(tomaster[myid()], z)
   p = similar(z)

   for i=1:cg_iters
      z[:] = take!( toworkers[myid()] )
      p[:] = take!( toworkers[myid()] )

      if isnan(z[1]) || i==cg_iters break end

      mul!(z, fim, p)
      put!(tomaster[myid()], z)
   end
   return p # p is NPG
end
function scatter(tomaster, val::Vector{T}) where T<:AbstractFloat
   @sync for w in workers()
      @async put!(tomaster[w], copy(val))
   end
end
function gather!(toworkers, val::Vector{T}) where T<:AbstractFloat
   val[:] .= T(0.0)
   for w in workers()
      val[:] += take!(toworkers[w])
   end
   val[:] /= nworkers()
end
function par_cg_manager(nparam::Integer,
                        toworkers,
                        tomaster,
                        cg_iters::Integer=10,
                        reg::T=1e-4,
                        tol::Float64=1e-10) where T<:AbstractFloat
   vpg = Array{T}(undef, nparam)
   z = Array{T}(undef, nparam)
   x = zeros(T, nparam)

   gather!(tomaster, vpg)
   scatter(toworkers, vpg)
   rdr = dot(vpg, vpg)  # outputs a scalar
   p = copy(vpg)
   r = copy(vpg)

   count = 0
   for i=1:cg_iters
      count += 1
      gather!(tomaster, z) # get Z vector
      v = rdr/dot(p, z)      # scalar

      x += v*p
      r -= v*z

      rdr_new = dot(r, r)    # scalar
      ratio = rdr_new/rdr
      rdr = rdr_new

      if rdr < tol || i == cg_iters
         z[1:5] .= NaN
         scatter(toworkers, z)
         scatter(toworkers, x) # send NPG
         break
      end

      p[:] .*= ratio
      p[:] += r

      scatter(toworkers, z)
      scatter(toworkers, p)
   end
   println("\tconverged in $count iterations.")
   return vpg, x # x == npg
end

# GPU solver
#=
function gpu_cg_solve{T<:AbstractFloat,
                      DT<:AbstractFloat}(fim::Knet.KnetArray{T,2}, # ((n+1)*m + m) X T*K
                                         vpg::Vector{DT}, # ((n+1)*m + m)
                                         cg_iters::Integer=10,
                                         reg::Float64=1e-4,
                                         tol::Float64=1e-10)
   # gradll (fat): num_params x T*K
   # vpg : vector num_params

   # Initialize cg variables
   if T == DT
      vpg32 = reshape(vpg, length(vpg), 1)
   else
      vpg32 = Array{T,2}(reshape(vpg, length(vpg), 1))
   end
   r = KnetArray{T}(vpg32)
   p = KnetArray{T}(vpg32)
   x = KnetArray{T}(zeros(size(vpg32)))
   rdr = dot(vpg, vpg) #(p' * p)[1] #dot(p, p)  # outputs a scalar
   z = fim * p

   @inbounds for i=1:cg_iters
      #v = rdr/dot(p, z)      # scalar
      v = rdr/((p'*z)[1])      # scalar

      x += v*p
      r -= v*z

      #rdr_new = dot(r, r)    # scalar
      rdr_new = (r' * r)[1] #dot(r, r)    # scalar
      ratio = rdr_new/rdr
      rdr = rdr_new

      if rdr < tol
         break   # this should break the for loop
      end

      p = r + ratio*p

      z[:] = fim*p
   end
   return x
end
=#

end


text:
# Import Libraries
using Turing
using Turing: Variational
using Distributions
using AbstractGPs, KernelFunctions
using PyPlot
using StatsFuns
using JSON3
using Flux
using ProgressBars
import Random
import LinearAlgebra

# Define a kernel.
sekernel(alpha, rho) = alpha^2 * transform(SEKernel(), invsqrt2/rho)

function compute_f(kernel, X, eta, beta=0, jitter=0)
  K = kernelmatrix(kernel, X, obsdim=1) + LinearAlgebra.I * jitter
  return LinearAlgebra.cholesky(K).L * eta .+ beta
end

@model function GPClassify(y, X, jitter=1e-6)
    # Priors.
    alpha ~ LogNormal(0, 1)
    rho ~ LogNormal(0, 1)
    beta ~ Normal(0, 1)  # intercept.
    eta ~ filldist(Normal(0, 1), length(y))

    # Latent GP
    kernel = sekernel(alpha, rho)
    f = compute_f(sekernel(alpha, rho), X, eta, beta, jitter)

    # Sampling Distribution.
    y ~ arraydist(Bernoulli.(logistic.(f)))
end;

# To extract parameters from trained variational distribution
# (surrogate posterior).
function make_surrogate_sampler(m, q; nsamples=300)
    qsamples = rand(q, nsamples)
    _, sym2range = Variational.bijector(m; sym_to_ranges=Val(true))
    return sym -> qsamples[collect(sym2range[sym][1]), :]
end;

### Read data ###
# See data in Notebook.
data_path = joinpath(@__DIR__, "../data/gp-classify-data-N50.json")

# Load data in JSON format.
data = let
    x = open(f -> read(f, String), data_path)
    JSON3.read(x)
end

### Name the data (X, y) ###
X = [data["x1"] data["x2"]]
y = Int64.(data["y"])

### Create model.###
m =  GPClassify(Float64.(y), X);
kernel_params = [:alpha, :rho, :beta];

### Fit via ADVI ###
Random.seed!(7)
@time q = vi(m, ADVI(1, 1000))  # num_elbo_samples, num_iterations.

# Get posterior samples
surrogate_sampler = make_surrogate_sampler(m, q, nsamples=500)
advi_samples = Dict{Symbol, Any}(sym => vec(surrogate_sampler(sym))
                                 for sym in kernel_params)
advi_samples[:eta] = extract_gp(:eta);

### Fit via HMC ###
Random.seed!(0)
burn = 500
nsamples = 500
@time hmc_chain = sample(m, HMC(0.05, 20), burn + nsamples);

# Get posterior samples
hmc_samples = Dict{Symbol, Any}([
    sym => vec(group(hmc_chain, sym).value.data)[end-nsamples+1:end]
for sym in kernel_params])
hmc_samples[:eta] = Matrix(group(hmc_chain, :eta).
                           value.data[end-nsamples+1:end, :, 1]')

### Fit via NUTS ###
Random.seed!(0)
burn = 500
nsamples = 500
@time nuts_chain = sample(m, NUTS(burn, 0.8, max_depth=10), burn + nsamples)

# Get posterior samples
nuts_samples = Dict{Symbol, Any}([
    sym => vec(group(nuts_chain, sym).value.data)
for sym in kernel_params])
nuts_samples[:eta] = Matrix(group(nuts_chain, :eta).value.data[:, :, 1]')

max_stars_repo_name: luiarthur/TuringBnpBenchmarks

max_issues_repo_name: luiarthur/TuringBnpBenchmarks

max_forks_repo_name: luiarthur/TuringBnpBenchmarks


text:
using SampledVectors
using Test

@testset "SampledVectors.jl" begin
    #basic
    l = SampledVector{Int}(23)
    for i in 1:1234
        @test sampledindexes(l)|>collect == sampled(l)
        push!(l, i)
    end
    @test sampled(l) == [1, 65, 129, 193, 257, 321, 385, 449, 513, 577, 641, 705, 769, 833, 897, 961, 1025, 1089, 1153, 1217, 1234]

    #empty!
    @test empty!(l) === l
    @test length(l) == 0
    @test collect(l) |> isempty
    @test step(l) == 1
    @test sampledindexes(l)|>collect |> isempty
    @test sampled(l)|>collect |> isempty

    # edge case: capacity=2
    l = SampledVector{Int}(BigInt(2))
    for i in 1:27
        push!(l, i)
    end
    @test sampled(l) == [1, 27]
    @test step(l) == 32

    # vector init
    l = SampledVector([1,2,3])
    for i in 4:7
        push!(l, i)
    end
    @test sampled(l) == [1, 5, 7]
    @test capacity(l) == 3

    #setcapacity!
    setcapacity!(l, 5) #larger
    @test capacity(l) == 5
    foreach(x->push!(l, x), 8:16)
    l2 = SampledVector{Int}(5)
    foreach(x->push!(l2, x), 1:16)
    @test sampled(l) == sampled(l2)
    @test sampledindexes(l)|>collect == sampled(l)
    @test step(l) == step(l2)
    setcapacity!(l, 4) #smaller
    @test capacity(l) == 4
    @test sampled(l) == [1, 9, 16]
    @test sampledindexes(l)|>collect == sampled(l)
    push!(l, 17)
    push!(l, 18) #has remains; sampled(l) has even length
    @test sampled(l) == [1, 9, 17, 18]
    setcapacity!(l, 3)
    @test capacity(l) == 3
    @test sampled(l) == [1, 17, 18]
    @test sampledindexes(l)|>collect == sampled(l)
    setcapacity!(l, 2)
    @test capacity(l) == 2
    @test sampled(l) == [1, 18]
    @test sampledindexes(l)|>collect == sampled(l)

    #factor=3
    l = SampledVector{Int}(15, factor=3)
    for i in 1:48
        @test sampledindexes(l)|>collect == sampled(l)
        push!(l, i)
    end
    @test sampled(l) == [1, 10, 19, 28, 37, 46, 48]
    @test factor(l) == 3
    @test step(l) == 9

    #set factor
    setfactor!(l, 4)
    for i in 49:200
        @test sampledindexes(l)|>collect == sampled(l)
        push!(l, i)
    end
    @test sampled(l) == [1, 37, 73, 109, 145, 181, 200]
    @test factor(l) == 4
    @test step(l) == 4*9
end


text:
"""
    TranslationUnit(idx, source, args)
    TranslationUnit(idx, source, args, unsavedFiles, options)
Parse the given source file and the translation unit corresponding to that file.
"""
mutable struct TranslationUnit
    ptr::CXTranslationUnit
    idx::Index
    function TranslationUnit(
        idx::Index,
        source_filename,
        command_line_args,
        num_command_line_args,
        unsaved_files,
        num_unsaved_files,
        options,
    )
        ptr = clang_parseTranslationUnit(
            idx,
            source_filename,
            command_line_args,
            num_command_line_args,
            unsaved_files,
            num_unsaved_files,
            options,
        )
        @assert ptr != C_NULL "failed to parse file: $source_filename"
        obj = new(ptr, idx)
        finalizer(obj) do x
            if x.ptr != C_NULL
                clang_disposeTranslationUnit(x)
                x.ptr = C_NULL
            end
        end
        return obj
    end
end
function TranslationUnit(idx, source, args, unsavedFiles, options)
    return TranslationUnit(
        idx, source, args, length(args), unsavedFiles, length(unsavedFiles), options
    )
end
function TranslationUnit(idx, source, args, options)
    return TranslationUnit(idx, source, args, length(args), C_NULL, 0, options)
end
function TranslationUnit(idx, source, args)
    return TranslationUnit(
        idx, source, args, length(args), C_NULL, 0, CXTranslationUnit_None
    )
end

Base.unsafe_convert(::Type{CXTranslationUnit}, x::TranslationUnit) = x.ptr

"""
    spelling(tu::TranslationUnit) -> String
Return the original translation unit source file name.
"""
function spelling(tu::TranslationUnit)
    cxstr = clang_getTranslationUnitSpelling(tu)
    ptr = clang_getCString(cxstr)
    s = unsafe_string(ptr)
    clang_disposeString(cxstr)
    return s
end

"""
    getTranslationUnitCursor(tu::TranslationUnit) -> CLCursor
    getTranslationUnitCursor(tu::CXTranslationUnit) -> CXCursor
Return the cursor that represents the given translation unit.
"""
getTranslationUnitCursor(tu::CXTranslationUnit) = clang_getTranslationUnitCursor(tu)
getTranslationUnitCursor(tu::TranslationUnit)::CLCursor = clang_getTranslationUnitCursor(tu)

## TODO:
# clang_createTranslationUnit
# clang_createTranslationUnit2
# clang_defaultEditingTranslationUnitOptions
# clang_parseTranslationUnit2
# clang_parseTranslationUnit2FullArgv
# clang_defaultSaveOptions
# clang_saveTranslationUnit
# clang_suspendTranslationUnit
# clang_defaultReparseOptions
# clang_reparseTranslationUnit
# clang_getTUResourceUsageName
# clang_getCXTUResourceUsage
# clang_disposeCXTUResourceUsage
# clang_getTranslationUnitTargetInfo
# clang_TargetInfo_dispose
# clang_TargetInfo_getTriple
# clang_TargetInfo_getPointerWidth

# helper
"""
    parse_header(
        index::Index,
        header::AbstractString,
        args::Vector{String}=[],
        flags=CXTranslationUnit_None,
    ) -> TranslationUnit
Return the TranslationUnit for a given header. This is the main entry point for parsing.

# Arguments
- `header::AbstractString`: the header file to parse.
- `index::Index`: CXIndex pointer (pass to avoid re-allocation).
- `args::Vector{String}`: compiler switches as string array, eg: ["-x", "c++", "-fno-elide-type"].
- `flags`: bitwise OR of CXTranslationUnit_Flags.

See also [`parse_headers`](@ref).
"""
function parse_header(
    index::Index,
    header::AbstractString,
    args::Vector{String}=[],
    flags=CXTranslationUnit_None,
)
    # TODO: support parsing in-memory with CXUnsavedFile arg
    #       to _parseTranslationUnit.jj
    #unsaved_file = CXUnsavedFile()
    !isfile(header) && throw(ArgumentError("$header not found."))

    return TranslationUnit(index, header, args, flags)
end

"""
    parse_headers(
        index::Index,
        headers::Vector{String},
        args::Vector{String}=[],
        flags=CXTranslationUnit_None,
    ) -> Vector{TranslationUnit}
Return a [`TranslationUnit`](@ref) vector for the given headers.

See also [`parse_header`](@ref).
"""
function parse_headers(
    index::Index,
    headers::Vector{String},
    args::Vector{String}=[],
    flags=CXTranslationUnit_None,
)
    map(unique(headers)) do header
        parse_header(index, header, args, flags)
    end
end


text:
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("RvSpectML")   end
using Pkg
 Pkg.activate(".")

verbose = true
 if verbose   println("# Loading RvSpecML")    end
 using RvSpectML
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics, Dates

all_spectra = include(joinpath(pkgdir(EchelleInstruments),"examples/read_neid_solar_data_20190918.jl"))

order_list_timeseries = extract_orders(all_spectra,pipeline_plan)

linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks","G2.espresso.mas")
line_list_df = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan,  v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = 21e3)

mask_scale_factors = [ 1.0,  2, 4, 6, 8, 10 ] # , 12, 14, 16 ]
 rms_rvs = zeros(length(mask_scale_factors))
 rms_binned_rvs = zeros(length(mask_scale_factors))
 println("Starting Tophat CCFs (w/o var)...")
 for (i,mask_scale_factor) in enumerate(mask_scale_factors)
   println("# mask_scale_factor = ", mask_scale_factor)
   local (ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_scale_factor=mask_scale_factor, mask_type=:tophat,ccf_mid_velocity=ccf_mid_velocity, recalc=true)
   alg_fit_rv = EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=2, init_guess_ccf_σ=ccf_mid_velocity)
   local rvs_ccf = calc_rvs_from_ccf_total(ccfs, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false, alg_fit_rv=alg_fit_rv)
   rms_rvs[i] = std(rvs_ccf.-mean(rvs_ccf))
   rms_binned_rvs[i] = std(RvSpectML.bin_rvs_consecutive(rvs_ccf.-mean(rvs_ccf), 4))
 end

if make_plot(pipeline_plan,:ccf_total)
   using Plots
   plt1 = plot(mask_scale_factors, rms_rvs, color=1, label="Tophat w/o var")
   plt2 = plot(mask_scale_factors, rms_binned_rvs, color=1, label=:none)
   xlabel!("Mask width scale parameter")
   plot(plt1,plt2,layout=(2,1) )
   ylabel!("RMS (m/s)")
end

mask_scale_factors2 = [ 1.0,  2, 4, 6, 8, 10, 12 ]
 rms_rvs2 = zeros(length(mask_scale_factors2))
 rms_binned_rvs2 = zeros(length(mask_scale_factors2))
 println("Starting Tophat CCFs (w/ vars)...")
 for (i,mask_scale_factor) in enumerate(mask_scale_factors2)
   println("# mask_scale_factor = ", mask_scale_factor)
   ((ccfs_expr, ccf_vars_expr), v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:tophat, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true, recalc=true)
   rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, ccf_vars_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
   rms_rvs2[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
   rms_binned_rvs2[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
 end
 if make_plot(pipeline_plan,:ccf_total)
   scatter!(plt1,mask_scale_factors2, rms_rvs2, color=2, label="Tophat w/ var")
   scatter!(plt2,mask_scale_factors2, rms_binned_rvs2, color=2, label=:none)
   plot(plt1,plt2,layout=(2,1) )
   xlabel!("Mask width scale parameter")
   ylabel!("RMS (m/s)")
 end

#=
((ccfs, ccf_vars), v_grid ) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:tophat, mask_scale_factor=4.0,  ccf_mid_velocity=ccf_mid_velocity, calc_ccf_var=true, recalc=true)
(ccfs_norm, ccf_vars_norm) = EchelleCCFs.calc_normalized_ccfs(ccfs, ccf_vars)
ccf_template = EchelleCCFs.calc_ccf_template(ccfs_norm, ccf_vars_norm)
ccfs_sample_covar = 1.0*EchelleCCFs.calc_ccf_sample_covar(ccfs_norm, ccf_vars_norm, assume_normalized=true )
mean_covar_on_diag = EchelleCCFs.RVFromCCF.compute_mean_off_diag(ccfs_sample_covar, 0)
mrvt = EchelleCCFs.MeasureRvFromCCFTemplate(v_grid=v_grid, template=ccf_template, measure_width_at_frac_depth=0.5, frac_of_width_to_fit=3.0, mean_var=mean_covar_on_diag)
(rvs_t, σ_rvs_t) = EchelleCCFs.measure_rvs_from_ccf(v_grid,ccfs_norm, ccf_vars_norm, alg=mrvt)
rms_rv_binned = bin_rvs_consecutive(rvs_t.-mean(rvs_t),4)
println("# RMS of RVs: ", std(rvs_t), "  binned RVs: ", std(rms_rv_binned), "  median(rvs_t): ", median(σ_rvs_t) )

function triangular_shape_near_diagonal(n::Integer, width::Real; norm::Real = 1)
  @. triangle_helper(x, p) = p[1]*max(1-abs(x)/p[2], 0)
  shape_model = collect( triangle_helper(i-j,[norm,width]) for i in 1:n, j in 1:n )
end
shape_model = triangular_shape_near_diagonal(length(v_grid), 8)

mrvt3 = EchelleCCFs.MeasureRvFromCCFTemplateNonDiagCovar(v_grid=v_grid, template=ccf_template, measure_width_at_frac_depth=0.5, frac_of_width_to_fit=2, mean_var=mean_covar_on_diag, near_diag_covar=shape_model )
(rvs_t3, σ_rvs_t3) = EchelleCCFs.measure_rvs_from_ccf(v_grid,ccfs_norm, ccf_vars_norm, ccfs_sample_covar, alg=mrvt3)
rms_rv_binned = bin_rvs_consecutive(rvs_t3.-mean(rvs_t3),4)
println("# RMS of RVs: ", std(rvs_t3), "  binned RVs: ", std(rms_rv_binned), "  median(σ_rvs_t): ", median(σ_rvs_t3) )
=#

mask_scale_factors3 = [ 1.0,  2, 4, 6, 8, 10, 12 ]
 rms_rvs3 = zeros(length(mask_scale_factors3))
 rms_binned_rvs3 = zeros(length(mask_scale_factors3))
  println("Starting Half Cos...")
  for (i,mask_scale_factor) in enumerate(mask_scale_factors3)
    println("# mask_scale_factor = ", mask_scale_factor)
    (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:halfcos, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true)
    rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
    rms_rvs3[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
    rms_binned_rvs3[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
  end
  if make_plot(pipeline_plan,:ccf_total)
    plot!(plt1,mask_scale_factors3, rms_rvs3, color=3, label="Half Cos")
    scatter!(plt2,mask_scale_factors3, rms_binned_rvs3, color=3, label=:none)
    plot(plt1,plt2,layout=(2,1))
  end

mask_scale_factors4 = [ 1.0,  2.0, 3.0, 4.0, 4.5, 5.0, 5.5, 6.0 ]
 rms_rvs4 = zeros(length(mask_scale_factors4))
   rms_binned_rvs4 = zeros(length(mask_scale_factors4))
   rms_rvs4_quad = zeros(length(mask_scale_factors4))
   rms_rvs4_cent = zeros(length(mask_scale_factors4))
   println("Starting Gaussian...")
    for (i,mask_scale_factor) in enumerate(mask_scale_factors4)
      println("# mask_scale_factor = ", mask_scale_factor)
      (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:gaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true)
      rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
      rms_rvs4[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
      rms_binned_rvs4[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
      rvs_ccf_quad = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false, alg_fit_rv= EchelleCCFs.MeasureRvFromCCFQuadratic())
      #rvs_ccf_cent = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false, alg_fit_rv= EchelleCCFs.MeasureRvFromCCFCentroid())
      rms_rvs4_quad[i] = std(rvs_ccf_quad.-mean(rvs_ccf_quad))
      #rms_rvs4_cent[i] = std(rvs_ccf_cent.-mean(rvs_ccf_cent))
    end

  if make_plot(pipeline_plan,:ccf_total)
    scatter!(plt1,mask_scale_factors4, rms_rvs4, color=4, label="Gaussian- Gaussian")
    scatter!(plt2,mask_scale_factors4, rms_binned_rvs4, color=4, label=:none)
    plot!(plt1,mask_scale_factors4, rms_rvs4_quad, color=5, label="Gaussian- Quadratic")
    #scatter!(plt1,mask_scale_factors4, rms_rvs4_cent, color=6, label="Gaussian- Centroid")
    plot(plt1,plt2,layout=(2,1) )
  end

#=
# Commented out since this is the EXPRES LSF model and slow, so people only run it for NEID when they really want to
mask_scale_factors5 = [ 4.0, 6, 8 ]
 rms_rvs5 = zeros(length(mask_scale_factors5))
   rms_binned_rvs5 = zeros(length(mask_scale_factors5))
   println("Starting Super Gaussian...")
    for (i,mask_scale_factor) in enumerate(mask_scale_factors5)
      println("# mask_scale_factor = ", mask_scale_factor)
      (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:supergaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true, use_old = false)
      rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
      rms_rvs5[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
      rms_binned_rvs5[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
    end
    if make_plot(pipeline_plan,:ccf_total)
      plot!(plt1,mask_scale_factors5, rms_rvs5, color=5, label="Super-Gaussian")
      plot!(plt2,mask_scale_factors5, rms_binned_rvs5, color=5, label=:none)
      plot(plt1,plt2,layout=(2,1) )
    end
=#

if make_plot(pipeline_plan,:ccf_total)
 plt1 = plot(mask_scale_factors, rms_rvs, color=1, label="Tophat w/o vars")
 plt2 = plot(mask_scale_factors, rms_binned_rvs, color=1, label=:none)
 plot(plt1,plt2,layout=(2,1) )
 xlabel!("Mask width scale parameter")
 ylabel!("RMS (m/s)")
 scatter!(plt1,mask_scale_factors2, rms_rvs2, color=2, label="Tophat w/ vars")
 scatter!(plt2,mask_scale_factors2, rms_binned_rvs2, color=2, label=:none)
 plot!(plt1,mask_scale_factors3, rms_rvs3, color=3, label="Half Cos")
 plot!(plt2,mask_scale_factors3, rms_binned_rvs3, color=3, label=:none)
 plot(plt1,plt2,layout=(2,1) )
 scatter!(plt1,mask_scale_factors4, rms_rvs4, color=4, label="Gaussian- Gaussian")
 scatter!(plt2,mask_scale_factors4, rms_binned_rvs4, color=4, label=:none)
 scatter!(plt1,mask_scale_factors4, rms_rvs4_quad, color=5, label="Gaussian- Quadratic")
 #scatter!(plt1,mask_scale_factors4, rms_rvs4_cent, color=6, label="Gaussian- Centroid")
 #plot!(plt1,mask_scale_factors5, rms_rvs5, color=5, label="Super-Gaussian")
 #plot!(plt2,mask_scale_factors5, rms_binned_rvs5, color=5, label=:none)

 ylims!(plt1,0.8,1.0)
 ylims!(plt2,0,0.5)
 plot(plt1,plt2,layout=(2,1) )
 if save_plot(pipeline_plan,:ccf_total)
   savefig("neid_mask_width_comp.png")
 end
end


text:
@testset "Search space" begin
    @testset "in()" begin
        for i in 1:NumTestRepetitions
            reps = rand(1:10)
            ss1 = symmetric_search_space(reps, (0.0, 1.0))
            ind = rand_individual(ss1)
            for j in 1:reps
                @test (mins(ss1)[j] <= ind[j] <= maxs(ss1)[j])
            end
        end
    end

    @testset "Symmetric search space with default range" begin
        ss1 = symmetric_search_space(1)
        @test numdims(ss1) == 1
        @test ranges(ss1) == [(0.0, 1.0)]
        @test range_for_dim(ss1,1) == (0.0, 1.0)

        for i in 1:NumTestRepetitions
            ind = rand_individual(ss1)
            @test size(ind) == (1,)
            @test in(ind, ss1)
        end

        ss3 = symmetric_search_space(3)
        @test numdims(ss3) == 3
        @test ranges(ss3) == [(0.0, 1.0), (0.0, 1.0), (0.0, 1.0)]
        @test range_for_dim(ss3,1) == (0.0, 1.0)
        @test range_for_dim(ss3,2) == (0.0, 1.0)
        @test range_for_dim(ss3,3) == (0.0, 1.0)

        for i in 1:NumTestRepetitions
            ind = rand_individual(ss3)
            @test size(ind) == (3,)
            @test in(ind, ss3)
        end
    end

    @testset "SymmetricSearchSpace with given range" begin
        ss1 = symmetric_search_space(1, (-1.0, 1.0))
        @test numdims(ss1) == 1
        @test ranges(ss1) == [(-1.0, 1.0)]
        @test range_for_dim(ss1,1) == (-1.0, 1.0)

        for i in 1:NumTestRepetitions
            reps = rand(1:100)
            range = (rand(), rand())
            ss = symmetric_search_space(reps, range)
            @test numdims(ss) == reps
            @test all([(dr == range) for dr in ranges(ss)])
        end
    end

    @testset "rand_individual is within the search space" begin
        for i in 1:NumTestRepetitions
            reps = rand(1:100)
            mm = sort(rand(2,1), dims=1)
            range = (mm[1], mm[2])
            ss = symmetric_search_space(reps, range)
            ind = rand_individual(ss)
            @test length(ind) == numdims(ss)
            @test in(ind, ss)
        end
    end

    @testset "rand_individuals creates many individuals and all are within the search space" begin
        for i in 1:NumTestRepetitions
            reps = rand(1:10)
            mm = sort(rand(2,1), dims=1)
            range = (mm[1], mm[2])
            ss = symmetric_search_space(reps, range)
            numinds = rand(1:10)
            inds = rand_individuals(ss, numinds)
            @test size(inds,1) == numdims(ss)
            @test size(inds,2) == numinds
            for j in 1:numinds
                @test in(inds[:,j], ss)
            end
        end
    end

    @testset "rand_individuals correctly handles column-wise generation in assymetric search spaces" begin
        for _ in 1:NumTestRepetitions÷10
            numdimensions = rand(1:13)
            minbounds = rand(numdimensions)
            ds = rand(1:10, numdimensions) .* rand(numdimensions)
            maxbounds = minbounds .+ ds
            parambounds = collect(zip(minbounds, maxbounds))
            ss = RangePerDimSearchSpace(parambounds)
            @test mins(ss) == minbounds
            @test maxs(ss) == maxbounds
            @test round.(deltas(ss), digits=6) == round.(ds, digits=6)

            # Now generate 100 individuals and make sure they are all within bounds
            inds = rand_individuals(ss, 100)
            @test size(inds, 2) == 100
            @inbounds for i in 1:size(inds, 2)
                for d in 1:numdimensions
                    @test (minbounds[d] <= inds[d,i] <= maxbounds[d])
                end
            end
        end
    end

    @testset "RangePerDimSearchSpace" begin
        ss = RangePerDimSearchSpace([(0.0, 1.0)])
        @test mins(ss) == [0.0]
        @test maxs(ss) == [1.0]
        @test deltas(ss) == [1.0]

        ss = RangePerDimSearchSpace([(0.0, 1.0), (0.5, 10.0)])
        @test mins(ss) == [0.0, 0.5]
        @test maxs(ss) == [1.0, 10.0]
        @test deltas(ss) == [1.0, 9.5]
    end

    @testset "rand_individuals_lhs samples in LHS intervals" begin
        ss = RangePerDimSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])

        inds = rand_individuals_lhs(ss, 2)
        @test size(inds, 1) == 3
        @test size(inds, 2) == 2

        sorted = sort(inds, dims=2) # Sort per row --> in their ordered intervals
        @test (0.0 <= sorted[1,1] <= 0.5)
        @test (0.5 <= sorted[1,2] <= 1.0)

        @test (2.0 <= sorted[2,1] <= 2.5)
        @test (2.5 <= sorted[2,2] <= 3.0)

        @test (4.0 <= sorted[3,1] <= 4.5)
        @test (4.5 <= sorted[3,2] <= 5.0)
    end

    @testset "feasible finds feasible points in the search space" begin
        ss = RangePerDimSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])

        # We use the double transpose below to ensure the actual and expected
        # values have the same type (matrices, not vectors).
        @test BlackBoxOptim.feasible([1.1, 2.0, 4.0], ss) == [1.0, 2.0, 4.0]
        @test BlackBoxOptim.feasible([1.1, 3.0, 4.0], ss) == [1.0, 3.0, 4.0]
        @test BlackBoxOptim.feasible([1.1, 2.0, 5.0], ss) == [1.0, 2.0, 5.0]
        @test BlackBoxOptim.feasible([1.1, 3.0, 5.0], ss) == [1.0, 3.0, 5.0]

        @test BlackBoxOptim.feasible([-0.1, 2.0, 4.0], ss) == [0.0, 2.0, 4.0]
        @test BlackBoxOptim.feasible([-0.1, 3.0, 4.0], ss) == [0.0, 3.0, 4.0]
        @test BlackBoxOptim.feasible([-0.1, 2.0, 5.0], ss) == [0.0, 2.0, 5.0]
        @test BlackBoxOptim.feasible([-0.1, 3.0, 5.0], ss) == [0.0, 3.0, 5.0]

        @test BlackBoxOptim.feasible([0.0, 1.9, 4.0], ss) == [0.0, 2.0, 4.0]
        @test BlackBoxOptim.feasible([0.0, 1.9, 4.0], ss) == [0.0, 2.0, 4.0]
        @test BlackBoxOptim.feasible([1.0, 1.9, 5.0], ss) == [1.0, 2.0, 5.0]
        @test BlackBoxOptim.feasible([1.0, 1.9, 5.0], ss) == [1.0, 2.0, 5.0]

        @test BlackBoxOptim.feasible([0.0, 3.3, 4.0], ss) == [0.0, 3.0, 4.0]
        @test BlackBoxOptim.feasible([0.0, 3.2, 4.0], ss) == [0.0, 3.0, 4.0]
        @test BlackBoxOptim.feasible([1.0, 3.1, 5.0], ss) == [1.0, 3.0, 5.0]
        @test BlackBoxOptim.feasible([1.0, 3.9, 5.0], ss) == [1.0, 3.0, 5.0]

        @test BlackBoxOptim.feasible([-0.4, 3.3, 14.5], ss) == [0.0, 3.0, 5.0]
    end

    @testset "diameters" begin
        ss = RangePerDimSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])
        diams = diameters(ss)

        @test length(diams) == 3
        @test diams[:] == [1.0, 1.0, 1.0]
    end

    @testset "concat(ss1, ss2)" begin
        ss1 = RangePerDimSearchSpace([(0.0, 1.0), (2.0, 3.0), (4.0, 5.0)])
        ss2 = RangePerDimSearchSpace([(6.0, 7.0), (8.0, 9.0)])

        sscat = vcat(ss1, ss2)
        @test numdims(sscat) == 5
        @test mins(sscat) == [0.0, 2.0, 4.0, 6.0, 8.0]
        @test maxs(sscat) == [1.0, 3.0, 5.0, 7.0, 9.0]
        @test deltas(sscat) == fill(1.0, 5)
    end
end


text:
function dotproduct!(y::Array{Float64,1}, a::Array{Float64,1}, b::Array{Float64,1})
    n = length(y)
    @inbounds for i = 1 : n
       y[i] = a[i] * b[i]
    end
    return nothing
end

function dotproduct1!(a::Array{Float64,1}, b::Array{Float64,1})
    n = length(a)
    @inbounds for i = 1 : n
        a[i] = a[i] * b[i]
    end
    return nothing
end

function mysum1!(y::Array{Float64,1}, a::Array{Float64,1}, b::Array{Float64,1})
    n = length(y)
    @inbounds for i = 1 : n
       y[i] = a[i] + b[i]
    end
    return nothing
end

function mysum2!(a::Array{Float64,1}, b::Array{Float64,1})
    n = length(a)
    @inbounds for i = 1 : n
       a[i] = a[i] + b[i]
    end
    return nothing
end

function OneStepForward!(snapShot2::SnapShot, snapShot1::SnapShot, fidMtx::FidMtx, tmp::Array{Float64,1}, tmp1::Array{Float64,1})
    snapShot2.it  = snapShot1.it + 1
    dotproduct!(snapShot2.Vxx, fidMtx.MVxxBVxx, snapShot1.Vxx); mysum1!(tmp, snapShot1.Txxx, snapShot1.Txxz);
    AmulB!(tmp1, fidMtx.MVxxBTxx, tmp);                         mysum2!(snapShot2.Vxx, tmp1);

    mysum1!(tmp, snapShot1.Txzx, snapShot1.Txzz);
    dotproduct!(snapShot2.Vxz, fidMtx.MVxzBVxz, snapShot1.Vxz); AmulB!(tmp1, fidMtx.MVxzBTxz, tmp); mysum2!(snapShot2.Vxz, tmp1);
    dotproduct!(snapShot2.Vzx, fidMtx.MVzxBVzx, snapShot1.Vzx); AmulB!(tmp1, fidMtx.MVzxBTxz, tmp); mysum2!(snapShot2.Vzx, tmp1);

    dotproduct!(snapShot2.Vzz, fidMtx.MVzzBVzz, snapShot1.Vzz); mysum1!(tmp, snapShot1.Tzzx, snapShot1.Tzzz);
    AmulB!(tmp1, fidMtx.MVzzBTzz, tmp);                         mysum2!(snapShot2.Vzz, tmp1);

    mysum1!(tmp, snapShot2.Vxx, snapShot2.Vxz);
    dotproduct!(snapShot2.Txxx, fidMtx.MTxxxBTxxx, snapShot1.Txxx); AmulB!(tmp1, fidMtx.MTxxxBVx, tmp); mysum2!(snapShot2.Txxx, tmp1);
    dotproduct!(snapShot2.Tzzx, fidMtx.MTzzxBTzzx, snapShot1.Tzzx); AmulB!(tmp1, fidMtx.MTzzxBVx, tmp); mysum2!(snapShot2.Tzzx, tmp1);
    dotproduct!(snapShot2.Txzz, fidMtx.MTxzzBTxzz, snapShot1.Txzz); AmulB!(tmp1, fidMtx.MTxzzBVx, tmp); mysum2!(snapShot2.Txzz, tmp1);

    mysum1!(tmp, snapShot2.Vzx, snapShot2.Vzz);
    dotproduct!(snapShot2.Txxz, fidMtx.MTxxzBTxxz, snapShot1.Txxz); AmulB!(tmp1, fidMtx.MTxxzBVz, tmp); mysum2!(snapShot2.Txxz, tmp1);
    dotproduct!(snapShot2.Tzzz, fidMtx.MTzzzBTzzz, snapShot1.Tzzz); AmulB!(tmp1, fidMtx.MTzzzBVz, tmp); mysum2!(snapShot2.Tzzz, tmp1);
    dotproduct!(snapShot2.Txzx, fidMtx.MTxzxBTxzx, snapShot1.Txzx); AmulB!(tmp1, fidMtx.MTxzxBVz, tmp); mysum2!(snapShot2.Txzx, tmp1);
    return nothing
end

function OneStepForward!(snapShot2::SnapShot, snapShot1::SnapShot, fidMtx::FidMtx)
    snapShot2.it  = snapShot1.it + 1
    snapShot2.Vxx[:] = fidMtx.MVxxBVxx .* snapShot1.Vxx + fidMtx.MVxxBTxx*(snapShot1.Txxx+snapShot1.Txxz)
    snapShot2.Vxz[:] = fidMtx.MVxzBVxz .* snapShot1.Vxz + fidMtx.MVxzBTxz*(snapShot1.Txzx+snapShot1.Txzz)
    snapShot2.Vzx[:] = fidMtx.MVzxBVzx .* snapShot1.Vzx + fidMtx.MVzxBTxz*(snapShot1.Txzx+snapShot1.Txzz)
    snapShot2.Vzz[:] = fidMtx.MVzzBVzz .* snapShot1.Vzz + fidMtx.MVzzBTzz*(snapShot1.Tzzx+snapShot1.Tzzz)

    snapShot2.Txxx[:] = fidMtx.MTxxxBTxxx .* snapShot1.Txxx + fidMtx.MTxxxBVx*(snapShot2.Vxx+snapShot2.Vxz)
    snapShot2.Txxz[:] = fidMtx.MTxxzBTxxz .* snapShot1.Txxz + fidMtx.MTxxzBVz*(snapShot2.Vzx+snapShot2.Vzz)
    snapShot2.Tzzx[:] = fidMtx.MTzzxBTzzx .* snapShot1.Tzzx + fidMtx.MTzzxBVx*(snapShot2.Vxx+snapShot2.Vxz)
    snapShot2.Tzzz[:] = fidMtx.MTzzzBTzzz .* snapShot1.Tzzz + fidMtx.MTzzzBVz*(snapShot2.Vzx+snapShot2.Vzz)
    snapShot2.Txzx[:] = fidMtx.MTxzxBTxzx .* snapShot1.Txzx + fidMtx.MTxzxBVz*(snapShot2.Vzx+snapShot2.Vzz)
    snapShot2.Txzz[:] = fidMtx.MTxzzBTxzz .* snapShot1.Txzz + fidMtx.MTxzzBVx*(snapShot2.Vxx+snapShot2.Vxz)
    return nothing
end

function OneStepAdjoint!(SnapShot2::SnapShot, SnapShot1::SnapShot, fidMtx::FidMtx)
    SnapShot2.it = SnapShot1.it - 1
    tmp = (fidMtx.MTxxxBVx)'*SnapShot1.Txxx + (fidMtx.MTzzxBVx)'*SnapShot1.Tzzx + (fidMtx.MTxzzBVx)'*SnapShot1.Txzz
    SnapShot2.Vxx = SnapShot1.Vxx + tmp
    SnapShot2.Vxz = SnapShot1.Vxz + tmp
    tmp = (fidMtx.MTxxzBVz)'*SnapShot1.Txxz + (fidMtx.MTzzzBVz)'*SnapShot1.Tzzz + (fidMtx.MTxzxBVz)'*SnapShot1.Txzx
    SnapShot2.Vzx = SnapShot1.Vzx + tmp
    SnapShot2.Vzz = SnapShot1.Vzz + tmp
    SnapShot2.Txxx = fidMtx.MTxxxBTxxx .* SnapShot1.Txxx
    SnapShot2.Txxz = fidMtx.MTxxzBTxxz .* SnapShot1.Txxz
    SnapShot2.Tzzx = fidMtx.MTzzxBTzzx .* SnapShot1.Tzzx
    SnapShot2.Tzzz = fidMtx.MTzzzBTzzz .* SnapShot1.Tzzz
    SnapShot2.Txzx = fidMtx.MTxzxBTxzx .* SnapShot1.Txzx
    SnapShot2.Txzz = fidMtx.MTxzzBTxzz .* SnapShot1.Txzz

    tmp = (fidMtx.MVxxBTxx)'*SnapShot2.Vxx
    SnapShot2.Txxx = tmp + SnapShot2.Txxx
    SnapShot2.Txxz = tmp + SnapShot2.Txxz
    tmp = (fidMtx.MVzzBTzz)'*SnapShot2.Vzz
    SnapShot2.Tzzx = tmp + SnapShot2.Tzzx
    SnapShot2.Tzzz = tmp + SnapShot2.Tzzz
    tmp = (fidMtx.MVxzBTxz)'*SnapShot2.Vxz + (fidMtx.MVzxBTxz)'*SnapShot2.Vzx
    SnapShot2.Txzx = tmp + SnapShot2.Txzx
    SnapShot2.Txzz = tmp + SnapShot2.Txzz
    SnapShot2.Vxx = fidMtx.MVxxBVxx .* SnapShot2.Vxx
    SnapShot2.Vxz = fidMtx.MVxzBVxz .* SnapShot2.Vxz
    SnapShot2.Vzx = fidMtx.MVzxBVzx .* SnapShot2.Vzx
    SnapShot2.Vzz = fidMtx.MVzzBVzz .* SnapShot2.Vzz
    return nothing
end

function OneStepAdjoint!(SnapShot2::SnapShot, SnapShot1::SnapShot, fidMtx::FidMtx, tmp::Array{Float64,1}, tmp1::Array{Float64,1})
    SnapShot2.it = SnapShot1.it - 1
    AcmulB!(tmp , fidMtx.MTxxxBVx, SnapShot1.Txxx); AcmulB!(tmp1, fidMtx.MTzzxBVx, SnapShot1.Tzzx); mysum2!(tmp, tmp1);
    AcmulB!(tmp1, fidMtx.MTxzzBVx, SnapShot1.Txzz); mysum2!(tmp , tmp1); mysum1!(SnapShot2.Vxx, SnapShot1.Vxx, tmp); mysum1!(SnapShot2.Vxz, SnapShot1.Vxz, tmp);
    AcmulB!(tmp , fidMtx.MTxxzBVz, SnapShot1.Txxz); AcmulB!(tmp1, fidMtx.MTzzzBVz, SnapShot1.Tzzz); mysum2!(tmp, tmp1);
    AcmulB!(tmp1, fidMtx.MTxzxBVz, SnapShot1.Txzx); mysum2!(tmp , tmp1); mysum1!(SnapShot2.Vzx, SnapShot1.Vzx, tmp); mysum1!(SnapShot2.Vzz, SnapShot1.Vzz, tmp);
    dotproduct!(SnapShot2.Txxx, fidMtx.MTxxxBTxxx, SnapShot1.Txxx)
    dotproduct!(SnapShot2.Txxz, fidMtx.MTxxzBTxxz, SnapShot1.Txxz)
    dotproduct!(SnapShot2.Tzzx, fidMtx.MTzzxBTzzx, SnapShot1.Tzzx)
    dotproduct!(SnapShot2.Tzzz, fidMtx.MTzzzBTzzz, SnapShot1.Tzzz)
    dotproduct!(SnapShot2.Txzx, fidMtx.MTxzxBTxzx, SnapShot1.Txzx)
    dotproduct!(SnapShot2.Txzz, fidMtx.MTxzzBTxzz, SnapShot1.Txzz)

    AcmulB!(tmp, fidMtx.MVxxBTxx, SnapShot2.Vxx); mysum2!(SnapShot2.Txxx, tmp); mysum2!(SnapShot2.Txxz, tmp);
    AcmulB!(tmp, fidMtx.MVzzBTzz, SnapShot2.Vzz); mysum2!(SnapShot2.Tzzx, tmp); mysum2!(SnapShot2.Tzzz, tmp);
    AcmulB!(tmp, fidMtx.MVxzBTxz, SnapShot2.Vxz);
    AcmulB!(tmp1,fidMtx.MVzxBTxz, SnapShot2.Vzx); mysum2!(tmp, tmp1);
    mysum2!(SnapShot2.Txzx, tmp); mysum2!(SnapShot2.Txzz, tmp);
    dotproduct1!(SnapShot2.Vxx, fidMtx.MVxxBVxx);
    dotproduct1!(SnapShot2.Vxz, fidMtx.MVxzBVxz);
    dotproduct1!(SnapShot2.Vzx, fidMtx.MVzxBVzx);
    dotproduct1!(SnapShot2.Vzz, fidMtx.MVzzBVzz);
    return nothing
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
export RegularMesh, getRegularMesh
export getCellCenteredGrid, getNodalGrid, getEdgeGrids, getFaceGrids
export getCellCenteredAxes, getNodalAxes
export getVolume, getVolumeInv, getFaceArea, getFaceAreaInv, getLength, getLengthInv

"""
	mutable struct jInv.Mesh.RegularMesh <: AbstractTensorMesh

	Regular mesh in 1D, 2D, and 3D

	Fields:

		domain::Vector{Float64}  - physical domain [min(x1) max(x1) min(x2) max(x2)]
		h::Vector{Float64}       - cell size
		x0::Vector{Float64}      - origin
		dim::Int                 - dimension of mesh
		n::Vector{Int64}         - number of cells in each dimension
		nc::Int                  - total number of cells
		nf::Vector{Int64}        - number of faces in each dimension
		ne::Vector{Int64}        - number of edges in each dimension

		Persistent Operators:

		Operators should not be accessed directly. They will be built, if needed,
		when accessing them using specified method. clear!(M) will release all
		memory.

			Div::SparseMatrixCSC    - divergence (faces -> cell-centers)
									  Access via: getDivergenceMatrix(M)
			Grad::SparseMatrixCSC   - gradient (nodal -> edges)
									  Access via: getNodalGradientMatrix(M)
			Curl::SparseMatrixCSC   - curl (edges -> faces)
									  Access via: getCurlMatrix(M)
			Af::SparseMatrixCSC     - face average (faces -> cell-centers)
									  Access via: getFaceAverageMatrix(M)
			Ae::SparseMatrixCSC     - edge average (edges -> cell-centers)
									  Access via: getEdgeAverageMatrix(M)
			An::SparseMatrixCSC     - nodal average (nodes -> cell-centers)
									  Access via: getNodalAverageMatrix(M)
			V::SparseMatrixCSC      - cell volumes (diagonal matrix)
									  Access via: getVolume(M)
			F::SparseMatrixCSC      - face area (diagonal matrix)
									  Access via: getFaceArea(M)
			L::SparseMatrixCSC      - edge length (diagonal matrix)
									  Access via: getLength(M)
			Vi::SparseMatrixCSC     - inverse cell volumes (diagonal matrix)
									  Access via: getVolumeInv(M)
			Fi::SparseMatrixCSC     - inverse face area (diagonal matrix)
									  Access via: getFaceAreaInv(M)
			Li::SparseMatrixCSC     - inverse edge length (diagonal matrix)
									  Access via: getLengthAreaInv(M)
			nLap::SparseMatrixCSC   - nodal Laplacian
									  Access via: getNodalLaplacian(M)

	Examples:
	M2D  = getRegularMesh([1.2 2.4 2.2 5.0],[3,4])
	M3D  = getRegularMesh([1.2 2.4 2.2 5.0 0 1],[3,4,7])

"""
mutable struct RegularMesh <: AbstractTensorMesh
	domain::Vector{Float64}
	h::Vector{Float64}
	x0::Vector{Float64}
	dim::Int
	n::Vector{Int64}
	nc::Int
	nf::Vector{Int64}
	ne::Vector{Int64}
	Div::SparseMatrixCSC
	Grad::SparseMatrixCSC
	Curl::SparseMatrixCSC
	Af::SparseMatrixCSC
	Ae::SparseMatrixCSC
	An::SparseMatrixCSC
	V::SparseMatrixCSC
	F::SparseMatrixCSC
	L::SparseMatrixCSC
	Vi::SparseMatrixCSC
	Fi::SparseMatrixCSC
	Li::SparseMatrixCSC
	nLap::SparseMatrixCSC
end

"""
	function jInv.Mesh.getRegularMesh

	Constructs regular mesh

	Input:
		domain - physical domain rectangular
		n      - number of cells in each dimension

	Examples:
	M2D  = getRegularMesh([1.2 2.4 2.2 5.0],[3,4])
	M3D  = getRegularMesh([1.2 2.4 2.2 5.0 0 1],[3,4,7])

"""
function getRegularMesh(domain,n)
	domain = vec(float(domain))
	n     = vec(n)
	nc    = prod(n)
	h     = vec((domain[2:2:end]-domain[1:2:end])./n)
	dim   = round(Int64,length(domain)/2)
	if dim==1
		nf = [(n[1]+1)]
		ne = [(n[1]+1)]
	elseif dim==2
		nf = [(n[1]+1)*n[2]; n[1]*(n[2]+1);]
		ne = [n[1]*(n[2]+1); (n[1]+1)*n[2];  ]
	elseif dim==3
		nf = [(n[1]+1)*n[2]*n[3]; n[1]*(n[2]+1)*n[3]; n[1]*n[2]*(n[3]+1); ]
		ne = [n[1]*(n[2]+1)*(n[3]+1); (n[1]+1)*n[2]*(n[3]+1); (n[1]+1)*(n[2]+1)*n[3]; ]
	end
	x0 = domain[1:2:end]
	empt = spzeros(0,0);
return RegularMesh(domain,h,x0,dim,n,nc,nf,ne,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt,empt)
end

import Base.==
function ==(M1::RegularMesh,M2::RegularMesh)
	isEqual = fill(true,20)

	# check mandatory fields
	isEqual[1] =  (M1.domain == M2.domain)
	isEqual[2] =  (M1.h     == M2.h)
	isEqual[3] =  (M1.x0    == M2.x0)
	isEqual[4] =  (M1.dim   == M2.dim)
	isEqual[5] =  (M1.n     == M2.n)
	isEqual[6] =  (M1.nc    == M2.nc)
	isEqual[7] =  (M1.nf    == M2.nf)
	isEqual[8] =  (M1.ne    == M2.ne)

	# check fields that might be empty
	if !(isempty(M1.Div)) && !(isempty(M2.Div))
		isEqual[9] = (M1.Div == M2.Div)
	end
	if !(isempty(M1.Grad)) && !(isempty(M2.Grad))
		isEqual[10] = (M1.Grad == M2.Grad)
	end
	if !(isempty(M1.Curl)) && !(isempty(M2.Curl))
		isEqual[11] = (M1.Curl == M2.Curl)
	end
	if !(isempty(M1.Af)) && !(isempty(M2.Af))
		isEqual[12] = (M1.Af == M2.Af)
	end
	if !(isempty(M1.Ae)) && !(isempty(M2.Ae))
		isEqual[13] = (M1.Ae == M2.Ae)
	end
	if !(isempty(M1.An)) && !(isempty(M2.An))
		isEqual[14] = (M1.An == M2.An)
	end
	if !(isempty(M1.V)) && !(isempty(M2.V))
		isEqual[15] = (M1.V == M2.V)
	end
	if !(isempty(M1.F)) && !(isempty(M2.F))
		isEqual[16] = (M1.F == M2.F)
	end
	if !(isempty(M1.L)) && !(isempty(M2.L))
		isEqual[17] = (M1.L == M2.L)
	end
	if !(isempty(M1.Vi)) && !(isempty(M2.Vi))
		isEqual[18] = (M1.Vi == M2.Vi)
	end
	if !(isempty(M1.Fi)) && !(isempty(M2.Fi))
		isEqual[19] = (M1.Fi == M2.Fi)
	end
	if !(isempty(M1.Li)) && !(isempty(M2.Li))
		isEqual[20] = (M1.Li == M2.Li)
	end
	return all(isEqual)
end

# --- grid constructor
function getCellCenteredGrid(Mesh::RegularMesh)
# X = getCellCenteredGrid(Mesh::RegularMesh)
	return getCellCenteredGrid(Mesh.domain,Mesh.n)
end

function getNodalGrid(Mesh::RegularMesh)
# X = getNodalGrid(Mesh::RegularMesh)
	return getNodalGrid(Mesh.domain,Mesh.n)
end

function getEdgeGrids(Mesh::RegularMesh)
	return getEdgeGrids(Mesh.domain,Mesh.n)
end

function getFaceGrids(Mesh::RegularMesh)
	return getFaceGrids(Mesh.domain,Mesh.n)
end

function getCellCenteredGrid(domain,n)
# X = getCellCenteredGrid(domain,n)
	dim = round(Int64,length(domain)/2)
	if dim==1
		xc = getCellCenteredAxes(domain,n)
	elseif dim==2
		x1,x2 = getCellCenteredAxes(domain,n)
		X1,X2 = ndgrid(x1,x2)
		xc = [vec(X1) vec(X2)]
	elseif dim==3
		x1,x2,x3 = getCellCenteredAxes(domain,n)
		X1,X2,X3 = ndgrid(x1,x2,x3)
		xc = [vec(X1) vec(X2) vec(X3)]
	end
	return xc
end

function getNodalGrid(domain,n)
# X = getNodalGrid(domain,nc)
	dim = round(Int64,length(domain)/2)
	if dim==1
		xc = getNodalAxes(domain,n)
	elseif dim==2
		x1,x2 = getNodalAxes(domain,n)
		X1,X2 = ndgrid(x1,x2)
		xc = [vec(X1) vec(X2)]
	elseif dim==3
		x1,x2,x3 = getNodalAxes(domain,n)
		X1,X2,X3 = ndgrid(x1,x2,x3)
		xc = [vec(X1) vec(X2) vec(X3)]
	end
	return xc
end

function getEdgeGrids(domain,nc)
# X = getEdgeGrids(domain,nc)
	dim = round(Int64,length(domain)/2)
	h   = (domain[2:2:end]-domain[1:2:end])./nc
	if dim==2
		x1n,x2n = getNodalAxes(domain,nc)
		x1c,x2c = getCellCenteredAxes(domain,nc)
		# edge-1 grid
		X1t,X2t = ndgrid(x1c,x2n)
		x1 = [vec(X1t) vec(X2t)]
		# edge-2 grid
		X1t,X2t = ndgrid(x1n,x2c)
		x2 = [vec(X1t) vec(X2t)]
		return (x1,x2)
	elseif dim==3
		x1n,x2n,x3n = getNodalAxes(domain,nc)
		x1c,x2c,x3c = getCellCenteredAxes(domain,nc)
		# edge-1 grid
		X1,X2,X3 = ndgrid(x1c,x2n,x3n)
		x1 = [vec(X1) vec(X2) vec(X3)]
		# edge-2 grid
		X1,X2,X3 = ndgrid(x1n,x2c,x3n)
		x2 = [vec(X1) vec(X2) vec(X3)]
		# edge-3 grid
		X1,X2,X3 = ndgrid(x1n,x2n,x3c)
		x3 = [vec(X1) vec(X2) vec(X3)]
		return (x1,x2,x3)
	end
end

function getFaceGrids(domain,nc)
# X = getFaceGrids(domain,nc)
	dim = round(Int64,length(domain)/2)
	h   = (domain[2:2:end]-domain[1:2:end])./nc
	if dim==2
		x1n,x2n = getNodalAxes(domain,nc)
		x1c,x2c = getCellCenteredAxes(domain,nc)
		# face-1 grid
		X1t,X2t = ndgrid(x1n,x2c)
		x1 = [vec(X1t) vec(X2t)]
		# face-2 grid
		X1t,X2t = ndgrid(x1c,x2n)
		x2 = [vec(X1t) vec(X2t)]
		return (x1,x2)
	elseif dim==3
		x1n,x2n,x3n = getNodalAxes(domain,nc)
		x1c,x2c,x3c = getCellCenteredAxes(domain,nc)
		# face-1 grid
		X1,X2,X3 = ndgrid(x1n,x2c,x3c)
		x1 = [vec(X1) vec(X2) vec(X3)]
		# face-2 grid
		X1,X2,X3 = ndgrid(x1c,x2n,x3c)
		x2 = [vec(X1) vec(X2) vec(X3)]
		# face-3 grid
		X1,X2,X3 = ndgrid(x1c,x2c,x3n)
		x3 = [vec(X1) vec(X2) vec(X3)]
		return (x1,x2,x3)
	end
end

function getNodalAxes(Mesh::RegularMesh)
	return getNodalAxes(Mesh.domain,Mesh.n)
end

function getNodalAxes(domain,nc)
	dim = round(Int64,length(domain)/2)

	if dim==1
		x1 = collect(range(domain[1], step=(domain[2]-domain[1])/nc,length=nc+1))
		return x1
	elseif dim==2
		x1 = collect(range(domain[1], step=(domain[2]-domain[1])/nc[1],length=nc[1]+1))
		x2 = collect(range(domain[3], step=(domain[4]-domain[3])/nc[2],length=nc[2]+1))
		return (x1,x2)
	elseif dim==3
		x1 = collect(range(domain[1], step=(domain[2]-domain[1])/nc[1],length=nc[1]+1))
		x2 = collect(range(domain[3], step=(domain[4]-domain[3])/nc[2],length=nc[2]+1))
		x3 = collect(range(domain[5], step=(domain[6]-domain[5])/nc[3],length=nc[3]+1))
		return (x1,x2,x3)
	end
end

function getCellCenteredAxes(Mesh::RegularMesh)
	return getCellCenteredAxes(Mesh.domain,Mesh.n)
end

function getCellCenteredAxes(domain,nc)
	dim = round(Int64,length(domain)/2)
	h   = vec(domain[2:2:end]-domain[1:2:end])./vec(nc)
	if dim==1
		x1 = collect(range(domain[1]+h[1]/2, step=h[1], length=nc[1]))
		return x1
	elseif dim==2
		x1 = collect(range(domain[1]+h[1]/2, step=h[1], length=nc[1]))
		x2 = collect(range(domain[3]+h[2]/2, step=h[2], length=nc[2]))
		return (x1,x2)
	elseif dim==3
		x1 = collect(range(domain[1]+h[1]/2, step=h[1], length=nc[1]))
		x2 = collect(range(domain[3]+h[2]/2, step=h[2], length=nc[2]))
		x3 = collect(range(domain[5]+h[3]/2, step=h[3], length=nc[3]))
		return (x1,x2,x3)
	end
end

# --- linear operators for tensor mesh
function getVolume(Mesh::RegularMesh;saveMat::Bool=true)
# Mesh.V = getVolume(Mesh::RegularMesh) computes volumes v, returns diag(v)
	if isempty(Mesh.V)
		V = sparse(prod(Mesh.h)*I, Mesh.nc,Mesh.nc)
		if saveMat
			Mesh.V = V;
		end
		return V;
	else
		return Mesh.V;
	end
end

function getVolumeInv(Mesh::RegularMesh; saveMat::Bool = true)
# Mesh.Vi = getVolumeInv(Mesh::RegularMesh) returns sdiag(1 ./v)
	if isempty(Mesh.Vi)
		Vi = sparse((1/prod(Mesh.h))*I, Mesh.nc,Mesh.nc)
		if saveMat
			Mesh.Vi = Vi;
		end
		return Vi;
	else
		return Mesh.Vi;
	end
end

function getFaceArea(Mesh::RegularMesh; saveMat = true)
# Mesh.F = getFaceArea(Mesh::RegularMesh) computes face areas a, returns  sdiag(a)
	if isempty(Mesh.F)
		if Mesh.dim==2
			f1  = Mesh.h[2]*sparse(1.0I,Mesh.nf[1],Mesh.nf[1])
			f2  = Mesh.h[1]*sparse(1.0I,Mesh.nf[2],Mesh.nf[2])
			F = blockdiag(f1,f2)
		elseif Mesh.dim==3
			f1  = (Mesh.h[3]*Mesh.h[2])*sparse(1.0I,Mesh.nf[1],Mesh.nf[1])
			f2  = (Mesh.h[3]*Mesh.h[1])*sparse(1.0I,Mesh.nf[2],Mesh.nf[2])
			f3  = (Mesh.h[2]*Mesh.h[1])*sparse(1.0I,Mesh.nf[3],Mesh.nf[3])
			F = blockdiag(blockdiag(f1,f2),f3)
		end
		if saveMat
			Mesh.F = F;
		end
		return F;
	else
		return Mesh.F;
	end
end
function getFaceAreaInv(Mesh::RegularMesh)
# Mesh.Fi = getFaceAreaInv(Mesh::RegularMesh) computes inverse of face areas, returns sdiag(1 ./a)
	if isempty(Mesh.Fi)
		if Mesh.dim==2
			f1i  = (1/Mesh.h[2])*sparse(1.0I,Mesh.nf[1],Mesh.nf[1])
			f2i  = (1/Mesh.h[1])*sparse(1.0I,Mesh.nf[2],Mesh.nf[2])
			Mesh.Fi = blockdiag(f1i,f2i)
		elseif Mesh.dim==3
			f1i  = (1/(Mesh.h[3]*Mesh.h[2]))*sparse(1.0I,Mesh.nf[1],Mesh.nf[1])
			f2i  = (1/(Mesh.h[3]*Mesh.h[1]))*sparse(1.0I,Mesh.nf[2],Mesh.nf[2])
			f3i  = (1/(Mesh.h[2]*Mesh.h[1]))*sparse(1.0I,Mesh.nf[3],Mesh.nf[3])
			Mesh.Fi = blockdiag(blockdiag(f1i,f2i),f3i)
		end
	end
	return Mesh.Fi
end

function getLength(Mesh::RegularMesh)
# Mesh.L = getLength(Mesh::RegularMesh) computes edge lengths l, returns sdiag(l)
	if isempty(Mesh.L)
		if Mesh.dim==2
			l1  = Mesh.h[1]*sparse(1.0I,Mesh.ne[1],Mesh.ne[1])
			l2  = Mesh.h[2]*sparse(1.0I,Mesh.ne[2],Mesh.ne[2])
			Mesh.L   = blockdiag(l1,l2)
		elseif Mesh.dim==3
			l1  = Mesh.h[1]*sparse(1.0I,Mesh.ne[1],Mesh.ne[1])
			l2  = Mesh.h[2]*sparse(1.0I,Mesh.ne[2],Mesh.ne[2])
			l3  = Mesh.h[3]*sparse(1.0I,Mesh.ne[3],Mesh.ne[3])
			Mesh.L   = blockdiag(blockdiag(l1,l2),l3)
		end
	end
	return Mesh.L
end

function getLengthInv(Mesh::RegularMesh)
# Mesh.L = getLength(Mesh::RegularMesh) computes inverse of edge lengths l, returns sdiag(1 ./l)
	if isempty(Mesh.Li)
		if Mesh.dim==2
			l1i  = (1/Mesh.h[1])*sparse(1.0I,Mesh.ne[1],Mesh.ne[1])
			l2i  = (1/Mesh.h[2])*sparse(1.0I,Mesh.ne[2],Mesh.ne[2])
			Mesh.Li   = blockdiag(l1i,l2i)
		elseif Mesh.dim==3
			l1i  = (1/Mesh.h[1])*sparse(1.0I,Mesh.ne[1],Mesh.ne[1])
			l2i  = (1/Mesh.h[2])*sparse(1.0I,Mesh.ne[2],Mesh.ne[2])
			l3i  = (1/Mesh.h[3])*sparse(1.0I,Mesh.ne[3],Mesh.ne[3])
			Mesh.Li  = blockdiag(blockdiag(l1i,l2i),l3i)
		end
	end
	return Mesh.Li
end


text:
using NeXLUncertainties

"""
    StandardLabel

A ReferenceLabel that represents a reference spectrum or reference properties associated with a set of
characteristic x-rays (CharXRay) objects over a contiguous range of spectrum channels that is to be
used as a standard.
"""
struct StandardLabel <: ReferenceLabel
    spectrum::SpectrumOrProperties # The spectrum used as the reference for fitting...
    # roi::UnitRange{Int} # Not available in k-ratios
    xrays::Vector{CharXRay}
    standard::Material
    hash::UInt
    function StandardLabel(cxrl::CharXRayLabel, stdMat::Material)
        new(cxrl.spectrum, cxrl.xrays, stdMat, hash(cxrl.spectrum, hash(cxrl.xrays, hash(stdMat, UInt(0xBADF00D5)))))
    end
    function StandardLabel(props::Dict{Symbol,Any}, xrays::Vector{CharXRay}, stdMat::Material)
        new(props, xrays, stdMat, hash(props, hash(xrays, hash(stdMat, UInt(0xBADF00D5)))))
    end

end

"""
    matches(cxrl::CharXRayLabel, std::KRatio)
    matches(cxrl::CharXRayLabel, std::StandardLabel)

Does the measured `CharXRayLabel` match the standard?
"""
function matches(cxrl::CharXRayLabel, std::KRatio)::Bool
    cp, sp = properties(cxrl), std.stdProps
    return issetequal(cxrl.xrays, std.xrays) && # Same X-rays
        isequal(cp[:Composition], sp[:Composition]) && # Same reference material
        isapprox(cp[:BeamEnergy], sp[:BeamEnergy], rtol=0.001) && # Same beam energy
        isapprox(cp[:TakeOffAngle], sp[:TakeOffAngle], atol=deg2rad(0.1)) # Same take-off-angle
end
function matches(cxrl::CharXRayLabel, std::StandardLabel)::Bool
    cp, sp = properties(cxrl), properties(std)
    return issetequal(cxrl.xrays, std.xrays) && # Same X-rays
        isequal(cp[:Composition], sp[:Composition]) && # Same reference material
        isapprox(cp[:BeamEnergy], sp[:BeamEnergy], rtol=0.001) && # Same beam energy
        isapprox(cp[:TakeOffAngle], sp[:TakeOffAngle], atol=deg2rad(0.1)) # Same take-off-angle
end

function Base.show(io::IO, cl::StandardLabel)
    comp = composition(cl)
    compname = isnothing(comp) ? "Unspecified" : name(comp)
    print(io,"k_std[$(name(cl.xrays)), $(name(cl.standard)) fit by $compname]")
end

"""
`StandardizeModel` is designed to combine k-ratio measurements relative to common reference spectra
to re-standardize a microanalysis measurement.  As a `MeasurementModel` it combines k-ratio data
while maintaining covariance relationships between the measured values.
"""
struct StandardizeModel <: MeasurementModel
    measured::Vector{CharXRayLabel}
    standards::Vector{StandardLabel}

    StandardizeModel(measured::AbstractVector{CharXRayLabel}, standards::AbstractVector{StandardLabel}) = new(measured, standards)
end

function NeXLUncertainties.compute(
    sm::StandardizeModel,
    inputs::LabeledValues,
    withJac::Bool,
)::MMResult
    outputs, results = Array{CharXRayLabel}(undef, length(sm.measured)), zeros(Float64, length(sm.measured))
    for (i, meas) in enumerate(sm.measured)
        j = findfirst(std->matches(meas, std), sm.standards)
        if (!isnothing(j)) && (inputs[sm.standards[j]] > 0.0) # Re-standardized
            std = sm.standards[j]
            sp = copy(properties(std))
            sp[:Composition] = std.standard
            outputs[i] = CharXRayLabel(sp, meas.roi, meas.xrays) # restandardized label
            results[i] = inputs[meas] / inputs[std]
        else
            outputs[i] = meas
            results[i] = inputs[meas]
        end
    end
    if withJac
        jac = zeros(Float64, length(sm.measured), length(sm.measured)+length(sm.standards))
        for (i, meas) in enumerate(sm.measured)
            j = findfirst(std->matches(meas, std), sm.standards)
            if (!isnothing(j)) &&  (inputs[sm.standards[j]] > 0.0) # Re-standardized
                std = sm.standards[j]
                jac[i, indexin(inputs, meas)] = 1.0 / inputs[std]
                jac[i, indexin(inputs, std)] = -inputs[meas]/(inputs[std]^2)
            else
                jac[i, indexin(inputs, meas)] = 1.0
            end
        end
    end
    return (LabeledValues(outputs, results), jac)
end

function __remap_peaktoback(pbs::Dict{<:ReferenceLabel,NTuple{2,Float64}}, sm::StandardizeModel)::Dict{ReferenceLabel,NTuple{2,Float64}}
    res = Dict{ReferenceLabel,NTuple{2,Float64}}()
    for (meas, pb) in pbs
        i = findfirst(std->matches(meas, std), sm.standards)
        if !isnothing(i)
            std = sm.standards[i]
            sp = copy(properties(std))
            sp[:Composition] = std.standard
            res[CharXRayLabel(sp, meas.roi, meas.xrays)] = pb
        else
            res[meas] = pb
        end
    end
    return res
end

"""
    NeXLCore.standardize(ffr::FilterFitResult, standard::FilterFitResult, material::Material, els=elms(material))::FilterFitResult
    NeXLCore.standardize(ffrs::Vector{FilterFitResult}, standard::FilterFitResult, material::Material, els=elms(material))::Vector{FilterFitResult}
    NeXLCore.standardize(ffr::FilterFitResult, standards::AbstractArray{KRatio})::FilterFitResult
    NeXLCore.standardize(ffrs::Vector{FilterFitResult}, standards::AbstractArray{KRatio})::FilterFitResult

Apply the standard `KRatio`s to the `FilterFitResult` producing a re-standardized `FilterFitResult`.
"""
function NeXLCore.standardize(ffr::FilterFitResult, standard::FilterFitResult, material::Material, els=elms(material))::FilterFitResult
    stds = filter(labels(standard.kratios)) do lbl
        (lbl isa CharXRayLabel) && (element(lbl) in els)
    end
    ext = extract(standard.kratios, stds)
    rext = uvs([ StandardLabel(l, material) for l in labels(ext)], ext.values, ext.covariance)
    # Problem here with equivalent CharXRayLabel(s) for unknown and standard...
    inp = cat(ffr.kratios, rext)
    stdize = StandardizeModel(Vector{CharXRayLabel}(labels(ffr.kratios)), Vector{StandardLabel}(labels(rext)))
    return FilterFitResult(
            ffr.label, # Same
            stdize(inp), # remapped
            ffr.roi, # Same
            ffr.raw, # Same
            ffr.residual, # Same - residual is due to references
            __remap_peaktoback(ffr.peakback, stdize) # Remapped to standard labels
        )
end

function NeXLCore.standardize(ffrs::AbstractVector{FilterFitResult}, standard::FilterFitResult, material::Material, els=elms(material))::Vector{FilterFitResult}
    map(ffrs) do ffr
        standardize(ffr, standard, material, els)
    end
end

function NeXLCore.standardize(ffr::FilterFitResult, standards::AbstractArray{KRatio})::FilterFitResult
    @assert all(std->isstandard(std), standards) "One or more required property is missing from a standard (for details, see the above warnings.)"
    stdis = map(cxrl->findfirst(std->matches(cxrl, std), standards), labels(ffr.kratios))
    stds = uvs(map(filter(i->!isnothing(i), stdis)) do i
        std = standards[i]
        StandardLabel(std.stdProps, std.xrays, std.unkProps[:Composition])=>uv(std.kratio)
    end...)
    return if length(stds)>0
        stdize = StandardizeModel(Vector{CharXRayLabel}(labels(ffr.kratios)), Vector{StandardLabel}(labels(stds)))
        inp = cat(ffr.kratios, stds)
        FilterFitResult(
            ffr.label,
            stdize(inp),
            ffr.roi,
            ffr.raw,
            ffr.residual,
            __remap_peaktoback(ffr.peakback, stdize)
        )
    else
        @warn "No suitable standards were provided to standardize this measurement."
        ffr
    end
end

function NeXLCore.standardize(ffrs::AbstractVector{FilterFitResult}, standards::AbstractArray{KRatio})::Vector{FilterFitResult}
    map(ffrs) do ffr
        standardize(ffr, standards)
    end
end


text:
using EEG, Lexicon

save("doc/API.md", EEG)


text:
"""
vectormesh(x::AbstractRange, y::AbstractRange)
Returns a square/rectangular grid of points generated from a x and y ranges.
"""
function vectormesh(x::AbstractRange, y::AbstractRange)

	[collect(i) for i in collect(Iterators.product(x, y))]

end

"""
transform_steps(transform_matrix::Matrix; n_steps::Int64 = 240)
Breaks a 2x2 transformation matrix into incremental n_steps. Useful for interactive/animated transforms.
"""
function transform_steps(transform_matrix::Matrix; n_steps::Int64 = 240)

	basis = [1 0; 0 1]
	transforms = []

	for j in 1:n_steps+1
		push!(transforms, basis .+ (j/n_steps).*(transform_matrix .- basis))
	end

	transforms
end


text:
export bcf, vcf

"""
    bcf(infile::String; ; rename_snp::Bool, silent::Bool, allow_monomorphic::Bool)
Load a BCF file into memory as a PopData object. Population information needs to be provided separately.
- `infile` : path to BCF file (can be gzipped)

**Keyword Arguments**
- `rename_loci` : true/false of whether to simplify loci names to "snp_#" (default: `false`)
- `allow_monomorphic` : true/false of whether to keep monomorphic loci (default: `false`)
- `silent`: true/false of whether to print extra file information (default: `false`).
Alleles are recoded according to the following schema:

| **Base**   |  A   |  T   |  C   |  G   |
| :--------  | :--: | :--: | :--: | :--: |
| **Allele** |  1   |  2   |  3   |  4   |

### Mixed-ploidy data
If importing mixed-ploidy data (such as poolseq), you will need to perform an additional
step to convert the genotype column into the correct `GenoArray` type:
```julia
julia> mydata = bcf("path/to/file.bcf", silent = true, rename_loci = true) ;

julia> mydata.loci.genotype =  mydata.loci.genotype |> Array{Union{Missing, NTuple}}
```
"""
function bcf(infile::String; rename_loci::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false)
    error("Please load in GeneticVariation.jl with \`using GeneticVariation\`")
end

### VCF parsing ###

"""
    vcf(infile::String; ; rename_snp::Bool, silent::Bool, allow_monomorphic::Bool)
Load a VCF file into memory as a PopData object. Population information needs to be provided separately.
- `infile` : path to VCF file (can be gzipped)

**Keyword Arguments**
- `rename_loci` : true/false of whether to simplify loci names to "snp_#" (default: `false`)
- `allow_monomorphic` : true/false of whether to keep monomorphic loci (default: `false`)
- `silent`: true/false of whether to print extra file information (default: `false`).
Alleles are recoded according to the following schema:

| **Base**   |  A   |  T   |  C   |  G   |
| :--------  | :--: | :--: | :--: | :--: |
| **Allele** |  1   |  2   |  3   |  4   |

### Mixed-ploidy data
If importing mixed-ploidy data (such as poolseq), you will need to perform an additional
step to convert the genotype column into the correct `GenoArray` type:
```julia
julia> mydata = vcf("path/to/file.vcf", silent = true, rename_loci = true) ;

julia> mydata.loci.genotype =  mydata.loci.genotype |> Array{Union{Missing, NTuple}}

```
"""
function vcf(infile::String; rename_snp::Bool = false, silent::Bool = false, allow_monomorphic::Bool = false)
    error("Please load in GeneticVariation.jl with \`using GeneticVariation\`")
end