text:
"""
    vec!(x::Vector, X::Matrix)

Perform `copy!(x, vec(X))` without any intermediate arrays.
"""
function vec!{T <: Float}(
    x :: DenseVector{T},
    X :: DenseMatrix{T};
    k :: Int = length(x),
    n :: Int = size(X,1),
    p :: Int = size(X,2)
)
    # check for conformable dimensions
    k == n*p || throw(DimensionMismatch("Arguments x and X must have same number of elements"))

    # will copy X into x in column-major order
    @inbounds for j = 1:p
        @inbounds for i = 1:n
            x[n*(j-1) + i] = X[i,j]
        end
    end

    return nothing
end

"""
    issymmetric(x) -> Bool

Check if a dense `n` x `n` matrix `x` is symmetric. The worst-case complexity occurs for symmetric `x`, where `issymetric` performs `n^2 - n` comparisons and returns `true`. If `x` is not symmetric, then `issymetric returns early with `false`.
"""
function issymetric{T <: Float}(
    x :: DenseMatrix{T}
)
    # size of x?
    (m,n) = size(x)

    # can only work with square matrices
    m == n || throw(ArgumentError("Argument x must be square"))

    # loop through x
    @inbounds for i = 1:n
        @inbounds for j = 1:n
            i == j && continue
            x[i,j] == x[j,i] || return false
        end
    end

    return true
end

"""
    vecnorm(x,y) -> T

Compute the Euclidean distance between two matrices `x` and `y` without introducing any intermediate arrays.
"""
function vecnorm{T <: Float}(
    x :: DenseMatrix{T},
    y :: DenseMatrix{T}
)
    @assert size(x) == size(y)
    s = zero(T)
    @inbounds for j in eachindex(x)
            s += (x[i] - y[i])^2
    end
    return sqrt(s)
end

"""
    mask!(x, v, val, mask_val [, n=length(x)])

A subroutine to mask entries of a vector `x`.

Arguments:

- `x` is the vector to mask.
- `v` is the `Int` vector used in masking.
- `val` is the value of `v` to use in masking.
- `mask_val` is the actual value of the mask, i.e. if `v[i] = val` then `x[i] = mask_val`.

Optional Arguments:

- `n = length(x)`.
"""
function mask!{T <: Float}(
    x        :: DenseVector{T},
    v        :: DenseVector{Int},
    val      :: Int,
    mask_val :: T;
    n        :: Int = length(x)
)
    n == length(v) || throw(BoundsError("Vector x and its mask must have same length"))
    @inbounds for i = 1:n
        if v[i] == val
            x[i] = mask_val
        end
    end
    return nothing
end

"Can also be called with an `Int` argument `n` instead of the data vector `y`."
function cv_get_folds(n::Int, q::Int)
    m, r = divrem(n, q)
    shuffle!([repmat(1:q, m); 1:r])
end

"""
CREATE UNSTRATIFIED CROSSVALIDATION PARTITION

    cv_get_folds(y,q) -> Vector{Int}

This function will partition the `n` components of `y` into `q` disjoint sets for unstratified `q`-fold crossvalidation.

Arguments:

- `y` is the `n`-vector to partition.
- `q` is the number of disjoint sets in the partition.
"""
function cv_get_folds(y::DenseVector, q::Int)
    n, r = divrem(length(y), q)
    shuffle!([repmat(1:q, n); 1:r])
end

"""
    standardize!(j, x [, return_values=false, n=size(x,1), p=size(x,2)])

For the `j`-th column of an `n` x `p` matrix `x`,
compute both the mean `m` and standard deviation `s` of `x[:,j]`
and overwrite `x[:,j]` with its standardized variant.
This function assumes all finite values, so any `NA` or `NaN` will poison the computations!

Arguments:

- `j` is the index of the column of `x` to use in calculations
- `x` is the matrix to use in calculations

Optional Arguments:

- `return_values` is a `Bool` to indicate whether or not to return `m` and `s`. Defaults to `false`.
- `n` and `p` are the dimensions of `x`.

Output:

- if `return_values` is true, then `standardize!` returns the mean `m` and standard deviation `s` of `x[:,j]`.
"""
function standardize!{T <: Float}(
    j :: Int,
    x :: DenseMatrix{T};
    return_values :: Bool = false,
    n :: Int = size(x,1),
    p :: Int = size(x,2)
)

    j <= p || throw(ArgumentError("Column index $i exceeds number of columns $p in x"))

    # accumulation variables
    m = zero(T)
    s = zero(T)

    # temporary float
    t = zero(T)

    # accumulate mean
    @inbounds for i = 1:n
        m += x[i,j]
    end

    # normalize mean
    m /= n

    # calculate variance
    @inbounds for i = 1:n
        t = x[i,j] - m
        t *= t
        x[i,j] = t
        s += t
    end

    # get (unbiased) std
    s = sqrt(s / (n-1))

    # normalize column
    @inbounds for i = 1:n
        x[i,j] /= s
    end

    return_values && return m, s
    return nothing
end

"""
    difference!(Z, X, Y [, a=1.0, b=1.0])

Compute the matrix difference `Z = a*Y - b*Z`, overwriting `Z`.
"""
function difference!{T <: Float}(
    Z :: DenseMatrix{T},
    X :: DenseMatrix{T},
    Y :: DenseMatrix{T};
    a :: T = one(T),
    b :: T = one(T),
)
    m,n = size(Z)
    (m,n) == size(X) == size(Y) || throw(DimensionMismatch("Arguments, Z, X, and Y must have same size"))
    @inbounds for j = 1:n
        @inbounds for i = 1:m
            Z[i,j] = a*X[i,j] - b*Y[i,j]
        end
    end
    return nothing
end

"""
    difference!(x, y, z [, a=1.0, b=1.0, n=length(x)])

Compute the difference `x = a*y - b*z`, overwriting `x`.
"""
function difference!{T <: Float}(
    x :: DenseVector{T},
    y :: DenseVector{T},
    z :: DenseVector{T};
    a :: T = one(T),
    b :: T = one(T),
    n :: Int = length(x)
)
    @inbounds for i = 1:n
        x[i] = a*y[i] - b*z[i]
    end
    return nothing
end

"""
    ypatzmw!(x, y, a, z, w[, n=length(x)])

Compute x = y + a(z - w) in one pass, overwriting x. Equivalent to

    copy!(x,w)
    axpy!(n, 1.0, z, 1, -1.0, x)
    axpy!(n, 1.0, y, 1, a, x)

The unusual name represents *Y* *P*lus *A* *T*imes *Z* *m*inus *W*.
Note that this function discards bounds checks for speed,
so the onous is on the user to supply vectors of the same dimension!
"""
function ypatzmw!{T <: Float}(
    x :: DenseVector{T},
    y :: DenseVector{T},
    a :: T,
    z :: DenseVector{T},
    w :: DenseVector{T};
    n :: Int = length(x))
    @inbounds for i = 1:n
        x[i] = y[i] + a*(z[i] - w[i])
    end
    return nothing
end

"""
    project_k!(x, k)

This function projects a vector `x` onto the set S_k = { y in R^p : || y ||_0 <= k }.
It does so by first finding the pivot `a` of the `k` largest components of `x` in magnitude.
`project_k!` then thresholds `x` by `abs(a)`, sending small components to 0.

Arguments:

- `b` is the vector to project.
- `k` is the number of components of `b` to preserve.
"""
function project_k!{T <: Float64}(
    x    :: DenseVector{T},
    k    :: Int;
)
    a = select(x, k, by = abs, rev = true) :: T
    threshold!(x,abs(a))
    return nothing
end

"""
    project_k!(idx::BitArray{1}, x, k)

This function computes the indices that project a vector `x` onto the set S_k = { y in R^p : || y ||_0 <= k }.
It does so by first finding the pivot `a` of the `k` largest components of `x` in magnitude.
`project_k!` then fills `idx` with the result of `abs(x) .> a`.

Arguments:

- `idx` is a `BitArray` to fill with the projection of `x`.
- `x` is the vector to project.
- `k` is the number of components of `b` to preserve.
"""
function project_k!{T <: Float64}(
    idx  :: BitArray{1},
    x    :: DenseVector{T},
    k    :: Int;
)
    a = select(x, k, by = abs, rev = true) :: T
    threshold!(idx,x,abs(a))
    return nothing
end

"""
    project_k!(b, bk, perm, k)

This function projects a vector `b` onto the set S_k = { x in R^p : || x ||_0 <= k }.
It does so by first finding the indices of the `k` largest components of `x` in magnitude.
Those components are saved in the array `bk`, and then `b` is filled with zeros.
The nonzero components from `bk` are then returned to their proper indices in `b`.

Arguments:

- `b` is the vector to project.
- `bk` is a vector to store the values of the largest `k` components of `b` in magnitude.
- `perm` is an `Int` array that indexes `b`.
- `k` is the number of components of `b` to preserve.
"""
function project_k!{T <: Float}(
    b    :: DenseVector{T},
    bk   :: DenseVector{T},
    perm :: DenseVector{Int},
    k    :: Int;
)
    kk = k == 1 ? 1 : 1:k
    select!(perm, kk, by = (i)->abs(b[i]), rev = true)
    fill_perm!(bk, b, perm, k=k)    # bk = b[perm[kk]]
    fill!(b,zero(T))
    @inbounds for i = 1:k
        b[perm[i]] = bk[i]
    end
    return nothing
end

"""
    project_k!(X, k [, n=size(X,1), p=size(X,2), x=zeros(n), xk=zeros(k), perm=collect(1:n)])

Apply `project_k!` onto each of the columns of a matrix `X`. This enforces the *same* sparsity level on each column.
"""
function project_k!{T <: Float}(
    X    :: DenseMatrix{T},
    k    :: Int;
    n    :: Int = size(X,1),
    p    :: Int = size(X,2),
    x    :: DenseVector{T}   = zeros(n),
)
    length(x)    = n || throw(DimensionMismatch("Arguments X and x must have same row dimension"))
    @inbounds for i = 1:p
        update_col!(x, X, i, n=n, p=p, a=one(T))
        project_k!(x, k)
        update_col!(X, x, i, n=n, p=p, a=one(T))
    end
    return nothing
end

"""
    project_k!(X, K [, n=size(X,1), p=size(X,2), x=zeros(n), xk=zeros(k), perm=collect(1:n)])

Apply `project_k!` onto each of the columns of a matrix `X`, where the `i`th column has sparsity level `K[i]`. This function permits a *different* sparsity level for each column.
"""
function project_k!{T <: Float}(
    X    :: DenseMatrix{T},
    K    :: DenseVector{Int};
    n    :: Int = size(X,1),
    p    :: Int = size(X,2),
    x    :: DenseVector{T} = zeros(eltype(X), n),
)
    length(x) == n || throw(DimensionMismatch("Arguments X and x must have same row dimension"))
    length(K) == p || throw(DimensionMismatch("Argument K must have one entry per column of x"))
    @inbounds for i = 1:p
        k  = K[i]
        update_col!(x, X, i, n=n, p=p, a=one(T))
        project_k!(x, k)
        update_col!(X, x, i, n=n, p=p, a=one(T))
    end
    return nothing
end

"""
    project_k!(x, X, k [, n=size(X,1), p=size(X,2), xk=zeros(k), perm=collect(1:p*n)])

Apply `project_k!` onto the matrix `X` as if it were a vector. The argument `x` facilitates the projection.
Sparsity is enforced on the matrix *as a whole*, so columns may vary in their sparsity.
"""
function project_k!{T <: Float}(
    x    :: DenseVector{T},
    X    :: DenseMatrix{T},
    k    :: Int;
    n    :: Int = size(X,1),
    p    :: Int = size(X,2),
)
    length(x) == n*p || throw(DimensionMismatch("Arguments X and x must have same number of elements"))
    vec!(x, X, k=n*p, n=n, p=p)
    a = select(x, k, by=abs, rev=true) :: T
    threshold!(X, abs(a))
    return nothing
end

"""
    selectpermk!(z, x, k [, p=length(x)])

This subroutine extends `sortperm` to get the top `k` components of a vector in magnitude.
By performing only a partial sort, it saves in compute time and memory.
This function is superceded in Julia v0.4 with `selectperm!`, so it is included only for backward compatibility with Julia v0.3-.

Arguments:

- `z` is the vector of indices to use in the actual sorting process.
- `x` is the vector to use for sorting. Note that `x` itself will not be sorted.
- `k` is the number of elements to choose from the partial sort.

Optional Arguments:

- `p` is the number of elements to sort. Defaults to `length(x)`.
"""
function selectpermk!{T <: Float}(
    z :: DenseVector{Int},
    x :: DenseVector{T},
    k :: Int;
    p :: Int = length(x)
)
    k <= p                 || throw(ArgumentError("selectperm: k cannot exceed length of x!"))
    length(z) == length(x) || throw(DimensionMismatch("Arguments z and x do not have the same length"))
    kk = k == 1 ? 1 : 1:k
    select!(z, kk, by = (i)->abs(x[i]), rev = true)
    return nothing
end

"""
    selectpermk(x, k [, p=length(x), z=collect(1:p(])

This function produces a `p`-vector `z` whose top `k` entries partially sort `x` by magnitude.
"""
function selectpermk{T <: Float}(
    x :: DenseVector{T},
    k :: Int;
    p :: Int              = length(x),
    z :: DenseVector{Int} = collect(1:p)
)
    k <= p || throw(ArgumentError("selectperm: k cannot exceed length of x!"))
    kk = k == 1 ? 1 : 1:k
    return select!(z, kk, by = (i)->abs(x[i]), rev = true)
end

"""
    threshold(X::SparseMatrixCSC, tol)

Send to zero all nonzero values of a sparse vector `X` whose absolute value is less than `tol`.
`threshold!` will throw an error if `X` contains more than one column.
"""
function threshold{T <: Float}(
    x   :: SparseMatrixCSC{T,Int},
    tol :: T
)
    x.n == 1 || throw(ArgumentError("x must be a sparse vector"))
    y = sparsevec(x.rowval,max(abs(x.nzval),tol),x.m)
    return y
end

"""
    threshold!(X::Matrix, tol)

Send to zero all values of a matrix `X` below tolerance `tol` in absolute value.
"""
function threshold!{T <: Float}(
    x   :: DenseMatrix{T},
    tol :: T
)
    m,n = size(x)
    @inbounds for j = 1:n
        @inbounds for i = 1:m
            a = x[i,j]
            if abs(a) < tol
                x[i,j] = zero(T)
            end
        end
    end
    return nothing
end

"""
    threshold!(x, a, tol)

If fed a floating point number `a` in addition to vector `x` and tolerance `tol`,
then `threshold!` will send to zero all components of `x` where `abs(x - a) < tol`.
"""
function threshold!{T <: Float}(
    x   :: DenseVector{T},
    a   :: T,
    tol :: T;
)
    @inbounds for i in eachindex(x)
        x[i] = abs(x[i] - a) < tol ? a : x[i]
    end
    return nothing
end

"""
    threshold!(x, tol)

This subroutine compares the absolute values of the components of a vector `x`
against a thresholding tolerance `tol`. All elements below `tol` are sent to zero.

Arguments:

- `x` is the vector to threshold.
- `tol` is the thresholding tolerance
"""
function threshold!{T <: Float}(
    x   :: DenseVector{T},
    tol :: T;
)
    @inbounds for i in eachindex(x)
        x[i] = abs(x[i]) < tol ? zero(T) : x[i]
    end
    return nothing
end

"""
threshold!(idx::BitArray{1}, x, tol)

This subroutine compares the absolute values of the components of a vector `x`
against a thresholding tolerance `tol`. It then fills `idx` with the result of `abs(x) .> tol`.

Arguments:

- `idx` is the `BitArray` to fill with the thresholding of `x`.
- `x` is the vector to threshold.
- `tol` is the thresholding tolerance
"""
function threshold!{T <: Float}(
    idx :: BitArray{1},
    x   :: DenseVector{T},
    tol :: T;
)
    @assert length(idx) == length(x)
    @inbounds for i in eachindex(x)
        idx[i] = abs(x[i]) >= tol
    end
    return nothing
end

"""
    update_indices!(idx, x)

Calculate the vector of booleans indicating nonzero status of betas.

Arguments:

- `idx` is a `BitArray` of `p` `Bool`s. It contains `true` for each nonzero component of `x` and `false` otherwise.
- `x` is the `p`-vector to index.
"""
function update_indices!{T <: Float}(
    idx :: BitArray{1},
    x   :: DenseVector{T};
)
    @assert length(idx) == length(x)
    @inbounds for i in eachindex(x)
        idx[i] = x[i] != zero(T)# ? true : false
    end
    return nothing
end

"""
    update_col!(z, x, j [, n=size(x,1), p=size(x,2), a=1.0])

This subroutine overwrites an `n`-vector `z` with the `j`th column of a matrix `x` scaled by `a`.
It is more efficient than `z = a*x[:,j]`.

Arguments:

- `z` is the `n`-vector to fill with `x[:,j]`.
- `x` is the `n` x `p` matrix to use in filling `z`.
- `j` indexes the column of `x` to use in filling `z`.

Optional Arguments:

- `n` is the leading dimension of `x`. Defaults to `size(x,1)`.
- `p` is the trailing dimension of `x`. Defaults to `size(x,2)`.
- `a` scales the entries of `z`. Defaults to `1.0` (no scaling).
"""
function update_col!{T <: Float}(
    z :: DenseVector{T},
    x :: DenseMatrix{T},
    j :: Int;
    n :: Int = size(x,1),
    p :: Int = size(x,2),
    a :: T = one(T)
)
    length(z) == n || throw(DimensionMismatch("update_col!: arguments z and X must have same number of rows!"))
    j <= p || throw(DimensionMismatch("update_col!: index j must not exceed number of columns p!"))
    @inbounds for i = 1:n
        z[i] = a*x[i,j]
    end
    return nothing
end

"""
    update_col!(x, z, j [, n=size(x,1), p=size(x,2), a=1.0])

Fills the `j`th column of an `n` x `p` matrix `x` with the entries in the `n`-vector `z`.
"""
function update_col!{T <: Float}(
    x :: DenseMatrix{T},
    z :: DenseVector{T},
    j :: Int;
    n :: Int = size(x,1),
    p :: Int = size(x,2),
    a :: T = one(T)
)
    length(z) == n || throw(DimensionMismatch("update_col!: arguments z and X must have same number of rows!"))
    j <= p || throw(DimensionMismatch("update_col!: index j must not exceed number of columns p!"))
    @inbounds for i = 1:n
        x[i,j] = a*z[i]
    end

    return nothing
end

"""
    update_col!(z, x, j, q [, n=size(x,1), p=size(x,2), a=1.0])

Updates the `q`th column of a matrix `z` with the `j`th column of a matrix `x` scaled by `a`. Both `x` and `z` must have the same leading dimension `n`. Both `j` and `q` must not exceed the trailing dimension `p` of `x`.
"""
function update_col!{T <: Float}(
    z :: DenseMatrix{T},
    x :: DenseMatrix{T},
    j :: Int,
    q :: Int;
    n :: Int = size(x,1),
    p :: Int = size(x,2),
    a :: T = one(T)
)
    size(z,1) == n || throw(DimensionMismatch("update_col!: arguments z and X must have same number of rows!"))
    j <= p         || throw(DimensionMismatch("update_col!: index j must not exceed number of columns p!"))
    @inbounds for i = 1:n
        z[i,q] = a*x[i,j]
    end
    return nothing
end

"""
    update_xk!(xk, x, idx [, k=size(xk,2), p=length(idx), n=size(x,1)])

This subroutine updates a submatrix `xk` of `k` columns of `x` as selected by `idx`.
It efficiently performs `xk = x[:,idx]`.

Arguments:

- `xk` is the `n` x `k` array to fill from `x`.
- `x` is the `n` x `p` array to use in filling `xk`.
- `idx` is a BitArray with a `true` for every column of `x` to write into `xk` and `false` otherwise.

Optional Arguments:

- `k` is the trailing dimension of `xk`. Defaults to `size(xk,2)`.
- `p` is the trailing dimension of `x`. Defaults to `length(idx)` and is used to test conformability with `x`.
- `n` is the leading dimension of `x` and `xk`. Defaults to `size(x,1)`.
"""
function update_xk!{T <: Float}(
    xk  :: DenseMatrix{T},
    x   :: DenseMatrix{T},
    idx :: BitArray{1};
    k   :: Int = size(xk,2),
    p   :: Int = length(idx),
    n   :: Int = size(x,1)
)
    k >= sum(idx) || throw(DimensionMismatch("update_xk!: number of trues in idx must not exceed the number of columns in xk!"))
    p == size(x,2)   || throw(DimensionMismatch("update_xk!: number of components in idx must equal the number of columns in x!"))

    # counter j is used to track the number of trues in idxvec
    j = 0

    # loop over entire vector idxvec
    @inbounds for i = 1:p

        # if current component of idxvec is a true, then increment j and fill column
        if idx[i]
            j += 1
            update_col!(xk, x, i, j, n=n, p=p)
        end

        j == k && return nothing
    end
    return nothing
end

"""
    fill_perm!(x, y, idx)

This subroutine fills a `k`-vector `x` from a `p`-vector `y` via an index vector `idx`.
This variant admits BitArray index vectors.

Arguments:

- `x` is the `k`-vector to fill.
- `y` is the `p`-vector to use in filling `x`.
- `idx` is either a `BitArray` or `Int` vector` that indexes the components of `y` to put into `x`. If `idx` contains `Int`s, then only the first `k` indices are used. Otherwise, `fill_perm!()` traverses `idx` until it encounters `k` `true`s.
"""
function fill_perm!{T <: Float}(
    x   :: DenseVector{T},
    y   :: DenseVector{T},
    idx :: BitArray{1}
)
    # x should have one element per "true" in idx
    k = length(x)
    #@assert k == sum(idx)

    # counter j is used to track the number of trues in idx
    j = 0

    # loop over entire vector idx
    @inbounds for i in eachindex(idx)

        # if current component of idx is a true, then increment j and fill x from y
        if idx[i]
            j += 1
            x[j] = y[i]
        end

        # once x has k components, then it is completely filled and we return it
        j == k && return nothing
    end

    return nothing
end

function fill_perm!{T <: Float}(
    x   :: DenseVector{T},
    y   :: DenseVector{T},
    idx :: DenseVector{Int}
)
    @assert length(x) <= length(idx)
    @inbounds for i in eachindex(x)
            x[i] = y[idx[i]]
    end
    return nothing
end

"""
    update_xb!(xβ, x, β, idx::Vector{Int}, k)

This function efficiently performs the "sparse" matrix-vector product `x*β`, of an `n` x `p` matrix `x` and a `p`-vector `β` with `k` nonzeroes.
The nonzeroes are encoded in the first `k` elements of the `Int vector `indices`.

Arguments:

- `xβ` is the array to overwrite with `x*β`.
- `x` is the `n` x `p` design matrix.
- `β` is the `p`-dimensional parameter vector.
- `idx` is a vector of integers that indexes `β`. The first `k` elements of `idx` should correspond to the `k` nonzeroes of `β`.
- `k` is the number of nonzeroes in `β`.
"""
function update_xb!{T <: Float}(
    xβ  :: DenseVector{T},
    x   :: DenseMatrix{T},
    β   :: DenseVector{T},
    idx :: DenseVector{Int},
    k   :: Int;
)
    k == length(idx) || throw(ArgumentError("Length of `idx` should equal $k"))
    fill!(xβ, zero(T))
    @inbounds for i in idx
        @inbounds for j in eachindex(xβ)
            xβ[j] += β[i]*x[j,i]
        end
    end

    return nothing
end

"""
    update_xb!(xβ, x, b, idx::BitArray{1}, k)

This function efficiently performs the "sparse" matrix-vector product `x*β`, of an `n` x `p` matrix `x` and a `p`-vector `β` with `k` nonzeroes.
The nonzeroes are encoded in the first `k` elements of the `Int vector `indices`.

Arguments:

- `xβ` is the array to overwrite with `x*β`.
- `x` is the `n` x `p` design matrix.
- `β` is the `p`-dimensional parameter vector.
- `idx` is a `BitArray` that indexes `β`. It must have `k` instances of `true`.
- `k` is the number of nonzeroes in `β`.
"""
function update_xb!{T <: Float}(
    xβ  :: DenseVector{T},
    x   :: DenseMatrix{T},
    β   :: DenseVector{T},
    idx :: BitArray{1},
    k   :: Int;
)
    sum(idx) <= k || throw(ArgumentError("Argument indices with $(sum(indices)) trues should have at most $k of them"))
    fill!(xβ, zero(T))
    numtrue = 0
    @inbounds for j in eachindex(idx)
        if idx[j]
            numtrue += 1
            @inbounds for i in eachindex(xβ)
                xβ[i] += β[j]*x[i,j]
            end
        end
        numtrue >= k && break
    end
    return nothing
end

"""
    update_xb(x, b, indices, k [, p=length(b), n=size(x,1)])

This function efficiently performs the "sparse" matrix-vector product `x*b`, of an `n` x `p` matrix `x` and a `p`-vector `b` with `k` nonzeroes.
The nonzeroes are encoded in the first `k` elements of the `Int vector `indices`.
"""
function update_xb{T <: Float}(
    x       :: Matrix{T},
    b       :: Vector{T},
    indices :: DenseVector{Int},
    k       :: Int;
    p       :: Int = length(b),
    n       :: Int = size(x,1)
)
    Xb = zeros(T, n)
    update_xb!(Xb, x, b, indices, k, p=p, n=n)
    return Xb
end

function update_xb{T <: Float}(
    x       :: SharedMatrix{T},
    b       :: SharedVector{T},
    indices :: DenseVector{Int},
    k       :: Int;
    pids    :: DenseVector = procs(),
    p       :: Int = length(b),
    n       :: Int = size(x,1)
)
    Xb = SharedArray(T, n, init = S -> S[localindexes(S)] = zero(T), pids=procs())
    update_xb!(Xb, x, b, indices, k, p=p, n=n)
    return Xb
end

"""
    count_partialnz(x, k0, k)

This subroutine counts nonzeroes in `x` from `k0` to `k` (inclusive).
It is more efficient than countnz(x[k0,k]).
"""
function count_partialnz{T <: Float}(
    x  :: DenseVector{T},
    k0 :: Int,
    k  ::Int
)
    k0 <= k || throw(ArgumentError("Start index must not exceed end index"))
    s = 0
    for i = k0:k
        if x[i] != zero(T)
            s += 1
        end
    end
    return s
end

"""
    fill_partial!(x, a, k0, k)

This subroutine efficiently performs `x[k0:k] = a`.
For more complicated partial fills, consider `fill_perm!()`.
"""
function fill_partial!{T <: Float}(
    x  :: DenseVector{T},
    a  :: T,
    k0 :: Int,
    k  :: Int
)
    1 <= k0 <= k || throw(ArgumentError("fill_partial!: Start index must lie between 1 and end index"))
    @inbounds for i = k0:k
        x[i] = a
    end
    return x
end


text:
"""
$(SIGNATURES)

Get codes of GNSS system as a Matrix where each column
represents a PRN.
```julia-repl
julia> get_code(GPSL1())
```
"""
function get_codes(gnss::AbstractGNSS)
    gnss.codes
end

"""
$(SIGNATURES)

Generate the code signal inplace for PRN-Number `prn` of system `gnss` at chip rate
`code_frequency`, sampled at sampling rate `sampling_frequency`. Make sure, that
`sampling_frequency` is larger than `code_frequency` to avoid overflows with the
modulo calculation.
"""
function gen_code!(
    code::AbstractVector,
    gnss::AbstractGNSS,
    prn::Integer,
    sampling_frequency::Frequency,
    code_frequency::Frequency = get_code_frequency(gnss),
    start_phase = 0.0,
    start_index::Integer = 0
)
    code_frequency > sampling_frequency && error("The code freqeuncy must not be larger than the sampling frequency.")
    num_samples = length(code)
    fixed_point = sizeof(Int) * 8 - 1 - min_bits_for_code_length(gnss)
    FP = Fixed{Int, fixed_point}
    total_code_length = get_code_length(gnss) * get_secondary_code_length(gnss)
    fp_total_code_length = FP(total_code_length)
    delta = FP(code_frequency / sampling_frequency)
    code_phase = FP(mod(FP(mod(start_phase, total_code_length)) + start_index * delta, total_code_length))
    @inbounds for i ∈ 1:num_samples
        code[i] = get_code_unsafe(gnss, code_phase, prn)
        code_phase += delta
        code_phase -= (code_phase >= fp_total_code_length) * fp_total_code_length
    end
    return code
end

"""
$(SIGNATURES)

Generate the code signal for PRN-Number `prn` of system `gnss` at chip rate
`code_frequency`, sampled at sampling rate `sampling_frequency`. Make sure, that
`sampling_frequency` is larger than `code_frequency` to avoid overflows with the
modulo calculation.
"""
function gen_code(
    num_samples::Integer,
    gnss::AbstractGNSS,
    prn::Integer,
    sampling_frequency::Frequency,
    code_frequency::Frequency = get_code_frequency(gnss),
    start_phase = 0.0,
    start_index::Integer = 0
)
    code = zeros(Int16, num_samples)
    gen_code!(
        code,
        gnss,
        prn,
        sampling_frequency,
        code_frequency,
        start_phase,
        start_index
    )
end

"""
$(SIGNATURES)

Get code of type <: `AbstractGNSS` at phase `phase` of PRN `prn`.
```julia-repl
julia> get_code(GPSL1(), 1200.3, 1)
```
"""
Base.@propagate_inbounds function get_code(
    gnss::AbstractGNSS,
    phase,
    prn::Integer
)
    floored_phase = floor(Int, phase)
    get_code_unsafe(
        gnss,
        mod(floored_phase, get_code_length(gnss) * get_secondary_code_length(gnss)),
        prn
    )
end

"""
$(SIGNATURES)

Get code of type <: `AbstractGNSS` at phase `phase` of PRN `prn`.
The phase will not be wrapped by the code length. The phase has to be smaller
than the code length incl. secondary code.
```julia-repl
julia> get_code_unsafe(GPSL1(), 10.3, 1)
```
"""
Base.@propagate_inbounds function get_code_unsafe(
    gnss::AbstractGNSS,
    phase,
    prn::Integer
)
    gnss.codes[floor(Int, phase) + 1, prn]
end

"""
$(SIGNATURES)

Get code to center frequency ratio
```julia-repl
julia> get_code_center_frequency_ratio(GPSL1())
```
"""
@inline function get_code_center_frequency_ratio(gnss::AbstractGNSS)
    get_code_frequency(gnss) / get_center_frequency(gnss)
end

"""
$(SIGNATURES)

Get the minimum number of bits that are needed to represent the code length
"""
@inline function min_bits_for_code_length(gnss::AbstractGNSS)
    ndigits(get_code_length(gnss) * get_secondary_code_length(gnss); base=2)
end

"""
$(SIGNATURES)

Calculate the spectral power of a BPSK modulated signal with chiprate `fc`
at baseband frequency `f`
"""
function get_code_spectrum_BPSK(fc::Frequency, f)
    return get_code_spectrum_BPSK(fc/1Hz, f)
end
function get_code_spectrum_BPSK(fc, f::Frequency)
    return get_code_spectrum_BPSK(fc, f/1Hz)
end
function get_code_spectrum_BPSK(fc::Frequency, f::Frequency)
    return get_code_spectrum_BPSK(fc/1Hz, f/1Hz)
end
function get_code_spectrum_BPSK(fc, f)
    return sinc(f/fc)^2 / fc
end

"""
$(SIGNATURES)
Calculate the spectral power of a sine phased BOC modulated signal with chiprate
`fc` and subcarrier frequency `fs` at baseband frequency `f`
"""
function get_code_spectrum_BOCsin(fc::Frequency, fs::Frequency, f)
    return get_code_spectrum_BOCsin(fc/1Hz, fs/1Hz, f)
end
function get_code_spectrum_BOCsin(fc, fs, f::Frequency)
    return get_code_spectrum_BOCsin(fc, fs, f/1Hz)
end
function get_code_spectrum_BOCsin(fc::Frequency, fs::Frequency, f::Frequency)
    return get_code_spectrum_BOCsin(fc/1Hz, fs/1Hz, f/1Hz)
end
function get_code_spectrum_BOCsin(fc, fs, f)
    return ((sinc(f/fc) * tan(pi*f/(2*fs)))^2/ fc)
end

"""
$(SIGNATURES)
Calculate the spectral power of a cosine phased BOC modulated signal with chiprate
`fc` and subcarrier frequency `fs` at baseband frequency `f`
"""
function get_code_spectrum_BOCcos(fc::Frequency, fs::Frequency, f)
    return get_code_spectrum_BOCcos(fc/1Hz, fs/1Hz, f)
end
function get_code_spectrum_BOCcos(fc, fs, f::Frequency)
    return get_code_spectrum_BOCcos(fc, fs, f/1Hz)
end
function get_code_spectrum_BOCcos(fc::Frequency, fs::Frequency, f::Frequency)
    return get_code_spectrum_BOCcos(fc/1Hz, fs/1Hz, f/1Hz)
end
function get_code_spectrum_BOCcos(fc, fs, f)
    return (2 * sinc(f/fc) * sinpi(f/4fs)^2 / cospi(f/2fs))^2 / fc
end

ext: jl


text:
__precompile__()

module AlgebraicDiffEq

  using DiffEqBase
  import DiffEqBase: solve
  using Parameters

end # module


text:
using LinearAlgebra

function compute_spectrum(A, k, params, Nmodes)

λ_in, λvec_in = eigen(A);
        index = sortperm(imag(λ_in), rev=true);
         λvec = λvec_in[:,index];
            λ = λ_in[index];

       growth = imag(λ[1:Nmodes]);
         freq = real(λ[1:Nmodes]);
        modes = λvec[:,1:Nmodes];

    return growth, freq, modes
end


text:
# HS problem 30 in NLS format
#
#   Source:
#   W. Hock and K. Schittkowski,
#   Test examples for nonlinear programming codes,
#   Lecture Notes in Economics and Mathematical Systems 187,
#   Springer Verlag Berlin Heidelberg, 1981
#   10.1007/978-3-642-48320-2
#
# A. S. Siqueira, Curitiba/BR, 04/2018.

export hs30

"Hock-Schittkowski problem 30 in NLS format"
function hs30()

  model = Model()
  lvar = [1.0; -10.0; -10.0]
  @variable(model, lvar[i] <= x[i=1:3] <= 10, start=1.0)
  @NLexpression(model, F[i=1:3], x[i] + 0.0)
  @NLconstraint(model, x[1]^2 + x[2]^2 >= 1.0)

  return MathProgNLSModel(model, F, name="hs30")
end


text:
using LinearAlgebra: norm, mul!
using Distances: evaluate, PeriodicEuclidean

export AbstractCell
export PeriodicCell
export InfiniteCell
export set_periodicity!
export set_vectors!
export apply_cell_boundaries!
export evaluate_periodic_distance
export check_atoms_in_cell

const periodic_distance = PeriodicEuclidean([1, 1, 1])

abstract type AbstractCell end

struct InfiniteCell <: AbstractCell end

"""
    PeriodicCell{T<:AbstractFloat} <: AbstractCell

Optionally periodic cell
"""
struct PeriodicCell{T<:AbstractFloat} <: AbstractCell
    vectors::Matrix{T}
    inverse::Matrix{T}
    periodicity::Vector{Bool}
    tmp_vector1::Vector{T}
    tmp_vector2::Vector{T}
    tmp_bools::Vector{Bool}
    function PeriodicCell{T}(vectors::AbstractMatrix, periodicity::Vector{Bool}) where {T}
        new{T}(vectors, inv(vectors), periodicity,
            zeros(size(vectors)[1]), zeros(size(vectors)[1]), zeros(Bool, size(vectors)[1]))
    end
end

Base.eltype(::PeriodicCell{T}) where {T} = T

function PeriodicCell(vectors::AbstractMatrix{T}) where {T<:AbstractFloat}
    PeriodicCell{T}(vectors, [true, true, true])
end

function PeriodicCell(vectors::AbstractMatrix{<:Integer})
    PeriodicCell{Float64}(vectors, [true, true, true])
end

function set_periodicity!(cell::PeriodicCell, periodicity::Vector{Bool})
    cell.periodicity .= periodicity
end

function set_vectors!(cell::PeriodicCell, vectors::Matrix)
    cell.vectors .= vectors
    cell.inverse .= inv(cell.vectors)
end

function apply_cell_boundaries!(cell::PeriodicCell, R::AbstractMatrix)
    @views for i in axes(R, 2) # atoms
        apply_cell_boundaries!(cell, R[:,i])
    end
end
apply_cell_boundaries!(::InfiniteCell, ::AbstractArray) = nothing

function apply_cell_boundaries!(cell::PeriodicCell, R::AbstractVector)
    mul!(cell.tmp_vector1, cell.inverse, R)
    for j in axes(R, 1) # DoFs
        if cell.periodicity[j]
            cell.tmp_vector1[j] = mod(cell.tmp_vector1[j], 1)
        end
    end
    mul!(R, cell.vectors, cell.tmp_vector1)
end

"""
    check_atoms_in_cell(cell::PeriodicCell, R::AbstractMatrix)::Bool

True if all atoms are inside the cell, false otherwise.
"""
function check_atoms_in_cell(cell::PeriodicCell, R::AbstractMatrix)::Bool
    @views for i in axes(R, 2) # atoms
        mul!(cell.tmp_vector1, cell.inverse, R[:,i])
        @. cell.tmp_bools = (cell.tmp_vector1 > 1) | (cell.tmp_vector1 < 0)
        any(cell.tmp_bools) && return false
    end
    true
end

function evaluate_periodic_distance(cell::PeriodicCell, r1::AbstractVector, r2::AbstractVector)
    mul!(cell.tmp_vector1, cell.inverse, r1)
    mul!(cell.tmp_vector2, cell.inverse, r2)
    evaluate(periodic_distance, cell.tmp_vector1, cell.tmp_vector2)
end


text:
using Random

struct SingleVertexHeuristic <: NearOptimalMethod
    batch_size::Int
end

struct RandomizedSingleVertexHeuristic <: NearOptimalMethod
    batch_size::Int
end

function solve_near_optimal(bp::BilevelProblem, mt::SingleVertexHeuristic, δ::Real, optimizer; silent=true, sublimit=500)
    eta = mt.batch_size
    model = Model(optimizer)
    silent && MOI.set(model, MOI.Silent(), true)
    (model, x, v, λ, σ, s, _, _, _, _, _) = bilevel_optimality(model, bp, upperlevel = true)
    (mu, _) = size(bp.G)
    add_upper_level(model, bp, x, v)
    optimize!(model)
    st = termination_status(model)
    if st != MOI.OPTIMAL && st != MOI.DUAL_INFEASIBLE
        return (st, model, nothing, nothing)
    end
    dual_problems = build_dual_adversarial(bp, optimizer)
    for (m, _, _) in dual_problems
        optimize!(m)
        if termination_status(m) != MOI.OPTIMAL
            return (MOI.INFEASIBLE, k)
        end
    end
    C=BitSet()
    count = 1 # number of subproblems added at outer iteration
    while count > 0 # outer iteration
        k = 1
        while k in C
            k += 1
        end
        count = 0
        xv = JuMP.value.(x)
        vv = JuMP.value.(v)
        while count <= eta && k <= mu
            (m, α, β) = dual_problems[k]
            MOI.set(m, MOI.TimeLimitSec(), sublimit)
            @objective(m, Min, dot(α, bp.b - bp.A * xv) + β * (dot(bp.d, vv) + δ))
            optimize!(m)
            if termination_status(m) != MOI.OPTIMAL
                @info "Unexpected dual status $k status: $(termination_status(m))"
                return (termination_status(m), :DualError, model, k)
            end
            # new vertex to add
            if objective_value(m) > bp.q[k] - dot(bp.G[k,:], xv)
                push!(C, k)
                αv = JuMP.value.(α)
                βv = JuMP.value.(β)
                @constraint(model,
                    dot(αv, bp.b - bp.A * x) + βv * (dot(bp.d, v) + δ) <= bp.q[k] - dot(bp.G[k,:], x)
                )
                count += 1
            end
            k0 = k
            while k in C || k == k0
                k += 1
            end
        end
        optimize!(model)
        st = termination_status(model)
        if st != MOI.OPTIMAL
            return (st, :NotFound, model)
        end
    end
    return (MOI.OPTIMAL, model, x, v, C)
end

function solve_near_optimal(bp::BilevelProblem, mt::RandomizedSingleVertexHeuristic, δ::Real, optimizer; silent=true, sublimit=500, rng=Random.MersenneTwister(33))
    eta = mt.batch_size
    model = Model(optimizer)
    silent && MOI.set(model, MOI.Silent(), true)
    (model, x, v, λ, σ, s, _, _, _, _, _) = bilevel_optimality(model, bp, upperlevel = true)
    (mu, _) = size(bp.G)
    add_upper_level(model, bp, x, v)
    optimize!(model)
    st = termination_status(model)
    if st != MOI.OPTIMAL && st != MOI.DUAL_INFEASIBLE
        return (st, model, nothing, nothing)
    end
    dual_problems = build_dual_adversarial(bp, optimizer)
    for (m, _, _) in dual_problems
        optimize!(m)
        if termination_status(m) != MOI.OPTIMAL
            return (MOI.INFEASIBLE, k)
        end
    end
    C=BitSet()
    count = 1 # number of subproblems added at outer iteration
    while count > 0 # outer iteration
        korder = permute!(
            [k for k in 1:mu if k ∉ C],
            randperm(rng, mu - length(C)),
        )
        k = pop!(korder)
        count = 0
        xv = JuMP.value.(x)
        vv = JuMP.value.(v)
        while count <= eta
            @info "expanding sub-problem $k"
            (m, α, β) = dual_problems[k]
            MOI.set(m, MOI.Silent(), true)
            MOI.set(m, MOI.TimeLimitSec(), sublimit)
            @objective(m, Min, dot(α, bp.b - bp.A * xv) + β * (dot(bp.d, vv) + δ))
            optimize!(m)
            if termination_status(m) != MOI.OPTIMAL
                @info "Unexpected dual status $k status: $(termination_status(m))"
                return (termination_status(m), :DualError, model, k)
            end
            # new vertex to add
            if objective_value(m) > bp.q[k] - dot(bp.G[k,:], xv)
                push!(C, k)
                αv = JuMP.value.(α)
                βv = JuMP.value.(β)
                @constraint(model,
                    dot(αv, bp.b - bp.A * x) + βv * (dot(bp.d, v) + δ) <= bp.q[k] - dot(bp.G[k,:], x)
                )
                count += 1
            end
            if isempty(korder)
                break
            end
            k = pop!(korder)
        end
        optimize!(model)
        st = termination_status(model)
        if st != MOI.OPTIMAL
            return (st, :NotFound, model)
        end
    end
    return (MOI.OPTIMAL, model, x, v, C)
end

function build_dual_adversarial(bp, optimizer)
    (mu, nl) = size(bp.H)
    ml = size(bp.B, 1)
    return map(1:mu) do k
        m = Model(optimizer)
        @variable(m, α[i=1:ml] >= 0)
        @variable(m, β >= 0)
        @constraint(m,
            cons[j=1:nl],
            sum(bp.B[i,j] * α[i] for i in 1:ml) + bp.d[j] * β ≥ bp.H[k,j]
        )
        (m, α, β)
    end
end


text:
using jInv.Utils
using Test

a = rand(1:300,300)
ua = unique(a)

# testing uniqueidx
b,ii,jj = uniqueidx(a)
b2      = sortunique(a)
t1 = sort(ua)
t2 = sortperm(ua)

@test all(b .== t1)
@test all(a[ii] .==ua[t2])
@test all(b[jj].==a)
@test b == b2


text:
import ClimateMachine.BalanceLaws:
# declaration
    vars_state,
# initialization
    init_state_prognostic!,
    init_state_auxiliary!,
    nodal_init_state_auxiliary!,
# rhs computation
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
# boundary conditions
    boundary_conditions,
    boundary_state!

import ClimateMachine.DGMethods.NumericalFluxes:
    CentralNumericalFluxGradient,
    CentralNumericalFluxSecondOrder,
    NumericalFluxFirstOrder,
    NumericalFluxSecondOrder,
    RusanovNumericalFlux,
    numerical_boundary_flux_second_order!,
    numerical_flux_second_order!,
    numerical_boundary_flux_first_order!,
    numerical_flux_first_order!,
    normal_boundary_flux_second_order!

##
#     Declaration of variables
##
"""
    vars_state
    - returns a NamedTuple of data types.
"""
function vars_state(model::ModelSetup, aux::Auxiliary, T)

    @vars begin
        x::T
        y::T
        z::T
        T_sfc::T  # land surface temperature
    end
end

function vars_state(::ModelSetup, ::Prognostic, T)
    @vars begin
        ρ::T
        ρu::SVector{3, T}
        ρθ::T
        F_ρθ_accum::T
    end
end

function vars_state(::ModelSetup, ::Gradient, T)
    @vars begin
        ∇ρ::T
        ∇u::SVector{3, T}
        ∇θ::T
    end
end

function vars_state(::ModelSetup, ::GradientFlux, T)
    @vars begin
        μ∇ρ::SVector{3, T}
        ν∇u::SMatrix{3, 3, T, 9}
        κ∇θ::SVector{3, T}
    end
end

##
#     Initialization of variables
##
"""
    init_state_xyz!
    - sets up the initial fields within our state variables
    (e.g., prognostic, auxiliary, etc.), however it seems to not initialized
    the gradient flux variables by default.
"""

function nodal_init_state_auxiliary!(m::ModelSetup, state_auxiliary, tmp, geom)
    init_state_auxiliary!(m, m.physics.orientation, state_auxiliary, geom)

    state_auxiliary.T_sfc = 0
end

function init_state_auxiliary!(::ModelSetup, ::Union{NoOrientation, FlatOrientation}, state_auxiliary, geom)
    FT = eltype(state_auxiliary)
    _grav = FT(grav(param_set))
    @inbounds r = geom.coord[3]
    state_auxiliary.x = geom.coord[1]
    state_auxiliary.y = geom.coord[2]
    state_auxiliary.z = geom.coord[3]
end

function init_state_prognostic!(model::ModelSetup, state::Vars, aux::Vars, localgeo, t)
    x = aux.x
    y = aux.y
    z = aux.z

    parameters = model.parameters
    ic = model.initial_conditions

    state.ρ = ic.ρ(parameters, x, y, z)
    state.ρu = ic.ρu(parameters, x, y, z)
    state.ρθ = ic.ρθ(parameters, x, y, z)

    state.F_ρθ_accum = 0

    return nothing
end

"""
    LHS computations
"""
@inline function flux_first_order!(model::ModelSetup, flux::Grad, state::Vars, aux::Vars, t::Real, direction)
    flux.ρu += calc_pressure(model.physics.eos, state) * I
    calc_advective_flux!(flux, model.physics.advection, state, aux, t)

    return nothing
end

@inline function compute_gradient_argument!(model::ModelSetup, grad::Vars, state::Vars, aux::Vars, t::Real)
    calc_diffusive_flux_argument!(grad, model.physics.diffusion, state, aux, t)

    return nothing
end

@inline function compute_gradient_flux!(model::ModelSetup, gradflux::Vars, grad::Grad, state::Vars, aux::Vars, t::Real)
    calc_diffusive_flux!(gradflux, model.physics.diffusion, grad, state, aux, t)

    return nothing
end

@inline function flux_second_order!(
    model::ModelSetup,
    flux::Grad,
    state::Vars,
    gradflux::Vars,
    ::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρ += gradflux.μ∇ρ
    flux.ρu += gradflux.ν∇u
    flux.ρθ += gradflux.κ∇θ

    return nothing
end

"""
    RHS computations
"""
@inline function source!(model::ModelSetup, source::Vars, state::Vars, gradflux::Vars, aux::Vars, t::Real, direction)
    source.F_ρθ_accum = calculate_land_sfc_fluxes(model, state, aux, t)

    return nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine

max_forks_repo_name: CliMA/CouplerMachine


text:
#RSDeltaSigmaPort: Synthesize NTFs
#-------------------------------------------------------------------------------

#==Warnings
===============================================================================#
function _synthesizeNTF_warn_Hinf()
	msg = string(@__FILE__, ": Unable to achieve specified Hinf.\n",
		"Setting all NTF poles to zero."
	)
	@warn(msg)
end
function _synthesizeNTF_warn_iter()
	msg = string(@__FILE__, ": Danger! Iteration limit exceeded.")
	@error(msg)
end
function _synthesizeNTF_warn_Hinf_iter()
	msg = string(@__FILE__, ": Danger! Hinf iteration limit exceeded.")
	@error(msg)
end

#==synthesizeNTF_setup()
===============================================================================#
function synthesizeNTF_setup(order::Int, OSR::Int, opt, f0::Float64, isNTF1::Bool)
	local z
	order = Float64(order) #Type stability
	dw = π/OSR

	if f0!=0		#Bandpass design-- halve the order temporarily.
		order = order/2
		dw = π/(2*OSR)
	end
	order_i = array_round(order)

	#Determine the zeros.
	if length(opt)==1
		if opt==0
			z = zeros(order_i)
		else
			opt_for_zeros = isNTF1 ? (1+mod(opt-1,2)) : opt
			z = dw*ds_optzeros(order, opt_for_zeros)
			if isempty(z)
				throw_unreachable()
				return;
			end
		end
		if f0!=0		# Bandpass design-- shift and replicate the zeros.
			order = order*2
			z = z .+ 2π*f0
			z = conj.(z)
			z = vcat(z, -z)
		end
		z = exp.(j*z)
	else
		z = opt[:] * complex(1.0) #Ensure is Vector{Complex{Float64}}
	end

	#Order might possibly be a N.5 value (float) here. Not sure if this is actually the intent.
	return (order, z) #Order potentially gets modified. Not sure this is actually the intent.
end

#==synthesizeNTF0()
===============================================================================#
function _synthesizeNTF0_z(order::Float64, OSR::Int, z::Vector{Complex{Float64}}, H_inf::Float64, f0::Float64)
	order_i = array_round(order)
	Ts = 1.0
	ntf = _zpk(z,zeros(order_i),1,Ts)
	itn_limit = 100
	#ntf_z, ntf_p, ntf_k, = zpkdata(ntf)

	# Iteratively determine the poles by finding the value of the x-parameter
	# which results in the desired H_inf.
	if f0 == 0			# Lowpass design
		HinfLimit = 2^order  # !!! The limit is actually lower for opt=1 and low OSR
		if H_inf >= HinfLimit
			_synthesizeNTF_warn_Hinf()
			ntf.p = zeros(order_i)
		else
			x=0.3^(order-1) # starting guess
			converged = 0
			local fprev = 0 #Initialize
			local delta_x = 0 #Initialize
			for itn=1:itn_limit
				me2 = -0.5*(x^(2/order))
				w = (2*collect(1:order_i) .- 1)*π/order #VERIFYME
				mb2 = 1 .+ me2*exp.(j*w)
				p = mb2 - sqrt.(mb2 .^ 2 .- 1)
				out = findall(abs.(p) .> 1)
				p[out] = 1 ./ p[out] # reflect poles to be inside the unit circle.
				ntf.z = z; ntf.p = cplxpair(p)
				f = real(evalTF(ntf,-1))-H_inf
				# [ x f ]
				if itn==1
					delta_x = -f/100
				else
					delta_x = -f*delta_x/(f-fprev)
				end

				xplus = x+delta_x
				if xplus>0
					x = xplus
				else
					x = x*0.1
				end
				fprev = f;

				if abs(f)<1e-10 || abs(delta_x)<1e-10
					converged = 1
					break
				end
				if x>1e6
					_synthesizeNTF_warn_Hinf()
					ntf.z = z; ntf.p = zeros(order_i)
					break;
				end
				if itn == itn_limit
					_synthesizeNTF_warn_iter()
				end
			end
		end
	else # Bandpass design.
		x = 0.3^(order/2-1)	# starting guess (not very good for f0~0)
		if f0>0.25
			z_inf=1
		else
			z_inf=-1
		end
		c2pif0 = cos(2π*f0)
		local fprev = 0 #Initialize
		local delta_x = 0 #Initialize
		for itn=1:itn_limit
			e2 = 0.5*x^(2/order)
			w = (2*collect(1:order_i) .- 1)*π/order #VERIFYME

			mb2 = c2pif0 .+ e2*exp.(j*w)
			p = mb2 .- sqrt.(mb2.^2 .- 1)
			# reflect poles to be inside the unit circle.
			out = findall(abs.(p) .> 1)
			p[out] = 1 ./ p[out]
			ntf.z = z; ntf.p = cplxpair(p)
			f = real(evalTF(ntf,z_inf))-H_inf
			# [x f]
			if itn==1
				delta_x = -f/100
			else
				delta_x = -f*delta_x/(f-fprev)
			end

			xplus = x+delta_x
			if xplus > 0
				x = xplus
			else
				x = x*0.1
			end

			fprev = f
			if abs(f)<1e-10 || abs(delta_x)<1e-10
				break
			end
			if x>1e6
				_synthesizeNTF_warn_Hinf()
				ntf.p = zeros(order_i)
				break
			end
			if itn == itn_limit
				_synthesizeNTF_warn_iter()
			end
		end
	end

	ntf.z = cleancplxpair(ntf.z)
	ntf.p = cleancplxpair(ntf.p)
	return ntf
end

function synthesizeNTF0(order::Int, OSR::Int, opt, H_inf::Float64, f0::Float64)
	(order, z) = synthesizeNTF_setup(order, OSR, opt, f0, false)
	_synthesizeNTF0_z(order, OSR, z, H_inf, f0)
end

#==synthesizeNTF1()
===============================================================================#
function _synthesizeNTF1_z(order::Float64, osr::Int, opt, z::Vector{Complex{Float64}}, H_inf::Float64, f0::Float64)
	order_i = array_round(order)
	zp = z[findall(angle.(z) .> 0)]
	x0 = (angle.(zp) .- 2π*f0) * osr / π
	if opt==4 && f0 != 0
		#Do not optimize the zeros at f0
		deleteat!(x0, abs.(x0) .< 1e-10)
	end

	Ts = 1.0
	ntf = _zpk(z,zeros(order_i),1,Ts)
	Hinf_itn_limit = 100

	opt_iteration = 5 #Max number of zero-optimizing/Hinf iterations
	while opt_iteration > 0
		#Iteratively determine the poles by finding the value of the x-parameter
		#which results in the desired H_inf.
		ftol = 1e-10
		z_inf = f0>0.25 ? 1 : -1
		if f0 == 0 #Lowpass design
			HinfLimit = 2^order #!!! The limit is actually lower for opt=1 and low osr
			if H_inf >= HinfLimit
				_synthesizeNTF_warn_Hinf()
				ntf.p = zeros(order_i,1)
			else
				x=0.3^(order-1) #starting guess
				converged = false
				local delta_x
				local fprev
				for itn in 1:Hinf_itn_limit
					me2 = -0.5*(x^(2/order))
					w = (2*collect(1:order_i) .- 1)*π/order #VERIFYME
					mb2 = 1 .+ me2*exp.(j*w)
					p = mb2 .- sqrt.(mb2 .^ 2 .- 1)
					out = findall(abs.(p) .> 1)
					p[out] = 1/p[out] #reflect poles to be inside the unit circle.
					p = cplxpair(p)
					ntf.z = z;	ntf.p = p
					f = real(evalTF(ntf, z_inf))-H_inf
					#[ x f ]
					if itn==1
						delta_x = -f/100
					else
						delta_x = -f*delta_x/(f-fprev)
					end
					xplus = x+delta_x
					x = (xplus>0) ? xplus : x*0.1
					fprev = f
					if abs(f)<ftol || abs(delta_x)<1e-10
						converged = true
						break
					end
					if x>1e6
						_synthesizeNTF_warn_Hinf()
						ntf.z = z;	ntf.p = zeros(order_i,1)
						break
					end
					if itn == Hinf_itn_limit
						_synthesizeNTF_warn_iter()
					end
				end
			end
		else #Bandpass design.
			x = 0.3^(order/2-1) #starting guess (not very good for f0~0)
			c2pif0 = cos(2π*f0)
			local delta_x
			local fprev
			for itn=1:Hinf_itn_limit
				e2 = 0.5*x^(2/order)
				w = (2*collect(1:order_i) .- 1)*π/order
				mb2 = c2pif0 .+ e2*exp.(j*w)
				p = mb2 .- sqrt.(mb2 .^ 2 .- 1)
				#reflect poles to be inside the unit circle.
				out = findall(abs.(p) .> 1)
				p[out] = 1/p[out]
				p = cplxpair(p)
				ntf.z = z; ntf.p = p
				f = real(evalTF(ntf, z_inf)) - H_inf
				#[x f]
				if itn==1
					delta_x = -f/100
				else
					delta_x = -f*delta_x/(f-fprev)
				end
				xplus = x+delta_x
				x = (xplus>0) ? xplus : x*0.1
				fprev = f
				if abs(f)<ftol || abs(delta_x)<1e-10
					break
				end
				if x>1e6
					_synthesizeNTF_warn_Hinf()
					p = zeros(order_i,1)
					ntf.p = p
					break
				end
				if itn == Hinf_itn_limit
					_synthesizeNTF_warn_Hinf_iter()
				end
			end
		end
		########################################################
		if opt < 3 #Do not optimize the zeros
			opt_iteration = 0
		else
			if f0 == 0
				ub = ones(size(x0))
				lb = zeros(size(x0))
			else
				ub = 0.5*ones(size(x0))
				lb = -ub
			end
			opt = Optim.Options(x_tol=0.001, f_tol=0.01, iterations=100)
			f2min(x) = ds_synNTFobj1(x,p,osr,f0) #Function to minimize (optimum point)
			x = optimize(f2min, x0, lb, ub, Fminbox(GradientDescent())).minimum
			x0 = x
			z = exp.(2π*j*(f0 .+ 0.5/osr*x))
			if f0>0
				z = padt(z, array_round(length(p)/2), val=exp(2π*j*f0))
			end
			z = vcat(z, conj.(z))
			if f0==0
				z = padt(z, length(p), val=1.0)
			end
			ntf.z = z; ntf.p = p
			if abs( real(evalTF(ntf, z_inf)) - H_inf ) < ftol
				opt_iteration = 0
			else
				opt_iteration = opt_iteration - 1
			end
		end
	end

	ntf.z = cleancplxpair(ntf.z)
	ntf.p = cleancplxpair(ntf.p)
	return ntf
end

function synthesizeNTF1(order::Int, OSR::Int, opt, H_inf::Float64, f0::Float64)
	(order, z) = synthesizeNTF_setup(order, OSR, opt, f0, true)
	_synthesizeNTF1_z(order, OSR, opt, z, H_inf, f0)
end

#==synthesizeNTF()
===============================================================================#
"""`NTF = synthesizeNTF(order=3, OSR=64, opt=0, H_inf=1.5, f0=0)`

Synthesize a noise transfer function for a delta-sigma modulator.

# Inputs
 - `order`: order of the modulator
 - `OSR`: oversampling ratio
 - `opt`: flag for optimized zeros\r
 - `--> opt=0` -> not optimized,
 - `--> opt=1` -> optimized,
 - `--> opt=2` -> optimized with at least one zero at band-center,
 - `--> opt=3` -> optimized zeros (Requires MATLAB6 and Optimization Toolbox),
 - `--> opt=[z]` -> zero locations in complex form (TODO?)
 - `H_inf`: maximum NTF gain
 - `f0`: center frequency (1->fs)

# Output
`NTF` is a `_zpk` object containing the zeros and poles of the NTF
(See [`_zpk`](@ref)).

# Note

This is actually a wrapper function which calls either the
appropriate version of `synthesizeNTF` based on the availability
of the `fmincon` function from the Optimization Toolbox
"""
function synthesizeNTF(order::Int=3, osr::Int=64; opt=0, H_inf::Float64=1.5, f0::Float64=0.0)
	if f0 > 0.5
		throw(ErrorException("f0 must be less than 0.5."))
	end
	if f0 != 0 && f0 < 0.25/osr
		@warn(string("(", @__FILE__, ") Creating a lowpass ntf."))
		f0 = 0
	end
	if f0 != 0 && mod(order,2) != 0
		throw(ErrorException("order must be even for a bandpass modulator."))
	end

	if length(opt)>1 && length(opt)!=order
		throw(ErrorException("The opt vector must be of length $order(=order)."))
	end

	if true #Use version that applies optimization
		#throw_notimplemented()
		return synthesizeNTF1(order,osr,opt,H_inf,f0)
	else
		return synthesizeNTF0(order,osr,opt,H_inf,f0)
	end
end

"""`NTF = synthesizeNTF(dsm)`

Synthesize a noise transfer function for a ΔΣ modulator.

# Output
`NTF` is a `_zpk` object containing the zeros and poles of the NTF
(See [`_zpk`](@ref)).

# See also:
 - `dsm` types: [`RealDSM`](@ref), [`QuadratureDSM`](@ref)
 - [`clans()`](@ref): "Closed-loop analysis of noise-shaper." An alternative
   method for selecting NTFs based on the 1-norm of the impulse
   response of the NTF.
 - [`synthesizeChebyshevNTF`](@ref): Select a type-2 highpass Chebyshev NTF.
   This function does a better job than `synthesizeNTF` if `OSR` or `H_inf`
   is low.
 - [`_zpk`](@ref)
""" synthesizeNTF

synthesizeNTF(dsm::RealDSM) =
	synthesizeNTF(dsm.order, dsm.OSR, opt=dsm.opt, H_inf=dsm.Hinf, f0=dsm.f0)

#synthesizeNTF(dsm::QuadratureDSM) =
#	synthesizeQNTF(dsm.order, dsm.OSR, dsm.f0, dsm.NG, dsm.ING)

#Last line


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

name = "NativeFileDialog"
version = v"1.1.6"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/mlabbe/nativefiledialog.git", "2850c97af02e6b4cce6ffb91a0445acb5eb9b241"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/nativefiledialog/src/
if [[ "${target}" == *-mingw* ]]; then
    c++ nfd_common.c nfd_win.cpp -DNDEBUG -DUNICODE -D_UNICODE -lole32 -luuid -Iinclude -O2 -Wall -Wextra -fno-exceptions -shared -o "${libdir}/libnfd.${dlext}"
elif [[ "${target}" == *-apple* ]]; then
    cc nfd_cocoa.m nfd_common.c -DNDEBUG -framework Foundation -framework AppKit -Iinclude -O2 -Wall -Wextra -fno-exceptions -fPIC -shared -o "${libdir}/libnfd.${dlext}"
elif [[ "${target}" == *-linux-* && "${nbits}" == 32 ]]; then
    # patch for old gcc in 32bit linux and other for musl
    atomic_patch -p2 ../../patches/32bit-linux-fix.diff
    cc nfd_common.c nfd_gtk.c -D_FILE_OFFSET_BITS=64 -DNDEBUG -Iinclude `pkg-config --cflags --libs gtk+-3.0` -O2 -Wall -Wextra -fno-exceptions -fPIC -shared -o "${libdir}/libnfd.${dlext}"
else
    cc nfd_common.c nfd_gtk.c -DNDEBUG -Iinclude `pkg-config --cflags --libs gtk+-3.0` -O2 -Wall -Wextra -fno-exceptions -fPIC -shared -o "${libdir}/libnfd.${dlext}"
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> arch(p) != "armv6l", supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct("libnfd", :libnfd)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GTK3_jll", uuid="77ec8976-b24b-556a-a1bf-49a033a670a6"))
    BuildDependency("Xorg_xorgproto_jll")
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")


text:
import OldOnly2
import ScientificTypes
import Test

Test.@testset "OldOnly2.jl" begin
    Test.@test OldOnly2.f(1) == 2
end


text:
using Test
using LogicCircuits

@testset "PSDD circuit loader test" begin

   circuit = load_logic_circuit(zoo_psdd_file("nltcs.10split.psdd"))

   @test circuit isa LogicCircuit
   @test isdecomposable(circuit)

   circuit = load_smooth_logic_circuit(zoo_psdd_file("nltcs.10split.psdd"))

   @test circuit isa LogicCircuit
   @test isdecomposable(circuit)

   circuit, vtree = load_struct_smooth_logic_circuit(
                     zoo_psdd_file("nltcs.10split.psdd"),
                     zoo_vtree_file("nltcs.vtree")) # not sure if this is the right vtree...

   @test circuit isa StructLogicCircuit
   @test vtree isa PlainVtree
   @test isdecomposable(circuit)
   @test issmooth(circuit)

end


text:
#=using Pkg
import Pkg; Pkg.add("MathOptInterface")
import Pkg; Pkg.add("Statistics")
import Pkg; Pkg.add("QuantEcon")
import Pkg; Pkg.add("Interpolations")
import Pkg; Pkg.add("NLsolve")
import Pkg; Pkg.add("Optim")
import Pkg; Pkg.add("Random")
import Pkg; Pkg.add("IterTools")
import Pkg; Pkg.add("JuMP")
import Pkg; Pkg.add("Ipopt")
import Pkg; Pkg.add("BenchmarkTools")
import Pkg; Pkg.add("GaussianProcesses")
=#

using
    MathOptInterface,
    BenchmarkTools,
    LinearAlgebra,
    Statistics,
    QuantEcon,
    Interpolations,
    NLsolve,
    Optim,
    Random,
    IterTools,
    JuMP,
    Ipopt,
    GaussianProcesses,
    Plots
    Random;

δ = 0.1
β = 0.86
ϕ = 0.6
ξ = 0.5
μ = 0.5
ψ = 1/ϕ
b = 1/(ϕ^ϕ * (1-ϕ)^(1-ϕ))

numSectors = 2
numPoints1D = 100
gridMax = 50
gridMin = 10

u(x,y) = (log(x) + log(y))
wInit(x) = (log(x[1]) + log(x[2]))/(1 - β)
function gradwInit(g, x, y)
    g[1] = 1/x/(1-β)
    g[2] = 1/y/(1-β)
end
function hesswInit(h, x, y)
    h[1,1] = -1/x^2/(1-β)
    h[1,2] = 0
    h[2,1] = 0
    h[2,2] = -1/y^2/(1-β)
end
grid = Vector{Float64}[]
yTarg = Array{Float64,2}[]
sTarg = Array{Float64,2}[]
mTarg = Array{Float64,3}[]
xTarg = Array{Float64,3}[]
wTarg = Float64[]
rng = MersenneTwister(1234)
#randvec = 20 .+ 30*rand(rng,numPoints1D)
#for p in product(randvec, randvec)
#    push!(grid, collect(p))
#end
for p in product(LinRange(gridMin,gridMax,numPoints1D),LinRange(gridMin,gridMax,numPoints1D))
    push!(grid, collect(p))
end
#wval = w.(grid)

#the following function generates target values (of the value function) for the ML process
##note that (unlike Sargent and Stachursky) the grid is written in terms of kapital
#this simplified the problem substantially as otherwise intermediate variables became dynamic
##also note that this accepts the smooth function
function TargComb(w, grid, β ; compute_policy = false)

    modTrial = Model(Ipopt.Optimizer)
    set_silent(modTrial)
    @variable(modTrial,  c[1:numSectors, 1:length(grid)] >= 0.0001)
    @variable(modTrial, y[1:numSectors, 1:length(grid)] >= 0.0001)
    @variable(modTrial, kp[1:numSectors, 1:length(grid)]>= 0.0001)
    @variable(modTrial, x[1:numSectors, 1:numSectors, 1:length(grid)]>=0.0001)
    @variable(modTrial, xComb[1:numSectors, 1:length(grid)]>=0.0001)
    @variable(modTrial, m[1:numSectors, 1:numSectors, 1:length(grid)]>=0.0001)
    @variable(modTrial, mComb[1:numSectors, 1:length(grid)]>=0.0001)
    @variable(modTrial, u[1:length(grid)])
    @variable(modTrial, w[1:length(grid)])
    k = grid
    for n in 1:length(grid)
        for i in 1:numSectors
            #output
            @NLconstraint(modTrial, y[i,n] == b * k[n][i]^ϕ * mComb[i,n]^(1-ϕ))
            #future capital
            @constraint(modTrial, kp[i,n] == xComb[i,n] + (1-δ) * k[n][i])
            #supply equals demand (eq'm constraint)
            @constraint(modTrial, 0 == y[i,n] - c[i,n] - sum(m[i,:,n]) - sum(x[i,:,n]))
        end
        @NLconstraint(modTrial, xComb[1,n] == x[1,1,n]^ξ * x[2,1,n]^(1-ξ))
        @NLconstraint(modTrial, xComb[2,n] == x[1,2,n]^ξ * x[2,2,n]^(1-ξ))
        @NLconstraint(modTrial, mComb[1,n] == m[1,1,n]^μ * m[2,1,n]^(1-μ))
        @NLconstraint(modTrial, mComb[2,n] == m[1,2,n]^μ * m[2,2,n]^(1-μ))
        @NLconstraint(modTrial, u[n] == log(c[1,n]) + log(c[2,n]))
        @NLconstraint(modTrial, w[n] == (log(kp[1,n]) + log(kp[2,n]))/(1-β))
#=             register(modTrial, :wInit, 2, wInit, gradwInit)
            @NLconstraint(modTrial, w == wInit(kp[1],kp[2]))
            register(modTrial, :u, 2, u; autodiff = true) =#
    end
    @NLobjective(modTrial, Max, sum(u[i] + β * w[i] for i in 1:length(grid)))
    JuMP.optimize!(modTrial)
    push!(wTarg, JuMP.objective_value(modTrial))
    if compute_policy
        push!(xTarg, value.(x))
        push!(sTarg, value.(c))
        push!(mTarg, value.(m))
        push!(yTarg, value.(y))
    end
    #global intK[n] = value.(f)
    #global prevK = intK
    if compute_policy
        return wTarg
    end
    return
end

function evalBellComb()
 wTargC  = TargComb(wInit, grid, β; compute_policy=true)
    return wTargC
end

#=XTarg = Vector{Float64}[]
       for j = 1:numPoints1D^numSectors
           push!(XTarg, sum(xTarg[j][:,i] for i in 1:numSectors))
       end
MTarg = Vector{Float64}[]
       for j = 1:numPoints1D^numSectors
           push!(MTarg, sum(mTarg[j][:,i] for i in 1:numSectors))
       end
fTarg = sTarg + MTarg + XTarg

MPCf = Vector{Float64}[]
       for j = 1:numPoints1D^numSectors
           push!(MPCf, sTarg[j]./fTarg[j])
       end
MPC = Vector{Float64}[]
       for j = 1:numPoints1D^numSectors
           push!(MPC, sTarg[j]./yTarg[j])
       end
MPM = Vector{Float64}[]
       for j = 1:numPoints1D^numSectors
           push!(MPM, MTarg[j]./yTarg[j])
       end
MPX = Vector{Float64}[]
       for j = 1:numPoints1D^numSectors
           push!(MPX, XTarg[j]./yTarg[j])
       end
# a function for turning vectors of vectors into matrices
#function vvm(x)
#           dim1 = length(x)
#           dim2 = length(x[1])
#           matrix = zeros(Float64, dim1, dim2)
#           @inbounds @fastmath for i in 1:dim1, for j in 1:dim2
#                   matrix[i, j] = x[i][j]
#               end
#           end
#           return matrix
#end
#alternatively
gridM = reduce(hcat, grid)

mz = MeanZero()
ml = MeanLin([2., 3.])
mp = MeanPoly([2. 3.; 1. 2.])
kern = SE(0.0, 0.0)
gp = GP(gridM, wTarg, mz, kern)
#gp = GP(gridM, wTarg, MeanLin([2., 3.]), kern)
GaussianProcesses.optimize!(gp)
predict_f(gp, gridM)[1]
fstar(x) = predict_f(gp,hcat(x))[1]

plot!(contour(gp), heatmap(gp); fmt=:png)
#wTargp(x, y) = predict_f(gp,hcat([x, y]))[1]=#

max_stars_repo_name: uq-aibe/aim

max_issues_repo_name: uq-aibe/aim

max_forks_repo_name: uq-aibe/aim


text:
using TensorOperations
using KrylovKit
using LinearAlgebra

# MPS A-matrix is a 3-index tensor, A[i,s,j]
#    s
#    |
# i -A- j
#
# [s] acts on the local Hilbert space
# [i,j] act on the virtual bonds

# MPO W-matrix is a 4-index tensor, W[i,j,s,t]
#     s
#     |
#  i -W- j
#     |
#     t
#
# [s,t] act on the local Hilbert space,
# [i,j] act on the virtual bonds

function contract_from_left_MPS(A,E,B)
    @tensor temp_1[:] := E[-1,1] * A[1,-2,-3]
    @tensor temp[:] := temp_1[1,2,-1] * conj( B[1,2,-2] )
    return temp
end

## initial E and F matrices for the left and right vacuum states
function initial_L(W,N)
    sW=size(W)
    L = ones(sW[1],1,1)#im*zeros(sW[1],1,1)
    #L[1] = 1
    #println(size(L))
    return L
end

function initial_R(W,N)
    sW=size(W)
    R = ones(sW[2],1,1)#im*zeros(sW[2],1,1)
    #R[end] = 1
    #println(size(R))
    return R
end

function construct_R(MPS, MPO, Blist, N)
    R = Dict()
    R[N] = initial_R(MPO[N],N)
    for i in N:-1:2
        R[i-1] = contract_from_right(R[i], MPS[i], MPO[i])
    end
    return R
end

function construct_L(Alist, MPO, Blist,N)
    L = Dict()
    L[1] = initial_L(MPO[1],N)
    return L
end

## tensor contraction from the left hand side
## +-    +--A-
## |     |  |
## L' =  L--R-
## |     |  |
## +-    +--B-

function contract_from_left(E_L,X,W)
    @tensor temp_1[:] := E_L[1,-3,-4]  * X[1,-2,-1]
    @tensor temp_2[:] := temp_1[-1, -2, -5, 1 ]* conj( X[1, -4,-3] )
    @tensor temp[:] := temp_2[ -1, 2, -3, 3, 1 ] * W[1,-2, 2, 3]
    return temp
end

## tensor contraction from the right hand side
##  -+     -A--+
##   |      |  |
##  -R' =  -W--R
##   |      |  |
##  -+     -B--+
function contract_from_right(E_R,Y,W)
    @tensor temp_1[:] :=  Y[-1,-2,1]  * E_R[ 1, -3,-4]
    @tensor temp_2[:] := temp_1[ -1, -2, -5,  1] * conj( Y[-3,-4,1] )
    @tensor temp[:] := temp_2[ -1, 1, -3, 2, 6] * W[ -2 6 1 2 ]
    return temp
end

function one_site_dmrg( MPS::Dict,
                        MPO::Dict,
                        chi::Int64,
                        sweeps::Int64,
                        N::Int64)

    d = 2

    #MPS=init_random_MPS(d,chi,N)
    MPS=Normalize(MPS)

    L = construct_L(MPS, MPO, MPS, N)
    R = construct_R(MPS, MPO, MPS, N)

    global MPS,L,R

    for sweep in 1:Int(sweeps/2)
        for i in 1:N-1
            println(i)
            Energy,MPS[i],MPS[i+1] = optimize_one_site( MPS[i], MPS[i+1], MPO[i], L[i], R[i], "right", chi)
            println("Sweep ",sweep*2," Sites ",i," ",i+1," Energy ",Energy)
            L[i+1] = contract_from_left(L[i], MPS[i], MPO[i])
        end

        for i in N:-1:2
            Energy,MPS[i-1],MPS[i] = optimize_one_site( MPS[i-1], MPS[i], MPO[i], L[i], R[i], "left", chi)
            println("Sweep ",sweep*2,"Sites ",i," ",i+1,"Energy %-16.15f\n",Energy)
            R[i-1] = contract_from_right(R[i], MPS[i], MPO[i])
        end

    end

    return MPS

end

function optimize_one_site(A, B, W, E, F, dir, chi )
    function H_lin(v)
        ##     +--A--+
        ##     |  |  |
        ##     L--W--R
        ##     |  |  |
        #println("E ",size(E))
        #println("MPS ",size(v))
        #println("F ",size(F))
        #println("W ",size(W))
        @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ 1, -1, -4]
        @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
        @tensor temp[:] := temp_2[1, 2, -1, 3,-3] * W[ 2 3 ; 1 -2 ]
        return temp
    end

    if (dir == "right" )
        sA=size(A)

        en,V,info = eigsolve( H_lin , A ,  1  , :SR ; issymmetric = true  )
        en = en[1]
        V = V[1]
        sV = size(V)

        V=reshape(V,(sV[1]*sV[2],sV[3]))#s i j -> i*s, j
        A,S,V = svd(V,full=false)
        V=V'
        A = reshape( A , ( sA[1], sA[2], :) )
        #"ij,jl,slk->sik"
        S=diagm(S)
        @tensor B[:] := S[-1,1 ] * V[ 1,2 ] * B[2,-2,-3]
    else
        sB=size(B)
        print(sB)
        en,V,info = eigsolve( H_lin , B ,  1  , :SR ; issymmetric = true  )
        en = en[1]
        V = V[1]

        sV=size(V)
        V=reshape(V,(sV[1],sV[2]*sV[3]))  #s i j -> i s j -> i, s*j

        U,S,B = svd(V,full=false)
        B=B'
        B = reshape(B,(:,sB[2],sB[3]))
        #"sij,jk,kl->sil"
        S=diagm(S)
        @tensor A[:] := A[-1,-2,2] * U[ 2,3 ] * S[ 3,-3 ]
    end

    return en, A, B

end

function init_random_MPS(d,m,N)

    MPS=Dict()
    MPS[1] = im*rand(1,d,m)
    for i in 2:N-1
        MPS[i]= im*rand(m,d,m)
    end
    MPS[N] = im*rand(m,d,1)
    return MPS
end

function Construct_Ising_MPO(h,J,N)
    d=2
    D=3
    id = [ 1 0 ; 0 1 ]
    sp = [ 0 1 ; 0 0 ]
    sm = [ 0 0 ; 1 0 ]
    sz = [ 1 0 ; 0 -1 ]
    sx = [ 0 1 ; 1 0 ]
    id = [ 1 0 ; 0 1 ]
    W = im *  zeros(D,D,d,d)
    W1 = im *  zeros(1,D,d,d)
    W2 = im *  zeros(D,1,d,d)
    W[1,1,:,:]=id
    W[2,1,:,:]=sz
    W[3,1,:,:]=-h*sx
    W[3,2,:,:]=-J*sz
    W[3,3,:,:]=id;

    W1[1,1,:,:]=-h*sx
    W1[1,2,:,:]=-J*sz
    W1[1,3,:,:]=id

    W2[1,1,:,:]=id
    W2[2,1,:,:]=sz
    W2[3,1,:,:]=-h*sx

    MPO = Dict()
    MPO[1] = W1
    for i in 2:(N-1)
        MPO[i] = W
    end
    MPO[N] = W2

    return MPO
end

function Normalize(MPS)
    norm = MPS_dot(MPS,MPS)
    for i in 1:length(MPS)
        MPS[i] = MPS[i]/norm^(1/(2*length(MPS)))
    end
    return MPS
end

function MPS_dot(A,B)
    E=ones(1,1)
    for i in 1:length(A)
        E = contract_from_left_MPS( A[i] , E , B[i] )
    end
    return E[1]
end

end

function initial_R(W,N)
    sW=size(W)
    R = ones(sW[2],1,1)#im*zeros(sW[2],1,1)
    #R[end] = 1
    #println(size(R))
    return R
end

function construct_R(MPS, MPO, Blist, N)
    R = Dict()
    R[N] = initial_R(MPO[N],N)
    for i in N:-1:2
        R[i-1] = contract_from_right(R[i], MPS[i], MPO[i])
    end
    return R
end

function construct_L(Alist, MPO, Blist,N)
    L = Dict()
    L[1] = initial_L(MPO[1],N)
    return L
end

## tensor contraction from the left hand side
## +-    +--A-
## |     |  |
## L' =  L--R-
## |     |  |
## +-    +--B-

function contract_from_left(E_L,X,W)
    @tensor temp_1[:] := E_L[1,-3,-4]  * X[-2,1,-1]
    @tensor temp_2[:] := temp_1[-1, -2, -5, 1 ]* conj( X[-4, 1,-3] )
    @tensor temp[:] := temp_2[ -1, 2, -3, 3, 1 ] * W[1,-2, 2, 3]
    return temp
end

## tensor contraction from the right hand side
##  -+     -A--+
##   |      |  |
##  -R' =  -W--R
##   |      |  |
##  -+     -B--+
function contract_from_right(E_R,Y,W)
    @tensor temp_1[:] :=  Y[-2,-1,1]  * E_R[ 1, -3,-4]
    @tensor temp_2[:] := temp_1[ -1, -2, -5,  1] * conj( Y[-4,-3,1] )
    @tensor temp[:] := temp_2[ -1, 1, -3, 2, 6] * W[ -2 6 1 2 ]
    return temp
end

function one_site_dmrg( MPS::Dict,
                        MPO::Dict,
                        chi::Int64,
                        sweeps::Int64,
                        N::Int64)

    d = 2

    #MPS=init_random_MPS(d,chi,N)
    MPS=Normalize(MPS)

    L = construct_L(MPS, MPO, MPS, N)
    R = construct_R(MPS, MPO, MPS, N)

    global MPS,L,R

    for sweep in 1:Int(sweeps/2)
        for i in 1:N-1
            println(i)
            Energy,MPS[i],MPS[i+1] = optimize_one_site( MPS[i], MPS[i+1], MPO[i], L[i], R[i], "right", chi)
            println("Sweep ",sweep*2," Sites ",i," ",i+1," Energy ",Energy)
            L[i+1] = contract_from_left(L[i], MPS[i], MPO[i])
        end

        for i in N:-1:2
            Energy,MPS[i-1],MPS[i] = optimize_one_site( MPS[i-1], MPS[i], MPO[i], L[i], R[i], "left", chi)
            println("Sweep ",sweep*2,"Sites ",i," ",i+1,"Energy %-16.15f\n",Energy)
            R[i-1] = contract_from_right(R[i], MPS[i], MPO[i])
        end

    end

    return MPS

end

function optimize_one_site(A, B, W, E, F, dir, chi )
    function H_lin(v)
        ##     +--A--+
        ##     |  |  |
        ##     L--W--R
        ##     |  |  |
        #println("E ",size(E))
        #println("MPS ",size(v))
        #println("F ",size(F))
        #println("W ",size(W))
        @tensor temp_1[ :] := E[ 1, -2, -3] *  v[ -1, 1, -4]
        @tensor temp_2[:] :=  temp_1[-1,-2,-3,1] * F[1,-4,-5]
        @tensor temp[:] := temp_2[1, 2, -2, 3,-3] * W[ 2 3 ; 1 -1 ]
        return temp
    end

    if (dir == "right" )
        sA=size(A)

        en,V,info = eigsolve( H_lin , A ,  1  , :SR ; issymmetric = true  )
        en = en[1]
        V = V[1]
        sV = size(V)

        V=reshape(V,(sV[1]*sV[2],sV[3]))#s i j -> i*s, j
        A,S,V = svd(V,full=false)
        V=V'
        A = reshape( A , ( sA[1], sA[2], :) )
        #"ij,jl,slk->sik"
        S=diagm(S)
        @tensor B[:] := S[-2,1 ] * V[ 1,2 ] * B[-1,2,-3]
    else
        sB=size(B)
        print(sB)
        en,V,info = eigsolve( H_lin , B ,  1  , :SR ; issymmetric = true  )
        en = en[1]
        V = V[1]

        V=permutedims(V,(2,1,3))
        sV=size(V)
        V=reshape(V,(sV[1],sV[2]*sV[3]))  #s i j -> i s j -> i, s*j

        U,S,B = svd(V,full=false)
        B=B'
        B = permutedims(reshape(B,(:,sB[1],sB[3])),(2,1,3))
        #"sij,jk,kl->sil"
        S=diagm(S)
        @tensor A[:] := A[-1,-2,2] * U[ 2,3 ] * S[ 3,-3 ]
    end

    return en, A, B

end

function init_random_MPS(d,m,N)

    MPS=Dict()
    MPS[1] = im*rand(d,1,m)
    for i in 2:N-1
        MPS[i]= im*rand(d,m,m)
    end
    MPS[N] = im*rand(d,m,1)
    return MPS
end

function Construct_Ising_MPO(h,J,N)
    d=2
    D=3
    id = [ 1 0 ; 0 1 ]
    sp = [ 0 1 ; 0 0 ]
    sm = [ 0 0 ; 1 0 ]
    sz = [ 1 0 ; 0 -1 ]
    sx = [ 0 1 ; 1 0 ]
    id = [ 1 0 ; 0 1 ]
    W = im *  zeros(D,D,d,d)
    W1 = im *  zeros(1,D,d,d)
    W2 = im *  zeros(D,1,d,d)
    W[1,1,:,:]=id
    W[2,1,:,:]=sz
    W[3,1,:,:]=-h*sx
    W[3,2,:,:]=-J*sz
    W[3,3,:,:]=id;

    W1[1,1,:,:]=-h*sx
    W1[1,2,:,:]=-J*sz
    W1[1,3,:,:]=id

    W2[1,1,:,:]=id
    W2[2,1,:,:]=sz
    W2[3,1,:,:]=-h*sx

    MPO = Dict()
    MPO[1] = W1
    for i in 2:(N-1)
        MPO[i] = W
    end
    MPO[N] = W2

    return MPO
end

function Normalize(MPS)
    norm = MPS_dot(MPS,MPS)
    for i in 1:length(MPS)
        MPS[i] = MPS[i]/norm^(1/(2*length(MPS)))
    end
    return MPS
end

function MPS_dot(A,B)
    E=ones(1,1)
    for i in 1:length(A)
        E = contract_from_left_MPS( A[i] , E , B[i] )
    end
    return E[1]
end

max_stars_repo_name: aanniikkeettr/Tensor-Networks

max_issues_repo_name: aanniikkeettr/Tensor-Networks

max_forks_repo_name: aanniikkeettr/Tensor-Networks


text:
const USE_GPU = false
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, Statistics
# pyplot()
do_viz = false

@parallel function check_Err!(ErrH::Data.Array, H::Data.Array)
    @all(ErrH) = @all(ErrH) - @all(H)
    return
end

@parallel function compute_Err!(ErrH::Data.Array, H::Data.Array)
    @all(ErrH) = @all(H)
    return
end

@parallel function compute_flux!(qx::Data.Array, qy::Data.Array, H::Data.Array, n::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(qx) = -@av_xi(H)^n*@d_xi(H)/dx
    @all(qy) = -@av_yi(H)^n*@d_yi(H)/dy
    return
end

@parallel function compute_ResH!(dtau::Data.Array, ResH::Data.Array, H::Data.Array, qx::Data.Array, qy::Data.Array, b::Data.Array, n::Data.Number, dmp::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(dtau) = min(dx*dx, dy*dy)/(1.0 + @inn(H)^n)/4.1./4.0
    @all(ResH) = (-@d_xa(qx)/dx  - @d_ya(qy)/dy + @inn(b)) + dmp*@all(ResH)
    return
end

@parallel function update_H!(H::Data.Array, dtau::Data.Array, ResH::Data.Array)
    @inn(H) = max(0.0, @inn(H) + @all(dtau)*@all(ResH))
    return
end

@views function sia_2D_xpu()
    # Physics
    lx, ly = 10.0, 10.0
    n      = 3.0
    # Numerics
    nx, ny = 128, 128
    niter  = 10000
    nout   = 100
    dmp    = 0.96
    ε      = 1e-8
    # Derived numerics
    dx, dy = lx/nx, ly/ny
    xc     = LinRange(dx/2, lx-dx/2, nx)
    yc     = LinRange(dy/2, ly-dy/2, ny)
    # Array allocation
    qx     = @zeros(nx-1,ny-2)
    qy     = @zeros(nx-2,ny-1)
    dtau   = @zeros(nx-2,ny-2)
    ResH   = @zeros(nx-2,ny-2)
    ErrH   = @zeros(nx  ,ny  )
    # Initial condition
    b      = 0.5*ones(nx,ny)
    H      = Data.Array( exp.(.-(xc.-lx./2.0).^2 .-(yc.-ly./2.0)'.^2) )
    rad    = (xc.-lx./2.0).^2 .+(yc.-ly./2.0)'.^2; b[rad.>lx/4] .= -0.5
    b      = Data.Array(b) # move data back to XPU
    # Time loop
    for iter = 1:niter
        if (iter==11) global t0 = Base.time() end
        @parallel compute_Err!(ErrH, H)
        @parallel compute_flux!(qx, qy, H, n, dx, dy)
        @parallel compute_ResH!(dtau, ResH, H, qx, qy, b, n, dmp, dx, dy)
        @parallel update_H!(H, dtau, ResH)
        @parallel check_Err!(ErrH, H)
        # Visualise
        if mod(iter,nout)==0 && do_viz
            p1 = heatmap(xc, yc, Array(H)', xlabel="lx", ylabel="ly", title="shallow ice, iter=$iter", aspect_ratio=1, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), c=:viridis)
            p2 = plot(xc, Array(H)[:,Int(round(ny/2))], xlabel="lx", ylabel="height", xlims=(xc[1], xc[end]), ylims=(0., 1.2), legend=false, framestyle=:box, aspect_ratio=2.5)
            l = @layout [ a{0.8h} ; b{0.2h} ]
            display(plot(p1, p2, layout=l ))
        end
        if mod(iter,nout)==0 maxErr=maximum(abs.(ErrH)); @printf("iter=%d, max(err)=%1.2e \n", iter, maxErr); if maxErr<ε global itg=iter; break; end end
    end
    time_s = (Base.time()-t0)
    @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", time_s, round((2/1e9*nx*ny*sizeof(lx))/(time_s/(itg-10)), sigdigits=2))
    return
end

@time sia_2D_xpu()

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
__precompile__()

module MarkovTransitionMatrices

using GenGlobal
using Distributions
using Optim
using StatsFuns
using SharedArrays
using Base.Iterators
using Distributed
using LinearAlgebra
using SparseArrays

export markov_transition, markovswitching_transition, markov_transition_moment_matching_parallel, markov_transition_moment_matching_serial

myDist(μ::Real, σ::Real) = Normal(μ, σ)
myDist(μ::Vector, Σ::Matrix) = MvNormal(μ, Σ)
mypdf(dist::UnivariateDistribution, s2) = pdf(dist, s2...)
mypdf(dist::MultivariateDistribution, s2) = pdf(dist, [s2...])

whichP(P::AbstractMatrix, ::Type{Val{true}})  = sparse(P)
whichP(P::SharedMatrix,   ::Type{Val{false}}) = sdata(P)
whichP(P::Matrix,         ::Type{Val{false}}) = P

function sparsify!(P::AbstractMatrix{T}, minp::Real) where {T<:AbstractFloat}
  P ./= sum(P, dims=2)
  P .*= (P .> minp)
  P ./= sum(P, dims=2)
  makesparse = minp > 0.0
  return whichP(P, Val{makesparse})
end

include("simple_no_matching.jl")
include("moment_matching.jl")
include("transition_from_data.jl")
# include("Farmer_Toda_VAR.jl")
include("normal_moments.jl")
include("new-moment-matching.jl")

# module end
end


text:
import Base.(./)

#=
 ./ : Elementwise divide for Assoc
=#

function ./(A::Assoc,B::Assoc)
#First, create the row and col of the intersection
ABrow = sortedintersect(A.row,B.row)
ABcol = sortedintersect(A.col,B.col)
#Filling the sparse matrix with

AA = spzeros(size(ABrow,1), size(ABcol,1))
rowMapping = searchsortedmapping(ABrow,A.row)
colMapping = searchsortedmapping(ABcol,A.col)
AA = A.A[rowMapping,colMapping]
#AA = round(Int64,AA)

BB = spzeros(size(ABrow,1), size(ABcol,1))
rowMapping = searchsortedmapping(ABrow,B.row)
colMapping = searchsortedmapping(ABcol,B.col)
BB = B.A[rowMapping,colMapping]
#BB = round(Int64,BB)

ABA = AA ./ BB
ABA = sparse(ABA)
return Assoc(ABrow,ABcol,promote([1.0],A.val)[1],ABA)
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
import Base.writemime
export writemime

# methods to translate

"""return the frequency ticks of a spectrogram, rounded to the nearest integers"""
function yticklabels(tfr::DSP.Periodograms.Spectrogram, yticks::Array)
    map(Int, map(round, yticks)), "frequency (Hz)"
end

"""return Hz values corresponding to given mel frequencies"""
function yticklabels(tfr::MelSpectrogram, yticks::Array)
    map(Int, map(round, mel_to_hz(yticks))), "frequency (Hz)"
end

"""return MFCC coefficient numbers"""
function yticklabels(tfr::MFCC, yticks::Array)
    map(Int, map(round, yticks)), "MFCC number"
end

heatmap(tfr::DSP.Periodograms.TFR) = log10(power(tfr))

function draw_heatmap(tfr::DSP.Periodograms.TFR)
    X = time(tfr)
    Y = freq(tfr)
    Z = heatmap(tfr)

    PyPlot.pcolormesh(X, Y, Z)
    PyPlot.xlim(first(X), last(X))
    PyPlot.ylim(first(Y), last(Y))

    (yticks, ylabel) = yticklabels(tfr, PyPlot.yticks()[1])
    PyPlot.gca()[:set_yticklabels](yticks)
    PyPlot.gca()[:spines]["top"][:set_visible](false)
    PyPlot.gca()[:spines]["right"][:set_visible](false)

    PyPlot.xlabel("time (seconds)")
    PyPlot.ylabel(ylabel)
end

"""Display a spectrogram"""
function writemime{R <: DSP.Periodograms.TFR}(io::IO, mime::MIME"image/png", tfr::R)
    @eval import PyPlot
    PyPlot.ioff()
    PyPlot.figure(figsize=(8, 4))
    draw_heatmap(tfr)
    writemime(io, mime, PyPlot.gcf())
end

"""Display multichannel spectrogram"""
function writemime{R <: DSP.Periodograms.TFR}(io::IO, mime::MIME"image/png", tfrs::Array{R, 1})
    nchannels = length(tfrs)

    @eval import PyPlot
    PyPlot.ioff()
    PyPlot.figure(figsize=(8, 8))

    for i = 1:nchannels
        PyPlot.subplot(nchannels, 1, i)
        draw_heatmap(tfrs[i])

        if i != nchannels
            PyPlot.gca()[:get_xaxis]()[:set_visible](false)
            PyPlot.gca()[:spines]["bottom"][:set_visible](false)
        end
    end
    writemime(io, mime, PyPlot.gcf())
end


text:
if haskey(ENV, "MOCHA_USE_CUDA")
  const test_gpu = true
else
  const test_gpu = false
end
const test_cpu   = true

using Mocha
using Base.Test

if test_cpu
  backend_cpu = CPUBackend()
  init(backend_cpu)
end

if test_gpu
  backend_gpu = GPUBackend()
  init(backend_gpu)
end

# run test in the whole directory, latest modified files
# are run first, this makes waiting time shorter when writing
# or modifying unit-tests
function test_dir(dir)
  map(reverse(Mocha.glob(dir, r".*\.jl$", sort_by=:mtime))) do file
    include("$dir/$file")
  end
end

test_dir("layers")

if test_gpu
  shutdown(backend_gpu)
end
if test_cpu
  shutdown(backend_cpu)
end

max_stars_repo_name: benmoran/MochaTheano

max_issues_repo_name: benmoran/MochaTheano

max_forks_repo_name: benmoran/MochaTheano


text:
export AbstractAgent, start!, step!, get_action

abstract type AbstractAgent end

"""
    start!(agent::AbstractAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)

Function for restarting the agent.
"""
function start!(agent::AbstractAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)
    throw("Implement start! function for agent $(typeof(agent))")
end

"""
    step!(agent::AbstractAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)

Function to take a step with an agent
"""
function step!(agent::AbstractAgent, env_s_tp1, r, terminal; rng=Random.GLOBAL_RNG, kwargs...)
    throw("Implement step! function for agent $(typeof(agent))")
end

"""
    get_action(agent::AbstractAgent, state; rng=Random.GLOBAL_RNG, kwargs...)

Helper function for developing agents. Returns an action given a state.
"""
function get_action(agent::AbstractAgent, state)
    throw("Implement get Action for agent")
end

export LinearAgents
include("agents/LinearQLearning.jl")
# abstract type AbstractAgent end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
# package for statistical distributions
using Distributions
# initializing pseudorandom number generation
srand(123)

# hidden coefficient of a linear model: y = mx + q
m = 3.0
q = 7.0
samples = 100
error_sigma = 0.5

# explanatory variable X
x = Uniform(0.0, 1000.0)
# there will be some error in measuring it, normally distributed
e = Normal(0.0, error_sigma)
# real values for X
x_latent = rand(x, samples)
# observed values, affected by noise
x_observed = x_latent + rand(e, samples)

# real values for X generated by the latent process
y_latent = x_latent * m + q
# observed values, affected by noise
y_observed = y_latent + rand(e, samples)

# the data available to us
println("x_observed:")
println(repr(x_observed))
println("y_observed:")
println(repr(y_observed))

# (ordinary) least squares estimation
# matrix with 2 columns: x_observed and (1, 1, 1, ...) for q estimation
X = [x_observed ones(length(x_observed))]
pseudoinverse = inv(X'*X)*X'
parameters = pseudoinverse*y_observed

println("actual parameters:")
println(m, ", ", q)
println("estimated parameters:")
println(repr(parameters))

max_stars_repo_name: nhausmna/julia-hello-world

max_issues_repo_name: nhausmna/julia-hello-world

max_forks_repo_name: nhausmna/julia-hello-world


text:
module PkgWithUndefinedExports

export undefined_name

end  # module


text:
using Random
using LinearAlgebra

kde_names(m::HybridNB) = collect(keys(m.c_kdes[m.classes[1]]))
discrete_names(m::HybridNB) = collect(keys(m.c_discrete[m.classes[1]]))

function compare_models!(m3::HybridNB, m4::HybridNB)
    @test m3.classes == m4.classes
    @test m3.priors == m4.priors
    @test kde_names(m3) == kde_names(m4)
    @test discrete_names(m3) == discrete_names(m4)

    for c in m3.classes
        for (p1, p2) = zip(m3.c_discrete[c], m4.c_discrete[c])
            @test p1.second.pairs == p2.second.pairs
            @test p1.first == p2.first
        end
        for (p1, p2) in zip(m3.c_kdes[c], m4.c_kdes[c])
            @test p1.first == p2.first
            @test p1.second.kde.x == p2.second.kde.x
            @test p1.second.kde.density == p2.second.kde.density
        end
    end

end

@testset "Core Functions" begin
    # 6 data samples with 2 variables belonging to 2 classes
    X = [-1.0 -2.0 -3.0 1.0 2.0 3.0;
         -1.0  -1.0 -2.0 1.0 1.0 2.0]
    y = [1, 1, 1, 2, 2, 2]

    @testset "Multinomial NB" begin
        m = MultinomialNB([:a, :b, :c], 5)
        X1 = [1 2 5 2;
             5 3 -2 1;
             0 2 1 11;
             6 -1 3 3;
             5 7 7 1]
        y1 = [:a, :b, :a, :c]
        fit(m, X1, y1)
        @test predict(m, X1) == y1
    end

    @testset "Gaussian NB" begin
        m = GaussianNB(unique(y), 2)
        fit(m, X, y)
        @test predict(m, X) == y
    end

    @testset "Hybrid NB" begin

        N1 = 100000
        N2 = 160000
        Np = 1000

        Random.seed!(0)

        # test with names as Symbols
        perm = Random.randperm(N1+N2)
        labels = [ones(Int, N1); zeros(Int, N2)][perm]
        f_c1 = [0.35randn(N1); 3.0 .+ 0.2randn(N2)][perm]
        f_c2 = [-4.0 .+ 0.35randn(N1); -3.0 .+ 0.2randn(N2)][perm]
        f_d = [rand(1:10, N1); rand(12:25, N2)][perm]

        N = AbstractString
        training_c = Dict{N, Vector{Float64}}("c1" => f_c1[1:end-Np],   "c2" => f_c2[1:end-Np])
        predict_c  = Dict{N, Vector{Float64}}("c1" => f_c1[end-Np:end], "c2" => f_c2[end-Np:end])
        training_d = Dict{N, Vector{Int}}("d1" => f_d[1:end-Np])
        predict_d  = Dict{N, Vector{Int}}("d1" => f_d[end-Np:end])

        model = train(HybridNB, training_c, training_d, labels[1:end-Np])
        y_h = predict(model, predict_c, predict_d)
        @test all(y_h .== labels[end-Np:end])

        mktempdir() do dir
            write_model(model, joinpath(dir, "test.h5"))
            m2 = load_model(joinpath(dir, "test.h5"))
            compare_models!(model, m2)
        end

        #testing reading and writing the model file with Symbols
        m3 = HybridNB(y)
        fit(m3, X, y)
        @test all(predict(m3, X) .== y)

        mktempdir() do dir
            write_model(m3, joinpath(dir, "test.h5"))
            m4 = load_model(joinpath(dir, "test.h5"))
            compare_models!(m3, m4)
        end
    end

    @testset "Restructure features" begin
        M = rand(3, 4)
        V = restructure_matrix(M)
        Mp = to_matrix(V)
        @test all(M .== Mp)
    end

    @testset "Multinomial NB - probabistic predictions" begin

        # some word counts in children's books about colours:
        red = [2 0 1 0 1]
        blue = [4 1 2 3 2]
        green = [0 2 0 6 1]
        X = vcat(red, blue, green)

        # gender of author:
        y = [:m, :f, :m, :f, :m]

        # Lagrangian smoothing replaces above data with
        # red = [2, 0, 1, 0, 1, 1, 1]'
        # blue = [4, 1, 2, 3, 2, 1, 1]'
        # green = [0, 2, 0, 6, 1, 1, 1]'
        # y = [:m, :f, :m, :f, :m, :m, :f]
        # which gives total class counts of :f => 3, :m => :4

        # working out, by hand, estimates of p(color=red|male), etc,
        # with Lagrangian smoothing:
        red_given_m = 5/16
        blue_given_m = 9/16
        green_given_m = 2/16
        red_given_f = 1/15
        blue_given_f = 5/15
        green_given_f = 9/15

        # let `m(r, b, g)` be Naive Bayes prediction of probablity of
        # class `:m`, given counts `red=r`, `blue=b` and
        # `green=g`. Similar for `f(r, b, g)`:
        m_(red, blue, green) =
            4/7*(red_given_m^red)*(blue_given_m^blue)*(green_given_m^green)
        f_(red, blue, green) =
            3/7*(red_given_f^red)*(blue_given_f^blue)*(green_given_f^green)
        normalizer(red, blue, green) =
            m_(red, blue, green) + f_(red, blue, green)
        m(a...) = m_(a...)/normalizer(a...)
        f(a...) = f_(a...)/normalizer(a...)

        # new data:
        red = [1 1]
        blue = [1 2]
        green = [1 3]
        Xnew = vcat(red, blue, green)

        # now get NaiveBayes.jl  predictions:
        model = MultinomialNB([:m, :f], 3, alpha=1)
        fit(model, X, y)
        classes, logprobs = predict_logprobs(model, Xnew)
        @test classes == [:f, :m] # implementation changes might
                                  # change order here?
        probs = exp.(logprobs)
        # NaiveBayes does not normalize probabilities, so:
        col_sums = sum(probs, dims=1)
        probs = probs ./ col_sums
        probs_f = probs[1,:]
        probs_m = probs[2,:]

        # compare with above:
        @test m(Xnew[:,1]...) ≈ probs_m[1]
        @test m(Xnew[:,2]...) ≈ probs_m[2]
        @test f(Xnew[:,1]...) ≈ probs_f[1]
        @test f(Xnew[:,2]...) ≈ probs_f[2]

    end

end