text:
# gwas.jl
# this file should ultimately be moved to a different package

## Reference implementation for no missing values, compare output with Base.linreg
# function univariate_ols{T}(x::Vector{T}, y::Vector{T})
#     xbar = mean(x)
#     ybar = mean(y)
#     n = length(x)
#
#     numerator = dot(x, y) - n*xbar*ybar
#     denominator = dot(x, x) - n*abs2(xbar)
#     b = numerator/denominator
#
#     R2 = dot(y, y) - n*abs2(ybar) - b*numerator
#     t = b*sqrt((n-2)*denominator/R2)
#     pval = tTestpvalue(t, n - 2)
#
#     b, ybar - b*xbar, t, pval
# end

import Distributions

# two-sided test
function tTestpvalue(t::Real, dof::Real)
    dof <= zero(dof) && return 1.0

    # Compute using F-distribution test for t^2
    return Distributions.ccdf(Distributions.FDist(1, dof), abs2(t))
end

immutable UnivariateOLSFit
    column::Int
    n::Int
    y0::Float64
    β::Float64

    se::Float64
    t::Float64
    pval::Float64
end

colno(fit::UnivariateOLSFit) = fit.column
npoints(fit::UnivariateOLSFit) = fit.n
ndof(fit::UnivariateOLSFit) = fit.n - 2
intercept(fit::UnivariateOLSFit) = fit.y0
coef(fit::UnivariateOLSFit) = fit.β
stderr(fit::UnivariateOLSFit) = fit.se
tvalue(fit::UnivariateOLSFit) = fit.t
pvalue(fit::UnivariateOLSFit) = fit.pval

function Base.show(io::IO, fit::BEDMatrices.UnivariateOLSFit)
    # print(io, "(", fit.varname, ", ", fit.n, ", ", fit.y0, ", ", fit.β, ", ", fit.pval, ")")
    print(io, "col. ", colno(fit), ": ", "y = ", intercept(fit), ifelse(coef(fit) > 0, " + ", " - "), abs(coef(fit)), "*x", "; p = ", pvalue(fit))
end

function Base.show(io::IO, ::MIME"text/plain", fit::BEDMatrices.UnivariateOLSFit)
    print(io, "BEDMatrices.UnivariateOLSFit:\n")
    label = string("column ", colno(fit), " with ", npoints(fit), " entries\n")
    print(io, "  ", label)
    print(io, "  intercept: ", fit.y0, "\tslope: ", fit.β, "\n")
    print(io, "  std err: ", fit.se, "\tp-value: ", fit.pval)
end

function gwas_writecsv(fout::IO, fits::Vector{UnivariateOLSFit}; header::Bool=true, labels::Vector{String}=Vector{String}(0))
    if header
        if length(labels) > 0
            write(fout, "\"snp_id\",")
        end
        write(fout, "\"coef\",", "\"std err\",", "\"t-value\",", "\"Pr(>|t|)\"\n")
    end
    if length(labels) > 0
        for fit in fits
            @inbounds write(fout, "\"", labels[colno(fit)], "\",", string(coef(fit)), ",", string(stderr(fit)), ",", string(tvalue(fit)), ",", string(pvalue(fit)), "\n")
        end
    else
        for fit in fits
            write(fout, string(coef(fit)), ",", string(stderr(fit)), ",", string(tvalue(fit)), ",", string(pvalue(fit)), "\n")
        end
    end
end

"""
    gwas_writecsv(file, fits::Vector{UnivariateOLSFit}; header::Bool=true, labels::Vector{String}=Vector{String}(0))

Writes `fits` to `file` in CSV format. `file` may either be a filename
or a writeable `IO`. `header` determine whether a header is written to
`file`. `labels` should be a list of labels for the fits; otherwise
omitted.

"""
function gwas_writecsv(filename::AbstractString, fits::Vector{UnivariateOLSFit}; header::Bool=true, labels::Vector{String}=Vector{String}(0))
    open(filename, "w") do fout
        gwas_writecsv(fout, fits; header=header, labels=labels)
    end
end

sub_mean(y::AbstractArray, ::Colon) = mean(y)

function sub_mean{R <: AbstractUnitRange}(y::AbstractArray, rranges::Vector{R})
    n = 0
    acc = 0.0
    for rrange in rranges
        @simd for r in rrange
            @inbounds acc += y[r]
        end
        n += length(rrange)
    end

    acc/n
end
function sub_mean(y::AbstractArray, rows)
    acc = 0.0
    @simd for r in rows
        @inbounds acc += y[r]
    end

    acc/length(rows)
end

sub_dot(x, y, ::Colon) = dot(x, y)

function sub_dot{R <: AbstractUnitRange}(x::AbstractArray, y::AbstractArray, rranges::Vector{R})
    # @boundscheck
    acc = 0.0
    for rrange in rranges
        @simd for r in rrange
            @inbounds acc += x[r]*y[r]
        end
    end

    acc
end

function sub_dot(x, y, rows)
    # @boundscheck
    acc = 0.0
    @simd for r in rows
        @inbounds acc += x[r]*y[r]
    end

    acc
end

"""
    _column_olsfit{T}(B::BEDMatrix, col::Integer, rows, y::Vector{T}, ybar::T, y2::T)

Performs OLS fit of `y` over `B[:, col]` and returns `(β, y0, n, ssr,
varx)`: the univariate OLS slope, `β`; the intercept, `y0`; the number
of non-missing entries, `n`; the sum of squared residuals, `ssr`; and
the variance of `B[:, col]`, `varx`.

"""
function _column_olsfit{T}(B::BEDMatrix, col::Integer, rows, y::Vector{T}, ybar::T, y2::T)
    N_zeros, N_ones, N_twos, N_nas, x_y, nasup_y, nasup_y2 = BEDMatrices.column_dist_dot(B, col, y, rows)
    ntot = N_zeros + N_ones + N_twos + N_nas

    @inbounds begin
        n = ntot
        if N_nas > 0  # there are missing values
            # adjust n
            n -= N_nas

            # adjust ybar and y2 to the sample
            ybar = (ntot*ybar - nasup_y)/n
            y2 = y2 - nasup_y2

            # insufficient data for a fit
            n < 2 && return UnivariateOLSFit(col, n, ybar, 0.0, 0.0, 0.0, 1.0)
        end

        xsum = N_ones + 2*N_twos
        x2 = N_ones + 4*N_twos
        numerator = x_y - xsum*ybar

        # I keep the denominator::Integer by writing the formula this
        # way.
        denominator = n*x2 - abs2(xsum)

        # Handle some singular cases in a defensible way
        if denominator === 0
            if numerator !== zero(numerator)
                # all the points have same x value but different y values
                return UnivariateOLSFit(col, n, -sign(xsum)*sign(numerator)*Inf,
                                        sign(numerator)*Inf, Inf, 0.0, 1.0)
            else
                # all the points have same x AND y values!
                return UnivariateOLSFit(col, n, NaN, NaN, NaN, NaN, NaN)
            end
        end

        β = n*numerator/denominator
        y0 = ybar - β*xsum/n

        R2 = y2 - n*abs2(ybar) - β*numerator

        se = sqrt(n*R2/((n - 2)*denominator))
        t = β/se
        pval = tTestpvalue(t, n - 2)
    end

    return UnivariateOLSFit(col, n, y0, β, se, t, pval)
end

function column_olsfit(B, col::Integer, y::AbstractArray, rows=(:))
    @boundscheck begin
        length(y) == size(B, 1) || throw(DimensionMismatch("vector must have same length as size(B, 1)"))
        checkbounds(B, :, col)
    end

    ybar = sub_mean(y, rows)
    y2 = sub_dot(y, y, rows)

    _column_olsfit(B, col, rows, y, ybar, y2)
end

st_column_olsfit(B, y::AbstractArray, rows=(:)) = st_column_olsfit(B, y, rows, indices(B, 2))
function st_column_olsfit(B, y::AbstractArray, rows, cols)
    ybar = sub_mean(y, rows)
    y2 = sub_dot(y, y, rows)

    gwas_results = Vector{UnivariateOLSFit}(length(cols))

    for col in cols
        @inbounds gwas_results[col] = _column_olsfit(B, col, rows, y, ybar, y2)
    end

    return gwas_results
end

mt_column_olsfit(B, y::AbstractArray, rows=(:)) = mt_column_olsfit(B, y, rows, indices(B, 2))
function mt_column_olsfit(B, y::AbstractArray, rows, cols)
    ybar::Float64 = sub_mean(y, rows)
    y2 = convert(Float64, sub_dot(y, y, rows))

    gwas_results = Vector{UnivariateOLSFit}(length(cols))

    _mt_fill_gwas!(gwas_results, B, y, ybar, y2, rows, cols)

    return gwas_results
end

# workaround for https://github.com/JuliaLang/julia/issues/15276, and https://github.com/JuliaLang/julia/issues/17395, https://github.com/yuyichao/explore/blob/8d52fb6caa745a658f2c9bbffd3b0f0fe4a2cc48/julia/issue-17395/scale.jl#L21
@noinline function _mt_fill_gwas!(results, B, y::AbstractArray, ybar, y2, rows, cols)
    @inbounds begin
        Threads.@threads for col in cols
            results[col] = _column_olsfit(B, col, rows, y, ybar, y2)
        end
    end
    results
end

mp_column_olsfit(B::BEDMatrix, y::AbstractArray, rows=(:)) = mp_column_olsfit(B::BEDMatrix, y::AbstractArray, rows, indices(B, 2))
function mp_column_olsfit(B::BEDMatrix, y::AbstractArray, rows, cols)
    @boundscheck begin
        length(y) == size(B, 1) || throw(DimensionMismatch("vector must have same length as size(B, 1)"))
    end
    ybar = sub_mean(y, rows)
    y2 = sub_dot(y, y, rows)

    gwas_results = pmap(col -> _column_olsfit(B, col, rows, y, ybar, y2), cols)
end

# function mt_column_olsfit(B::BEDSubMatrix, y::AbstractArray)
# end
#
# function mp_column_olsfit(B::BEDSubMatrix, y::AbstractArray)
# end

"""
    GWAS(B::BEDMatrix, y::AbstractArray; mode::Symbol=:multithreaded, outfile::IO=DevNull, verbose::Bool=false, method::Symbol=:ols)

Returns vector of `UnivariateOLSFit` containers with fields `(n, y0, β, se, t, pval)`
 for fits of `y` over _each_ column of `B`.

## Optional keyword arguments
* `mode` determines parallelism; can be one of `:single`, `:multithreaded`, or [Not implemented] `:multiprocessed`
* `method` determines type of regression; currently only `:ols`, ordinary least squares, supported
* `outfile::IO` gives IO stream for writing the results
* `verbose::Bool` determines whether informative messages are printed

"""
function GWAS(B::BEDMatrix, y::AbstractArray; mode::Symbol=:multithreaded, outfile::IO=DevNull, verbose::Bool=true, method::Symbol=:ols)
    mode in (:single, :multithreaded, :multiprocessed) || error("Invalid mode, $mode")
    mode === :multiprocessed && error("Multiprocessing not currently implemented.")

    method in (:ols,) || error("Unsupported method, $method")

    if mode === :multithreaded
        if Threads.nthreads() === 1
            verbose && info("Multithreaded requested, but only a single thread available.",
                            " Number of threads fixed by environment variable `JULIA_NUM_THREADS`,",
                            " see julia documentation for details.")
            mode = :single
        elseif verbose
            info("Using ",Threads.nthreads(), " threads.")
        end
    end

    results = mode === :single ? st_column_olsfit(B, y) : (mode === :multithreaded ? mt_column_olsfit(B, y) : mp_column_olsfit(B, y))

    if outfile !== DevNull
        gwas_writecsv(outfile, results, labels=colnames(B))
    end

    results
end

##################################################
# For general ::Matrix (not ::BEDMatrix)
function mean_impute!{T}(M::AbstractMatrix{T}, na::T)
    n, p = size(M)

    z = zero(T)
    μ = z/one(T)

    na_inds = fill(0, n)

    @inbounds begin
        for col in indices(M, 2)
            # set up
            total = z
            na_count = 0

            # find NAs and total on non-NAs
            for row in indices(M, 1)
                x = M[row, col]
                if x === na
                    na_count += 1
                    na_inds[na_count] = row
                else
                    total += x
                end
            end

            # replace NAs with mean
            μ = total/(n - na_count)
            for na_idx in 1:na_count
                M[na_inds[na_idx], col] = μ

                # clean up as we go
                na_inds[na_idx] = 0
            end
        end
    end

    M
end

function _matrix_column_covvar{T}(M::AbstractMatrix{T}, rows, col::Integer, y::AbstractArray, na::T)
    @boundscheck checkbounds(M, rows, col)

    @inbounds begin
        na_count = 0
        xsum = zero(M[1, col])
        x2 = xsum
        ysum = zero(y[1])
        y2 = ysum
        x_y = zero(M[1, col]*y[1])
        for row in rows
            x = M[row, col]

            if x === na
                na_count += 1
            else
                xsum += x
                x2 += abs2(x)
                ysum += y[row]
                y2 += abs2(y[row])
                x_y += x*y[row]
            end
        end
    end
    return na_count, xsum, ysum, x2, y2, x_y
end

function _matrix_column_olsfit{T}(M::AbstractMatrix{T}, rows, col::Integer, y::AbstractArray, na::T)
    @boundscheck checkbounds(M, rows, col)

    @inbounds begin
        na_count, xsum, ysum, x2, y2, x_y = _matrix_column_covvar(M, rows, col, y, na)
        n = length(rows) - na_count
        n < 2 && return UnivariateOLSFit(col, n, ysum/n, 0.0, 0.0, 0.0, 1.0)

        numerator = n*x_y - xsum*ysum
        denominator = n*x2 - abs2(xsum)

        # handle singular case(s)?

        β = numerator/denominator
        y0 = (ysum - β*xsum)/n

        R2 = y2 - (abs2(ysum) + β*numerator)/n
        se = sqrt(n*R2/((n-2)*denominator))
        t = β/se
        pval = tTestpvalue(t, n - 2)
    end

    return UnivariateOLSFit(col, n, y0, β, se, t, pval)
end

st_matrix_column_olsfit{T}(M::AbstractMatrix{T}, y::AbstractArray, na::T) = st_matrix_column_olsfit(M, y, indices(M, 1), indices(M, 2), na)
st_matrix_column_olsfit{T}(M::AbstractMatrix{T}, y::AbstractArray, rows, na::T) = st_matrix_column_olsfit(M, y, rows, indices(M, 2), na)

function st_matrix_column_olsfit{T}(M::AbstractMatrix{T}, y::AbstractArray, rows, cols, na::T)
    @boundscheck checkbounds(M, rows, cols)

    gwas_results = Vector{UnivariateOLSFit}(length(cols))

    @inbounds for (colidx, col) in enumerate(cols)
        gwas_results[colidx] = _matrix_column_olsfit(M, rows, col, y, na)
    end

    return gwas_results
end

mt_matrix_column_olsfit{T}(M::AbstractMatrix{T}, y::AbstractArray, na::T) = mt_matrix_column_olsfit(M, y, indices(M, 1), indices(M, 2), na)
mt_matrix_column_olsfit{T}(M::AbstractMatrix{T}, y::AbstractArray, rows, na::T) = mt_matrix_column_olsfit(M, y, rows, indices(M, 2), na)
function mt_matrix_column_olsfit{T}(M::AbstractMatrix{T}, y::AbstractArray, rows, cols, na::T)
    @boundscheck checkbounds(M, rows, cols)

    gwas_results = Vector{UnivariateOLSFit}(length(cols))

    _mt_matrix_fill_gwas!(gwas_results, M, y, rows, cols, na)
end

@noinline function _mt_matrix_fill_gwas!(results, M, y, rows, cols, na)
    @inbounds begin
        Threads.@threads for colidx in eachindex(cols)
            results[colidx] = _matrix_column_olsfit(M, rows, cols[colidx], y, na)
        end
    end

    results
end


text:
function GroundAtNightDetector(arg0::TopocentricFrame, arg1::PVCoordinatesProvider, arg2::jdouble, arg3::AtmosphericRefractionModel)
    return GroundAtNightDetector((TopocentricFrame, PVCoordinatesProvider, jdouble, AtmosphericRefractionModel), arg0, arg1, arg2, arg3)
end

function g(obj::GroundAtNightDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

ext: jl


text:
module ApproxBayes

using Distributions
using ProgressMeter
using StatsBase
using Compat
using DataFrames
using Colors
using Plots
#using PyPlot
using Distances

import Base.show

export
  # types
  ABCtype,
  Prior,
  Particle,
  abctype,
  ParticleRejection,
  SimData,
  ABCRejection,
  ABCSMC,
  ABCRejectionModel,
  ABCSMCModel,

  #functions
  ksdist,
  runabc,
  plotresults,
  plotparameterposterior,
  plotmodelposterior,
  writeoutput

### source files
include("types.jl")
include("ABCalgorithm.jl")
include("sampling.jl")
include("util.jl")
include("plots.jl")

end


text:
#using  NearestNeighbors: Metric
#import Distances: evaluate
using StatsBase: mode

#evaluate(dist::dwt, a, b) = dtw_distance(a, b, dist.w)

"""
   `dtw_distance(a, b, w)` is the basic dynamic time wraping function.
where `a` & `b` are the time series matrices and `w` is the percentage
of window for warping.
"""
function dtw_distance(a, b, w)
    l_a, l_b = length(a), length(b)
    FloatMax = maxintfloat(Float64)
    if w <= 0
        band = max(l_b, l_a)
    else
        band = floor(Int, w*max(l_b, l_a))
    end
    M = zeros(l_a+1,l_b+1)
    M[1, 2:end] .= FloatMax
    M[2:end, 1] .= FloatMax
    for k=2:l_a
        M[k, 2:end] =  (a .-  b[k]).^2
    end
    for i=2:l_a+1
        jstart = max(2, i-band)
        jstop = min(l_b+1, i+band+1)
        idx_inf_left = i-band-1

        if idx_inf_left >= 1
            M[i, idx_inf_left] = FloatMax
        end
        for j=jstart:jstop
            im = i-1
            jm = j-1
            M[i,j] = M[i,j] + min(min(M[im,j], M[i,jm]), M[im,jm])
        end
        if jstop < l_b + 1
            M[i,jstop] = FloatMax
        end
    end
    return M[l_a+1,l_b+1]
end

function Predict_new(m, X::Array, Y::Array, yplane::Array)
    k = m.n_neighbors
    n_train, serieslength = size(X)
    n_test = size(Y, 1)
    y_pred = zeros(n_test)
    DistanceMatrix = zeros(n_test, n_train)
    for i=1:n_test
        for j=1:n_train
            DistanceMatrix[i, j] = dtw_distance(X[j,:],Y[i,:], m.metric_params[1])
        end
    end
    index = select_sort(DistanceMatrix, k)
    y_index = yplane[Int.(index)]
    for i=1:n_test
        y_pred[i] = mode(y_index[i, :])
    end
    return y_pred, DistanceMatrix
end

function select_sort(A, k) # some times we get consicative indx
    n_test, n_train = size(A)
    index = zeros(n_test, k)
    for l=1:n_test
        for i=1:n_train
            if i > k
                break
            end
            min_idx = i
            for j=i+1:n_train
                if A[l, min_idx] > A[l, j]
                    min_idx = j
                end
            end
            A[l, i], A[l, min_idx] = A[l, min_idx], A[l, i]
            index[l, i] = min_idx
        end
    end
    return index
end

# Add tuning and fit!, predic interface


text:
"""
```
augment_states{T<:AbstractFloat}(m::AbstractModel, TTT::Matrix{T}, RRR::Matrix{T}, CCC::Matrix{T})
```

### Arguments

-`m`: the m object
-`TTT`, `RRR`, and `CCC`: matrices of the state transition equation

### Return values

- `TTT_aug`, `RRR_aug`, and `CCC_aug`: extend the corresponding input matrices to include
  jobservables which are growth rates.

### Description

Some observables in the model are growth rates, which are calculated as a linear combination
of a present and lagged state (which is not yet accounted for in the `TTT`, `RRR`,and `CCC`
matrices). To improve the performance of `gensys`, these additional states are added after
the model is solved. `augment_states` assigns an index to each lagged state, and extends the
input `TTT`, `RRR`, and `CCC` matrices to accommodate the additional states and capture the
lagged state value in the current state vector. `RRR` and `CCC` are mostly augmented with
zeros.

The diagram below shows how `TTT` is extended to `TTT_aug`.

                TTT_aug
     (m.endogenous_states_additional
                  x
     m.endogenous_states_additional)
     _________________________________
    |                     |           |
    |          TTT        | endog_    |
    | (endogenous_states  | states_   |
    |          x          | augmented |
    |  endogenous_states) |           |
    |_____________________|           |
    |                                 |
    |    endogenous_states_augmented  |
    |_________________________________|

"""
function augment_states{T<:AbstractFloat}(m::Model1002, TTT::Matrix{T}, RRR::Matrix{T}, CCC::Vector{T})
    endo     = m.endogenous_states
    endo_new = m.endogenous_states_augmented
    exo      = m.exogenous_shocks

    n_endo = n_states(m)
    n_exo  = n_shocks_exogenous(m)
    @assert (n_endo, n_endo) == size(TTT)
    @assert (n_endo, n_exo)  == size(RRR)
    @assert (n_endo,)        == size(CCC)

    # Initialize augmented matrices
    n_states_add = length(endo_new)
    TTT_aug = zeros(n_endo + n_states_add, n_endo + n_states_add)
    TTT_aug[1:n_endo, 1:n_endo] = TTT
    RRR_aug = [RRR; zeros(n_states_add, n_exo)]
    CCC_aug = [CCC; zeros(n_states_add)]

    ### TTT modifications

    # Track Lags
    TTT_aug[endo_new[:y_t1],     endo[:y_t]] = 1.0
    TTT_aug[endo_new[:c_t1],     endo[:c_t]] = 1.0
    TTT_aug[endo_new[:i_t1],     endo[:i_t]] = 1.0
    TTT_aug[endo_new[:w_t1],     endo[:w_t]] = 1.0
    TTT_aug[endo_new[:π_t1_dup], endo[:π_t]] = 1.0
    TTT_aug[endo_new[:L_t1],     endo[:L_t]] = 1.0
    TTT_aug[endo_new[:u_t1],     endo[:u_t]] = 1.0
    TTT_aug[endo_new[:e_gdp_t1], endo_new[:e_gdp_t]] = 1.0
    TTT_aug[endo_new[:e_gdi_t1], endo_new[:e_gdi_t]] = 1.0

    # Expected inflation
    TTT_aug[endo_new[:Et_π_t], 1:n_endo] = (TTT^2)[endo[:π_t], :]

    # The 8th column of the addition to TTT corresponds to "v_lr" which is set equal to
    # e_lr – measurement errors for the two real wage observables built in
    # as exogenous structural shocks.
    TTT_aug[endo_new[:lr_t], endo_new[:lr_t]]               = m[:ρ_lr]
    TTT_aug[endo_new[:tfp_t], endo_new[:tfp_t]]             = m[:ρ_tfp]
    TTT_aug[endo_new[:e_gdpdef_t], endo_new[:e_gdpdef_t]]   = m[:ρ_gdpdef]
    TTT_aug[endo_new[:e_corepce_t], endo_new[:e_corepce_t]] = m[:ρ_corepce]
    TTT_aug[endo_new[:e_gdp_t], endo_new[:e_gdp_t]]         = m[:ρ_gdp]
    TTT_aug[endo_new[:e_gdi_t], endo_new[:e_gdi_t]]         = m[:ρ_gdi]

    ### RRR modfications

    # Expected inflation
    RRR_aug[endo_new[:Et_π_t], :] = (TTT*RRR)[endo[:π_t], :]

    # Measurement Error on long rate
    RRR_aug[endo_new[:lr_t], exo[:lr_sh]] = 1.0

    # Measurement Error on TFP
    RRR_aug[endo_new[:tfp_t], exo[:tfp_sh]] = 1.0

    # Measurement Error on GDP Deflator
    RRR_aug[endo_new[:e_gdpdef_t], exo[:gdpdef_sh]] = 1.0

    # Measurement Error on Core PCE
    RRR_aug[endo_new[:e_corepce_t], exo[:corepce_sh]] = 1.0

    # Measurement Error on GDP Deflator
    RRR_aug[endo_new[:e_gdp_t], exo[:gdp_sh]] = 1.0
    RRR_aug[endo_new[:e_gdp_t], exo[:gdi_sh]] = m[:ρ_gdpvar] * m[:σ_gdp] ^ 2

    RRR_aug[endo_new[:e_gdi_t], exo[:gdi_sh]] = 1.0

    ### CCC Modifications

    # Expected inflation
    CCC_aug[endo_new[:Et_π_t]] = (CCC + TTT*CCC)[endo[:π_t]]

    return TTT_aug, RRR_aug, CCC_aug
end


text:
#------------------------------------------------------------------------------
"""
    OlsGMFn(Y,X)

LS of Y on X; for one dependent variable, Gauss-Markov assumptions

# Usage
(b,u,Yhat,V,R2a) = OlsGMFn(Y,X)

# Input
- `Y::Array`:     Tx1, the dependent variable
- `X::Array`:     Txk matrix of regressors (including deterministic ones)

# Output
- `b::Array`:     kx1, regression coefficients
- `u::Array`:     Tx1, residuals Y - yhat
- `Yhat::Array`:  Tx1, fitted values X*b
- `V::Array`:     kxk matrix, covariance matrix of b
- `R2a::Number`:  scalar, R2 value

"""
function OlsGMFn(Y,X)
    T    = size(Y,1)
    b    = X\Y
    Yhat = X*b
    u    = Y - Yhat
    σ²   = var(u)
    V    = inv(X'X)*σ²
    R2a  = 1 - σ²/var(Y)
    return b,u,Yhat,V,R2a
end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    OlsFn(Y,X,m=0)

LS of Y on X; for one dependent variable, using Newey-West covariance matrix

# Usage
(b,u,Yhat,V,R2a) = OlsFn(Y,X,m)

# Input
- `Y::Array`:     Tx1, the dependent variable
- `X::Array`:     Txk matrix of regressors (including deterministic ones)
- `m::Int`:       scalar, bandwidth in Newey-West

# Output
- `b::Array`:     kx1, regression coefficients
- `u::Array`:     Tx1, residuals Y - Yhat
- `Yhat::Array`:  Tx1, fitted values X*b
- `V::Array`:     kxk matrix, covariance matrix of b
- `R2a::Number`:  scalar, R2 value

"""
function OlsFn(Y,X,m=0)
    T    = size(Y,1)
    b    = X\Y
    Yhat = X*b
    u    = Y - Yhat
    S0   = NWFn(X.*u,m)*T          #Newey-West covariance matrix
    Sxx  = X'X
    V    = inv(Sxx)'S0*inv(Sxx)
    R2a  = 1 - var(u)/var(Y)
    return b,u,Yhat,V,R2a
end
#------------------------------------------------------------------------------

max_stars_repo_name: barunik/EmpiricalFinancePhD

max_issues_repo_name: barunik/EmpiricalFinancePhD

max_forks_repo_name: barunik/EmpiricalFinancePhD


text:
push!(LOAD_PATH,"../src/")
module Reporstat
include("Inegi.jl")
end

max_stars_repo_name: yollotltamayo/Covid

max_issues_repo_name: yollotltamayo/Covid

max_forks_repo_name: yollotltamayo/Covid


text:
function dis2spt(w::Array{Float64,1}, dis::Array{Float64,1}, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, dt::Float64, it::Int64)
    if length(w) < it
       error("out of source time range")
    end
    if length(dis) != nz*nx
       error("length of dis does not much model size")
    end
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    spt = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, it)
    tmp = zeros(Nz, Nx); dis_tmp = reshape(dis, nz, nx)
    tmp[zupper+1: zupper+nz, ext+1:ext+nx] = 1/2 * w[it] * dis_tmp
    tmp = vec(tmp)
    spt.pz[:] = tmp[:]
    spt.px[:] = tmp[:]
    return spt
end

function spt2dis(w::Array{Float64,1}, spt::SnapShot)
    nz = spt.nz ; nx = spt.nx;
    ext= spt.ext; iflag = spt.iflag;
    it = spt.it
    if length(w) < spt.it
       error("out of source time range")
    end
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    pz = reshape(spt.pz, Nz, Nx)[zupper+1:zupper+nz, ext+1:ext+nx]
    px = reshape(spt.px, Nz, Nx)[zupper+1:zupper+nz, ext+1:ext+nx]
    dis = 1/2 * w[it] * pz + 1/2 * w[it] * px
    return vec(dis)
end

function spt2wlet(w::Array{Float64,1}, dis::Array{Float64,1}, spt::SnapShot)
    nz = spt.nz ; nx = spt.nx;
    ext= spt.ext; iflag = spt.iflag;
    it = spt.it
    if length(w) < it
       error("out of source time range")
    end
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    p = reshape(spt.pz, Nz, Nx)[zupper+1:zupper+nz, ext+1:ext+nx] + reshape(spt.px, Nz, Nx)[zupper+1:zupper+nz, ext+1:ext+nx]
    p = vec(p)
    w[it] = 1/2 * dot(dis, p)
    return nothing
end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
using Base: Test
using IncrementalInference, Distributions
# going to introduce two new constraint types
import Base: convert
import IncrementalInference: getSample

mutable struct DevelopPrior <: FunctorSingleton
  x::Distribution
end
getSample(dpl::DevelopPrior, N::Int=1) = (reshape(rand(dpl.x, N),1,N), )

mutable struct DevelopLikelihood <: FunctorPairwise
  x::Distribution
end
getSample(dpl::DevelopLikelihood, N::Int=1) = (reshape(rand(dpl.x, N),1,N), )
function (vv::DevelopLikelihood)(res::Array{Float64},
            userdata,
            idx::Int,
            meas::Tuple,
            wXi::Array{Float64,2},
            wXj::Array{Float64,2}  )::Void
  #
  res[1] = meas[1][idx] - (wXj[1,idx] - wXi[1,idx])
  nothing
end

N  = 100
fg = emptyFactorGraph()

v1 = addNode!(fg, :x1, ContinuousScalar, N=N)

pr = DevelopPrior(Normal(10.0,1.0))
f1 = addFactor!(fg,[:x1],pr)

ensureAllInitialized!(fg)

# Juno.breakpoint("/home/dehann/.julia/v0.5/IncrementalInference/src/ApproxConv.jl",121)

pts = evalFactor2(fg, f1, v1.index, N=N)

@test sum(abs.(pts - 1.0) .< 5) < 30
@test sum(abs.(pts - 10.0) .< 5) > 30

v2 = addNode!(fg, :x2, ContinuousScalar, N=N)
pp = DevelopLikelihood(Normal(100.0,1.0))
f2 = addFactor!(fg, [:x1;:x2], pp)

ensureAllInitialized!(fg)

@test abs(Base.mean(getVal(fg, :x2))-110.0) < 10.0

v3 = addNode!(fg, :x3, ContinuousScalar, N=N)
v4 = addNode!(fg, :x4, ContinuousScalar, N=N)

ppMH = DevelopLikelihood(Normal(90.0,1.0))
f3 = addFactor!(fg, [:x2;:x3;:x4], ppMH, multihypo=(1.0,0.5,0.5))

# @test getData(f3).fnc.hypoverts == [:x3, :x4]
@test sum(abs.(getData(f3).fnc.hypotheses.p[1] .- 0.0)) < 0.1  # 1.0 becomes 0.0 for computational convenience
@test sum(abs.(getData(f3).fnc.hypotheses.p[2:3] .- 0.5)) < 0.1

setVal!(v2,1*ones(1,100))
setVal!(v3,2*ones(1,100))
setVal!(v4,3*ones(1,100))

@testset "Test multi-hypothesis factor convolution exploration" begin

  pts = evalFactor2(fg, f3, v2.index, N=N)

  @test 99 < sum(pts .<= -70.0)

  pts = evalFactor2(fg, f3, v3.index, N=N)

  @test 25 < sum(pts .== 3.0) < 75

  pts = evalFactor2(fg, f3, v4.index, N=N)

  @test 25 < sum(pts .== 2.0) < 75

end

println("Packing converters")

mutable struct PackedDevelopPrior <: PackedInferenceType
  x::String
  PackedDevelopPrior() = new()
  PackedDevelopPrior(x) = new(x)
end
function convert(::Type{PackedDevelopPrior}, d::DevelopPrior)
  PackedDevelopPrior(string(d.x))
end
function convert(::Type{DevelopPrior}, d::PackedDevelopPrior)
  DevelopPrior(IncrementalInference.normalfromstring(d.x))
end

mutable struct PackedDevelopLikelihood <: PackedInferenceType
  x::String
  PackedDevelopLikelihood() = new()
  PackedDevelopLikelihood(x) = new(x)
end
function convert(::Type{PackedDevelopLikelihood}, d::DevelopLikelihood)
  PackedDevelopLikelihood(string(d.x))
end
function convert(::Type{DevelopLikelihood}, d::PackedDevelopLikelihood)
  DevelopLikelihood(IncrementalInference.extractdistribution(d.x))
end

@testset "test packing and unpacking the data structure" begin

  topack = getData(f1)
  dd = convert(PackedFunctionNodeData{PackedDevelopPrior},topack)
  unpacked = convert(FunctionNodeData{CommonConvWrapper{DevelopPrior}},dd)

  @test abs(unpacked.fnc.usrfnc!.x.μ - 10.0) < 1e-10
  @test abs(unpacked.fnc.usrfnc!.x.σ - 1.0) < 1e-10

  topack = getData(f3)
  dd = convert(PackedFunctionNodeData{PackedDevelopLikelihood},topack)
  unpacked = convert(FunctionNodeData{CommonConvWrapper{DevelopLikelihood}},dd)

  # @test unpacked.fnc.hypoverts == Symbol[:x3; :x4]
  @test sum(abs.(unpacked.fnc.hypotheses.p[1] .- 0.0)) < 0.1
  @test sum(abs.(unpacked.fnc.hypotheses.p[2:3] .- 0.5)) < 0.1
  # str = "Symbol[:x3, :x4];[0.5, 0.5]"
  # IncrementalInference.parsemultihypostr(str)

end

# using KernelDensityEstimatePlotting
#
#
#
# plot(kde!(pts))

#


text:
#!/usr/bin/env julia

scriptDir=@__DIR__

if VERSION >= v"0.7.0-DEV.2005"
    using Pkg
end

ENV["PYTHON"] = string(ENV["CONDA_PREFIX"], "/bin/python")

@info "Installing the packages from" scriptDir * "/julia_requirements.txt"
pkgName=in
open(scriptDir * "/julia_requirements.txt") do f
    for i in enumerate(eachline(f))

      pkgId=i[1]
      pkgName=i[2]
      @info "Installing" pkgName
      Pkg.add(pkgName)

      @info "Testing: ", pkgName
      # load the package this needs to be called at top-level
      Expr(:toplevel, :(module ($pkgName) end))

    end
end

# Show the package list
@info "Installed Julia packages:"
Pkg.installed()
Pkg.status()

max_stars_repo_name: chrisdane/ESMValTool

max_issues_repo_name: chrisdane/ESMValTool

max_forks_repo_name: chrisdane/ESMValTool


text:
# Plane model segmentation
# http://pointclouds.org/documentation/tutorials/planar_segmentation.php

using PCLCommon
using PCLSampleConsensus
using PCLSegmentation
using Cxx

function showxyz(c, i)
    x,y,z = c[i,:x],c[i,:y],c[i,:z]
    println("x,y,z: $x,$y,$z")
end

cloud = PointCloud{PointXYZ}(15, 1)
for i = 0:length(cloud)-1
    cloud[i,:x] = randn()
    cloud[i,:y] = randn()
    cloud[i,:z] = 1.0
end

# Set a few outliers
cloud[0,:z] = 2.0
cloud[3,:z] = -2.0
cloud[6,:z] = 4.0

println("Point cloud data: $(length(cloud)) points")
for i = 0:length(cloud)-1
    showxyz(cloud, i)
end

coefficients = ModelCoefficients()
inliers = PointIndices()

seg = SACSegmentation{PointXYZ}()
setOptimizeCoefficients(seg, true)
setModelType(seg, SACMODEL_PLANE)
setMethodType(seg, SAC_RANSAC)
setDistanceThreshold(seg, 0.01)

setInputCloud(seg, cloud)
segment(seg, inliers, coefficients)

coefvalues = icxx"$(coefficients.handle)->values;"
indices = icxx"$(inliers.handle)->indices;"

if length(inliers) == 0
    error("Could not estimate a planar model for the given dataset.")
end

print("Model Coefficients: ")
for i in 0:length(coefvalues)-1
    print(Float32(coefvalues[i]), " ")
end
println("")

println("Model inliers: $(length(indices))")
for i in 0:length(indices)-1
    print(i, "  ")
    showxyz(cloud, Int32(indices[i]))
end


text:
using SeisIO.UW
import SeisIO.UW: readuwevt, uwdf, uwdf!, uwpf, uwpf!

printstyled("  UW\n", color=:light_green)
uwf0 = joinpath(path, "SampleFiles/UW/00*W")
uwf1 = joinpath(path, "SampleFiles/UW/99011116541")
uwf2 = joinpath(path, "SampleFiles/UW/94100613522o")
uwf3 = joinpath(path, "SampleFiles/UW/02062915175o")
uwf4 = joinpath(path, "SampleFiles/UW/00012502123W")
uwf5 = joinpath(path, "SampleFiles/UW/02062915205o")

S = read_data("uw", uwf0)
S1 = SeisData()
uwdf!(S1, joinpath(path, "SampleFiles/UW/00012502123W"))
@test S == S1

printstyled("  readuwevt\n", color=:light_green)

# Can we read from pickfile only? datafile only?
printstyled("    read (pickfile, datafile) from pickfile name\n", color=:light_green)
W = readuwevt(uwf1*"o")

printstyled("    read (pickfile, datafile) from datafile name\n", color=:light_green)

redirect_stdout(out) do
  W = readuwevt(uwf1*"W", v=3)
end

for i in ["UW.WWVB..TIM","UW.TCG..TIM","UW.TDH..EHZ","UW.VLM..EHZ"]
  @test !isempty(findall(W.data.id.==i))
  n = findfirst(W.data.id.==i)
  @test ≈(W.data.fs[n], 100.0)
end

# Can we read from filename stub?
printstyled("    read (pickfile, datafile) from filename stub\n", color=:light_green)
W = readuwevt(uwf1)
@test W.hdr.id == W.source.eid == "99011116541"
@test W.hdr.mag.val == 3.0f0
@test occursin("99011116541o", W.hdr.src)
@test W.hdr.ot == DateTime("1999-01-11T16:54:11.96")

S = breaking_seis()
n = S.n
S += convert(SeisData, W.data)
@test S.n == n + W.data.n

δt = 1.0e-6*(rem(W.hdr.ot.instant.periods.value*1000 - SeisIO.dtconst, 60000000))

i = findfirst(W.data.id.=="UW.TDH..EHZ")
pha = W.data[i].pha["P"].tt
@test ≈(pha + δt, 14.506)

i = findfirst(W.data.id.=="UW.VLM..EHZ")
pha = W.data[i].pha["S"].tt
@test ≈(pha + δt, 24.236)

i = findfirst(W.data.id.=="UW.VFP..EHZ")
@test W.data.misc[i]["dur"] == 116.0

printstyled("    pickfile handling\n", color=:light_green)

# What about when there is no data file?
W = readuwevt(uwf2)
@test W.hdr.mag.val == 0.9f0
@test occursin("94100613522o", W.hdr.src)
@test W.hdr.ot == DateTime("1994-10-06T13:52:39.02")

W = readuwevt(uwf3)
@test W.hdr.id == "041568"

printstyled("    data file with a time correction structure\n", color=:light_green)
redirect_stdout(out) do
  W = readuwevt(uwf4, v=2)
end

printstyled("    pick file with nonnumeric error info\n", color=:light_green)
redirect_stdout(out) do
  W = readuwevt(uwf5, v=2)
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id$
#*------------------------------------------------------------------
type DoublePipe_Basic
	DoublePipe_Basic()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components",
				:Hidden=>true
			)),
			fill(molweight(Dict{Symbol,Any}(
				:Brief=>"Component Mol Weight",
				:Hidden=>true
			)),(NComp)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Flag for Fluid Alocation ",
				:Valid=>["outer","inner"],
				:Default=>"outer",
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Inner Flow Regime ",
				:Valid=>["laminar","transition","turbulent"],
				:Default=>"laminar",
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Outer Flow Regime ",
				:Valid=>["laminar","transition","turbulent"],
				:Default=>"laminar",
				:Hidden=>true
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Laminar Flow for the Inner Side",
				:Valid=>["Hausen","Schlunder"],
				:Default=>"Hausen"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Transition Flow for the Inner Side",
				:Valid=>["Gnielinski","Hausen"],
				:Default=>"Gnielinski"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Turbulent Flow for the Inner Side",
				:Valid=>["Petukhov","SiederTate"],
				:Default=>"Petukhov"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Laminar Flow for the Outer Side",
				:Valid=>["Hausen","Schlunder"],
				:Default=>"Hausen"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Transition Flow for the OuterSide",
				:Valid=>["Gnielinski","Hausen"],
				:Default=>"Gnielinski"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Heat Transfer Correlation in Turbulent Flow for the Outer Side",
				:Valid=>["Petukhov","SiederTate"],
				:Default=>"Petukhov"
			)),
			DanaSwitcher(Dict{Symbol,Any}(
				:Brief=>"Options for convergence Calculations ",
				:Valid=>["Simplified","Full"],
				:Default=>"Full"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Estimated Duty",
				:Default=>70,
				:Lower=>1e-6,
				:Upper=>1e10
			)),
			DoublePipe_Geometry(Dict{Symbol,Any}(
				:Brief=>"Double pipe geometry",
				:Symbol=>" "
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Inner Stream",
				:PosX=>0,
				:PosY=>0.5225,
				:Symbol=>"_{inInner}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Outer Stream",
				:PosX=>0.2805,
				:PosY=>0,
				:Symbol=>"_{inOuter}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Inner Stream",
				:PosX=>1,
				:PosY=>0.5225,
				:Symbol=>"_{outInner}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Outlet Outer Stream",
				:PosX=>0.7264,
				:PosY=>1,
				:Symbol=>"_{outOuter}"
			)),
			Details_Main(Dict{Symbol,Any}(
				:Brief=>"Some Details in the Heat Exchanger",
				:Symbol=>" "
			)),
			Main_DoublePipe(Dict{Symbol,Any}(
				:Brief=>"Inner Side of the Heat Exchanger",
				:Symbol=>"_{Inner}"
			)),
			Main_DoublePipe(Dict{Symbol,Any}(
				:Brief=>"Outer Side of the Heat Exchanger",
				:Symbol=>"_{Outer}"
			)),
			[
				:(Outer.Properties.Average.T = 0.5*InletOuter.T + 0.5*OutletOuter.T),
				:(Inner.Properties.Average.T = 0.5*InletInner.T + 0.5*OutletInner.T),
				:(Outer.Properties.Average.P = 0.5*InletOuter.P+0.5*OutletOuter.P),
				:(Inner.Properties.Average.P = 0.5*InletInner.P+0.5*OutletInner.P),
				:(Inner.Properties.Wall.Twall = 0.5*Outer.Properties.Average.T + 0.5*Inner.Properties.Average.T),
				:(Outer.Properties.Wall.Twall = 0.5*Outer.Properties.Average.T + 0.5*Inner.Properties.Average.T),
				:(Outer.Properties.Average.Mw = sum(M*InletOuter.z)),
				:(Inner.Properties.Average.Mw = sum(M*InletInner.z)),
				:(Inner.Properties.Inlet.Fw = sum(M*InletInner.z)*InletInner.F),
				:(Inner.Properties.Outlet.Fw = sum(M*OutletInner.z)*OutletInner.F),
				:(Outer.Properties.Inlet.Fw = sum(M*InletOuter.z)*InletOuter.F),
				:(Outer.Properties.Outlet.Fw = sum(M*OutletOuter.z)*OutletOuter.F),
				:(OutletOuter.F = InletOuter.F),
				:(OutletInner.F = InletInner.F),
				:(OutletOuter.z=InletOuter.z),
				:(OutletInner.z=InletInner.z),
				:(Details.A=Geometry.Pi*Geometry.DoInner*Geometry.Lpipe),
				:(Inner.Properties.Average.Cp = PP.LiquidCp(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Average.rho = PP.LiquidDensity(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Inlet.rho = PP.LiquidDensity(InletInner.T,InletInner.P,InletInner.z)),
				:(Inner.Properties.Outlet.rho = PP.LiquidDensity(OutletInner.T,OutletInner.P,OutletInner.z)),
				:(Inner.Properties.Average.Mu = PP.LiquidViscosity(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Average.K = PP.LiquidThermalConductivity(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Wall.Mu = PP.LiquidViscosity(Inner.Properties.Wall.Twall,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Average.Cp = PP.VapourCp(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Average.rho = PP.VapourDensity(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Inlet.rho = PP.VapourDensity(InletInner.T,InletInner.P,InletInner.z)),
				:(Inner.Properties.Outlet.rho = PP.VapourDensity(OutletInner.T,OutletInner.P,OutletInner.z)),
				:(Inner.Properties.Average.Mu = PP.VapourViscosity(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Average.K = PP.VapourThermalConductivity(Inner.Properties.Average.T,Inner.Properties.Average.P,InletInner.z)),
				:(Inner.Properties.Wall.Mu = PP.VapourViscosity(Inner.Properties.Wall.Twall,Inner.Properties.Average.P,InletInner.z)),
				:(Outer.Properties.Average.Cp = PP.LiquidCp(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Average.rho = PP.LiquidDensity(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Inlet.rho = PP.LiquidDensity(InletOuter.T,InletOuter.P,InletOuter.z)),
				:(Outer.Properties.Outlet.rho = PP.LiquidDensity(OutletOuter.T,OutletOuter.P,OutletOuter.z)),
				:(Outer.Properties.Average.Mu = PP.LiquidViscosity(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Average.K = PP.LiquidThermalConductivity(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Wall.Mu = PP.LiquidViscosity(Outer.Properties.Wall.Twall,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Average.Cp = PP.VapourCp(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Average.rho = PP.VapourDensity(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Inlet.rho = PP.VapourDensity(InletOuter.T,InletOuter.P,InletOuter.z)),
				:(Outer.Properties.Outlet.rho = PP.VapourDensity(OutletOuter.T,OutletOuter.P,OutletOuter.z)),
				:(Outer.Properties.Average.Mu = PP.VapourViscosity(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Average.K = PP.VapourThermalConductivity(Outer.Properties.Average.T,Outer.Properties.Average.P,InletOuter.z)),
				:(Outer.Properties.Wall.Mu = PP.VapourViscosity(Outer.Properties.Wall.Twall,Outer.Properties.Average.P,InletOuter.z)),
				:(Details.Q = InletOuter.F*(InletOuter.h-OutletOuter.h)),
				:(Details.Q = InletInner.F*(OutletInner.h-InletInner.h)),
				:(Details.Q = InletInner.F*(InletInner.h-OutletInner.h)),
				:(Details.Q = InletOuter.F*(OutletOuter.h - InletOuter.h)),
				:(Inner.PressureDrop.fi*Inner.PressureDrop.Re = 16),
				:((Inner.PressureDrop.fi-0.0035)*(Inner.PressureDrop.Re^0.42) = 0.264),
				:((Inner.PressureDrop.fi-0.0035)*(Inner.PressureDrop.Re^0.42) = 0.264),
				:(Outer.PressureDrop.fi*Outer.PressureDrop.Re = 16),
				:((Outer.PressureDrop.fi-0.0035)*(Outer.PressureDrop.Re^0.42) = 0.264),
				:((Outer.PressureDrop.fi-0.0035)*(Outer.PressureDrop.Re^0.42) = 0.264),
				:(Inner.HeatTransfer.fi = 1/(0.79*ln(Inner.HeatTransfer.Re)-1.64)^2),
				:(Inner.HeatTransfer.Nu = 3.665 + ((0.19*((Geometry.DiInner/Geometry.Lpipe)*Inner.HeatTransfer.Re*Inner.HeatTransfer.PR)^0.8)/(1+0.117*((Geometry.DiInner/Geometry.Lpipe)*Inner.HeatTransfer.Re*Inner.HeatTransfer.PR)^0.467))),
				:(Inner.HeatTransfer.Nu = (49.027896+4.173281*Inner.HeatTransfer.Re*Inner.HeatTransfer.PR*(Geometry.DiInner/Geometry.Lpipe))^(1/3)),
				:(Inner.HeatTransfer.fi = 1/(0.79*ln(Inner.HeatTransfer.Re)-1.64)^2),
				:(Inner.HeatTransfer.Nu*(1+(12.7*sqrt(0.125*Inner.HeatTransfer.fi)*((Inner.HeatTransfer.PR)^(2/3) -1))) = 0.125*Inner.HeatTransfer.fi*(Inner.HeatTransfer.Re-1000)*Inner.HeatTransfer.PR),
				:(Inner.HeatTransfer.Nu =0.116*(Inner.HeatTransfer.Re^(0.667)-125)*Inner.HeatTransfer.PR^(0.333)*(1+(Geometry.DiInner/Geometry.Lpipe)^0.667)),
				:(Inner.HeatTransfer.fi = 1/(1.82*log(Inner.HeatTransfer.Re)-1.64)^2),
				:(Inner.HeatTransfer.Nu*(1.07+(12.7*sqrt(0.125*Inner.HeatTransfer.fi)*((Inner.HeatTransfer.PR)^(2/3) -1))) = 0.125*Inner.HeatTransfer.fi*Inner.HeatTransfer.Re*Inner.HeatTransfer.PR),
				:(Inner.HeatTransfer.Nu = 0.027*(Inner.HeatTransfer.PR)^(1/3)*(Inner.HeatTransfer.Re)^(4/5)),
				:(Inner.HeatTransfer.fi = 1/(1.82*log(Inner.HeatTransfer.Re)-1.64)^2),
				:(Outer.HeatTransfer.fi = 1/(0.79*ln(Outer.HeatTransfer.Re)-1.64)^2),
				:(Outer.HeatTransfer.Nu = 3.665 + ((0.19*((Outer.HeatTransfer.Dh/Geometry.Lpipe)*Outer.HeatTransfer.Re*Outer.HeatTransfer.PR)^0.8)/(1+0.117*((Outer.HeatTransfer.Dh/Geometry.Lpipe)*Outer.HeatTransfer.Re*Outer.HeatTransfer.PR)^0.467))),
				:(Outer.HeatTransfer.Nu = (49.027896+4.173281*Outer.HeatTransfer.Re*Outer.HeatTransfer.PR*(Outer.HeatTransfer.Dh/Geometry.Lpipe))^(1/3)),
				:(Outer.HeatTransfer.fi = 1/(0.79*ln(Outer.HeatTransfer.Re)-1.64)^2),
				:(Outer.HeatTransfer.Nu*(1+(12.7*sqrt(0.125*Outer.HeatTransfer.fi)*((Outer.HeatTransfer.PR)^(2/3) -1))) = 0.125*Outer.HeatTransfer.fi*(Outer.HeatTransfer.Re-1000)*Outer.HeatTransfer.PR),
				:(Outer.HeatTransfer.Nu = 0.116*(Outer.HeatTransfer.Re^(0.667)-125)*Outer.HeatTransfer.PR^(0.333)*(1+(Outer.HeatTransfer.Dh/Geometry.Lpipe)^0.667)),
				:(Outer.HeatTransfer.fi = 1/(0.79*ln(Outer.HeatTransfer.Re)-1.64)^2),
				:(Outer.HeatTransfer.fi = 1/(1.82*log(Outer.HeatTransfer.Re)-1.64)^2),
				:(Outer.HeatTransfer.Nu*(1.07+(12.7*sqrt(0.125*Outer.HeatTransfer.fi)*((Outer.HeatTransfer.PR)^(2/3) -1))) = 0.125*Outer.HeatTransfer.fi*Outer.HeatTransfer.Re*Outer.HeatTransfer.PR),
				:(Outer.HeatTransfer.Nu = 0.027*(Outer.HeatTransfer.PR)^(1/3)*(Outer.HeatTransfer.Re)^(4/5)),
				:(Outer.HeatTransfer.fi = 1/(1.82*log(Outer.HeatTransfer.Re)-1.64)^2),
				:(Outer.PressureDrop.Pdrop = Outer.PressureDrop.Pd_fric),
				:(Inner.PressureDrop.Pdrop = Inner.PressureDrop.Pd_fric),
				:(OutletOuter.P = InletOuter.P - Outer.PressureDrop.Pdrop),
				:(OutletInner.P = InletInner.P - Inner.PressureDrop.Pdrop),
				:(Outer.PressureDrop.Pd_fric = (2*Outer.PressureDrop.fi*Geometry.Lpipe*Outer.Properties.Average.rho*Outer.HeatTransfer.Vmean^2)/(Outer.PressureDrop.Dh*Outer.HeatTransfer.Phi)),
				:(Inner.PressureDrop.Pd_fric = (2*Inner.PressureDrop.fi*Geometry.Lpipe*Inner.Properties.Average.rho*Inner.HeatTransfer.Vmean^2)/(Geometry.DiInner*Inner.HeatTransfer.Phi)),
				:(Outer.PressureDrop.Pdrop = Outer.PressureDrop.Pd_fric),
				:(Inner.PressureDrop.Pdrop = Inner.PressureDrop.Pd_fric),
				:(OutletOuter.P = InletOuter.P - Outer.PressureDrop.Pdrop),
				:(OutletInner.P = InletInner.P - Inner.PressureDrop.Pdrop),
				:(Outer.PressureDrop.Pd_fric = 0.01*InletOuter.P),
				:(Inner.PressureDrop.Pd_fric = 0.01*InletInner.P),
				:(Inner.HeatTransfer.hcoeff = (Inner.HeatTransfer.Nu*Inner.Properties.Average.K/Geometry.DiInner)*Inner.HeatTransfer.Phi),
				:(Outer.HeatTransfer.hcoeff= (Outer.HeatTransfer.Nu*Outer.Properties.Average.K/Outer.HeatTransfer.Dh)*Outer.HeatTransfer.Phi),
				:(Outer.PressureDrop.Pd_ret = 0*"kPa"),
				:(Inner.PressureDrop.Pd_ret = 0*"kPa"),
				:(Outer.HeatTransfer.Phi = (Outer.Properties.Average.Mu/Outer.Properties.Wall.Mu)^0.14),
				:(Inner.HeatTransfer.Phi = (Inner.Properties.Average.Mu/Inner.Properties.Wall.Mu)^0.14),
				:(Outer.HeatTransfer.PR = ((Outer.Properties.Average.Cp/Outer.Properties.Average.Mw)*Outer.Properties.Average.Mu)/Outer.Properties.Average.K),
				:(Inner.HeatTransfer.PR = ((Inner.Properties.Average.Cp/Inner.Properties.Average.Mw)*Inner.Properties.Average.Mu)/Inner.Properties.Average.K),
				:(Outer.HeatTransfer.Re = (Outer.Properties.Average.rho*Outer.HeatTransfer.Vmean*Outer.HeatTransfer.Dh)/Outer.Properties.Average.Mu),
				:(Outer.PressureDrop.Re = (Outer.Properties.Average.rho*Outer.HeatTransfer.Vmean*Outer.PressureDrop.Dh)/Outer.Properties.Average.Mu),
				:(Inner.HeatTransfer.Re = (Inner.Properties.Average.rho*Inner.HeatTransfer.Vmean*Inner.HeatTransfer.Dh)/Inner.Properties.Average.Mu),
				:(Inner.PressureDrop.Re = Inner.HeatTransfer.Re),
				:(Outer.HeatTransfer.Vmean*(Outer.HeatTransfer.As*Outer.Properties.Average.rho) = Outer.Properties.Inlet.Fw),
				:(Inner.HeatTransfer.Vmean*(Inner.HeatTransfer.As*Inner.Properties.Average.rho) = Inner.Properties.Inlet.Fw),
				:(Details.Uc*((Geometry.DoInner/(Inner.HeatTransfer.hcoeff*Geometry.DiInner) )+(Geometry.DoInner*ln(Geometry.DoInner/Geometry.DiInner)/(2*Geometry.Kwall))+(1/(Outer.HeatTransfer.hcoeff)))=1),
				:(Details.Ud*(Geometry.Rfi*(Geometry.DoInner/Geometry.DiInner) + Geometry.Rfo + (Geometry.DoInner/(Inner.HeatTransfer.hcoeff*Geometry.DiInner) )+(Geometry.DoInner*ln(Geometry.DoInner/Geometry.DiInner)/(2*Geometry.Kwall))+(1/(Outer.HeatTransfer.hcoeff)))=1),
			],
			[
				"Outer	Stream Average Temperature","Inner Stream Average Temperature","Outer Stream Average Pressure","Inner Stream Average Pressure","Inner Stream Wall Temperature","Outer Stream Wall Temperature","Outer Stream Average Molecular Weight","Inner Stream Average Molecular Weight","Flow Mass Inlet Inner Stream","Flow Mass Outlet Inner Stream","Flow Mass Inlet Outer Stream","Flow Mass Outlet Outer Stream","Molar Balance Outer Stream","Molar Balance Inner Stream","Outer Stream Molar Fraction Constraint","Inner Stream Molar Fraction Constraint","Exchange Surface Area for one segment of pipe","Average Heat Capacity Inner Stream","Average Mass Density Inner Stream","Inlet Mass Density Inner Stream","Outlet Mass Density Inner Stream","Average Viscosity Inner Stream","Average	Conductivity Inner Stream","Viscosity Inner Stream at wall temperature","Average Heat Capacity InnerStream","Average Mass Density Inner Stream","Inlet Mass Density Inner Stream","Outlet Mass Density Inner Stream","Average Viscosity Inner Stream","Average Conductivity Inner Stream","Viscosity Inner Stream at wall temperature","Average Heat Capacity Outer Stream","Average Mass Density Outer Stream","Inlet Mass Density Outer Stream","Outlet Mass Density Outer Stream","Average Viscosity Outer Stream","Average Conductivity Outer Stream","Viscosity Outer Stream at wall temperature","Average Heat Capacity Outer Stream","Average Mass Density Outer Stream","Inlet Mass Density Outer Stream","Outlet Mass Density Outer Stream","Average Viscosity Outer Stream","Average Conductivity Outer Stream","Viscosity Outer Stream at wall temperature","Energy Balance Outer Stream","Energy Balance Inner Stream","Energy Balance Hot Stream","Energy Balance Cold Stream","Inner Side Friction Factor for Pressure Drop - laminar Flow","using Turbulent Flow","Inner Side Friction Factor","Outer Side Friction Factor - laminar Flow","using Turbulent Flow","Outer Side Friction Factor","Inner Side Friction Factor for Heat Transfer - laminar Flow","Nusselt Number","Nusselt Number","Inner Side Friction Factor for Heat Transfer - transition Flow","Nusselt Number","Nusselt Number","Inner Side Friction Factor for Heat Transfer - turbulent Flow","Nusselt Number","Nusselt Number","Inner Side Friction Factor for Heat Transfer - turbulent Flow","Outer Side Friction Factor for Heat Transfer - laminar Flow","Nusselt Number","Nusselt Number","Outer Side Friction Factor for Heat Transfer - transition Flow","Nusselt Number","Nusselt Number","Outer Side Friction Factor for Heat Transfer - transition Flow","Outer Side Friction Factor for Heat Transfer - turbulent Flow","Nusselt Number","Nusselt Number","Outer Side Friction Factor for Heat Transfer - turbulent Flow","Total Pressure Drop Outer Stream","Total Pressure Drop Inner Stream","Pressure Drop Outer Stream","Pressure Drop Inner Stream","Outer Pipe Pressure Drop for friction","Inner Pipe Pressure Drop for friction","Total Pressure Drop Outer Stream","Total Pressure Drop Inner Stream","Pressure Drop Outer Stream","Pressure Drop Inner Stream","Outer Pipe Pressure Drop for friction","Inner Pipe Pressure Drop for friction","Inner Pipe Film Coefficient","Outer Pipe Film Coefficient","Outer Pipe Pressure Drop due to return","Inner Pipe Pressure Drop due to return","Outer Pipe Phi correction","Inner Pipe Phi correction","Outer Pipe Prandtl Number","Inner Pipe Prandtl Number","Outer Pipe Reynolds Number for Heat Transfer","Outer Pipe Reynolds Number for Pressure Drop","Inner Pipe Reynolds Number for Heat Transfer","Inner Pipe Reynolds Number for Pressure Drop","Outer Pipe Velocity","Inner Pipe Velocity","Overall Heat Transfer Coefficient Clean","Overall Heat Transfer Coefficient Dirty",
			],
			[:PP,:NComp,:M,:HotSide,:innerFlowRegime,:outerFlowRegime,:InnerLaminarCorrelation,:InnerTransitionCorrelation,:InnerTurbulentCorrelation,:OuterLaminarCorrelation,:OuterTransitionCorrelation,:OuterTurbulentCorrelation,:CalculationApproach,:Qestimated,],
			[:Geometry,:InletInner,:InletOuter,:OutletInner,:OutletOuter,:Details,:Inner,:Outer,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	M::Array{molweight}
	HotSide::DanaSwitcher
	innerFlowRegime::DanaSwitcher
	outerFlowRegime::DanaSwitcher
	InnerLaminarCorrelation::DanaSwitcher
	InnerTransitionCorrelation::DanaSwitcher
	InnerTurbulentCorrelation::DanaSwitcher
	OuterLaminarCorrelation::DanaSwitcher
	OuterTransitionCorrelation::DanaSwitcher
	OuterTurbulentCorrelation::DanaSwitcher
	CalculationApproach::DanaSwitcher
	Qestimated::power
	Geometry::DoublePipe_Geometry
	InletInner::stream
	InletOuter::stream
	OutletInner::streamPH
	OutletOuter::streamPH
	Details::Details_Main
	Inner::Main_DoublePipe
	Outer::Main_DoublePipe
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export DoublePipe_Basic
function set(in::DoublePipe_Basic)
	#"Inner Pipe Cross Sectional Area for Flow"
	Inner.HeatTransfer.As = 0.25*Geometry.Pi*Geometry.DiInner*Geometry.DiInner
	 #"Outer Pipe Cross Sectional Area for Flow"
	Outer.HeatTransfer.As = 0.25*Geometry.Pi*(Geometry.DiOuter*Geometry.DiOuter - Geometry.DoInner*Geometry.DoInner)
	 #"Inner Pipe Hydraulic Diameter for Heat Transfer"
	Inner.HeatTransfer.Dh = Geometry.DiInner
	 #"Outer Pipe Hydraulic Diameter for Heat Transfer"
	Outer.HeatTransfer.Dh = (Geometry.DiOuter*Geometry.DiOuter-Geometry.DoInner*Geometry.DoInner)/Geometry.DoInner
	 #"Inner Pipe Hydraulic Diameter for Pressure Drop"
	Inner.PressureDrop.Dh = Geometry.DiInner
	 #"Outer Pipe Hydraulic Diameter for Pressure Drop"
	Outer.PressureDrop.Dh=Geometry.DiOuter-Geometry.DoInner

end
function setEquationFlow(in::DoublePipe_Basic)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	if InletInner.v == 0
		addEquation(18)
		addEquation(19)
		addEquation(20)
		addEquation(21)
		addEquation(22)
		addEquation(23)
		addEquation(24)
	else
		addEquation(25)
		addEquation(26)
		addEquation(27)
		addEquation(28)
		addEquation(29)
		addEquation(30)
		addEquation(31)
	end
	if InletOuter.v == 0
		addEquation(32)
		addEquation(33)
		addEquation(34)
		addEquation(35)
		addEquation(36)
		addEquation(37)
		addEquation(38)
	else
		addEquation(39)
		addEquation(40)
		addEquation(41)
		addEquation(42)
		addEquation(43)
		addEquation(44)
		addEquation(45)
	end
	let switch=HotSide
		if HotSide==InletInner.T > InletOuter.T
			set(switch,"inner")
		end
		if HotSide==InletInner.T < InletOuter.T
			set(switch,"outer")
		end
		if switch=="outer"
			addEquation(46)
			addEquation(47)
		elseif switch=="inner"
			addEquation(48)
			addEquation(49)
		end
	end
	let switch=innerFlowRegime
		if innerFlowRegime==Inner.PressureDrop.Re > 2300
			set(switch,"transition")
		end
		if innerFlowRegime==Inner.PressureDrop.Re < 2300
			set(switch,"laminar")
		end
		if innerFlowRegime==Inner.PressureDrop.Re > 10000
			set(switch,"turbulent")
		end
		if innerFlowRegime==Inner.PressureDrop.Re < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(50)
		elseif switch=="transition"
			addEquation(51)
		elseif switch=="turbulent"
			addEquation(52)
		end
	end
	let switch=outerFlowRegime
		if outerFlowRegime==Outer.PressureDrop.Re > 2300
			set(switch,"transition")
		end
		if outerFlowRegime==Outer.PressureDrop.Re < 2300
			set(switch,"laminar")
		end
		if outerFlowRegime==Outer.PressureDrop.Re > 10000
			set(switch,"turbulent")
		end
		if outerFlowRegime==Outer.PressureDrop.Re < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(53)
		elseif switch=="transition"
			addEquation(54)
		elseif switch=="turbulent"
			addEquation(55)
		end
	end
	let switch=innerFlowRegime
		if innerFlowRegime==Inner.HeatTransfer.Re > 2300
			set(switch,"transition")
		end
		if innerFlowRegime==Inner.HeatTransfer.Re < 2300
			set(switch,"laminar")
		end
		if innerFlowRegime==Inner.HeatTransfer.Re > 10000
			set(switch,"turbulent")
		end
		if innerFlowRegime==Inner.HeatTransfer.Re < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(56)
			let switch=InnerLaminarCorrelation
				if switch=="Hausen"
					addEquation(57)
				elseif switch=="Schlunder"
					addEquation(58)
				end
			end
		elseif switch=="transition"
			addEquation(59)
			let switch=InnerTransitionCorrelation
				if switch=="Gnielinski"
					addEquation(60)
				elseif switch=="Hausen"
					addEquation(61)
				end
			end
		elseif switch=="turbulent"
			let switch=InnerTurbulentCorrelation
				if switch=="Petukhov"
					addEquation(62)
					addEquation(63)
				elseif switch=="SiederTate"
					addEquation(64)
					addEquation(65)
				end
			end
		end
	end
	let switch=outerFlowRegime
		if outerFlowRegime==Outer.HeatTransfer.Re > 2300
			set(switch,"transition")
		end
		if outerFlowRegime==Outer.HeatTransfer.Re < 2300
			set(switch,"laminar")
		end
		if outerFlowRegime==Outer.HeatTransfer.Re > 10000
			set(switch,"turbulent")
		end
		if outerFlowRegime==Outer.HeatTransfer.Re < 10000
			set(switch,"transition")
		end
		if switch=="laminar"
			addEquation(66)
			let switch=OuterLaminarCorrelation
				if switch=="Hausen"
					addEquation(67)
				elseif switch=="Schlunder"
					addEquation(68)
				end
			end
		elseif switch=="transition"
			let switch=OuterTransitionCorrelation
				if switch=="Gnielinski"
					addEquation(69)
					addEquation(70)
				elseif switch=="Hausen"
					addEquation(71)
					addEquation(72)
				end
			end
		elseif switch=="turbulent"
			let switch=OuterTurbulentCorrelation
				if switch=="Petukhov"
					addEquation(73)
					addEquation(74)
				elseif switch=="SiederTate"
					addEquation(75)
					addEquation(76)
				end
			end
		end
	end
	let switch=CalculationApproach
		if switch=="Full"
			addEquation(77)
			addEquation(78)
			addEquation(79)
			addEquation(80)
			addEquation(81)
			addEquation(82)
		elseif switch=="Simplified"
			addEquation(83)
			addEquation(84)
			addEquation(85)
			addEquation(86)
			addEquation(87)
			addEquation(88)
		end
	end
	addEquation(89)
	addEquation(90)
	addEquation(91)
	addEquation(92)
	addEquation(93)
	addEquation(94)
	addEquation(95)
	addEquation(96)
	addEquation(97)
	addEquation(98)
	addEquation(99)
	addEquation(100)
	addEquation(101)
	addEquation(102)
	addEquation(103)
	addEquation(104)
end
function atributes(in::DoublePipe_Basic,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Basic Equations for rigorous double pipe heat exchanger model."
	fields[:Info]="Thermal analysis of double pipe heat exchanger using the NTU or LMTD Method.

== References ==

[1] E.A.D. Saunders, Heat Exchangers: Selection, Design and
 Construction, Longman, Harlow, 1988.

[2] Serth, Robert W., Process Heat Transfer: Principles and Applications, Elsevier, 2007.

[3] Gnielinski, V., Forced convection in ducts, in Heat Exchanger Design Handbook, Vol. 2
 Hemisphere Publishing Corp., New York, 1988."
	drive!(fields,_)
	return fields
end
DoublePipe_Basic(_::Dict{Symbol,Any})=begin
	newModel=DoublePipe_Basic()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
const v = 6.0  # velocity of seismic waves [km/s]

function earthquake_epicenter(x1, y1, t1, x2, y2, t2, x3, y3, t3)
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
abstract type ConstantApproach end;
struct FixConstant      <: ConstantApproach end;
struct DoNotFixConstant <: ConstantApproach end;

#"""
#    struct FESpaceWithConstantFixed{CS,<:ConstantApproach} <: SingleFieldFESpace
#      space::SingleFieldFESpace
#      constraint_style::Val{CS}
#      dof_to_fix::Int
#    end
#"""
struct FESpaceWithConstantFixed{CS,CA<:ConstantApproach} <: SingleFieldFESpace
  space::SingleFieldFESpace
  constraint_style::Val{CS}
  dof_to_fix::Int
  @doc """
      FESpaceWithConstantFixed(space::SingleFieldFESpace, fix_constant::Bool,
      dof_to_fix::Int=num_free_dofs(space))
  """
  function FESpaceWithConstantFixed(space::SingleFieldFESpace,
     fix_constant::Bool,
     dof_to_fix::Int=num_free_dofs(space))
    cs = constraint_style(space)
    CS = get_val_parameter(cs)
    if (fix_constant && num_dirichlet_dofs(space)==0)
      new{CS,FixConstant}(space,cs,dof_to_fix)
    else
      new{CS,DoNotFixConstant}(space,cs,dof_to_fix)
    end
  end
end

const FESpaceWithLastDofRemoved{CS} = FESpaceWithConstantFixed{CS,FixConstant}

@deprecate(FESpaceWithLastDofRemoved(space::SingleFieldFESpace),
           FESpaceWithConstantFixed(space,true))

get_constant_approach(::FESpaceWithConstantFixed{CS,CA}) where {CS,CA}=CA()

# Genuine functions
function num_free_dofs(f::FESpaceWithConstantFixed{CS,FixConstant}) where {CS}
  num_free_dofs(f.space)-1
end

function num_free_dofs(f::FESpaceWithConstantFixed{CS,DoNotFixConstant}) where {CS}
  num_free_dofs(f.space)
end

function zero_free_values(f::FESpaceWithConstantFixed)
  zeros(num_free_dofs(f))
end

function get_cell_dofs(f::FESpaceWithConstantFixed{CS,FixConstant}) where {CS}
  cell_dofs = get_cell_dofs(f.space)
  CellDofsWithDofFixed(cell_dofs,f.dof_to_fix)
end

function get_cell_dofs(f::FESpaceWithConstantFixed{CS,DoNotFixConstant}) where {CS}
  get_cell_dofs(f.space)
end

num_dirichlet_dofs(f::FESpaceWithConstantFixed{CS,FixConstant}) where {CS} = 1

num_dirichlet_dofs(f::FESpaceWithConstantFixed{CS,DoNotFixConstant}) where {CS} = 0

function zero_dirichlet_values(f::FESpaceWithConstantFixed)
  T = Float64 # TODO
  zeros(T,num_dirichlet_dofs(f))
end

num_dirichlet_tags(f::FESpaceWithConstantFixed{CS,FixConstant}) where {CS} = 1

num_dirichlet_tags(f::FESpaceWithConstantFixed{CS,DoNotFixConstant}) where {CS} = 0

get_dirichlet_dof_tag(f::FESpaceWithConstantFixed{CS,FixConstant}) where {CS} = Int8[1,]

get_dirichlet_dof_tag(f::FESpaceWithConstantFixed{CS,DoNotFixConstant}) where {CS} = Int8[]

function scatter_free_and_dirichlet_values(
  f::FESpaceWithConstantFixed{CS,FixConstant},fv,dv) where {CS}
  @assert length(dv) == 1
  _dv = similar(dv,eltype(dv),0)
  _fv = VectorWithEntryInserted(fv,f.dof_to_fix,dv[1])
  scatter_free_and_dirichlet_values(f.space,_fv,_dv)
end

function scatter_free_and_dirichlet_values(
  f::FESpaceWithConstantFixed{CS,DoNotFixConstant},fv,dv) where {CS}
  @assert length(dv) == 0
 scatter_free_and_dirichlet_values(f.space,fv,dv)
end

function gather_free_and_dirichlet_values(
  f::FESpaceWithConstantFixed{CS,FixConstant},cv) where {CS}
  _fv, _dv = gather_free_and_dirichlet_values(f.space,cv)
  @assert length(_dv) == 0
  fv = VectorWithEntryRemoved(_fv,f.dof_to_fix)
  dv = _fv[f.dof_to_fix:f.dof_to_fix]
  (fv, dv)
end

function gather_free_and_dirichlet_values(
  f::FESpaceWithConstantFixed{CS,DoNotFixConstant},cv) where {CS}
  gather_free_and_dirichlet_values(f.space,cv)
end

function gather_free_and_dirichlet_values!(
  fv,dv,f::FESpaceWithConstantFixed{CS,FixConstant},cv) where {CS}
  _fv, _dv = gather_free_and_dirichlet_values(f.space,cv)
  @assert length(_dv) == 0
  fv    .= VectorWithEntryRemoved(_fv,f.dof_to_fix)
  dv[1]  = _fv[f.dof_to_fix]
  (fv, dv)
end

function gather_free_and_dirichlet_values!(
  fv,dv,f::FESpaceWithConstantFixed{CS,DoNotFixConstant},cv) where {CS}
  gather_free_and_dirichlet_values(f.space,cv)
end

function TrialFESpace(f::FESpaceWithConstantFixed{CS,CA}) where {CS,CA}
  U = TrialFESpace(f.space)
  FESpaceWithConstantFixed(U,CA==FixConstant,f.dof_to_fix)
end

# Delegated functions

function get_cell_basis(f::FESpaceWithConstantFixed)
  get_cell_basis(f.space)
end

function get_cell_dof_basis(f::FESpaceWithConstantFixed)
  get_cell_dof_basis(f.space)
end

get_cell_axes(t::FESpaceWithConstantFixed)= get_cell_axes(t.space)

get_cell_axes_with_constraints(t::FESpaceWithConstantFixed)= get_cell_axes_with_constraints(t.space)

CellData.CellField(t::FESpaceWithConstantFixed,cell_vals) = CellField(t.space,cell_vals)

constraint_style(::Type{<:FESpaceWithConstantFixed{CS}}) where CS = Val{CS}()

# Helpers

struct CellDofsWithDofFixed{A<:AbstractArray} <: AbstractVector{Vector{Int}}
  cell_dofs::A
  dof_to_fix::Int
end

Base.size(a::CellDofsWithDofFixed) = (length(a.cell_dofs),)

Base.IndexStyle(::Type{<:CellDofsWithDofFixed}) = IndexLinear()

function Base.getindex(a::CellDofsWithDofFixed,i::Integer)
  cache = array_cache(a)
  getindex!(cache,a,i)
end

function array_cache(a::CellDofsWithDofFixed)
  @assert eltype(a.cell_dofs) == Vector{Int}
  b = testitem(a.cell_dofs)
  c = CachedArray(b)
  cache = array_cache(a.cell_dofs)
  (c, cache)
end

@inline function getindex!(d,a::CellDofsWithDofFixed,i::Integer)
  c, cache = d
  b = getindex!(cache,a.cell_dofs,i)
  setsize!(c,size(b))
  r = c.array
  for j in 1:length(b)
    bj = b[j]
    if bj == a.dof_to_fix
      rj = -1
    elseif bj < a.dof_to_fix
      rj = bj
    else
      rj = bj-1
    end
    r[j] = rj
  end
  r
end


text:
struct FEOperatorFromTerms <: FEOperator
  trial::FESpace
  test::FESpace
  assem::Assembler
  terms
end

"""
"""
function FEOperator(trial::FESpace,test::FESpace,assem::Assembler,terms::FETerm...)
  FEOperatorFromTerms(trial,test,assem,terms)
end

"""
"""
function FEOperator(trial::FESpace,test::FESpace,terms::FETerm...)
  assem = SparseMatrixAssembler(trial,test)
  FEOperator(trial,test,assem,terms...)
end

function FEOperator(mat::Type{<:AbstractSparseMatrix},trial::FESpace,test::FESpace,terms::FETerm...)
  assem = SparseMatrixAssembler(mat,trial,test)
  FEOperator(trial,test,assem,terms...)
end

function FEOperator(mat::Type{<:AbstractSparseMatrix},vec::Type{<:AbstractVector},trial::FESpace,test::FESpace,terms::FETerm...)
  assem = SparseMatrixAssembler(mat,vec,trial,test)
  FEOperator(trial,test,assem,terms...)
end

function get_test(op::FEOperatorFromTerms)
  op.test
end

function get_trial(op::FEOperatorFromTerms)
  op.trial
end

function allocate_residual(op::FEOperatorFromTerms,uh)
  @assert is_a_fe_function(uh)
  v = get_cell_basis(op.test)
  vecdata = collect_cell_residual(uh,v,op.terms)
  allocate_vector(op.assem, vecdata)
end

function residual!(b::AbstractVector,op::FEOperatorFromTerms,uh)
  @assert is_a_fe_function(uh)
  v = get_cell_basis(op.test)
  vecdata = collect_cell_residual(uh,v,op.terms)
  assemble_vector!(b,op.assem, vecdata)
  b
end

function allocate_jacobian(op::FEOperatorFromTerms,uh)
  @assert is_a_fe_function(uh)
  du = get_cell_basis(op.trial)
  v = get_cell_basis(op.test)
  matdata = collect_cell_jacobian(uh,du,v,op.terms)
  allocate_matrix(op.assem, matdata)
end

function jacobian!(A::AbstractMatrix,op::FEOperatorFromTerms,uh)
  @assert is_a_fe_function(uh)
  du = get_cell_basis(op.trial)
  v = get_cell_basis(op.test)
  matdata = collect_cell_jacobian(uh,du,v,op.terms)
  assemble_matrix!(A,op.assem,matdata)
  A
end

function residual_and_jacobian!(b::AbstractVector,A::AbstractMatrix,op::FEOperatorFromTerms,uh)
  @assert is_a_fe_function(uh)
  du = get_cell_basis(op.trial)
  v = get_cell_basis(op.test)
  data = collect_cell_jacobian_and_residual(uh,du,v,op.terms)
  assemble_matrix_and_vector!(A, b, op.assem, data)
  (b,A)
end

function residual_and_jacobian(op::FEOperatorFromTerms,uh)
  @assert is_a_fe_function(uh)
  du = get_cell_basis(op.trial)
  v = get_cell_basis(op.test)
  data = collect_cell_jacobian_and_residual(uh,du,v,op.terms)
  A, b = assemble_matrix_and_vector(op.assem, data)
  (b, A)
end


text:
# JuliaData/Tables.jl integration

# const TimeArrayColIter = Tables.EachColumn{<:TimeArray}
module TablesIntegration

using ..TimeSeries

using Tables

# S: the column names, including the extra column name for time index
# N: number of rows
# M: number of columns
struct TableIter{T<:AbstractTimeSeries,S,N,M}
    x::T
end

function TableIter(ta::T) where {T<:TimeArray}
    N, M = size(ta, 1), size(ta, 2) + 1
    S = Tuple(TimeSeries.replace_dupes!([:timestamp; colnames(ta)]))
    TableIter{T,S,N,M}(ta)
end

data(i::TableIter) = getfield(i, :x)

TimeSeries.timestamp(i::TableIter) = timestamp(data(i))
TimeSeries.colnames(i::TableIter)  = colnames(data(i))

Base.getindex(x::TableIter{<:TimeArray}, i::Integer) =
    i == 1 ? timestamp(x) : values(data(x)[colnames(x)[i - 1]])
Base.getindex(x::TableIter{<:TimeArray}, j::Integer, i::Integer) =
    i == 1 ? timestamp(x)[j] : values(data(x)[colnames(x)[i - 1]])[j]

Base.length(::TableIter{<:TimeArray,S,N,M}) where {S,N,M}    = M
Base.lastindex(::TableIter{<:TimeArray,S,N,M}) where {S,N,M} = M
Base.lastindex(::TableIter{<:TimeArray,S,N,M}, d::Integer) where {S,N,M} =
    ifelse(d == 1, N, ifelse(d == 2, M, 1))
Base.size(::TableIter{<:TimeArray,S,N,M}) where {S,N,M} = (N, M)

Base.propertynames(x::TableIter{<:TimeArray,S}) where {S} = S
Base.getproperty(x::TableIter{<:TimeArray,S}, c::Symbol) where {S} =
    c == S[1] ? timestamp(x) : values(getproperty(data(x), c))

function Base.iterate(x::TableIter, i::Integer = 1)
    i > length(x) && return nothing
    x[i], i + 1
end

Tables.istable(::Type{<:AbstractTimeSeries}) = true
Tables.rowaccess(::Type{<:TimeArray}) = true
Tables.rows(ta::TimeArray) = Tables.rows(Tables.columntable(ta))
Tables.columnaccess(::Type{<:TimeArray}) = true
Tables.columns(ta::TimeArray) = TableIter(ta)
Tables.eachcolumn(i::TableIter) = i
Tables.schema(ta::AbstractTimeSeries{T,N,D}) where {T,N,D} = Tables.schema(TableIter(ta))
Tables.schema(i::TableIter{T,S}) where {T,S} = Tables.Schema(S, coltypes(data(i)))

coltypes(x::AbstractTimeSeries{T,N,D}) where {T,N,D} = (D, (T for _ ∈ 1:size(x, 2))...)

function TimeSeries.TimeArray(x; timestamp::Symbol)
    Tables.istable(x) || throw(ArgumentError("TimeArray requires a table input"))

    sch = Tables.schema(x)
    names = sch.names
    (timestamp ∉ names) && throw(ArgumentError("time index `$timestamp` not found"))
    names′ = filter(!isequal(timestamp), collect(sch.names))

    cols = Tables.columns(x)
    val = mapreduce(n -> collect(getproperty(cols, n)), hcat, names′)
    TimeArray(collect(getproperty(cols, timestamp)), val, names′, x)
end

end  # TablesIntegration


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # [Benders Decomposition (Lazy Constraints)](@id benders_decomposition_lazy)

# **Originally Contributed by**: Mathieu Besançon

# This notebook describes how to implement the Benders decomposition in JuMP
# using lazy constraints. We keep the same notation and problem form as the
# first notebook *Benders decomposition*.

# ```math
# \begin{aligned}
# & \text{maximize} \quad &&c_1^T x+c_2^T v \\
# & \text{subject to} \quad &&A_1 x+ A_2 v \preceq b \\
# & &&x \succeq 0, x \in \mathbb{Z}^n \\
# & &&v \succeq 0, v \in \mathbb{R}^p \\
# \end{aligned}
# ```

# where $b \in \mathbb{R}^m$, $A_1 \in \mathbb{R}^{m \times n}$,
# $A_2 \in \mathbb{R}^{m \times p}$ and $\mathbb{Z}$ is the set of integers.
#
# Here the symbol $\succeq$ ($\preceq$) stands for element-wise greater (less)
# than or equal to. Any mixed integer programming problem can be written in the
# form above.

# For a detailed explanation on the Benders decomposition algorithm, see the
# introduction notebook.

# ### Lazy constraints

# Some optimization solvers allow users to interact with them during the
# solution process by providing user-defined functions which are triggered under
# certain conditions.

# The generic term for these functions is **callback**.

# In integer optimization, the main callback types are lazy constraints,
# user-cuts and heuristic solutions. See the [Callbacks](@ref callbacks_manual)
# section for an introduction on how to use them.
#
# Some callbacks define a new constraint which is only activated when necessary,
# i.e., when a current solution does not respect them. It can avoid building an
# optimization model with too many constraints up-front.
#
# This is the case for Benders decomposition, since the sub-problem defines an
# exponential number of primal vertices and therefore dual cuts.
#
# A detailed explanation on the distinction between user-cuts and lazy
# constraints is also available on [Paul Rubin's blog](https://orinanobworld.blogspot.com/2012/08/user-cuts-versus-lazy-constraints.html).
# He also [describes this approach to Benders Decomposition](http://orinanobworld.blogspot.ca/2011/10/benders-decomposition-then-and-now.html).

# We use the data from the original notebook and change the solution algorithm
# to leverage lazy constraints:

# * Step 1 (Initialization)
# * Step 2 (defining the subproblem model)
# * Step 3 (registering the lazy constraint of the subproblem)

# ### Data for the problem

# The input data is from page 139, Integer programming by Garfinkel and
# Nemhauser[[1]](#c1).

c1 = [-1; -4]
c2 = [-2; -3]

dim_x = length(c1)
dim_u = length(c2)

b = [-2; -3]

A1 = [1 -3;
     -1 -3]
A2 = [1 -2;
     -1 -1]

M = 1000;

# ### Loading the necessary packages

using JuMP
import GLPK
import Test

# Subproblem creation

function build_subproblem()
    sub_problem_model = Model(GLPK.Optimizer)
    @variable(sub_problem_model, u[1:dim_u] >= 0)
    @constraint(
        sub_problem_model,
        constr_ref_subproblem[j = 1:size(A2, 2)],
        A2[:, j]' * u >= c2[j],
    )
    return (sub_problem_model, u)
end

# ### Master Problem Description

master_problem_model = Model(GLPK.Optimizer);

(sub_problem_model, u) = build_subproblem();

# Variable Definition

@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int)
@variable(master_problem_model, t <= 1e6)

# Objective Setting

@objective(master_problem_model, Max, t)

print(master_problem_model)

# Track the calls to the callback

iter_num = 0

# Define lazy constraints

function benders_lazy_constraint_callback(cb_data)
    global iter_num
    iter_num += 1
    println("Iteration number = ", iter_num)

    x_current = callback_value.(Ref(cb_data), x)
    fm_current = callback_value(cb_data, t)

    c_sub = b - A1 * x_current
    @objective(sub_problem_model, Min, c1' * x_current + c_sub' * u)
    optimize!(sub_problem_model)

    print("\nThe current subproblem model is \n")
    print(sub_problem_model)

    t_status_sub = termination_status(sub_problem_model)
    p_status_sub = primal_status(sub_problem_model)

    fs_x_current = objective_value(sub_problem_model)

    u_current = value.(u)

    γ = b' * u_current

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current  ≈  fm_current # we are done
        @info("No additional constraint from the subproblem")
    end

    if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current < fm_current
        println("\nThere is a suboptimal vertex, add the corresponding constraint")
        cv = A1' * u_current - c1
        new_optimality_cons = @build_constraint(t + cv' * x <= γ)
        MOI.submit(
            master_problem_model,
            MOI.LazyConstraint(cb_data),
            new_optimality_cons,
        )
    end

    if t_status_sub == MOI.INFEASIBLE_OR_UNBOUNDED
        println("\nThere is an  extreme ray, adding the corresponding constraint")
        ce = A1' * u_current
        new_feasibility_cons = @build_constraint(dot(ce, x) <= γ)
        MOI.submit(
            master_problem_model,
            MOI.LazyConstraint(cb_data),
            new_feasibility_cons,
        )
    end
end

MOI.set(
    master_problem_model,
    MOI.LazyConstraintCallback(),
    benders_lazy_constraint_callback,
)

optimize!(master_problem_model)

t_status = termination_status(master_problem_model)
p_status = primal_status(master_problem_model)

if p_status == MOI.INFEASIBLE_POINT
    println("The problem is infeasible :-(")
end

if t_status == MOI.INFEASIBLE_OR_UNBOUNDED
    fm_current = M
    x_current = M * ones(dim_x)
end

if p_status == MOI.FEASIBLE_POINT
    fm_current = value(t)
    x_current = value.(x)
end

Test.@test value(t) ≈ -4 #hide

println("Status of the master problem is ", t_status,
        "\nwith fm_current = ", fm_current,
        "\nx_current = ", x_current)

# ### References

# ```@raw html
# <a id='c1'></a>
# ```

# 1. Garfinkel, R. & Nemhauser, G. L. Integer programming. (Wiley, 1972).

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
## generic linear algebra routines
CuMatOrAdj{T} = Union{CuMatrix, LinearAlgebra.Adjoint{T, <:CuMatrix{T}}, LinearAlgebra.Transpose{T, <:CuMatrix{T}}}
CuOrAdj{T} = Union{CuVecOrMat, LinearAlgebra.Adjoint{T, <:CuVecOrMat{T}}, LinearAlgebra.Transpose{T, <:CuVecOrMat{T}}}

function LinearAlgebra.tril!(A::CuMatrix{T}, d::Integer = 0) where T
  function kernel!(_A, _d)
    li = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    m, n = size(_A)
    if 0 < li <= m*n
      i, j = Tuple(CartesianIndices(_A)[li])
      if i < j - _d
        _A[i, j] = 0
      end
    end
    return nothing
  end

  blk, thr = cudims(A)
  @cuda blocks=blk threads=thr kernel!(A, d)
  return A
end

function LinearAlgebra.triu!(A::CuMatrix{T}, d::Integer = 0) where T
  function kernel!(_A, _d)
    li = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    m, n = size(_A)
    if 0 < li <= m*n
      i, j = Tuple(CartesianIndices(_A)[li])
      if j < i + _d
        _A[i, j] = 0
      end
    end
    return nothing
  end

  blk, thr = cudims(A)
  @cuda blocks=blk threads=thr kernel!(A, d)
  return A
end

function LinearAlgebra.copy_transpose!(dst::CuArray, src::CuArray)
  function kernel(dst, src)
    I = @cuindex dst
    dst[I...] = src[reverse(I)...]
    return
  end
  blk, thr = cudims(dst)
  @cuda blocks=blk threads=thr kernel(dst, src)
  return dst
end

# Matrix division

function Base.:\(_A::CuMatOrAdj, _B::CuOrAdj)
    A, B = copy(_A), copy(_B)
    A, ipiv = CuArrays.CUSOLVER.getrf!(A)
    return CuArrays.CUSOLVER.getrs!('N', A, ipiv, B)
end


text:
struct Incremental end
struct AllAtOnce end

mode(::Type) = Incremental()
mode(::Type{Scatter}) = AllAtOnce()
mode(::Type{BarPlot}) = AllAtOnce()

grouped_plot(P, args...; kwargs...) = _grouped_plot(mode(P), P, args...; kwargs...)

function _grouped_plot(::AllAtOnce, P, ax, df, group_dict, args, kws, group_pairs, style_pairs)
	xyz = [get(df, arg) for arg in args]

	pairs = lookup_symbols(df, group_pairs, style_pairs, group_dict)

	plt = plot!(P, ax, categorical_positions.(xyz)...; kws..., pairs...)

	categorical_ticks!(ax, xyz[1], xyz[2])

	(; plt)
end

function _grouped_plot(::Incremental, P, ax, gdf, group_dict, args, kws, group_pairs, style_pairs)

	if length(group_pairs) > 0
		grp = [var_key(last(p)) for p in pairs(group_pairs)]
		unique!(grp)
	else
		grp = Symbol[]
	end

	if haskey(group_pairs, :group)
		group_pairs = delete(group_pairs, :group)
	end

	out = combine(groupby(gdf, grp)) do df
		xyz = [get(df, arg) for arg in args]

		pairs = lookup_symbols(df, group_pairs, style_pairs, group_dict)

		plt = plot!(P, ax, categorical_positions.(xyz)...; kws..., pairs...)

		categorical_ticks!(ax, xyz[1], xyz[2])

		(; plt)
	end

	(; out.plt)
end

max_stars_repo_name: greimel/TabularMakie

max_issues_repo_name: greimel/TabularMakie

max_forks_repo_name: greimel/TabularMakie


text:
"""
    module NonLinearReactionAdvectionDiffusionWithFrontData

Пакет реализует исследование в рамках гранта РФФИ #20-31-70016 "Численные методы
решения обратных задач для нелинейных сингулярно возмущённых уравнений типа
реакция-диффузия-адвекция с данными о положении фронта реакции".

Автор исходного кода пакета:
- [Андрей Борзунов](mailto:aborzunov@physics.msu.ru), Кафедра математики физического факультета МГУ им. Ломоносова.

Исполнители гранта:
 - Мельникова Алина Александровна
 - Левашова Наталия Тимуровна
 - Лукьяненко Дмитрий Витальевич (Руководитель)
 - Быцюра Светлана Владимировна
 - Борзунов Андрей Анатольевич
 - Исаев Темур Фуркатович
 - Аргун Рауль Ларикович
 - Горбачев Александр Викторович

Документация: https://github.com/aborzunov/NonLinearReactionAdvectionDiffusionWithFrontData.jl
"""
module NonLinearReactionAdvectionDiffusionWithFrontData

using ForwardDiff
using LinearAlgebra
using Dierckx
using LaTeXStrings
using Plots
using Printf
using Missings
using ProgressMeter

export u_init, dparams, dparams_nonuniform;
export make_plot, make_gif, make_minimzation_gif;
export minimization_draft, draft;

export solve, solve_adjoint;
export phidetermination, Φ;
export apply_on_dynamic_mesh;
export f1, f2;
export initial_guess;
export generate_obs_data;

export delta;
export J, J_q, minimize;

ENV["GKS_ENCODING"] = "utf-8"

include("utils.jl")
include("direct.jl")
include("degenerated.jl")
include("initial_guess.jl")
include("adjoint.jl")
include("functional.jl")
include("plotting.jl")

end # module


text:
struct Callback{P, F, S}
    func :: F
    schedule :: S
    parameters :: P
end

@inline (callback::Callback)(sim) = callback.func(sim, callback.parameters)
@inline (callback::Callback{<:Nothing})(sim) = callback.func(sim)

"""
    Callback(func, schedule=IterationInterval(1); parameters=nothing)

Return `Callback` that executes `func` on `schedule`
with optional `parameters`. `schedule = IterationInterval(1)` by default.

If `isnothing(parameters)`, `func(sim::Simulation)` is called.
Otherwise, `func` is called via `func(sim::Simulation, parameteres)`.
"""
Callback(func, schedule=IterationInterval(1); parameters=nothing) =
    Callback(func, schedule, parameters)


text:
function randnumtable(weights::Array{Float64}, table::Array{Int})

    numtable = zeros(Int, size(table))

    B = unique( table )
    w = log( weights )

    for i in 1:length(B)

        max = B[i]
        if max > 0

            m = 1:max

            stirnums = map( x -> abs(stirlings1(max, x)), m)
            stirnums /= maximum(stirnums)

            for (idx, j) in enumerate(find(table .== max))

                clike = m .* w[idx]
                clike = cumsum(stirnums .* exp(clike - maximum(clike)))

                numtable[j] = 1+sum(rand() * clike[max] .> clike)
            end

        end

    end

    return numtable

end


text:
#!/usr/bin/env julia
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Libdl

# Build a system image binary at sysimg_path.dlext. Allow insertion of a userimg via
# userimg_path.  If sysimg_path.dlext is currently loaded into memory, don't continue
# unless force is set to true. Allow targeting of a CPU architecture via cpu_target.
function default_sysimg_path(debug=false)
    if Sys.isunix()
        splitext(Libdl.dlpath(debug ? "sys-debug" : "sys"))[1]
    else
        joinpath(dirname(Sys.BINDIR), "lib", "julia", debug ? "sys-debug" : "sys")
    end
end

"""
    build_sysimg(sysimg_path=default_sysimg_path(), cpu_target="native", userimg_path=nothing; force=false)

Rebuild the system image. Store it in `sysimg_path`, which defaults to a file named `sys.ji`
that sits in the same folder as `libjulia.{so,dylib}`, except on Windows where it defaults
to `Sys.BINDIR/../lib/julia/sys.ji`.  Use the cpu instruction set given by `cpu_target`.
Valid CPU targets are the same as for the `-C` option to `julia`, or the `-march` option to
`gcc`.  Defaults to `native`, which means to use all CPU instructions available on the
current processor. Include the user image file given by `userimg_path`, which should contain
directives such as `using MyPackage` to include that package in the new system image. New
system image will not replace an older image unless `force` is set to true.
"""
function build_sysimg(sysimg_path=nothing, cpu_target="native", userimg_path=nothing; force=false, debug=false)
    if sysimg_path === nothing
        sysimg_path = default_sysimg_path(debug)
    end

    # Quit out if a sysimg is already loaded and is in the same spot as sysimg_path, unless forcing
    sysimg = Libdl.dlopen_e("sys")
    if sysimg != C_NULL
        if !force && Base.samefile(Libdl.dlpath(sysimg), "$(sysimg_path).$(Libdl.dlext)")
            @info("System image already loaded at $(Libdl.dlpath(sysimg)), set force=true to override.")
            return nothing
        end
    end

    # Canonicalize userimg_path before we enter the base_dir
    if userimg_path !== nothing
        userimg_path = abspath(userimg_path)
    end

    # Enter base and setup some useful paths
    base_dir = dirname(Base.find_source_file("sysimg.jl"))
    cd(base_dir) do
        julia = joinpath(Sys.BINDIR, debug ? "julia-debug" : "julia")
        cc, warn_msg = find_system_compiler()

        # Ensure we have write-permissions to wherever we're trying to write to
        try
            touch("$sysimg_path.ji")
        catch
            err_msg =  "Unable to modify $sysimg_path.ji, ensure parent directory exists "
            err_msg *= "and is writable; absolute paths work best.)"
            error(err_msg)
        end

        # Copy in userimg.jl if it exists
        if userimg_path !== nothing
            if !isfile(userimg_path)
                error("$userimg_path is not found, ensure it is an absolute path.")
            end
            if isfile("userimg.jl")
                error("$(joinpath(base_dir, "userimg.jl")) already exists, delete manually to continue.")
            end
            cp(userimg_path, "userimg.jl")
        end
        try
            # Start by building basecompiler.{ji,o}
            basecompiler_path = joinpath(dirname(sysimg_path), "basecompiler")
            @info("Building basecompiler.o")
            @info("$julia -C $cpu_target --output-ji $basecompiler_path.ji --output-o $basecompiler_path.o compiler/compiler.jl")
            run(`$julia -C $cpu_target --output-ji $basecompiler_path.ji --output-o $basecompiler_path.o compiler/compiler.jl`)

            # Bootstrap off of that to create sys.{ji,o}
            @info("Building sys.o")
            @info("$julia -C $cpu_target --output-ji $sysimg_path.ji --output-o $sysimg_path.o -J $basecompiler_path.ji --startup-file=no sysimg.jl")
            run(`$julia -C $cpu_target --output-ji $sysimg_path.ji --output-o $sysimg_path.o -J $basecompiler_path.ji --startup-file=no sysimg.jl`)

            if cc !== nothing
                link_sysimg(sysimg_path, cc, debug)
            end

            # Spit out all our warning messages
            for w in warn_msg
                @warn w
            end

            if cc == nothing
                @info("System image successfully built at $sysimg_path.ji.")
            end

            if !Base.samefile("$(default_sysimg_path(debug)).ji", "$sysimg_path.ji")
                if isfile("$sysimg_path.$(Libdl.dlext)")
                    @info("To run Julia with this image loaded, run: `julia -J $sysimg_path.$(Libdl.dlext)`.")
                else
                    @info("To run Julia with this image loaded, run: `julia -J $sysimg_path.ji`.")
                end
            else
                @info("Julia will automatically load this system image at next startup.")
            end
        finally
            # Cleanup userimg.jl
            if userimg_path !== nothing && isfile("userimg.jl")
                rm("userimg.jl")
            end
        end
    end
end

# Search for a compiler to link sys.o into sys.dl_ext. Honor CC environment variable.
function find_system_compiler()
    warn_msg = String[] # save warning messages into an array

    if haskey(ENV, "CC")
        if !success(`$(ENV["CC"]) -v`)
            push!(warn_msg, "Using compiler override $(ENV["CC"]), but unable to run `$(ENV["CC"]) -v`.")
        end
        return ENV["CC"], warn_msg
    end

    # On Windows, check to see if WinRPM is installed, and if so, see if gcc is installed
    if Sys.iswindows()
        try
            Core.eval(Main, :(using WinRPM))
            winrpmgcc = joinpath(WinRPM.installdir, "usr", "$(Sys.ARCH)-w64-mingw32",
                "sys-root", "mingw", "bin", "gcc.exe")
            if success(`$winrpmgcc --version`)
                return winrpmgcc, warn_msg
            else
                throw()
            end
        catch
            push!(warn_msg, "Install GCC via `Pkg.add(\"WinRPM\"); WinRPM.install(\"gcc\")` to generate sys.dll for faster startup times.")
        end
    end

    # See if `cc` exists
    try
        if success(`cc -v`)
            return "cc", warn_msg
        end
    catch
        push!(warn_msg, "No supported compiler found; startup times will be longer.")
        return nothing, warn_msg
    end
end

# Link sys.o into sys.$(dlext)
function link_sysimg(sysimg_path=nothing, cc=find_system_compiler(), debug=false)
    if sysimg_path === nothing
        sysimg_path = default_sysimg_path(debug)
    end
    julia_libdir = dirname(Libdl.dlpath(debug ? "libjulia-debug" : "libjulia"))

    FLAGS = ["-L$julia_libdir"]

    push!(FLAGS, "-shared")
    push!(FLAGS, debug ? "-ljulia-debug" : "-ljulia")
    if Sys.iswindows()
        push!(FLAGS, "-lssp")
    end

    if Sys.isapple()
        whole_archive = "-Wl,-all_load"
        no_whole_archive = ""
    else
        whole_archive = "-Wl,--whole-archive"
        no_whole_archive = "-Wl,--no-whole-archive"
    end

    sysimg_file = "$sysimg_path.$(Libdl.dlext)"
    cc_cmd = `$cc $FLAGS -o $sysimg_path.tmp $whole_archive $sysimg_path.o $no_whole_archive`
    @info("Linking sys.$(Libdl.dlext) with $cc_cmd")
    # Windows has difficulties overwriting a file in use so we first link to a temp file
    if success(pipeline(cc_cmd; stdout=stdout, stderr=stderr))
        mv("$sysimg_path.tmp", sysimg_file; force=true)
    end
    @info("System image successfully built at $sysimg_path.$(Libdl.dlext)")
    return
end

# When running this file as a script, try to do so with default values. If arguments are passed
# in, use them as the arguments to build_sysimg above.
# Also check whether we are running `genstdlib.jl`, in which case we don't want to build a
# system image and instead only need `build_sysimg`'s docstring to be available.
if !isdefined(Main, :GenStdLib) && !isinteractive()
    if length(ARGS) > 5 || ("--help" in ARGS || "-h" in ARGS)
        println("Usage: build_sysimg.jl <sysimg_path> <cpu_target> <usrimg_path.jl> [--force] [--debug] [--help]")
        println("   <sysimg_path>    is an absolute, extensionless path to store the system image at")
        println("   <cpu_target>     is an LLVM cpu target to build the system image against")
        println("   <usrimg_path.jl> is the path to a user image to be baked into the system image")
        println("   --debug          Using julia-debug instead of julia to build the system image")
        println("   --force          Set if you wish to overwrite the default system image")
        println("   --help           Print out this help text and exit")
        println()
        println(" Example:")
        println("   build_sysimg.jl /usr/local/lib/julia/sys core2 ~/my_usrimg.jl --force")
        println()
        println(" Running this script with no arguments is equivalent to:")
        println("   build_sysimg.jl $(default_sysimg_path()) native")
        return 0
    end

    debug_flag = "--debug" in ARGS
    filter!(x -> x != "--debug", ARGS)
    force_flag = "--force" in ARGS
    filter!(x -> x != "--force", ARGS)
    build_sysimg(ARGS...; force=force_flag, debug=debug_flag)
end