text:
"""
# Summary:
function dgbfa{T}(abd::Array{T,2}, lda::Int64, n::Int64, ml::Int64, mu::Int64)

This function is similar to the function gbtrf! provided in Base.LinAlg.LAPACK.
Julia documentation mentions the following:

*Note that the LAPACK API provided by Julia can and will change in the future.*
*Since this API is not user-facing, there is no commitment to support/deprecate*
*this specific set of functions in future releases.*

Due to this reason dgbfa function has been written as per the one provided in the
original version of these solvers at: [NewtonLib](http://elib.zib.de/pub/elib/codelib/NewtonLib/index.html)
based on LINPACK written for Matlab by Cleve Moler, University of New Mexico, Argonne National Lab.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable    | Description                                                  |
|:------------|:-------------------------------------------------------------|
| abd[lda,n]* | Contains the matrix in band storage. See the comments below. |
| lda         | Leading dimension of the matrix. lda >= 2*ml + mu + 1        |
| n           | Order of the original matrix                                 |
| ml          | Number of diagonals below the main diagonal. 0 <= ml < n     |
| mu          | Number of diagonals above the main diagonal. 0 <= mu < n     |

(* marks inout parameters)

## Output parameters
--------------------
| Variable    | Description                                                  |
|:------------|:-------------------------------------------------------------|
| abd[lda,n]* | An upper triangular matrix in band storage and the multipliers which were used to obtain it. Factorization can be written as a = l*u where l is a product of permutation and unit lower triangular matrices and u is upper triangular. |
| ipvt        | Vector containing the pivot indices                          |
| info = 0    | Normal value                                                 |
| info = k    | if u(k,k) = 0.0. This only indicates that dgbsl will divide by zero if called. |

## Band storage
If a is a band matrix, the following program segment will set up input matrix abd
```
ml = (bandwidth below the diagonal)
mu = (bandwidth above the diagonal)
m = ml + mu + 1
for j = 1:n
    i1 = max(1,j-mu)
    i2 = min(n,j+ml)
    for i = i1:i2
        k = i - j + m
        abd[k,j] = a[i,j]
    end
end
```
"""
function dgbfa{T}(abd::Array{T,2}, lda::Int64, n::Int64, ml::Int64, mu::Int64)
    # Begin
    # Copy the matrix so that it remains unchanged at exit
    ipvt = zeros(Int64,n)
    a = abd[:,:]
    m = ml + mu + 1
    info = 0
    # zero intial fill-in columns
    j0 = mu + 2
    j1 = min(n,m) - 1
    if j1 >= j0
        for jz = j0:j1
            i0 = m + 1 - jz
            a[i0:ml,jz:jz] = zeros(T,ml-i0+1,1)
        end
    end
    jz = j1
    ju = 0
    # Gaussian elimination with partial pivoting
    nm1 = n - 1
    if nm1 >= 1
        for k = 1:nm1
            kp1 = k + 1
            # zero next fill-in column
            jz += 1
            if jz <= n && ml >= 1
                a[1:ml,jz:jz] = zeros(T,ml)
            end
            # find l = pivot index
            lm = min(ml,n-k)
            l = maxabs(a[m:m+lm,k:k]) + m - 1
            ipvt[k] = l + k - m
            # zero pivot implies this column already triangularized
            if a[l,k] != 0.0
                # interchange if necessary
                if l != m
                    t = a[l,k]
                    a[l,k] = a[m,k]
                    a[m,k] = t
                end
                # compute multipliers
                t = -1.0/a[m,k]
                a[m+1:m+lm,k:k] *= t
                # row elimintaion with column indexing
                ju = min(max(ju,mu+ipvt[k]),n)
                mm = m
                if ju >= kp1
                    for j = kp1:ju
                        l = l - 1
                        mm = mm - 1
                        t = a[l,j]
                        if l!= mm
                            a[l,j] = a[mm,j]
                            a[mm,j] = t
                        end
                        a[mm+1:mm+lm,j:j] += t*a[m+1:m+lm,k:k]
                    end
                end
            else
                info = k
            end
        end
    end
    ipvt[n] = n
    if a[m,n] == 0.0
        info = n
    end
    return (a,ipvt,info)
end

"""
function dgbsl{T}(abd::Array{T,2}, lda::Int64, n::Int64, ml::Int64, mu::Int64,
    ipvt::Vector{Int64}, b::Vector{T}, flag::Int64)

Solves the double precision band system using the factors computed by dgbfa.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable    | Description                                                  |
|:------------|:-------------------------------------------------------------|
| abd[lda,n]  | Matrix output from dgbfa                                     |
| lda         | Leading dimension of the matrix abd                          |
| n           | Order of the original matrix                                 |
| ml          | Number of diagonals below the main diagonal                  |
| mu          | Number of diagonals above the main diagonal                  |
| ipvt        | Vector containing the pivot indices from dgbfa               |
| b           | Right hand side vector                                       |
| info = 0    | Solve a*x = b                                                |
| info != 0   | Solve trans(a)*x = b                                         |

## Output parameters
--------------------
| Variable    | Description                                                  |
|:------------|:-------------------------------------------------------------|
| x           | Solution vector                                              |
"""
function dgbsl{T}(abd::Array{T,2}, lda::Int64, n::Int64, ml::Int64, mu::Int64,
    ipvt::Vector{Int64}, b::Vector{T}, flag::Int64)
    # Begin
    x = b[:]
    m = mu + ml + 1
    nm1 = n - 1
    if flag == 0
        # solve a*x = b
        # first solve l*y = b
        if ml != 0 && nm1 >= 1
            for k = 1:nm1
                lm = min(ml,n-k)
                l = ipvt[k]
                t = x[l]
                if l != k
                    x[l] = x[k]
                    x[k] = t
                end
                x[k+1:k+lm] += t*abd[m+1:m+lm,k:k]
            end
        end
        # now solve u*x = y
        for kb = 1:n
            k = n + 1 - kb
            x[k] = x[k]/abd[m,k]
            lm = min(k,m) - 1
            la = m - lm
            lb = k - lm
            t = -x[k]
            x[lb:lb+lm-1] += t*abd[la:la+lm-1,k:k]
        end
    else
        # flag = nonzero => solve transpose(a)*x = b
        # first solve transpose(u)*y = b
        for k = 1:n
            lm = min(k,m) - 1
            la = m - lm
            lb = k - lm
            t = abd[la:la+lm-1,k:k].*b[lb:lb+lm-1]
            x[k] = (x[k]-t)/abd[m,k]
        end
        # now solve transpose(l)*x = y
        if ml != 0 && nm1 >= 1
            for kb = 1:n
                k = n - kb
                lm = min(ml,n-k)
                x[k] += lm,abd[m+1:m+lm,k:k].*b[k+1:k+lm]
                l = ipvt[k]
                if l != k
                    t = x[l]
                    x[l] = x[k]
                    x[k] = t
                end
            end
        end
    end
    return x
end

"""
function n2prjn{T}(n::Int64, iRank::Int64, u::Vector{T}, d::Vector{T},
    qe::Array{T,2}, p::Vector{Int64}, v::Vector{T})

Provides the projection to the appropriate subspace in case
    of rank - reduction. To be used in connection with DECCON/SOLCON.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description                                      |
|:---------|:-------------------------------------------------|
| n        | Number of parameters to be estimated             |
| iRank    | Pseudo rank of decomposed Jacobian matrix        |
| u[n]     | Scaled Newton correction                         |
| d[n]     | Diagonal elements of upper triangular matrix     |
| qe[n]    | Part of pseudoinverse Jacobian matrix            |
| p[n]     | Pivot vector resulting from matrix decomposition |
| v[n]     | Temporary work array                             |

## Output parameters
--------------------
| Variable | Description                 |
|:---------|:----------------------------|
| del      | Defect                      |
"""
function n2prjn{T}(n::Int64, iRank::Int64, u::Vector{T}, d::Vector{T},
    qe::Array{T,2}, p::Vector{Int64}, v::Vector{T})
    # Begin
    v[:] = u[p]
    iRk1 = iRank + 1
    del = 0.0
    for i = iRk1:n
        sh = sum(qe[1:i-1,i:i].*v[1:i-1])
        s = (v[i] - sh)/d[i]
        del += s*s
        v[i] = s
    end
    return del
end

"""
function deccon{T}(a::Array{T,2}, nRow::Int64, nCol::Int64, mCon::Int64,
    m::Int64, n::Int64, iRankC::Int64, iRank::Int64, cond::T,
    d::Vector{T}, pivot::Vector{Int64}, kRed::Int64, ah::Array{T,2})

Constrained QR-decomposition of (m,n)-system  with computation of pseudoinverse
in case of rank-defeciency. First mcon rows belong to equality constraints.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable      | Description                                                   |
|:--------------|:--------------------------------------------------------------|
| a[nRow,nCol]* | Array holding the (m,n) matrix to be decomposed               |
| nRow          | Declared number of rows of array a                            |
| nCol          | Declared number of columns of array a and rows and columns of ah |
| mCon          | Number of equality constraints. mCon <= n. Internally reduced if equality constraints are linearly dependent |
| m             | Current number of rows of array a                             |
| n             | Current number of columns of array a                          |
| iRankC*       | Prescribed maximum pseudo-rank of constrained part of a. iRankC <= mCon |
| iRank*        | Prescribed maximum pseudo-rank of matrix a. iRank <= n        |
| cond*         | Permitted upper bound for the subcondition of least squares part of a |
| kRed >= 0     | Householder triangularization. Build up pseudo-inverse if iRank < n |
| kRed < 0      | Reduction of pseudo-rank of matrix a, skipping Householder triangularization, build-up new pseudo-inverse |

(* marks inout parameters)

## Output parameters
--------------------
| Variable      | Description                                                   |
|:--------------|:--------------------------------------------------------------|
| a[nRow,nCol]* | Upper triangular part consisting of transformed matrix and lower triangular part consisting of Householder transformations |
| iRankC*       | New pseudo-rank of constrained part of a such that abs(d[1]/d[iRankC]) < 1/epMach |
| iRank*        | New pseudo-rank of matrix a such that abs(d[iRankC+1]/d[iRank]) < cond |
| cond*         | The subcondition number belonging to the least squares part of A. In case of rank reduction: subcondition number which led to rank reduction |
| d             | Diagonal elements of the upper triangular matrix              |
| pivot         | Index vector storing permutation of columns due to pivoting   |
| ah[nCol,nCol] | In case of rank-defect used to compute the psuedo-inverse     |
| iFail = 0     | deccon computation was successful                             |
| iFail = -2    | Numerically negative diagonal element encountered during computation of pseudo-inverse due to extremely bad conditioned matrix a. deccon is unable to continue rank-reduction |
"""
function deccon{T}(a::Array{T,2}, nRow::Int64, nCol::Int64, mCon::Int64,
    m::Int64, n::Int64, iRankC::Int64, iRank::Int64, cond::T,
    d::Vector{T}, pivot::Vector{Int64}, kRed::Int64, ah::Array{T,2})
    # Begin
    # --------------------------------------------------------------------------
    # 1 Initialization
    level = 1
    k     = 1
    dd    = 1
    jj    = 1
    i1    = 1
    hMax  = zero(T)
    iFail = 0
    v     = zeros(T,n)
    t     = zero(T)

    if iRank > n
        iRank = n
    end
    if iRank > m
        iRank = m
    end
    # --------------------------------------------------------------------------
    # 1.1 Special case m = 1 and n = 1
    if m == 1 && n == 1
        pivot[1] = 1
        d[1] = a[1,1]
        cond = one(T)
        return (iRankC, iRank, cond, cond, iFail)
    end
    if kRed >= 0
        # ----------------------------------------------------------------------
        # 1.2 Initialize pivot-array
        pivot[1:n] = collect(1:n)
        # ----------------------------------------------------------------------
        # 2. Constrained Householder triangularization
        jd = 1
        iRanC1 = iRankC + 1
        mh = mCon
        iRankH = iRankC
        iData = 0
        if mh == 0
            iRankH = iRank
            mh = m
            iData = 1
        end
        iRk1 = iRank
        for k = 1:iRk1
            level = 1
            if k != n
                k1 = k+1
                qLoop = true
                while qLoop
                    if jd != 0
                        for j = k:n
                            d[j] = sum(a[k:mh,j:j].^2)
                        end
                    end
                    # ----------------------------------------------------------
                    # 2.1 Column pivoting
                    s1 = d[k]
                    jj = k
                    for l1 = k:n
                        if d[l1] > s1
                            s1 = d[l1]
                            jj = l1
                        end
                    end
                    h = d[jj]
                    if jd == 1
                        hMax = h/max(1.0e1,cond*0.05)
                    end
                    jd = 0
                    if h < hMax
                        jd = 1
                    end
                    qLoop = !(h >= hMax)
                end
                if jj != k
                    # ----------------------------------------------------------
                    # Column interchange
                    i = pivot[k]
                    pivot[k] = pivot[jj]
                    pivot[jj] = i
                    d[jj] = d[k]
                    for l1 = 1:m
                        s1 = a[l1,jj]
                        a[l1,jj] = a[l1,k]
                        a[l1,k] = s1
                    end
                end
            end
            # end for k != n case
            h = sum(a[k:mh,k:k].^2)
            t = sqrt(h)
            # ------------------------------------------------------------------
            # 2.3.0 A-priori test on pseudo-rank
            if k == 1 || k == iRanC1
                dd = t/cond
            end
            if t <= dd || k > iRankH
                # --------------------------------------------------------------
                # Rank reduction
                iRankH = k-1
                if mh != mCon || iData == 1
                    iRank = iRankH
                    if iRankC == iRank
                        level = 4
                    else
                        level = 3
                    end
                else
                    iRankC = iRankH
                    if iRankC != mCon
                        mh = m
                        iRankH = iRank
                        jd = 1
                        iData = 1
                        continue
                    else
                        error("Internal Error of deccon")
                    end
                end
            end
            if level == 1
                # --------------------------------------------------------------
                # 2.4 Householder step
                s = a[k,k]
                t = -sign(s)*abs(t)
                d[k] = t
                a[k,k] = s-t
                if k != n
                    t = 1.0/(h-s*t)
                    for j = k1:n
                        s = sum(a[k:mh,k:k].*a[k:mh,j:j])
                        s *= t
                        s1 = -s
                        if s != 0.0
                            # Update the sub columns
                            a[k:m,j:j] = a[k:m,j:j]+a[k:m,k:k]*s1
                        end
                        # update sub column norms
                        d[j] = d[j]-a[k,j]^2
                    end
                    if k == iRankC
                        mh = m
                        jd = 1
                        iRankH = iRank
                    end
                    if k == iRk1
                        level = 3
                    end
                else
                    level = 4
                end
                # end for Householder step
            end
            if level > 1
                break
            end
        end
    else
        k = -1
        level = 3
    end
    # --------------------------------------------------------------------------
    # 3 Rank-deficient psuedo-inverse
    if level == 3
        iRk1 = iRank+1
        for j = iRk1:n
            for ii = 1:iRank
                i = iRk1 - ii
                s = a[i,j]
                if ii != 1
                    sh = sum(a[i:i,i1:iRank]'.*v[i1:iRank])
                    s = s - sh
                end
                i1 = i
                v[i] = s/d[i]
                ah[i,j] = v[i]
            end
            for i = iRk1:j
                s = sum(ah[1:i-1,i:i].*v[1:i-1])
                if i != j
                    v[i] = -s/d[i]
                    ah[i,j] = -v[i]
                end
            end
            if s > -1.0
                d[j] = sqrt(s+1.0)
            else
                iFail = -2
                return (iRankC, iRank, cond, v[1], iFail)
            end
        end
    end
    # --------------------------------------------------------------------------
    # 9 Exit
    if iRankC != 0
        sh = d[iRankC]
        if sh != 0.0
            sh = abs(d[1]/sh)
        end
    else
        sh = zero(T)
    end
    v[1] = sh
    if k == iRank
        t = d[iRank]
    end
    if iRankC+1 <= iRank && t != 0.0
        s = abs(d[iRankC+1]/t)
    else
        s = zero(T)
    end
    cond = s
    iFail = 0
    return (iRankC, iRank, cond, v[1], iFail)
end

"""
function solcon{T}(a::Array{T,2}, nRow::Int64, nCol::Int64, mCon::Int64, m::Int64,
    n::Int64, x::Vector{T}, b::Vector{T}, iRankC::Int64, iRank::Int64,
    d::Vector{T}, pivot::Vector{Int64}, kRed::Int64, ah::Array{T,2})

Best constrained linear least squares solution of (m,n)-system. First mcon rows
comprise mcon equality constraints. To be used in connection with subroutine deccon.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description             |
|:---------|:------------------------|
| a[m,n], nRow, nCol, m, n, mCon, iRankC, iRank, d[n], pivot[n], ah[n,n], kRed | Refer to input and output parameters of deccon     |
| b[m]*    | Right hand side of the linear system if kRed >= 0. Right hand side of the upper linear system if kRed < 0 |

(* marks inout parameters)

## Output parameters
--------------------
| Variable | Description                        |
|:---------|:-----------------------------------|
| x[n]     | Best LSQ-solution of linear system |
| b[m]*    | Right-hand side of upper triangular system (transformed right-hand side of linear system) |
"""
function solcon{T}(a::Array{T,2}, nRow::Int64, nCol::Int64, mCon::Int64, m::Int64,
    n::Int64, x::Vector{T}, b::Vector{T}, iRankC::Int64, iRank::Int64,
    d::Vector{T}, pivot::Vector{Int64}, kRed::Int64, ah::Array{T,2})
    # Begin
    v = zeros(T,n)
    s = zero(T)
    # --------------------------------------------------------------------------
    # 1 Solution for pseudo-rank zero
    if iRank == 0
        x[1:n] = zeros(n)
        return (iRankC, iRank)
    end
    if iRank <= iRankC && iRank != n
        iRanC1 = iRankC + 1
    end
    if kRed >= 0 && (m != 1 || n != 1)
        # ----------------------------------------------------------------------
        # 2 Constrained householder transformations of right-hand side
        mh = mCon
        if iRankC == 0
            mh = m
        end
        for j = 1:iRank
            s = sum(a[j:mh,j:j].*b[j:mh])
            s = s/(d[j]*a[j,j])
            b[j:m] += a[j:m,j:j]*s
            if j == iRankC
                mh = m
            end
        end
    end
    # --------------------------------------------------------------------------
    # 3 Solution of upper triangular system
    iRk1 = iRank + 1
    for ii = 1:iRank
        i = iRk1 - ii
        i1 = i + 1
        s = b[i]
        if ii != 1
            sh = sum(a[i:i,i1:iRank]'.*v[i1:iRank])
            s = s-sh
        end
        v[i] = s/d[i]
    end
    if iRank != n && iRank != iRankC
        # ----------------------------------------------------------------------
        # 3.2 Computation of the best constrained least squares-solution
        for j = iRk1:n
            s = sum(ah[1:j-1,j:j].*v[1:j-1])
            v[j] = -s/d[j]
        end
        j1 = 1
        for jj = 1:n
            j = n-jj+1
            s = zero(T)
            if jj != 1
                s = sum(ah[j:j,j1:n]'.*v[j1:n])
            end
            if jj != 1 && j <= iRank
                v[j] -= s
            else
                j1 = j
                v[j] = -(s+v[j])/d[j]
            end
        end
    end
    # --------------------------------------------------------------------------
    # 4 Back-permutation of solution components
    x[pivot] = v[:]
    return iRank
end

"""
function nFact{T}(n::Int64, lda::Int64, ml::Int64, mu::Int64, a::Array{T,2},
    mStor::Int64, l::Array{T,2}, u::Array{T,2}, p::Vector{Int64})

Call linear algebra subprogram for factorization of a (n,n)-matrix.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable  | Description                                             |
|:----------|:--------------------------------------------------------|
| n         | Order of the linear system                              |
| lda       | Leading dimension of the matrix a                       |
| ml        | Lower bandwidth of the matrix (only for banded systems) |
| mu        | Upper bandwidth of the matrix (only for banded systems) |
| a[lda,n]  | Matrix to be factorized                                 |
| mStor = 0 | Full storage mode for matrix                            |
| mStor = 1 | Band storage mode for matrix                            |

## Output parameters
--------------------
| Variable | Description                                        |
|:---------|:---------------------------------------------------|
| iFail    | Exit code in case of errors                        |
| = 0      | Matrix decomposition successful                    |
| = 1      | Decomposition failed - matrix numerically singular |
| l[lda,n] | Lower triangular part of decomposed matrix in case of full mode. The complete LU-decomposition in band mode |
| u[lda,n] | Upper triangular part of decomposed matrix in case of full mode. Unused in case of band mode. |
| p        | Vector of pivot indices                            |
"""
function nFact{T}(n::Int64, lda::Int64, ml::Int64, mu::Int64, a::Array{T,2},
    mStor::Int64, l::Array{T,2}, u::Array{T,2}, p::Vector{Int64})
    # Begin
    iFail = 0
    if mStor == 0
        try
            (l[:,:],u[:,:],p[:]) = lu(a)
        catch
            iFail = 1
        end
    elseif mStor == 1
        # Band mode: l holds the complete lu-factorization of p*a
        l[:,:] = a
        # Use LINPACK function to compute the LU in place
        # p stores the pivot vectors and not the permuted identity
        (l[:,:],p[:],iFail) = dgbfa(a,lda,n,ml,mu)
        if iFail != 0
            iFail = 1
        end
    end
    return iFail
end

"""
function nFact{T}(n::Int64, lda::Int64, ldaInv::Int64, ml::Int64, mu::Int64,
    a::Array{T,2}, aInv::Array{T,2}, cond::T, iRank::Int64,
    opt::OptionsNLEQ, p::Vector{Int64}, d::Vector{T}, iRepeat::Int64, iRankC::Int64)

Call linear algebra subprogram for factorization of a (n,n)-matrix
with rank decision and casual computation of the rank deficient
pseudo-inverse matrix.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable  | Description                                             |
|:----------|:--------------------------------------------------------|
| n         | Order of the linear system                              |
| lda       | Leading dimension of the matrix a                       |
| ldaInv    | Leading dimension of the matrix aInv                    |
| ml        | Lower bandwidth of the matrix (only for banded systems) |
| mu        | Upper bandwidth of the matrix (only for banded systems) |
| a[lda,n]  | Matrix to be factorized                                 |
| cond*     | Maximum permitted subcondition for the prescribed rank  |
| opt       | User prescribed options                                 |
| iRankC, iRank, cond | Refer to input and output parameters of deccon |
| iRank*    | Prescribed maximum pseudo-rank of matrix a              |
| cond*     | Permitted upper bound for the subcondition of leastsquares part of a |

## Output parameters
--------------------
| Variable       | Description                                        |
|:---------------|:---------------------------------------------------|
| aInv[ldaInv,n] | If matrix a is rank deficient this array holds the pseudo-inverse of a |
| iFail = 0      | Matrix decomposition successful                    |
| p              | Vector of pivot indices                            |
| d              | Refer to deccon                                    |
"""
function nFact{T}(n::Int64, lda::Int64, ldaInv::Int64, ml::Int64, mu::Int64,
    a::Array{T,2}, aInv::Array{T,2}, cond::T, iRank::Int64,
    opt::OptionsNLEQ, p::Vector{Int64}, d::Vector{T}, iRepeat::Int64, iRankC::Int64)
    # Begin
    mPrWarn = opt.options[OPT_PRINTWARNING]
    printIO = opt.options[OPT_PRINTIOWARN]
    mCon = 0
    iRepeat = -iRepeat
    if iRepeat == 0
        iRankC = mCon
    end
    (iRankC, iRank, cond, subCond, iFail) = deccon(a, lda, n, mCon, n, n,
                                            iRankC, iRank, cond, d, p,
                                            iRepeat, aInv)
    if iFail == -2 && mPrWarn == 1
        write(printIO,"\n","Deccon failed to compute Rank-deficient QR-Decomposition\n")
    end
    if iRank != 0
        cond = abs(d[1]/d[iRank])
        setOption!(wkNLEQ2, WK_SENS1, abs(d[1]))
    else
        cond = one(T)
        setOption!(wkNLEQ2, WK_SENS1, zero(T))
    end
    return (cond, iRankC, iFail)
end

"""
function nSolv{T}(n::Int64, lda::Int64, ml::Int64, mu::Int64, l::Array{T,2},
    u::Array{T,2}, p::Vector{Int64}, b::Vector{T}, mStor::Int64)

Call linear algebra subprogram for solution of the linear system a*z = b

T = Float64 or BigFloat

## Parameters
-------------------
| Variable | Description             |
|:---------|:------------------------|
| n, lda, ml, mu, l, u, p, b, mStor, iFail | Refer nFact |
"""
function nSolv{T}(n::Int64, lda::Int64, ml::Int64, mu::Int64, l::Array{T,2},
    u::Array{T,2}, p::Vector{Int64}, b::Vector{T}, mStor::Int64)
    # Begin
    if mStor == 0
        b[:] = b[p]
        b[:] = l\b
        b[:] = u\b
    elseif mStor == 1
        b[:] = dgbsl(l,lda,n,ml,mu,p,b,0)
    end
    iFail = 0
    return iFail
end

"""
function nSolv{T}(n::Int64, lda::Int64, ldaInv::Int64, ml::Int64, mu::Int64,
    a::Array{T,2}, aInv::Array{T,2}, b::Vector{T},
    z::Vector{T}, iRank::Int64, iRepeat::Int64, d::Vector{T},
    pivot::Vector{Int64}, iRankC::Int64)

Call linear algebra subprogram for solution of the linear system a*z = b.
To be used with the factorization routine of a (n,n)-matrix
with rank decision and casual computation of the rank deficient
pseudo-inverse matrix.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description                          |
|:---------|:-------------------------------------|
| b[n]*    | Right hand side of the linear system |
| n, lda, ldaInv, ml, mu, a, aInv, iRank, iRepeat, d, pivot, iRankC | Refer nFact |

(* marks inout parameters)

## Output parameters
--------------------
| Variable | Description                                                       |
|:---------|:------------------------------------------------------------------|
| b[n]*    | RHS transformed to the upper triangular part of the linear system |
"""
function nSolv{T}(n::Int64, lda::Int64, ldaInv::Int64, ml::Int64, mu::Int64,
    a::Array{T,2}, aInv::Array{T,2}, b::Vector{T},
    z::Vector{T}, iRank::Int64, iRepeat::Int64, d::Vector{T},
    pivot::Vector{Int64}, iRankC::Int64)
    # Begin
    mCon = 0
    iRepeat = -iRepeat
    iRank = solcon(a,lda,n,mCon,n,n,z,b,iRankC,iRank,d,pivot,
                iRepeat,aInv)
    iFail = 0
    return iFail
end


text:
"""
    `sort_data(data, file_type)`

Sorts data from wopwop text file outputs.

# ARGUMENTS
* `data::Array{Real, Real}`     : Name of the file to be read, path included
* `file_type::String`           : Specifies the type of output file from WOPWOP. Options:
                                    "pressure-time"
                                    "spl_spectrum"
                                    "oaspl"
"""
function sort_data(data, file_type)

    if file_type == "pressure-time"

        time = data[:,1]
        thickness = data[:,2]
        loading = data[:,3]
        total = data[:,4]
        thicknessX = data[:,5]
        thicknessY = data[:,6]
        thicknessZ = data[:,7]
        loadingX = data[:,8]
        loadingY = data[:,9]
        loadingZ = data[:,10]
        totalX = data[:,11]
        totalY = data[:,12]
        totalZ = data[:,13]

        return time, thickness, loading, total, thicknessX,
                thicknessY, thicknessZ, loadingX, loadingY,
                loadingZ, totalX, totalY, totalZ

    elseif file_type == "spl-spectrum"

        frequency = data[:,1]
        thickness_dB = data[:,2]
        loading_dB = data[:,3]
        total_dB = data[:,4]
        thickness_dBA = data[:,5]
        loading_dBA = data[:,6]
        total_dBA = data[:,7]

        return frequency, thickness_dB, loading_dB, total_dB,
                thickness_dBA, loading_dBA, total_dBA

    elseif file_type == "oaspl" || file_type == "oaspla"

        time, thickness_OASPL, loading_OASPL, total_OASPL =
            Float64[], Float64[], Float64[], Float64[]

        time = data[1]
        thickness_OASPL = data[2]
        loading_OASPL = data[3]
        total_OASPL = data[4]

        return time, thickness_OASPL, loading_OASPL, total_OASPL

    else
        println("ERROR: Invalid file_type input")
    end
end

"""
    `prepare_test(filepath, windowtype)`

Reads all output files from wopwop test cases and calculates the same results from the post_processing
    functions in FLOWNoise. Returns the needed variables to perform the required tests.

# ARGUMENTS
* `filepath::String`            : Name of the file to be read, path included
* `windowtype::String`          : Data windowing options:
                                    "None", "Hann"
"""
function prepare_test(filepath, windowtype)

    #Read time-pressure data
    filename = joinpath(filepath, "pressure")
    file_type = "pressure-time"
    header, data = FLOWNoise.read_wopwopoutput(filename, read_path="", verbose=true, v_lvl=0, tec=true)
    time, thickness, loading, total,
            thicknessX, thicknessY, thicknessZ,
            loadingX, loadingY, loadingZ,
            totalX, totalY, totalZ = sort_data(data, file_type)

    #Read Pressure Spectrum data
    filename = joinpath(filepath, "spl_spectrum")
    file_type = "spl-spectrum"
    header, data = FLOWNoise.read_wopwopoutput(filename, read_path="", verbose=true, v_lvl=0, tec=true)
    frequency, thickness_dB_actual, loading_dB_actual, total_dB_actual,
            thickness_dBA_actual, loading_dBA_actual, total_dBA_actual = sort_data(data, file_type)

    #Read non-weighted OASPL data
    filename = joinpath(filepath, "OASPLdB")
    file_type = "oaspl"
    header, data = FLOWNoise.read_wopwopoutput(filename, read_path="", verbose=true, v_lvl=0, tec=true)
    time_2, thickness_OASPL_actual, loading_OASPL_actual, total_OASPL_actual= sort_data(data, file_type)

    #Read A-weighted OASPL data
    filename = joinpath(filepath, "OASPLdBA")
    file_type = "oaspla"
    header, data = FLOWNoise.read_wopwopoutput(filename, read_path="", verbose=true, v_lvl=0, tec=true)
    time_3, thickness_OASPLA_actual, loading_OASPLA_actual, total_OASPLA_actual = sort_data(data, file_type)

    # Convert from time to frequency domain
    thickness_dB_calc = FLOWNoise.pressure2SPL(time, thickness)
    thickness_OASPL_calc = FLOWNoise.SPL2OASPL(thickness_dB_calc)
    loading_dB_calc = FLOWNoise.pressure2SPL(time, loading)
    loading_OASPL_calc = FLOWNoise.SPL2OASPL(loading_dB_calc)
    total_dB_calc = FLOWNoise.pressure2SPL(time, total)
    total_OASPL_calc = FLOWNoise.SPL2OASPL(total_dB_calc)

    thickness_dBA_calc = FLOWNoise.pressure2SPL(time, thickness, AweightingFlag=true)
    thickness_OASPLA_calc = FLOWNoise.SPL2OASPL(thickness_dBA_calc)
    loading_dBA_calc = FLOWNoise.pressure2SPL(time, loading, AweightingFlag=true)
    loading_OASPLA_calc = FLOWNoise.SPL2OASPL(loading_dBA_calc)
    total_dBA_calc = FLOWNoise.pressure2SPL(time, total, AweightingFlag=true)
    total_OASPLA_calc = FLOWNoise.SPL2OASPL(total_dBA_calc)

    println("Actual thickness OASPL = ", thickness_OASPL_actual)
    println("Calculated thickness OASPL = ", thickness_OASPL_calc)

    println("Actual loading OASPL = ", loading_OASPL_actual)
    println("Calculated loading OASPL = ", loading_OASPL_calc)

    println("Actual total OASPL = ", total_OASPL_actual)
    println("Calculated total OASPL = ", total_OASPL_calc)

    println("Actual thickness OASPLA = ", thickness_OASPLA_actual)
    println("Calculated thickness OASPLA = ", thickness_OASPLA_calc)

    println("Actual loading OASPLA = ", loading_OASPLA_actual)
    println("Calculated loading OASPLA = ", loading_OASPLA_calc)

    println("Actual total OASPLA = ", total_OASPLA_actual)
    println("Calculated total OASPLA = ", total_OASPLA_calc)

    println("")

    return thickness_dB_calc, loading_dB_calc, total_dB_calc,
            thickness_dBA_calc, loading_dBA_calc, total_dBA_calc,
            thickness_dB_actual, loading_dB_actual, total_dB_actual,
            thickness_dBA_actual, loading_dBA_actual, total_dBA_actual,
            thickness_OASPL_calc, loading_OASPL_calc, total_OASPL_calc,
            thickness_OASPLA_calc, loading_OASPLA_calc, total_OASPLA_calc,
            thickness_OASPL_actual, loading_OASPL_actual, total_OASPL_actual,
            thickness_OASPLA_actual, loading_OASPLA_actual, total_OASPLA_actual

end

"""
    `percent_error(x1, x2)`

Simple percent error calculation to clean up code.

"""
function percent_error(x1, x2)

    return abs((x1-x2)/x1)

end

ext: jl

max_stars_repo_name: byuflowlab/FLOWNoise

max_issues_repo_name: byuflowlab/FLOWNoise

max_forks_repo_name: byuflowlab/FLOWNoise


text:
using StochasticDiffEq, Test, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_linear, prob_sde_2Dlinear
Random.seed!(100)
prob = prob_sde_linear
f(u,p,t) = 2u
prob = SDEProblem{false}(f,prob.f,prob.u0,prob.tspan)
sol = solve(prob,EM(),dt=1//2^(4),saveat = 0.33)

sol.t == [0.0,.33,.33+.33,1]

sol = solve(prob,EM(),dt=1//2^(4),saveat = 0.33,save_start=false)

sol.t == [.33,.33+.33,1]
sol(0.4)

sol = solve(prob,EM(),dt=1//2^(4),saveat = 0.33,save_start=false,save_everystep=true)

@test 0.33 ∈ sol.t
@test 0.66 ∈ sol.t
sol.t != [.33,.33+.33,1]

prob = prob_sde_2Dlinear
f(du,u,p,t) = prob.f(du,u,p,t)
prob2 = SDEProblem(f,prob.g,vec(prob.u0),prob.tspan)
Random.seed!(200)
sol = solve(prob2,EM(),dt=1//2^(4),saveat = 0.33,save_start=false,save_everystep=true)
Random.seed!(200)
sol2= solve(prob2,EM(),dt=1//2^(4),saveat = 0.33,save_start=false,save_everystep=true,save_idxs=1:2:5)

@test sol(0.43)[1:2:5] == sol2(0.43)

# Test save_on switch
sol = solve(prob2, EM(), dt=1//2^(4), save_on=false, save_start=false, save_end=false)
@test isempty(sol.t) && isempty(sol.u)
sol = solve(prob2, EM(), dt=1//2^(4), saveat=0.33, save_on=false, save_start=false, save_end=false)
@test isempty(sol.t) && isempty(sol.u)

_saveat = [0.0,0.25,0.5,1.0]
integ = init(SDEProblem((u,p,t)->u,(u,p,t)->u,0.0,(0.0,1.0)),SOSRI(),saveat=_saveat)
add_tstop!(integ,2.0)
solve!(integ)
@test integ.sol.t == _saveat

integ = init(SDEProblem((u,p,t)->u,(u,p,t)->u,0.0,(0.0,1.0)),SOSRI(),saveat=_saveat,save_end = true)
add_tstop!(integ,2.0)
solve!(integ)
@test integ.sol.t == [0.0,0.25,0.5,1.0,2.0]

integ = init(SDEProblem((u,p,t)->u,(u,p,t)->u,0.0,(0.0,1.0)),SOSRI(),saveat=_saveat,save_end = false)
add_tstop!(integ,2.0)
solve!(integ)
@test integ.sol.t == _saveat

# Catch save for maxiters
ode = SDEProblem((u,p,t) -> u, (u,p,t) -> u, 1.0, (0.0, 1.0))
sol = solve(ode, SOSRI(), save_everystep=false) # okay, as expected
@test length(sol) == 2
@info "Warning Expected"
sol = solve(ode, SOSRI(), save_everystep=false, maxiters=3) # doesn't save the final solution anymore!
@test length(sol) == 2


text:
using Test
using GraphIO.NET

@testset "Pajek NET" begin
    for g in values(allgraphs)
        readback_test(NETFormat(), g)
    end
    fname = joinpath(testdir, "testdata", "kinship.net")
    @test length(loadgraphs(fname, NETFormat())) == 1
end


text:
#MDDatasets function tools
#-------------------------------------------------------------------------------

#==Basic tools
===============================================================================#

#Get the value of a particular keyword in the list of keyword arguments:
getkwarg(kwargs::Base.Iterators.Pairs, s::Symbol) = get(kwargs, s, nothing)

#==Ensure interface (similar to assert)
===============================================================================#
#=Similar to assert.  However, unlike assert, "ensure" is not meant for
debugging.  Thus, ensure is never meant to be compiled out.
=#
function ensure(cond::Bool, err)
	if !cond; throw(err); end
end

#Conditionnally generate error using "do" syntax:
function ensure(fn::Function, cond::Bool)
	if !cond; throw(fn()); end
end

#Last line


text:
using FEBOL
using Base.Test

# SearchDomain tests
m = SearchDomain(400, 170, 210)
@test m.theta == (170.0, 210.0, 0.0,0.0)
m = SearchDomain(400)
theta!(m)
theta!(m, 32, 40)

# BearingOnly
s = BearingOnly()
@test s.noise_sigma == 10.0
s = BearingOnly(5)
@test s.noise_sigma == 5.0

# bearing only simulation
m = SearchDomain(200)
s = FOV()
f = DF(m, 41, s, 0:1)
x = Vehicle(100, 100, 0, 5, s)
p = GreedyPolicy(x, 4)
simulate(m, x, f, p)
simulate(m, SimUnit(x,f,p))

# fov simulation
m = SearchDomain(200)
s = FOV()
f = DF(m, 41, s, 0:10:35)
x = Vehicle(100, 100, 0, 5, s)
p = GreedyPolicy(x, 4)
simulate(m, x, f, p)

# simunit way
tc = StepThreshold(10)
su = SimUnit(x, f, p, tc)
simulate(m, su)

# batch sims
simulate(m, su, 10)


text:
# Useful for LinearAlgebra tests
function generate_well_conditioned_matrix(rng, N)
    A = randn(rng, N, N)
    return A * A' + I
end

generate_well_conditioned_matrix(N) = generate_well_conditioned_matrix(Random.GLOBAL_RNG, N)


text:
"""
function addEventDeath_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions for division events.
"""
function addEventDeath_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

if "EventDeath" in keys(abm.declaredUpdates)

    condition = abm.declaredUpdates["EventDeath"].args[end]

    #Create function to check dead elements
    if platform == "cpu"
        up= :(removeList_Pos_ = Threads.atomic_add!(remV_,1) + 1)
    elseif platform == "gpu"
        up= :(removeList_Pos_ = CUDA.atomic_add!(pointer(remV_,1),Int(1)) + 1)
    end

    code = quote
        if $condition
            keepList_[ic1_] = 0
            $up
            removeList_[removeList_Pos_] = ic1_
        else
            keepList_[ic1_] = ic1_
        end
    end

    code = vectorize_(abm,code,p)

    f1 = simpleFirstLoopWrapInFunction_(platform,:checkDyingAgents_!,code)

    #Create function to put agents in their new position
    code = quote
        oldPos_ = keepList_[ic1_]
        newPos_ = removeList_[ic1_]
    end
    if !isempty(abm.declaredSymbols["Local"])
        push!(code.args,
            :(begin
                if oldPos_ > 0
                    for ic2_ in 1:$(length(abm.declaredSymbols["Local"]))
                        localV[newPos_, ic2_] = localV[oldPos_, ic2_]
                    end
                end
            end)
        )
    end
    if !isempty(abm.declaredSymbols["Identity"])
        push!(code.args,
            :(begin
                if oldPos_ > 0
                    for ic2_ in 1:$(length(abm.declaredSymbols["Identity"]))
                        identityV[newPos_, ic2_] = identityV[oldPos_, ic2_]
                    end
                end
            end)
        )
    end
    if !isempty(p.update["Local"])
        push!(code.args,
            :(begin
                if oldPos_ > 0
                    for ic2_ in 1:$(length(p.update["Local"]))
                        localVCopy[newPos_, ic2_] = localVCopy[oldPos_, ic2_]
                    end
                end
            end)
        )
    end
    if !isempty(p.update["Identity"])
        push!(code.args,
            :(begin
                if oldPos_ > 0
                    for ic2_ in 1:$(length(p.update["Identity"]))
                        identityVCopy[newPos_, ic2_] = identityVCopy[oldPos_, ic2_]
                    end
                end
            end)
        )
    end

    f2 = simpleFirstLoopWrapInFunction_(platform,:assignDeadAgentsSpaces_!,code)
    if platform == "cpu"
        f2 = postwalk(x->@capture(x,N) ? :(remV_[]) : x, f2)
    elseif platform == "gpu"
        f2 = postwalk(x->@capture(x,N) ? :(remV_[1]) : x, f2)
    end

    #Make wrap function of the algorithm
    if platform == "cpu"
        code = quote
            remV_[] = 0
            if N > 0
                @platformAdapt checkDyingAgents_!(ARGS_)
                sort!(keepList_,rev=true)
                @platformAdapt assignDeadAgentsSpaces_!(ARGS_)
            end
        end
    elseif platform == "gpu"
        code = quote
            remV_ .= 0
            if N > 0
                @platformAdapt checkDyingAgents_!(ARGS_)
                sort!(keepList_,rev=true)
                @platformAdapt assignDeadAgentsSpaces_!(ARGS_)
            end
        end
    end

    f3 = wrapInFunction_(:removeDeadAgents_!,code)

    # Add to program
    if platform == "cpu"
        push!(p.declareVar.args,
            :(begin
                remV_ = Threads.Atomic{Int}()
                keepList_ = zeros(Int,nMax)
                removeList_ = zeros(Int,nMax)
            end)
        )
    elseif platform == "gpu"
        push!(p.declareVar.args,
            :(begin
                remV_ = CUDA.zeros(Int,1)
                keepList_ = zeros(Int,nMax)
                removeList_ = zeros(Int,nMax)
            end)
        )
    end

    push!(p.args,:remV_,:keepList_,:removeList_)

    push!(p.declareF.args,f1,f2,f3)

    if platform == "cpu"
        push!(p.execInloop.args,:(removeDeadAgents_!(ARGS_); N -= remV_[]))
    else
        push!(p.execInloop.args,:(removeDeadAgents_!(ARGS_); N -= Core.Array(remV_)[1]))
    end
end

return nothing
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
# This is a crude replication of Carvalho and Vilela (2015)
## By Daniel Coutinho and Rafaela

####### WARNING############
## NOT DEBUGGED AS MUCH AS IT SHOULD!
#################

### TODO

## 1. Pi and Psi matrices
## 2. More comments
## 3. Coefficients as strings for latexify

#######

### ORDER
## c_t
##c_{t-1}
##i_{t-1}
##i_{t-2}
## pi_t
## gamma_t
## y_t
## tot_{t-1}
## y^{star}_t
## q_{t-1}
## s_t
## s_{t-1}
## pi^{star}_t
## pi_{D,t}
## pi_{D,t-1}
## pi_{I,t}
## pi_{I,t-1}
## a_t
## z_t
## i^{star}_t
## varepsilon_{cp,t}
## phi_t

## Order of the expectational errors

#c_t
#pi_t
#gamma_t
#pi^d_t
#pi^I_t
#s_t

## Order of fundamental shocks
#
#epsilon_a
#epsilon_gamma
#epsilon_cp
#epsion_phi
#epsilon_i
#epsilon_y*
#epsilon_pi*
#epsilon_i*

#Creating Matrices

G0 = zeros(23,23)
G1 = zeros(23,23)
Pi = zeros(23,6)
Psi = zeros(23,8)

kappa_d = (1-theta_D)*(1-theta_D*bet)/theta_D
kappa_i = (1-theta_I)*(1-theta_I*bet)/theta_I

#First equation

G0[1,1] = 1
G0[1,5] = (1-h)/sig
G0[1,6] = (1-h)/sig
G0[1,3] = -(1-h)/sig

G1[1,1] = (1-h)
G1[1,2] = h
G1[1,6] = (1-h)/sig

Pi[1,1] = -1
Pi[1,2] = -(1-h)/sig
Pi[1,3] = (1-h)/sig

#Second Equation

G1[2,7] = -1
G1[2,1] = 1-alfa
G1[2,8] = alfa*eta*(2-alfa)
G1[2,9] = alfa*eta
G1[2,10] = alfa

#Third Equation

G0[3,8] = 1

G1[3,8] = 1
G1[3,15] = -1
G1[3,17] = 1

#Fourth Equation

G1[4,11] = -1
G1[4,9] = 1
G1[4,8] = 1 -alfa

#Fifth equation

G0[5,11] = -1
G0[5,14] = -1

G1[5,11] = 1
G1[5,12] = 1
G1[5,13] = -1
G1[5,5] = -1

#Equation number 6

G0[6,15] = bet

G1[6,15] = 1+delta_D
G1[6,16] = bet*delta_D
G1[6,8] = -kappa_d*alfa
G1[6,7] = -kappa_d*phi
G1[6,1] = -kappa_d*sig/(1-h)
G1[6,2] = kappa_d*sig/(1-h)
G1[6,19] = kappa_d*(1+phi)

Pi[6,4] = -bet

#Equation number 7

G0[7,17] = bet

G1[7,17] = (1+bet*delta_I)
G1[7,18] = delta_I
G1[7,9] = kappa_i

Pi[7,5] = -bet

#Equation number 8

G0[8,8] = -alfa

G1[8,5] = -1
G1[8,15] = 1
G1[8,8] = alfa

#Equation number 9

G0[9,20] = 1
G0[9,8] = alfa

G1[9,20] = 1/bet
G1[9,9] = -alfa
G1[9,7] = 1
G1[9,1] = -1

# Equation number 10

G0[10,12] = 1
G0[10,21] = 1
G0[10,3] = -1

G1[10,20] = chi

Pi[10,6] = -1

# Equation number 11: Taylor rule

G0[11,3] = 1
G1[11,3] = rho_1
G1[11,4] = rho_2
G1[11,5] = (1-rho_1-rho_2)*lambda_pi
G1[11,7] = (1-rho_1-rho_2)*lambda_y
G1[11,12] = (1-rho_1-rho_2)*lambda_s
G1[11,13] = -(1-rho_1-rho_2)*lambda_s ## its for the delta tot

Psi[11,5] = sigma_i

###################
## Foreign Block ##
###################

# Equation 12

G0[12,10] = 1

G1[12,10] = a_1yy
G1[12,14] = a_1ypi
G1[12,21] = a_1yi

Psi[12,6] = sigma_y_star

#Equation 13

G0[13,14] = 1
G0[13,10] = a_0piy

G1[13,10] = a_1piy
G1[13,14] = a_1pipi
G1[13,21] = a_1pii

Psi[13,7] = sigma_pi_star

#Equation 14

G0[14,21] = 1
G0[14,14] = a_0ipi
G0[14,10] = a_0iy

G1[14,10] = a_1iy
G1[14,14] = a_1ipi
G1[14,21] = a_1ii

Psi[14,8] = sigma_i_star

#################
## Shocks block #
#################

# Equation 15

G0[15,19] = 1

G1[15,19] = rho_a

Psi[15,1] = sigma_a

# Equation 16

G0[16,6] = 1

G1[16,6] = rho_gamma

Psi[16,2] = sigma_gamma

# Equation 17

G0[17,22] = 1

G1[17,22] = rho_cp

Psi[17,3] = sigma_cp

# Equation 18

G0[18,23] = 1

G1[18,23] = rho_phi

Psi[18,4] = sigma_phi

####################
## Identity block ##
####################

# Equation 19

G0[19,1] = 1
G1[19,2] = 1

# Eqquation 20

G0[20,3] = 1
G1[20,4] = 1

# Equation 21

G0[21,12] = 1
G1[21,13] = 1

# Equation 22

G0[22,15] = 1
G1[22,16] = 1

# Equation 23

G0[23,17] = 1
G1[23,18] = 1

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using VariationalInequality
using JuMP, Ipopt

using Base.Test
using FactCheck

# include("../src/model.jl")
# include("../src/algorithms.jl")

# [21] L. Mathiesen, An algorithm based on a sequence of linear complementarity problems applied
# to a Walrasian equilibrium model: An example, Math. Programming, 37 (1987), pp. 1–18.

m = VIPModel()

@variable(m, x1>=0)
@variable(m, x2>=0)
@variable(m, x3>=0)

@NLconstraint(m, x1+x2+x3 <= 1)
@NLconstraint(m, x1-x2-x3 <= 0)

@NLexpression(m, F1, 0.9* (5*x2 + 3*x3)/x1 )
@NLexpression(m, F2, 0.1* (5*x2 +3*x3) /x2   -  5)
@NLexpression(m, F3, -3)

correspond(m, [F1, F2, F3], [x1, x2, x3])

solveVIP(m, algorithm=:fixed_point, max_iter=1000, step_size=0.01)
sol0, Fval0, gap0 = saveSolution(m)

solveVIP(m, algorithm=:extra_gradient, max_iter=1000, step_size=0.01)
sol1, Fval1, gap1 = saveSolution(m)

println(sol0)
println(sol1)

println(gap0)
println(gap1)

# (1, 0, 0)


text:
# test print with nicely printed diff on failures (from
# https://github.com/invenia/PkgTemplates.jl/blob/master/test/runtests.jl): access point
# is `test_show`.

using DeepDiffs: deepdiff
function print_diff(a, b)
    old = Base.have_color
    @eval Base have_color = true
    try
        println(deepdiff(a, b))
    finally
        @eval Base have_color = $old
    end
end

function test_text(expected::AbstractString, observed::AbstractString)
    if expected == observed
        @test true
    else
        print_diff(expected, observed)
        @test expected == observed
    end
end

function test_show(v, observed::AbstractString; mime=MIME"text/plain"())
    test_text(repr(mime, v), observed)
end


text:
module BigFloatX

using Compat

export BigFloatA, BigFloatB

import Base.MPFR: ROUNDING_MODE

import Base.GMP: ClongMax, CulongMax, CdoubleMax, Limb

import Base: convert, cconvert, unsafe_convert, string, precision,
    +, *, cmp, zero

type Mpfr_t
    prec::Clong
    sign::Cint
    exp::Clong
    d::Ptr{Limb}
end

include("a.jl")
include("b.jl")

end # module


text:
using CLIMA.ArtifactWrappers
using NCDatasets
using Test

# Get soundings dataset folder:
soundings_dataset = ArtifactWrapper(
    joinpath(@__DIR__, "Artifacts.toml"),
    "soundings",
    ArtifactFile[
        ArtifactFile(
            url = "https://caltech.box.com/shared/static/rjnvt2dlw7etm1c7mmdfrkw5gnfds5lx.nc",
            filename = "sounding_gabersek.nc",
        ),
        ArtifactFile(
            url = "https://caltech.box.com/shared/static/th1a0p5bduej2hx6lisphlpo7pm2wvys.nc",
            filename = "sounding_WKR88.nc",
        ),
        ArtifactFile(
            url = "https://caltech.box.com/shared/static/etz2f6qn45m5x4ow1nwx2mcams9fno9c.nc",
            filename = "sounding_JCP2013_with_pressure.nc",
        ),
        ArtifactFile(
            url = "https://caltech.box.com/shared/static/90um61xxjtxh45cwg2d2t5mucl8cbphe.nc",
            filename = "sounding_JCP2013.nc",
        ),
        ArtifactFile(
            url = "https://caltech.box.com/shared/static/satr1neltu8pedvukn125ywr5eks8m0r.nc",
            filename = "sounding_GC1991.nc",
        ),
        ArtifactFile(
            url = "https://caltech.box.com/shared/static/3kmfcqrlf9scx9166c2fx0g1c36i124l.nc",
            filename = "sounding_blend.nc",
        ),
    ],
)
data_folder = get_data_folder(soundings_dataset)

@testset "Get soundings data" begin
    @test isfile(joinpath(data_folder, "sounding_gabersek.nc"))
    @test isfile(joinpath(data_folder, "sounding_GC1991.nc"))
    @test isfile(joinpath(data_folder, "sounding_WKR88.nc"))
    @test isfile(joinpath(data_folder, "sounding_JCP2013.nc"))
    @test isfile(joinpath(data_folder, "sounding_blend.nc"))
    @test isfile(joinpath(data_folder, "sounding_JCP2013_with_pressure.nc"))
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
module ReplMaker

import REPL
import REPL.LineEdit
import REPL.LineEditREPL
import Base.display

export initrepl, enter_mode!

"""
```
         initrepl(parser::Function;
                  prompt_text = "myrepl> ",
                  prompt_color = :blue,
                  start_key = ')',
                  repl = Base.active_repl,
                  mode_name = :mylang,
                  show_function = nothing,
                  show_function_io = stdout,
                  valid_input_checker::Function = (s -> true),
                  keymap::Dict = REPL.LineEdit.default_keymap_dict,
                  completion_provider = REPL.REPLCompletionProvider(),
                  sticky = true,
                  startup_text=true
                  )
```
creates a custom repl mode which takes in code and parses it according to whatever parsing function you
provide in the argument `parser`. Choose which key initializes the repl mode with `start_key`, the name of
your repl mode with `mode_name` and optionally provide a function which checks if a given repl input is valid
before parsing with `valid_input_checker`. Autocompletion options are supplied through the argument
`completion_provider` which defaults to the standard julia REPL TAB completions.
"""
function initrepl(parser::Function;
                  prompt_text = "myrepl> ",
                  prompt_color = :blue,
                  start_key = ')',
                  repl = Base.active_repl,
                  mode_name = :mylang,
                  show_function = nothing,
                  show_function_io = stdout,
                  valid_input_checker::Function = (s -> true),
                  keymap::Dict = REPL.LineEdit.default_keymap_dict,
                  completion_provider = REPL.REPLCompletionProvider(),
                  sticky_mode = true,
                  startup_text=true
                  )

    color = Base.text_colors[prompt_color]

    if show_function != nothing
        repl = enablecustomdisplay(repl, show_function, show_function_io)
    end

    julia_mode = repl.interface.modes[1]
    prefix = repl.hascolor ? color : ""
    suffix = repl.hascolor ? (repl.envcolors ? Base.input_color : repl.input_color()) : ""

    lang_mode = LineEdit.Prompt(prompt_text;
                                prompt_prefix    = prefix,
                                prompt_suffix    = suffix,
                                keymap_dict      = keymap,
                                on_enter         = valid_input_checker,
                                complete         = completion_provider,
                                sticky           = sticky_mode
                                )

    lang_mode.on_done = REPL.respond(parser, repl, lang_mode)

    push!(repl.interface.modes, lang_mode)

    hp                     = julia_mode.hist
    hp.mode_mapping[mode_name |> Symbol] = lang_mode
    lang_mode.hist         = hp

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)

    mk = REPL.mode_keymap(julia_mode)

    if start_key in keys(julia_mode.keymap_dict)
        @warn "REPL key '$start_key' overwritten."
        alt = deepcopy(julia_mode.keymap_dict[start_key])
    else
        alt = (s, args...) -> LineEdit.edit_insert(s, start_key)
    end

    lang_keymap = Dict{Any,Any}(
    start_key => (s, args...) ->
      if isempty(s) || position(LineEdit.buffer(s)) == 0
        enter_mode!(s, lang_mode)
      else
        alt(s, args...)
      end
    )

    lang_mode.keymap_dict = LineEdit.keymap(Dict{Any,Any}[
        skeymap,
        mk,
        LineEdit.history_keymap,
        LineEdit.default_keymap,
        LineEdit.escape_defaults,
    ])

    julia_mode.keymap_dict = LineEdit.keymap_merge(julia_mode.keymap_dict, lang_keymap)

    startup_text && println("REPL mode $mode_name initialized. Press $start_key to enter and backspace to exit.")

  lang_mode
end

"""
an extension of `REPL.LineEdit.transition` to achieve handy use.

When
- giving 3 arguments(action, mistate, mode),
  it equals to `REPL.LineEdit.transition`.

- giving 2 arguments(mistate, mode),
  there is a default 'action' copying the
  current buffer to the target REPL mode.

- giving 1 argument(mode),
  the action is the same as the case of 2 arguments,
  and 'mistate' is default to be `Base.active_repl.mistate`
"""
function enter_mode!(f, s :: LineEdit.MIState, lang_mode)
  LineEdit.transition(f, s, lang_mode)
end

function enter_mode!(s :: LineEdit.MIState, lang_mode)
  function default_trans_action()
    buf = copy(LineEdit.buffer(s))
    LineEdit.state(s, lang_mode).input_buffer = buf
  end
  enter_mode!(default_trans_action, s, lang_mode)
end

function enter_mode!(lang_mode)
  enter_mode!(Base.active_repl.mistate, lang_mode)
end

"""
```
CustomREPLDisplay(io::IO, replshow::Function)
```
Returns a `CustomREPLDisplay <: AbstractDisplay`, which behaves similarly to (e.g.), `TextDisplay`, except with
a custom `replshow` method that will be called instead of `Base.show` for a REPL with a `CustomREPLDisplay`.
"""
struct CustomREPLDisplay <: AbstractDisplay
    io::IO
    replshow::Function
end

# Extend Base.display with new methods for CustomREPLDisplay that use replshow
display(d::CustomREPLDisplay, @nospecialize x) = display(d, MIME"text/plain"(), x)
display(d::CustomREPLDisplay, M::MIME"text/plain", @nospecialize x) = d.replshow(d.io, M, x)

"""
```
enablecustomdisplay(repl::LineEditREPL, replshow::Function=show, io::IO=stdout)
```
Make a new LineEditREPL that has all the properties of `repl`, except with `repl.specialdisplay` set to `CustomREPLDisplay(io)`.
A repl with a `CustomREPLDisplay` set will dispatch to `replshow`, instead of `show`, to allow for custom display for various data types in different REPLs.
The `replshow` function should support three argument `show`, with a default method of
```
replshow(io, M, x) = show(io, M, x)
```
where `io` is an `IO` object, `M` is a `MIME` type, and `x` is the object to be shown
"""
function enablecustomdisplay(repl::LineEditREPL, replshow::Function=show, io::IO=stdout)
    customrepl = LineEditREPL(
        repl.t,
        repl.hascolor,
        repl.prompt_color,
        repl.input_color,
        repl.answer_color,
        repl.shell_color,
        repl.help_color,
        repl.history_file,
        repl.in_shell,
        repl.in_help,
        repl.envcolors
    )
    customrepl.interface = repl.interface
    customrepl.backendref = repl.backendref
    customrepl.specialdisplay = CustomREPLDisplay(io, replshow)
    return customrepl
end

end


text:
# this function copies the getindex tool OpenMendel/SnpArrays.jl uses
# to pull out the byte values from the compressed hex genotypes
# Repo: https://github.com/OpenMendel/SnpArrays.jl
@inline function _plinkindex(s::Matrix{UInt8}, i::Integer, j::Integer)
    #@boundscheck checkbounds(s, i, j)
    ip3 = i + 3
    (s[ip3 >> 2, j] >> ((ip3 & 0x03) << 1)) & 0x03
end

@inline function _plinkindex(s::Matrix{UInt8})
    _dim1, _dim2 = size(s)
    @inbounds reshape([_plinkindex(s, i, j) for i in 1:(_dim1*4) for j in 1:_dim2], :, _dim2)
end

@inline function _SNP(genotype::UInt8)
    #=
    00	Homozygous for first allele (0x00)
    01	Missing genotype (0x01)
    10	Heterozygous  (0x02)
    11	Homozygous for second allele in .bim file (0x03)
    =#
    genotype == 0x00 ? NTuple{2,Int8}((Int8(1),Int8(1))) :
        genotype == 0x01 ? missing :
        genotype == 0x02 ? NTuple{2,Int8}((Int8(1),Int8(2))) : NTuple{2,Int8}((Int8(2),Int8(2)))
end

@inline function _SNP(genomatrix::AbstractArray{UInt8})
    Matrix{Union{Missing, NTuple{2,Int8}}}(_SNP.(genomatrix))
end

function _plinkped(infile::String, keepfields::Union{Symbol,Vector{Symbol}} = :all, silent::Bool = false)
    basefile = splitext(infile)[1]
    bimfound = isfile(basefile * ".bim")
    if bimfound
        bimfile = CSV.read(
            basefile * ".bim",
            DataFrame,
            header = [:chrom, :snp, :cM, :bp, :allele1, :allele2],
            missingstring = ["0"],
            drop = [3,4],
            types = Dict(:chrom => String, :allele1 => Char, :allele2 => Char),
            pool=0.3
        )
        locinames = bimfile.snp
    end
    pedfile = CSV.read(
        infile,
        DataFrame,
        header = false,
        missingstring = ["-9", "0"],
        typemap = Dict(Int64 => Int8)
        )
    n = size(pedfile)[2] - 6
    locinames = bimfound ? locinames : ["snp_" * "$i" for i in 1:n]
    rename!(pedfile, vcat([:population, :name, :sire, :dam, :sex, :phenotype], Symbol.(locinames)))
    return pedfile
    nsamples = length(pedfile.name)
    npopulations = length(unique(pedfile.population))
    sirecheck = !all(ismissing.(pedfile.sire)) ?  "(✔)" : "(𐄂)"
    damcheck = !all(ismissing.(pedfile.dam)) ?  "(✔)" : "(𐄂)"
    sexcheck = !all(ismissing.(pedfile.sex)) ?  "(✔)" : "(𐄂)"
    phenotypecheck = !all(ismissing.(pedfile.phenotype)) ?  "(✔)" : "(𐄂)"
    if !silent
        if !bimfound
            @info "\n $(truncatepath(abspath(infile)))\n data: loci = $(n), samples = $(nsamples), populations = $(npopulations)\n .fam: sire $(sirecheck), dam $(damcheck), sex $(sexcheck), phenotype $(phenotypecheck)\n $(basefile * ".bim") not found, generating new marker names"
        else
            @info "\n $(truncatepath(abspath(infile)))\n data: loci = $(n), samples = $(nsamples), populations = $(npopulations)\n .fam: sire $(sirecheck), dam $(damcheck), sex $(sexcheck), phenotype $(phenotypecheck)"
        end
        println()
    end
end

function _plinkbed(infile::String, famfields::Union{Symbol,Vector{Symbol}} = :all, bimfields::Union{Symbol,Vector{Symbol}} = :all, silent::Bool = false)
    basefile = splitext(infile)[1]
    !isfile(basefile * ".fam") && throw(ArgumentError("$(basefile * ".fam") is required but wasn't found"))
    famfile = CSV.read(
        basefile * ".fam", DataFrame,
        header = [:population, :name, :sire, :dam, :sex, :phenotype],
        missingstring = ["-9", "0"],
        types = Dict(:sex => Int8)
    )
    nsamples = length(famfile.name)
    npopulations = length(unique(famfile.population))
    sirecheck = !all(ismissing.(famfile.sire)) ?  "(✔)" : "(𐄂)"
    damcheck = !all(ismissing.(famfile.dam)) ?  "(✔)" : "(𐄂)"
    sexcheck = !all(ismissing.(famfile.sex)) ?  "(✔)" : "(𐄂)"
    phenotypecheck = !all(ismissing.(famfile.phenotype)) ?  "(✔)" : "(𐄂)"
    bimfound = isfile(basefile * ".bim")
    if bimfound
        bimfile = CSV.read(
            basefile * ".bim",
            DataFrame,
            header = [:chromosome, :snp, :cm, :bp, :allele1, :allele2],
            missingstring = ["0"],
            #drop = [3,4],
            types = Dict(:chromosome => String, :allele1 => Char, :allele2 => Char),
            pool=0.3
        )
        locinames = bimfile.snp
        chromcheck = !all(ismissing.(bimfile.chromosome)) ?  "(✔)" : "(𐄂)"
        cmcheck = !all(ismissing.(bimfile.cm)) ?  "(✔)" : "(𐄂)"
        bpcheck = !all(ismissing.(bimfile.bp)) ?  "(✔)" : "(𐄂)"
    end
    bedfile = basefile * ".bed"
    data = open(bedfile, "r") do io
        Base.read(io, UInt16) == 0x1b6c || throw(ArgumentError("Incorrect \"magic\" number in $bedfile. It should be 0x1b6c"))
        Base.read(io, UInt8) == 0x01 || throw(ArgumentError("$bedfile is not in the correct orientation"))
        return Base.read(io)
    end
    nrows = (nsamples + 3) >> 2   # the number of rows in the Matrix{UInt8}
    n, r = divrem(length(data), nrows)
    locinames = bimfound ? locinames : ["snp_" * "$i" for i in 1:n]
    iszero(r) || throw(ArgumentError("The filesize of $bedfile is not a multiple of $nrows (and it should be)"))
    if !silent
        if !bimfound
            @info "\n $(truncatepath(abspath(infile)))\n data: loci = $(n), samples = $(nsamples), populations = $(npopulations)\n .fam: sire $(sirecheck), dam $(damcheck), sex $(sexcheck), phenotype $(phenotypecheck)\n $(basefile * ".bim") not found, generating new marker names"
        else
            @info "\n $(truncatepath(abspath(infile)))\n data: loci = $(n), samples = $(nsamples), populations = $(npopulations)\n .fam: sire $(sirecheck), dam $(damcheck), sex $(sexcheck), phenotype $(phenotypecheck)\n .bim: chromosome $(chromcheck), position $(cmcheck), bp $(bpcheck)"
        end
        println()
    end
    locinames = bimfound ? locinames : ["snp_" * "$i" for i in 1:n]
    genodf = DataFrame(
        :name => PooledArray(repeat(famfile.name, n) , compress = true),
        :population => PooledArray(repeat(famfile.population, n) , compress = true),
        :locus => PooledArray(repeat(locinames, inner = nsamples) , compress = true),
        :genotype => vec(_SNP(_plinkindex(reshape(data, (nrows, n)))))
    )
    #TODO generate_meta is deprecated
    metadf = generate_meta(genodf)
    if famfields != :none
        if famfields == :all
            sirecheck == "(✔)" && insertcols!(metadf, :sire => famfile.sire)
            damcheck == "(✔)" && insertcols!(metadf, :dam => famfile.dam)
            sexcheck == "(✔)" && insertcols!(metadf, :sex => famfile.sex)
            phenotypecheck == "(✔)" && insertcols!(metadf, :phenotype => famfile.phenotype)
        else (typeof(famfields) == Symbol) | (typeof(famfields) == Vector{Symbol})
            _fields = typeof(famfields) == Symbol ? [famfields] : famfields
            checkfields = symdiff(_fields, intersect(_fields, [:sire, :dam, :sex, :phenotype]))
            length(checkfields) > 0 && throw(ArgumentError("Unknown .fam file field(s) requested: $checkfields"))
            [insertcols!(metadf, i => famfile[:, i]) for i in _fields]
        end
    end
    if bimfound && bimfields != :none
        if bimfields == :all
            chromcheck == "(✔)" && insertcols!(genodf, 3, :chromosome => PooledArray(repeat(coalesce.(bimfile.chromosome,0), inner = nsamples), compress = true))
            cmcheck == "(✔)" && insertcols!(genodf, :cm => PooledArray(repeat(Float64.(coalesce.(bimfile.cm, 0)), inner = nsamples), compress = true))
            bpcheck == "(✔)" &&  insertcols!(genodf, :bp => PooledArray(repeat(Int64.(coalesce.(bimfile.bp, 0)), inner = nsamples), compress = true))
        else (typeof(bimfields) == Symbol) | (typeof(bimfields) == Vector{Symbol})
            _fields = typeof(bimfields) == Symbol ? [bimfields] : bimfields
            checkfields = symdiff(_fields, intersect(_fields, [:chromosome, :cm, :bp]))
            length(checkfields) > 0 && throw(ArgumentError("Unknown .bim file field(s) requested: $checkfields"))
            :chromosome ∈ _fields && insertcols!(genodf, 3, :chromosome => PooledArray(repeat(coalesce.(bimfile.chromosome,0), inner = nsamples), compress = true))
            :cm ∈ _fields && insertcols!(genodf, :cm => PooledArray(repeat(Float64.(coalesce.(bimfile.cm, 0)), inner = nsamples), compress = true))
            :bp ∈ _fields &&  insertcols!(genodf, :bp => PooledArray(repeat(Int64.(coalesce.(bimfile.bp, 0)), inner = nsamples), compress = true))
            #[insertcols!(genodf, i => PooledArray(repeat(coalesce.(bimfile[:, i], 0), inner = nsamples), compress = true)) for i in _fields]
        end
    end
    return PopData(metadf, genodf)
end

"""
    plink(infile::String; keepfields::Symbol|Vector{Symbol}, silent::Bool)
Read a PLINK `.ped` or binary `.bed` file into memory as a `PopData` object.
Requires an accompanying `.fam` file in the same directory, but an accompanying `.bim` file is optional.
- `infile::String` : path to `.ped` or `.bed` file

### Keyword Arguments
- `famfields::Symbol|Vector{Symbol}` : which additional fields to import from the `.fam` file
    - `:all` [default]
    - `:none`
    - any one or combination of `[:sire, :dam, :sex, :phenotype]`
- `bimfields::Symbol|Vector{Symbol}` : which additional fields to import from the optional `.bim` file
    - `:all` [default]
    - `:none`
    - any one or combination of `[:chromosome, :cm, :bp]`
- `silent::Bool`   : whether to print file information during import (default: `false`)

## Example
```julia
para = plink("datadir/parakeet.ped", famfields = :sex)
parr = plink("datadir/parrot.bed", famfields = [:sire, :dam], bimfields = :chromosome)
```
"""
function plink(infile::String; famfields::Union{Symbol,Vector{Symbol}} = :all, bimfields::Union{Symbol,Vector{Symbol}} = :all, silent::Bool = false)
    isfile(infile) || throw(ArgumentError("$infile not found in working directory."))
    if endswith(infile, ".ped")
        _plinkped(infile, famfields, bimfields, silent)
    elseif endswith(infile, ".bed")
        _plinkbed(infile, famfields, bimfields, silent)
    else
        throw(ArgumentError("Filename $infile is not recognized as having a .ped or .bed extension"))
    end
end

### writing to PLINK format ###

@inline _genoconversion(genotype::T) where T<:Genotype = join(genotype, " ")
@inline _genoconversion(genotype::Missing) = "0 0"

"""
    plink(data::PopData; filename::String)

Write a biallelic `PopData` object to PLINK `.ped` format with an accompanying
`.fam` file. Genotypes are coded by the PLINK standard:
- Integers are the alleles
- `0` encodes missing
- After column 6, every two numbers indicate a diploid genotype.
## Example
```julia
sharks = dropmultiallelic(@gulfsharks) ;
plink(sharks, filename = "biallelic_sharks.ped")
```
"""
function plink(data::PopData; filename::String)
    data.metadata.biallelic != true && throw(ArgumentError("To write to PLINK format, data must be biallelic.\nThis can be done with dropmultiallelic(::PopData)\n"))
    basename = endswith(filename, r".bed|.ped") ? splitext(filename)[1] : filename
    # the .fam file
    tmp = data.metadata[:, r"population|name|sire|dam|sex|phenotype"]
    current_cols = propertynames(tmp)
    fill_cols = symdiff(current_cols,[:population, :name, :sire, :dam, :sex, :phenotype])
    for i in fill_cols
        tmp[!, i] .= Int8(0)
    end
    select!(tmp, :population, :name, :sire, :dam, :sex, :phenotype)
    tmp.sex = coalesce.(tmp.sex, 0)
    CSV.write(basename * ".fam", tmp, delim = " ", writeheader = false, quotestrings = false, missingstring = "0")
    # the .ped file
    open(filename, "w") do outped
        eachind = groupby(data.genodata, :name)
        eachindmeta = groupby(tmp, :name)
        for i in data.genodata.name.pool
            metarow = eachindmeta[(;name = i)]
            metainfo = [metarow.population[1], metarow.name[1], metarow.sire[1], metarow.dam[1], metarow.sex[1], metarow.phenotype[1]]
            genos = eachind[(;name = i)].genotype
            printrow = vcat(metainfo, _genoconversion.(genos))
            join(outped, printrow, " ")
            println(outped, "")
        end
    end
    # the bim file
    tmp = unique(data.genodata, :locus)[:, r"chrom|locus|cm|bp"]
    current_cols = propertynames(tmp)
    fill_cols = symdiff(current_cols,[:chrom, :locus, :cm, :bp])
    for i in fill_cols
        tmp[!, i] .= Int8(0)
    end
    select!(tmp, :chrom, :locus, :cm, :bp)
    CSV.write(basename * ".bim", tmp, delim = " ", writeheader = false, quotestrings = false)
    return
end

# BED FORMAT
"""
    plink(data::PopData; filename::String)

Write a biallelic `PopData` object to PLINK `.bed` format with an accompanying
`.fam` file. Genotypes are coded by the PLINK standard:
- `00` Homozygous for first allele (`0x00`)
- `01` Missing genotype (`0x01`)
- `10` Heterozygous  (`0x02`)
- `11` Homozygous for second allele (`0x03`)

## Example
```julia
sharks = dropmultiallelic(@gulfsharks) ;
plink(sharks, filename = "biallelic_sharks.ped")
```
"""


text:
export UnitGaussianNormalizer

mutable struct UnitGaussianNormalizer
    mean_
    std_
    eps_
end

function UnitGaussianNormalizer(x;eps_=1f-5)
    mean_ = mean(x,dims=length(size(x)))
    std_ = std(x,dims=length(size(x)))
    return UnitGaussianNormalizer(mean_,std_,eps_)
end

function encode(normalizer::UnitGaussianNormalizer,x)
    x = (x.-normalizer.mean_)./(normalizer.std_.+normalizer.eps_)
end

function decode(normalizer::UnitGaussianNormalizer,x;sample_idx=nothing)
    std_ = normalizer.std_ .+ normalizer.eps_
    mean_ = normalizer.mean_
    x = x .* std_ .+ mean_
    return x
end


text:
function getcartlocation(env, width) :: Int
    world_width = env.x_threshold * 2
    scale = width / world_width
    floor(env.state[1] * scale + width / 2.0)  # middle of cart
end

function showenv(env, req_width)
    rgb_arr = render!(env)
    n_channels, height, width = size(rgb_arr)

    # squash height
    lbound :: Int = floor(height * 0.4)
    ubound :: Int = floor(height * 0.8)
    rgb_arr = rgb_arr[:, lbound:ubound, :]

    # squash width
#     view_width :: Int = floor(width * 0.6)
    cart_loc = getcartlocation(env._env, width)
#     if cart_loc <= view_width ÷ 2
#         w_range = 1:view_width
#     elseif cart_loc >= width - view_width ÷ 2
#         w_range = (width - view_width):width
#     else
#         w_range = (cart_loc - view_width ÷ 2):(cart_loc + view_width ÷ 2)
#     end
#     rgb_arr = rgb_arr[:, :, w_range]

    # convert to floating point and resize
    img = Gray.(colorview(RGB{Float64}, rgb_arr))
    ratio = req_width / min(width, length(lbound:ubound))
    img = imresize(img, ratio=ratio)
    rgb_flt = convert(Array{Float32}, channelview(img))
    rgb_flt = permutedims(rgb_flt, [2, 1])
    width, height = size(rgb_flt)
#     rgb_batch = similar(rgb_flt, eltype(rgb_flt), (width, height, 1, 1))
    rgb_batch = zeros(width, height, 1, 1)
    rgb_batch[:, :, 1, 1] = rgb_flt
    (rgb_batch |> gpu, floor(cart_loc * ratio))
end

max_stars_repo_name: darsnack/rl-experiments

max_issues_repo_name: darsnack/rl-experiments

max_forks_repo_name: darsnack/rl-experiments


text:
using jlpkg, Test, Pkg, Pkg.TOML

const root = joinpath(dirname(dirname(pathof(jlpkg))))
const test_cmd = ```$(Base.julia_cmd()) $(jlpkg.default_julia_flags)
    --code-coverage=$(["none", "user", "all"][Base.JLOptions().code_coverage+1])
    $(joinpath(root, "src", "cli.jl"))```
const jlpkg_version = match(r"^version = \"(\d+.\d+.\d+)\"$"m,
        read(joinpath(root, "Project.toml"), String)).captures[1]

function download_release(v::VersionNumber)
    x, y, z = v.major, v.minor, v.patch
    julia_exec = cd(mktempdir()) do
        julia = "julia-$(x).$(y).$(z)-linux-x86_64"
        tarball = "$(julia).tar.gz"
        sha256 = "julia-$(x).$(y).$(z).sha256"
        run(`curl -o $(tarball) -L https://julialang-s3.julialang.org/bin/linux/x64/$(x).$(y)/$(tarball)`)
        run(`curl -o $(sha256) -L https://julialang-s3.julialang.org/bin/checksums/$(sha256)`)
        run(pipeline(`grep $(tarball) $(sha256)`, `sha256sum -c`))
        mkpath(julia)
        run(`tar -xzf $(tarball) -C $(julia) --strip-components 1`)
        return abspath(julia, "bin", "julia")
    end
    return julia_exec
end

mktempdir() do tmpdir; mktempdir() do depot
    withenv("PATH" => tmpdir * (Sys.iswindows() ? ';' : ':') * get(ENV, "PATH", ""),
            "JULIA_DEPOT_PATH" => depot, "JULIA_PKG_DEVDIR" => nothing) do
        # Installation
        jlpkg.install(destdir=tmpdir)
        @test_throws ErrorException jlpkg.install(destdir=tmpdir)
        jlpkg.install(destdir=tmpdir, force = true)
        jlpkg.install(command="pkg", destdir=tmpdir)
        @test realpath(Sys.which("jlpkg" * (Sys.iswindows() ? ".cmd" : ""))) ==
              realpath(joinpath(tmpdir, "jlpkg" * (Sys.iswindows() ? ".cmd" : "")))
        @test realpath(Sys.which("pkg" * (Sys.iswindows() ? ".cmd" : ""))) ==  realpath(joinpath(tmpdir, "pkg" * (Sys.iswindows() ? ".cmd" : "")))
        @test_logs((:warn, r"can not be found in PATH"),
            jlpkg.install(command="cmd-not-in-path", destdir=joinpath(tmpdir, "dir-not-in-path")))
        mktempdir(tmpdir) do tmpdir2; withenv("PATH" => get(ENV, "PATH", "") * (Sys.iswindows() ? ';' : ':') * tmpdir2) do
            @test_logs (:warn, r"points to a different program") jlpkg.install(destdir=tmpdir2)
        end end
        # Basic usage
        @test Sys.iswindows() ? success(`cmd /c jlpkg --update st`) : success(`jlpkg --update st`)
        @test Sys.iswindows() ? success(`cmd /c jlpkg --help`) : success(`jlpkg --help`)
        @test Sys.iswindows() ? success(`cmd /c pkg --update st`) : success(`pkg --update st`)
        @test Sys.iswindows() ? success(`cmd /c pkg --help`) : success(`pkg --help`)
        @test success(`$(test_cmd) --update --project=$tmpdir add Example=7876af07-990d-54b4-ab0e-23690620f79a`)
        withenv("JULIA_PROJECT" => tmpdir) do
            @test success(`$(test_cmd) add JSON=682c06a0-de6a-54ab-a142-c8b1cf79cde6`)
        end
        project = TOML.parsefile(joinpath(tmpdir, "Project.toml"))
        @test project["deps"]["Example"] == "7876af07-990d-54b4-ab0e-23690620f79a"
        @test project["deps"]["JSON"] == "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
        @test success(`$(test_cmd) --project=$tmpdir rm Example JSON`)
        project = TOML.parsefile(joinpath(tmpdir, "Project.toml"))
        @test isempty(get(project, "deps", []))
        cd(tmpdir) do
            @test success(`$(test_cmd) --project add Example=7876af07-990d-54b4-ab0e-23690620f79a`)
            withenv("JULIA_PROJECT" => "@.") do
                @test success(`$(test_cmd) add JSON=682c06a0-de6a-54ab-a142-c8b1cf79cde6`)
            end
            project = TOML.parsefile(joinpath(tmpdir, "Project.toml"))
            @test project["deps"]["Example"] == "7876af07-990d-54b4-ab0e-23690620f79a"
            @test project["deps"]["JSON"] == "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
            @test success(`$(test_cmd) --project rm Example JSON`)
            project = TOML.parsefile(joinpath(tmpdir, "Project.toml"))
            @test isempty(get(project, "deps", []))
        end
        @test success(`$(test_cmd) --update --project=$tmpdir add https://github.com/JuliaLang/Example.jl`)
        project = TOML.parsefile(joinpath(tmpdir, "Project.toml"))
        @test project["deps"]["Example"] == "7876af07-990d-54b4-ab0e-23690620f79a"
        withenv("JULIA_LOAD_PATH" => tmpdir) do # Should work even though Pkg is not in LOAD_PATH
            @test success(`$(test_cmd) --update st -m`)
        end
        # Test --julia flag
        if Sys.islinux() && get(ENV, "CI", nothing) == "true"
            julia10 = download_release(v"1.0.4")
            julia11 = download_release(v"1.1.1")
            stdout, stderr = joinpath.(tmpdir, ("stdout.txt", "stderr.txt"))
            @test success(pipeline(`$(test_cmd) --julia=$(julia11) --version`, stdout=stdout, stderr=stderr))
            @test occursin(", julia version 1.1.1", read(stdout, String))
            @test isempty(read(stderr, String))
            @test success(pipeline(`$(test_cmd) --julia=$(julia11) --julia=$(julia10) --version`, stdout=stdout, stderr=stderr))
            @test occursin(", julia version 1.0.4", read(stdout, String))
            @test isempty(read(stderr, String))
            @test !success(pipeline(`$(test_cmd) --julia=juliafoobar --version`, stdout=stdout, stderr=stderr))
            @test isempty(read(stdout, String))
            @test occursin("Error: IOError: could not spawn `juliafoobar", read(stderr, String))
        end
        # Smoke test all Pkg commands in interpreted mode
        @test success(`$(test_cmd) activate foo`)
        @test success(`$(test_cmd) add SpecialFunctions=276daf66-3868-5448-9aa4-cd146d93841b`)
        @test success(`$(test_cmd) build SpecialFunctions`)
        @test success(`$(test_cmd) develop Example`)
        @test success(`$(test_cmd) free Example`)
        @test success(`$(test_cmd) remove SpecialFunctions`)
        @test success(`$(test_cmd) gc`)
        @test success(`$(test_cmd) pin Example`)
        @test success(`$(test_cmd) resolve`)
        @test success(`$(test_cmd) instantiate`)
        @test success(`$(test_cmd) precompile`)
        @test success(`$(test_cmd) preview up`)
        @test success(`$(test_cmd) test Example`)
        @test success(`$(test_cmd) update`)
        @test success(`$(test_cmd) status`)
        @test success(`$(test_cmd) help`)
        cd(tmpdir) do
            @test success(`$(test_cmd) generate HelloWorld`)
        end
        if VERSION > v"1.1"
            @test success(`$(test_cmd) registry add https://github.com/fredrikekre/Staging`)
            @test success(`$(test_cmd) registry status`)
            @test success(`$(test_cmd) registry update`)
            @test success(`$(test_cmd) registry remove Staging`)
        end
        # Output
        stdout, stderr = joinpath.(tmpdir, ("stdout.txt", "stderr.txt"))
        @test success(pipeline(`$(test_cmd) --help`, stdout=stdout, stderr=stderr))
        @test !occursin("No input arguments, showing help:", read(stdout, String))
        @test occursin("jlpkg - command line interface", read(stdout, String))
        @test isempty(read(stderr, String))
        @test !success(pipeline(`$(test_cmd)`, stdout=stdout, stderr=stderr))
        @test isempty(read(stdout, String))
        @test occursin("No input arguments, showing help:", read(stderr, String))
        @test occursin("jlpkg - command line interface", read(stderr, String))
        @test success(pipeline(`$(test_cmd) --version`, stdout=stdout, stderr=stderr))
        @test occursin("jlpkg version $(jlpkg_version), julia version $(VERSION)", read(stdout, String))
        @test isempty(read(stderr, String))
        # Error paths
        @test !success(pipeline(`$(test_cmd) --project=$tmpdir rm`, stdout=stdout, stderr=stderr))
        @test isempty(read(stdout, String))
        @test occursin("PkgError:", read(stderr, String))
        @test !success(pipeline(`$(test_cmd) --project=$tmpdir st --help`, stdout=stdout, stderr=stderr))
        @test isempty(read(stdout, String))
        @test occursin("PkgError:", read(stderr, String))
        @test !success(pipeline(`$(test_cmd)`, stdout=stdout, stderr=stderr))
        @test isempty(read(stdout, String))
        @test occursin("No input arguments, showing help:", read(stderr, String))
        @test !success(pipeline(`$(test_cmd) --foobar`, stdout=stdout, stderr=stderr))
        @test isempty(read(stdout, String))
        @test occursin("Invalid argument `--foobar`, showing help:", read(stderr, String))
        @test !success(pipeline(`$(test_cmd) --foobar st`, stdout=stdout, stderr=stderr))
        @test isempty(read(stdout, String))
        @test occursin("Invalid argument `--foobar`, showing help:", read(stderr, String))
        @test !success(pipeline(`$(test_cmd) --update`, stdout=stdout, stderr=stderr))
        @test isempty(read(stdout, String))
        @test occursin("No Pkg REPL arguments, showing help:", read(stderr, String))
        # Test that --compile and --optimize options are default, see issue #1
        # by running jlpkg test suite with jlpkg --project test on CI
        @test Base.JLOptions().opt_level === Int8(2)
        @test Base.JLOptions().compile_enabled === Int8(1)
    end
end end


text:
function daycount_testcases()
    dateranges = [
        Date(2014, 6, 4)  Date(2014, 6, 4)
        Date(2014, 8, 5)  Date(2016, 5, 25)
        Date(2015, 8, 5)  Date(2017, 5, 25)
        Date(2016, 8, 5)  Date(2018, 5, 25)
        Date(2016, 8, 29) Date(2018, 5, 31)
        Date(2016, 8, 30) Date(2018, 5, 31)
        Date(2016, 8, 31) Date(2018, 5, 31)
        Date(2014, 2, 28) Date(2016, 2, 28)
        Date(2014, 2, 28) Date(2016, 2, 29)
        Date(2015, 2, 27) Date(2017, 2, 28)
        Date(2015, 2, 28) Date(2017, 2, 28)
        Date(2014, 2, 28) Date(2016, 2, 28)
        Date(2014, 2, 28) Date(2016, 2, 29)
        Date(2015, 2, 27) Date(2017, 2, 28)
        Date(2015, 2, 28) Date(2017, 2, 28)
    ]

    daycount_Thirty360Bond  = [0, 650, 650, 650, 632, 630, 630, 720, 721, 721, 720, 720, 721, 721, 720]
    daycount_Thirty360US    = [0, 650, 650, 650, 632, 630, 630, 720, 721, 721, 720, 718, 720, 721, 720]
    daycount_ThirtyE360     = [0, 650, 650, 650, 631, 630, 630, 720, 721, 721, 720, 720, 721, 721, 720]
    daycount_ThirtyE360ISDA = [0, 650, 650, 650, 631, 630, 630, 718, 720, 723, 720, 718, 719, 721, 718]
    daycount_Actual         = [0, 659, 659, 658, 640, 639, 638, 730, 731, 732, 731, 730, 731, 732, 731]

    maturity = [Date(), Date(), Date(), Date(), Date(), Date(), Date(), Date(2016, 3, 1), Date(2016, 3, 1), Date(2017, 3, 1), Date(2017, 3, 1), Date(2016, 2, 29), Date(2016, 2, 29), Date(2017, 2, 28), Date(2017, 2, 28)]
    eom = [false, false, false, false, false, false, false, false, false, false, false, true, true, true, true]

    dateranges, daycount_Thirty360Bond, daycount_Thirty360US, daycount_ThirtyE360, daycount_ThirtyE360ISDA, daycount_Actual, maturity, eom
end

function test_daycount_yearlength()
    test(daycounter, result) = begin
        @test [yearlength(daycounter, year) for year in [1900, 2000, 2015, 2016]] == result
    end

    test(Thirty360Bond(),        [360, 360, 360, 360])
    test(Thirty360US(false),     [360, 360, 360, 360])
    test(ThirtyE360(),           [360, 360, 360, 360])
    test(ThirtyE360ISDA(Date()), [360, 360, 360, 360])
    test(Actual365Fixed(),       [365, 365, 365, 365])
    test(Actual364(),            [364, 364, 364, 364])
    test(Actual360(),            [360, 360, 360, 360])
    test(ActualActualISDA(),     [365, 366, 365, 366])
end

function test_daycount_daycount()
    dateranges, daycount_Thirty360Bond, daycount_Thirty360US, daycount_ThirtyE360, daycount_ThirtyE360ISDA, daycount_Actual, maturity, eom = daycount_testcases()

    test(daycounter, result) = begin
        @test [daycount(daycounter(i), dateranges[i,:]...) for i in 1:size(dateranges, 1)] == result
    end

    test(i -> Thirty360Bond(),             daycount_Thirty360Bond)
    test(i -> Thirty360US(eom[i]),         daycount_Thirty360US)
    test(i -> ThirtyE360(),                daycount_ThirtyE360)
    test(i -> ThirtyE360ISDA(maturity[i]), daycount_ThirtyE360ISDA)
    test(i -> Actual365Fixed(),            daycount_Actual)
    test(i -> Actual364(),                 daycount_Actual)
    test(i -> Actual360(),                 daycount_Actual)
    test(i -> ActualActualISDA(),          daycount_Actual)
end

function test_daycount_yearfraction()
    dateranges, daycount_Thirty360Bond, daycount_Thirty360US, daycount_ThirtyE360, daycount_ThirtyE360ISDA, daycount_Actual, maturity, eom = daycount_testcases()

    test(daycounter, result) = begin
        @test [yearfraction(daycounter(i), dateranges[i,:]...) for i in 1:size(dateranges, 1)] == result
    end

    test(i -> Thirty360Bond(),             daycount_Thirty360Bond  ./ 360)
    test(i -> Thirty360US(eom[i]),         daycount_Thirty360US    ./ 360)
    test(i -> ThirtyE360(),                daycount_ThirtyE360     ./ 360)
    test(i -> ThirtyE360ISDA(maturity[i]), daycount_ThirtyE360ISDA ./ 360)
    test(i -> Actual365Fixed(),            daycount_Actual         ./ 365)
    test(i -> Actual364(),                 daycount_Actual         ./ 364)
    test(i -> Actual360(),                 daycount_Actual         ./ 360)
    test(i -> ActualActualISDA(), [
        0
        convert(Float64, (1 - 216 // 365) + 1 + 145 // 366)
        convert(Float64, (1 - 216 // 365) + 1 + 144 // 365)
        convert(Float64, (1 - 217 // 366) + 1 + 144 // 365)
        convert(Float64, (1 - 241 // 366) + 1 + 150 // 365)
        convert(Float64, (1 - 242 // 366) + 1 + 150 // 365)
        convert(Float64, (1 - 243 // 366) + 1 + 150 // 365)
        convert(Float64, (1 -  58 // 365) + 1 +  58 // 366)
        convert(Float64, (1 -  58 // 365) + 1 +  59 // 366)
        convert(Float64, (1 -  57 // 365) + 1 +  58 // 365)
        convert(Float64, (1 -  58 // 365) + 1 +  58 // 365)
        convert(Float64, (1 -  58 // 365) + 1 +  58 // 366)
        convert(Float64, (1 -  58 // 365) + 1 +  59 // 366)
        convert(Float64, (1 -  57 // 365) + 1 +  58 // 365)
        convert(Float64, (1 -  58 // 365) + 1 +  58 // 365)
    ])
end

function test_daycount_all()
    print(rpad("Finance.DayCount...", 50, ' '))

    test_daycount_yearlength()
    test_daycount_daycount()
    test_daycount_yearfraction()

    println("PASS")
end


text:
include("../imports.jl")

@testset "TabularClassificationSingle" begin
    df = DataFrame(A = 1:4, B = ["M", "F", "F", "M"], C = ["P", "F", "P", "F"])
    td = TableDataset(df)

    method = TabularClassificationSingle(["P", "F"], td; catcols=(:B,), contcols=(:A,), )
    testencoding(method.encodings, method.blocks)
    DLPipelines.checkmethod_core(method)
    @test_nowarn methodlossfn(method)
    @test_nowarn methodmodel(method)

    @testset "`encodeinput`" begin
        row = mockblock(method.blocks[1])

        xtrain = encodeinput(method, Training(), row)
        @test length(xtrain[1]) == length(method.blocks[1].catcols)
        @test length(xtrain[2]) == length(method.blocks[1].contcols)

        @test eltype(xtrain[1]) <: Number
    end

    @testset "`encodetarget`" begin
        category = "P"
        y = encodetarget(method, Training(), category)
        @test y ≈ [1, 0]
    end

    @test_nowarn TabularClassificationSingle(["P", "F"], td)
    @test TabularClassificationSingle(["P", "F"], td).blocks[1].catcols == (:B, :C)

    FastAI.test_method_show(method, ShowText(Base.DevNull()))
end


text:
"""
N is either one or two
* `x` coordinates before interpolation
* `xi` coordinates after interpolation
"""
mutable struct P_core{T<:Real}
	x::Array{Array{T,1},1}
	xi::Array{Array{T,1},1}
	np::Int
	nd::Int
	nx::Vector{Int}
	nxi::Vector{Int}
	ivecinterp::Vector{Matrix{Int64}}
	ivecspray::Vector{Matrix{Int64}}
	iximin::Vector{Int64}
	iximax::Vector{Int64}
	ixmin::Vector{Int64}
	ixmax::Vector{Int64}
	interp_func::Function
	spray_func::Function
	Battrib::Symbol # select spline
	y_x::Matrix{T}
	y_z::Matrix{T}
end

function P_core(x::T, xi::T, Battrib::Symbol=:B1) where T

	if(Battrib == :B1)
		np=2;
		interp_func = interp_B1_1D
		spray_func = spray_B1_1D
	elseif(Battrib == :B2)
		np=4;
		interp_func = interp_B2_1D!
		spray_func = spray_B2_1D!
	else
		error("invalid attrib")
	end
	!(length(x)==length(xi)) && error("incorrect dimensions x")
	nd=length(x)
	nx=[length(x[id]) for id in 1:nd]
	nxi=[length(xi[id]) for id in 1:nd]

	ivecinterp=[zeros(Int64,np,nxi[id]) for id in 1:nd]
	ivecspray=[zeros(Int64,np,nx[id]) for id in 1:nd]

	iximin=zeros(Int64,nd)
	iximax=zeros(Int64,nd)
	ixmin=zeros(Int64,nd)
	ixmax=zeros(Int64,nd)

	for id in 1:nd
		# yi is only updated "within" bounds of y
		# get two indices and choose the inner most one
		iximin[id], iximax[id] = indminn_inside(xi[id], x[id][1], x[id][end])
		# special case for NONZERO step
		if(iximin[id]==iximax[id])
			iximax[id]=indminn(xi[id], x[id][1]+step(x[id]),1)[1]
		end
		# spary values of y "within" bounds of yi
		# get 2 indices and choose the inner most one
		ixmin[id], ixmax[id] = indminn_inside(x[id],xi[id][1],xi[id][end])
		# special case for NONZERO step
		if(ixmin[id]==ixmax[id])
			ixmax[id]=indminn(x[id], xi[id][1]+step(xi[id]),1)[1]
		end

		# index interp
		for i in 1:nxi[id]
			ivec=ivecinterp[id]
			ivecc=view(ivec, :, i)
			xn=x[id]
			xit=xi[id]
			indminn!(ivecc,xn,xit[i])
		end
		# index spray
		for i in 1:nx[id]
			ivec=ivecspray[id]
			ivecc=view(ivec, :, i)
			xn=x[id]
			xit=xi[id]
			indminn!(ivecc,xit,xn[i])
		end

	end

	# some other allocations for 2D
	if(nd==2)
		y_x=zeros(Float64, length(x[2]), length(iximin[1]:iximax[1]))
		y_z = zeros(Float64, length(x[2]), length(iximin[1]:iximax[1]))
	else
		y_x=zeros(1,1)
		y_z=zeros(1,1)
	end
	return P_core(
	       Array.(x),Array.(xi),	np,	nd,	nx,	nxi,	ivecinterp,	ivecspray,	iximin,	iximax,	ixmin,
	ixmax,	interp_func,	spray_func,	Battrib, y_x, y_z)

end # P_core

function interp_spray!(y::AbstractVector{Float64}, yi::AbstractVector{Float64} , pa::P_core{Float64}, attrib)

	if(attrib == :interp)
		# yi is only updated "within" bounds of y
		for i in pa.iximin[1] : pa.iximax[1]
			yi[i] = zero(Float64)
		end
		ivec=pa.ivecinterp[1]
		@simd for i in pa.iximin[1]:pa.iximax[1]
			ivecc=view(ivec,:,i)
			pa.interp_func(i, ivecc, yi, pa.x[1], y, pa.xi[1][i])
		end
	elseif(attrib == :spray)
		for i in eachindex(y)
			y[i] = zero(Float64)
		end
		# spary values of yi "within" bounds of y
		ivec=pa.ivecinterp[1]
		@simd for i in pa.iximin[1]:pa.iximax[1]
			ivecc=view(ivec,:,i)
			pa.spray_func(ivecc, y, pa.x[1], pa.xi[1][i], yi[i])
		end
	end
end # interp_spray!

"""
This allocates memory, need to be fixed
"""
function interp_spray!(y::Matrix{Float64}, yi::Matrix{Float64} , pa::P_core{Float64}, attrib)

	if(attrib == :interp)
		for i1 in pa.iximin[2]:pa.iximax[2], i2 in pa.iximin[1]:pa.iximax[1]
			yi[i1, i2] = zero(Float64);
		end
		fill!(pa.y_x, zero(Float64))

		# first along x
		ivec=pa.ivecinterp[1]
		@simd for ix in pa.iximin[1]:pa.iximax[1]
			@simd for iz in eachindex(pa.x[2])
				ivecc=view(ivec,:,ix)
				iy=(ix-pa.iximin[1])*length(pa.x[2])+iz
				yy=view(y, iz, :)
				pa.interp_func(iy, ivecc, pa.y_x, pa.x[1], yy, pa.xi[1][ix])
			end
		end
		# then along z
		ivec=pa.ivecinterp[2]
		@simd for ix in pa.iximin[1]:pa.iximax[1]
			@simd for iz in pa.iximin[2]:pa.iximax[2]
				ivecc=view(ivec,:,iz)
				iy=iz+(ix-1)*length(pa.xi[2])
				yy_x=view(pa.y_x, :, ix-pa.iximin[1]+1)
				pa.interp_func(iy, ivecc, yi, pa.x[2], yy_x, pa.xi[2][iz])
			end
		end
	elseif(attrib == :spray)
		# spary values of y only within bounds of yi
		# get 2 indices and choose the inner most one
		for i in eachindex(y)
			y[i]=zero(Float64)
		end
		pa.y_z[:]=zero(Float64)
		# first along z
		ivec=pa.ivecinterp[2]
		@simd for ix in pa.iximin[1]:pa.iximax[1]
			@simd for iz in pa.iximin[2]:pa.iximax[2]
				ivecc=view(ivec,:,iz)
				yy_z=view(pa.y_z,:,ix-pa.iximin[1]+1)
				pa.spray_func(ivecc, yy_z, pa.x[2], pa.xi[2][iz], yi[iz, ix])
			end
		end
		# then along x
		ivec=pa.ivecinterp[1]
		@simd for ix in  pa.iximin[1]:pa.iximax[1]
			@simd for iz in eachindex(pa.x[2])
				ivecc=view(ivec,:,ix)
				yy=view(y,iz,:)
				pa.spray_func(ivecc, yy, pa.x[1],  pa.xi[1][ix], pa.y_z[iz,ix-pa.iximin[1]+1])
			end
		end
	end
end # interp_spray!

"""
this subroutine interpolates or sprays bilinearly [one is adjoint of another]
interpolation returns y using Y[ivec[1]], Y[ivec[2]]
sprayg returns Y[ivec[1]], Y[ivec[2]] using y

                        +
                        |
    Y[ivec[1]]= f(X[ivec[1]])           |      Y[ivec[2]]= f(X[ivec[2]])
      +-----------------x--------+
                 y=f(x) |
                        +

bilinear interpolation
Reference: http://www.ajdesigner.com/phpinterpolation/bilinear_interpolation_equation.php
"""
function interp_B1_1D(iy, ivec, yV, X, Y, x)
	denom=inv(X[ivec[2]]-X[ivec[1]])
	if(isinf(denom))
		@inbounds yV[iy] = Y[ivec[1]]
	else
		@inbounds yV[iy] = ((Y[ivec[1]]*(X[ivec[2]]-x) + Y[ivec[2]]*(x-X[ivec[1]])) * denom)
	end
end # interpolate_spray_B1_1D

function spray_B1_1D(ivec, yV, X, x, y)
	denom = inv(X[ivec[2]] - X[ivec[1]])
	if(isinf(denom))
		@inbounds yV[ivec[1]] += y
		@inbounds yV[ivec[2]] += y
	else
		@inbounds yV[ivec[1]] += (y * (X[ivec[2]]-x)*denom)
		@inbounds yV[ivec[2]] += (y * (x-X[ivec[1]])*denom)
	end
end # interpolate_spray_B1_1D

"""
this subroutine interpolates or sprays
using cubic bspline
interpolation returns y using Y[ivec[1]], y2, Y[ivec[3]], Y[ivec[4]]
sprayg returns Y[ivec[1]], y2, Y[ivec[3]], Y[ivec[4]] using y

                       +
                       |
   Y[ivec[1]]      y2          |    Y[ivec[3]]       Y[ivec[4]]
     +-------+---------x----+--------+
                y=f(x) |
                       +
"""
function interp_B2_1D!(iy, ivec, yV, X,Y,x)

	"some constants"
	h = (X[ivec[2]] - X[ivec[1]]);
	hcube = h*h*h;
	hcubeI = hcube^(-1.0)

	if(h == zero(Float64))
		yV[iy] =  (Float64(0.25) * (Y[ivec[1]] + Y[ivec[2]] + Y[ivec[3]] + Y[ivec[4]]))
	else

		if((((x-X[ivec[1]])*(x-X[ivec[2]])) < zero(Float64)) | (x==X[ivec[1]]))
			"left edge interpolate"
			c1 = (1.0 / 6.0 * (2.0 * h - (x - (X[ivec[1]]-h)))^3) * hcubeI
			c4 = (1.0 / 6.0 * (2.0 * h - ((X[ivec[4]]-h) - x))^3) * hcubeI
			c2 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[2]]-h) - x)^2.0 * (2.0 * h - (x - (X[ivec[2]]-h))) ) )) * hcubeI
			c3 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[3]]-h) - x)^2.0 * (2.0 * h + (x - (X[ivec[3]]-h))) ) )) * hcubeI
			@inbounds yV[iy] = (Y[ivec[1]] * (c2 + 2.0*c1) + Y[ivec[2]] * (c3 - c1) + Y[ivec[3]] * c4)
		elseif((((x-X[ivec[3]])*(x-X[ivec[4]])) < zero(Float64)) | (x==X[ivec[4]]))
			"right edge interpolate"
			c1 = (1.0 / 6.0 * (2.0 * h - (x - (X[ivec[1]]+h)))^3) * hcubeI
			c4 = (1.0 / 6.0 * (2.0 * h - ((X[ivec[4]]+h) - x))^3) * hcubeI
			c2 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[2]]+h) - x)^2.0 * (2.0 * h - (x - (X[ivec[2]]+h))) ) )) * hcubeI
			c3 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[3]]+h) - x)^2.0 * (2.0 * h + (x - (X[ivec[3]]+h))) ) )) * hcubeI
			@inbounds yV[iy] = (Y[ivec[2]] * c1 + Y[ivec[3]] * (c2 - c4) + Y[ivec[4]] * (c3 + 2.0*c4))
		else
			c1 = (1.0 / 6.0 * (2.0 * h - (x - X[ivec[1]]))^3) * hcubeI
			c4 = (1.0 / 6.0 * (2.0 * h - (X[ivec[4]] - x))^3) * hcubeI
			c2 = ((2.0 / 3.0 * hcube - (0.5 * (X[ivec[2]] - x)^2.0 * (2.0 * h - (x - X[ivec[2]])) ) )) * hcubeI
			c3 = ((2.0 / 3.0 * hcube - (0.5 * (X[ivec[3]] - x)^2.0 * (2.0 * h + (x - X[ivec[3]])) ) )) * hcubeI
			"center interpolate"
			@inbounds yV[iy] = ((Y[ivec[1]] * c1 + Y[ivec[2]] * c2 + Y[ivec[3]] * c3 + Y[ivec[4]] * c4))
		end
	end
	return yV[iy]

end # interp_B2_1D

function spray_B2_1D!(ivec, yV, X, x, y)

	"some constants"
	h = (X[ivec[2]] - X[ivec[1]]);
	hcube = h*h*h;
	hcubeI = hcube^(-1.0)

	if(h == zero(Float64)) then
		@inbounds yV[ivec[1]] += (y)
		@inbounds yV[ivec[2]] += (y)
		@inbounds yV[ivec[3]] += (y)
		@inbounds yV[ivec[4]] += (y)
	else
		if((((x-X[ivec[1]])*(x-X[ivec[2]])) < zero(Float64)) | (x==X[ivec[1]]))
			"left edge spray"
			c1 = (1.0 / 6.0 * (2.0 * h - (x - (X[ivec[1]]-h)))^3) * hcubeI
			c4 = (1.0 / 6.0 * (2.0 * h - ((X[ivec[4]]-h) - x))^3) * hcubeI
			c2 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[2]]-h) - x)^2.0 * (2.0 * h - (x - (X[ivec[2]]-h))) ) )) * hcubeI
			c3 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[3]]-h) - x)^2.0 * (2.0 * h + (x - (X[ivec[3]]-h))) ) )) * hcubeI

			@inbounds yV[ivec[1]] += (y * (c2 + 2.0*c1))
			@inbounds yV[ivec[2]] += (y * (c3 - c1))
			@inbounds yV[ivec[3]] += (y * c4)
		elseif((((x-X[ivec[3]])*(x-X[ivec[4]])) < zero(Float64)) | (x==X[ivec[4]]))
			"right edge spray"
			c1 = (1.0 / 6.0 * (2.0 * h - (x - (X[ivec[1]]+h)))^3) * hcubeI
			c4 = (1.0 / 6.0 * (2.0 * h - ((X[ivec[4]]+h) - x))^3) * hcubeI
			c2 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[2]]+h) - x)^2.0 * (2.0 * h - (x - (X[ivec[2]]+h))) ) )) * hcubeI
			c3 = ((2.0 / 3.0 * hcube - (0.5 * ((X[ivec[3]]+h) - x)^2.0 * (2.0 * h + (x - (X[ivec[3]]+h))) ) )) * hcubeI

			@inbounds yV[ivec[2]] += y * c1
			@inbounds yV[ivec[3]] += y * (c2 - c4)
			@inbounds yV[ivec[4]] += y * (c3+ 2.0 * c4)
		else
			"center spray"
			c1 = (1.0 / 6.0 * (2.0 * h - (x - X[ivec[1]]))^3) * hcubeI
			c4 = (1.0 / 6.0 * (2.0 * h - (X[ivec[4]] - x))^3) * hcubeI
			c2 = ((2.0 / 3.0 * hcube - (0.5 * (X[ivec[2]] - x)^2.0 * (2.0 * h - (x - X[ivec[2]])) ) )) * hcubeI
			c3 = ((2.0 / 3.0 * hcube - (0.5 * (X[ivec[3]] - x)^2.0 * (2.0 * h + (x - X[ivec[3]])) ) )) * hcubeI
			@inbounds yV[ivec[1]] += (y * c1)
			@inbounds yV[ivec[2]] += (y * c2)
			@inbounds yV[ivec[3]] += (y * c3)
			@inbounds yV[ivec[4]] += (y * c4)
		end
	end
end # interpolate_spray_B3_1D