text:
### species definition
macro species(typ, name, icon)
    esc(_species(typ, name, icon))
end

function _species(typ, name, icon)
    quote
        abstract type $name <: $(typ == :Animal ? AnimalSpecies : PlantSpecies) end
        Base.show(io::IO, ::Type{$name}) = print(io, $(QuoteNode(icon)))
        export $name
    end
end

macro plant(name, icon)
    return :(@species Plant $name $icon)
end

macro animal(name, icon)
    return :(@species Animal $name $icon)
end

### eating behavior
macro eats(species::Symbol, foodlist::Expr)
    return esc(_eats(species, foodlist))
end

function _generate_eat(eater::Type{<:AnimalSpecies}, food::Type{<:PlantSpecies}, multiplier)
    quote
        EcosystemCore.eats(::Animal{$(eater)}, ::Plant{$(food)}) = true
        function EcosystemCore.eat!(a::Animal{$(eater)}, p::Plant{$(food)}, w::World)
            if size(p)>0
                incr_energy!(a, $(multiplier)*size(p)*Δenergy(a))
                p.size = 0
            end
        end
    end
end

function _generate_eat(eater::Type{<:AnimalSpecies}, food::Type{<:AnimalSpecies}, multiplier)
    quote
        EcosystemCore.eats(::Animal{$(eater)}, ::Animal{$(food)}) = true
        function EcosystemCore.eat!(ae::Animal{$(eater)}, af::Animal{$(food)}, w::World)
            incr_energy!(ae, $(multiplier)*energy(af)*Δenergy(ae))
            kill_agent!(af, w)
        end
    end
end

_parse_eats(ex) = Dict(arg.args[2] => arg.args[3] for arg in ex.args if arg.head == :call && arg.args[1] == :(=>))

function _eats(species, foodlist)
    cfg = _parse_eats(foodlist)
    code = Expr(:block)
    for (k,v) in cfg
        push!(code.args, _generate_eat(eval(species), eval(k), v))
    end
    code
end


text:
__precompile__(true)

module BayesianNonparametricStatistics

using Distributions, PDMats, PyPlot

import Base: length, step

import Distributions: rand, ZeroVector

import PyPlot.plot

export AbstractModel, SDEModel, AbstractSDE, AbstractSamplePath,
AbstractGaussianProcess, SDE, SDEWithConstantVariance, SamplePathRange,
SamplePath, GaussianProcess, FaberSchauderExpansionWithGaussianCoefficients,
calculateposterior, fourier, faberschauder, faberschauderone,
fourierseriesprior, faberschauderprior,
calculatedependencystructurefaberschauderbasis, sumoffunctions

# model.jl implements AbstractModel and SDEModel.

# Hier beginnen met short-circuit evaluation invoegen!

include("model.jl")
# samplepath.jl implements SamplePath, SamplePathRange, AbstractSamplePath and
# extends Base.step and Base.length to the appropriate types.
include("samplepath.jl")
# basisfunctions.jl implements the Fourier and Faber-Schauder basis.
include("basisfunctions.jl")
# Implements AbstractSDE, SDE, SDEWithConstantVariance and extends
#Distribution.rand to sample an SDE sample path.
include("SDE.jl")
# Implements AbstractGaussianProcess, GaussianProcess,
# FaberSchauderExpansionWithGaussianCoefficients types and methods to calculate
# the posterior, with associated methods calculateLebesgueintegral,
# calculatestochasticintegral, calculatedependencystructurefaberschauderbasis.
include("gaussianprocess.jl")
end # End of module.

ext: jl


text:
"""
Voxel corner and edge indexing conventions

        Z
        |

        5------5------6              Extra edges not drawn
       /|            /|              -----------
      8 |           6 |              - face diagonals
     /  9          /  10                - 13: 1 to 3
    8------7------7   |                 - 14: 1 to 8
    |   |         |   |                 - 15: 1 to 6
    |   1------1--|---2  -- Y           - 16: 5 to 7
    12 /          11 /                  - 17: 2 to 7
    | 4           | 2                   - 18: 4 to 7
    |/            |/                 - body diagonal
    4------3------3                     - 19: 1 to 7

  /
 X
"""
immutable VoxelIndices{T <: Integer}
    voxCrnrPos::NTuple{8,NTuple{3,T}}
    voxEdgeCrnrs::NTuple{19, NTuple{2,T}}
    voxEdgeDir::NTuple{19,T}
    voxEdgeIx::NTuple{8,NTuple{8,T}}
    subTets::NTuple{6,NTuple{4,T}}
    tetEdgeCrnrs::NTuple{6,NTuple{2,T}}
    tetTri::NTuple{16,NTuple{6,T}}

    function VoxelIndices()
        voxCrnrPos = ((0, 0, 0),
                      (0, 1, 0),
                      (1, 1, 0),
                      (1, 0, 0),
                      (0, 0, 1),
                      (0, 1, 1),
                      (1, 1, 1),
                      (1, 0, 1))
        # the voxel IDs at either end of the tetrahedra edges, by edge ID
        voxEdgeCrnrs = ((1, 2),
                        (2, 3),
                        (4, 3),
                        (1, 4),
                        (5, 6),
                        (6, 7),
                        (8, 7),
                        (5, 8),
                        (1, 5),
                        (2, 6),
                        (3, 7),
                        (4, 8),
                        (1, 3),
                        (1, 8),
                        (1, 6),
                        (5, 7),
                        (2, 7),
                        (4, 7),
                        (1, 7))

        # direction codes:
        # 0 => +x, 1 => +y, 2 => +z,
        # 3 => +xy, 4 => +xz, 5 => +yz, 6 => +xyz
        voxEdgeDir = (1,0,1,0,1,0,1,0,2,2,2,2,3,4,5,3,4,5,6)

        # For a pair of corner IDs, the edge ID joining them
        # 0 denotes a pair with no edge
        voxEdgeIx = ((0,1,13,4,9,15,19,14),
                     (1,0,2,0,0,10,17,0),
                     (13,2,0,3,0,0,11,0),
                     (4,0,3,0,0,0,18,12),
                     (9,0,0,0,0,5,16,8),
                     (15,10,0,0,5,0,6,0),
                     (19,17,11,18,16,6,0,7),
                     (14,0,0,12,8,0,7,0))

        # voxel corners that comprise each of the six tetrahedra
        subTets = ((1,3,2,7),
                      (1,8,4,7),
                      (1,4,3,7),
                      (1,2,6,7),
                      (1,5,8,7),
                      (1,6,5,7))
        # tetrahedron corners for each edge (indices 1-4)
        tetEdgeCrnrs = ((1,2),
                        (2,3),
                        (1,3),
                        (1,4),
                        (2,4),
                        (3,4))

        # triangle cases for a given tetrahedron edge code
        tetTri = ((0,0,0,0,0,0),
                    (1,3,4,0,0,0),
                    (1,5,2,0,0,0),
                    (3,5,2,3,4,5),
                    (2,6,3,0,0,0),
                    (1,6,4,1,2,6),
                    (1,5,6,1,6,3),
                    (4,5,6,0,0,0),
                    (4,6,5,0,0,0),
                    (1,6,5,1,3,6),
                    (1,4,6,1,6,2),
                    (2,3,6,0,0,0),
                    (3,2,5,3,5,4),
                    (1,2,5,0,0,0),
                    (1,4,3,0,0,0),
                    (0,0,0,0,0,0))

        new(voxCrnrPos,
            voxEdgeCrnrs,
            voxEdgeDir,
            voxEdgeIx,
            subTets,
            tetEdgeCrnrs,
            tetTri)
    end
end
# (X,Y,Z)-coordinates for each voxel corner ID

"""
Checks if a voxel has faces. Should be false for most voxels.
This function should be made as fast as possible.
"""
function hasFaces{T<:Real}(vals::Vector{T}, iso::T)
    @inbounds v = vals[1]
    if v < iso
        @inbounds for i=2:8
            vals[i] >= iso && return true
        end
    else
        @inbounds for i=2:8
            vals[i] <  iso && return true
        end
    end
    false
end

"""
Determines which case in the triangle table we are dealing with
"""
function tetIx{T<:Real, IType <: Integer}(tIx::IType, vals::Vector{T}, iso::T, vxidx::VoxelIndices{IType})
    @inbounds v1 = vals[vxidx.subTets[tIx][1]]
    @inbounds v2 = vals[vxidx.subTets[tIx][2]]
    @inbounds v3 = vals[vxidx.subTets[tIx][3]]
    @inbounds v4 = vals[vxidx.subTets[tIx][4]]
    ifelse(v1 < iso, 1, 0) +
    ifelse(v2 < iso, 2, 0) +
    ifelse(v3 < iso, 4, 0) +
    ifelse(v4 < iso, 8, 0) + 1
end

"""
Determines a unique integer ID associated with the edge. This is used
as a key in the vertex dictionary. It needs to be both unambiguous (no
two edges get the same index) and unique (every edge gets the same ID
regardless of which of its neighboring voxels is asking for it) in order
for vertex sharing to be implemented properly.
"""
function vertId{IType <: Integer}(e::IType, x::IType, y::IType, z::IType,
                nx::IType, ny::IType, vxidx::VoxelIndices{IType})
    @inbounds dx, dy, dz = vxidx.voxCrnrPos[vxidx.voxEdgeCrnrs[e][1]]
    vxidx.voxEdgeDir[e]+7*(x-1+dx+nx*(y-1+dy+ny*(z-1+dz)))
end

"""
Assuming an edge crossing, determines the point in space at which it
occurs.
eps represents the "bump" factor to keep vertices away from voxel
corners (thereby preventing degeneracies).
"""
function vertPos{T<:Real, IType <: Integer}(e::IType, x::IType, y::IType, z::IType,
                          vals::Vector{T}, iso::T, eps::T, vxidx::VoxelIndices{IType})

    @inbounds ixs     = vxidx.voxEdgeCrnrs[e]
    @inbounds srcVal  = vals[ixs[1]]
    @inbounds tgtVal  = vals[ixs[2]]
    a       = (iso-srcVal)/(tgtVal-srcVal)
    a       = min(max(a, eps), one(T)-eps)
    b       = one(T)-a
    @inbounds c1x,c1y,c1z = vxidx.voxCrnrPos[ixs[1]]
    @inbounds c2x,c2y,c2z = vxidx.voxCrnrPos[ixs[2]]

    Point{3,T}(
          x+b*c1x+a*c2x,
          y+b*c1y+a*c2y,
          z+b*c1z+a*c2z
    )
end

"""
Gets the vertex ID, adding it to the vertex dictionary if not already
present.
"""
function getVertId{T<:Real, IType <: Integer}(e::IType, x::IType, y::IType, z::IType,
                            nx::IType, ny::IType,
                            vals::Vector{T}, iso::T,
                            vts::Dict{IType, Point{3,T}},
                            eps::T, vxidx::VoxelIndices{IType})

    vId = vertId(e, x, y, z, nx, ny, vxidx)
    if !haskey(vts, vId)
        vts[vId] = vertPos(e, x, y, z, vals, iso, eps, vxidx)
    end
    vId
end

"""
Given a sub-tetrahedron case and a tetrahedron edge ID, determines the
corresponding voxel edge ID.
"""
function voxEdgeId{IType <: Integer}(subTetIx::IType, tetEdgeIx::IType, vxidx::VoxelIndices{IType})
    @inbounds srcVoxCrnr::IType = vxidx.subTets[subTetIx][vxidx.tetEdgeCrnrs[tetEdgeIx][1]]
    @inbounds tgtVoxCrnr::IType = vxidx.subTets[subTetIx][vxidx.tetEdgeCrnrs[tetEdgeIx][2]]
    @inbounds v = vxidx.voxEdgeIx[srcVoxCrnr][tgtVoxCrnr]
    v
end

"""
Processes a voxel, adding any new vertices and faces to the given
containers as necessary.
"""
function procVox{T<:Real, IType <: Integer}(vals::Vector{T}, iso::T,
                          x::IType, y::IType, z::IType,
                          nx::IType, ny::IType,
                          vts::Dict{IType, Point{3,T}}, fcs::Vector{Face{3,IType,0}},
                          eps::T, vxidx::VoxelIndices{IType})

    # check each sub-tetrahedron in the voxel
    for i::IType = 1:6
        tIx = tetIx(i, vals, iso, vxidx)
        for j::IType in 1:3:4
            @inbounds e1 = vxidx.tetTri[tIx][j]
            # bail if there are no more faces
            e1 == 0 && break
            @inbounds e2 = vxidx.tetTri[tIx][j+1]
            @inbounds e3 = vxidx.tetTri[tIx][j+2]

            # add the face to the list
            push!(fcs, Face{3,IType,0}(
                      getVertId(voxEdgeId(i, e1, vxidx), x, y, z, nx, ny, vals, iso, vts, eps, vxidx),
                      getVertId(voxEdgeId(i, e2, vxidx), x, y, z, nx, ny, vals, iso, vts, eps, vxidx),
                      getVertId(voxEdgeId(i, e3, vxidx), x, y, z, nx, ny, vals, iso, vts, eps, vxidx)))
        end
    end
end

"""
Given a 3D array and an isovalue, extracts a mesh represention of the
an approximate isosurface by the method of marching tetrahedra.
"""
function marchingTetrahedra{T<:Real, IT <: Integer}(lsf::AbstractArray{T,3}, iso::T, eps::T, indextype::Type{IT})
    vts        = Dict{indextype, Point{3,T}}()
    fcs        = Array(Face{3,indextype,0}, 0)
    sizehint!(vts, div(length(lsf),8))
    sizehint!(fcs, div(length(lsf),4))
    const vxidx = VoxelIndices{indextype}()
    # process each voxel
    (nx::indextype,ny::indextype,nz::indextype) = size(lsf)
    vals = zeros(T, 8)
    for k::indextype = 1:nz-1, j::indextype = 1:ny-1, i::indextype = 1:nx-1
        for l::indextype=1:8
            @inbounds vals[l] = lsf[i+vxidx.voxCrnrPos[l][1], j+vxidx.voxCrnrPos[l][2], k+vxidx.voxCrnrPos[l][3]]
        end
        if hasFaces(vals,iso)
            procVox(vals, iso, i, j, k, nx, ny, vts, fcs, eps, vxidx)
        end
    end

    (vts,fcs)
end

function isosurface(lsf, isoval, eps, indextype=Int, index_start=one(Int))
    # get marching tetrahedra version of the mesh
    (vts, fcs) = marchingTetrahedra(lsf, isoval, eps, indextype)
    # normalize the mesh representation
    vtD = Dict{indextype,indextype}()
    sizehint!(vtD, length(vts))
    k = index_start
    for x in keys(vts)
        vtD[x] = k
        k += one(indextype)
    end
    fcAry = Face{3,indextype, index_start-1}[Face{3,indextype, index_start-1}(vtD[f[1]], vtD[f[2]], vtD[f[3]]) for f in fcs]
    vtAry = collect(values(vts))

    (vtAry, fcAry)
end

isosurface(lsf,isoval) = isosurface(lsf,isoval, convert(eltype(lsf), 0.001))

@compat function (::Type{MT}){MT <: AbstractMesh, T}(volume::Array{T, 3}, iso_val::Real, eps_val=0.001)
    iso_val = convert(T, iso_val)
    eps_val = convert(T, eps_val)
    vts, fcs = isosurface(volume, iso_val, eps_val)
    MT(vts, fcs)
end

@compat function (::Type{MT}){MT <: AbstractMesh}(df::SignedDistanceField, eps_val=0.001)
    vts, fcs = isosurface(df.data, 0.0, eps_val)
    MT(vts, fcs)
end


text:
# Types
struct SGD end
struct OJA end
struct CCIPCA end
struct GD end
struct RSGD end
struct SVRG end
struct RSVRG end
struct ALGORITHM971 end
struct HALKO end
struct ORTHITER end
struct RBKITER end
struct SINGLEPASS end
struct SINGLEPASS2 end
struct ARNOLDI end
struct LANCZOS end
struct TENXPCA end

struct ROBBINS_MONRO end
struct MOMENTUM end
struct NAG end
struct ADAGRAD end

# Row / Column
function tenxnm(tenxfile::AbstractString, group::AbstractString)
    N, M = h5open(tenxfile, "r") do file
        read(file[group*"/shape"])
    end
    return N[], M[]
end

# Index Pointer
function indptr(tenxfile::AbstractString, group::AbstractString)
    h5open(tenxfile, "r") do file
        read(file[group*"/indptr"]) .+ 1
    end
end

function loadchromium(tenxfile, group, idp, startp, endp, M, perm)
    # new objects
    newidp = ones(Int64, 1)
    newindices = zeros(Int64, 0)
    newdata = zeros(Int64, 0)
    file = h5open(tenxfile, "r")
    # Each column
    l = length(idp) - 1
    grp = file[group]
    progress = Progress(l)
    for i=1:l
        lo = idp[i]
        hi = idp[i+1]
        # Extract
        extindices = grp["indices"][lo:hi-1] .+ 1
        extdata = grp["data"][lo:hi-1]
        orderIndices = sortperm(extindices)
        extindices = extindices[orderIndices]
        extdata = extdata[orderIndices]
        # e.g. 325
        lower = searchsortedfirst(extindices, startp)
        # e.g. 663
        upper = searchsortedlast(extindices, endp)
        # update
        append!(newindices, extindices[lower:upper])
        append!(newdata, extdata[lower:upper])
        append!(newidp, length(newdata)+1)
        # Progress Bar
        next!(progress)
    end
    close(file)
    @assert minimum(newindices) >= startp
    @assert maximum(newindices) <= endp
    newindices = newindices .- (startp - 1)
    if perm
    counts = SparseMatrixCSC(endp - startp + 1, M, newidp, newindices, newdata[randperm(length(newdata))])
    else
    counts = SparseMatrixCSC(endp - startp + 1, M, newidp, newindices, newdata)
    end
    SparseArrays.sortSparseMatrixCSC!(counts)
    return counts
end

function sparseLog10(X::SparseMatrixCSC)
    m = X.m
    n = X.n
    colptr = X.colptr
    rowval = X.rowval
    nzval = log10.(X.nzval .+ 1)
    return SparseMatrixCSC(m, n, colptr, rowval, nzval)
end

# Check NaN value (only GD)
function checkNaN(W::AbstractArray, pca::GD)
    if any(isnan, W)
        error("NaN values are generated. Select other stepsize")
    end
end

# Check NaN value (other PCA)
function checkNaN(N::Number, s::Number, n::Number, W::AbstractArray, evalfreq::Number, pca::Union{OJA,SGD,CCIPCA,RSGD,SVRG,RSVRG})
    if mod((N*(s-1)+n), evalfreq) == 0
        if any(isnan, W)
            error("NaN values are generated. Select other stepsize")
        end
    end
end

# Output the result of PCA
function output(outdir::AbstractString, out::Tuple, expvar::Number)
    writecsv(joinpath(outdir, "Eigen_vectors.csv"), out[1])
    writecsv(joinpath(outdir, "Eigen_values.csv"), out[2])
    writecsv(joinpath(outdir, "Loadings.csv"), out[3])
    writecsv(joinpath(outdir, "Scores.csv"), out[4])
    writecsv(joinpath(outdir, "ExpVar.csv"), out[5])
    if out[5] > expvar && out[6] == 1
        touch(joinpath(outdir, "Converged"))
    end
end

writecsv(filename::AbstractString, data) = writedlm(filename, data, ',')
readcsv(filename::AbstractString) = readdlm(filename, ',')
readcsv(filename::AbstractString, ::Type{T}) where {T} = readdlm(filename, ',', T)

# Parse command line options (TENXPCA)
function parse_commandline(pca::TENXPCA)
    s = ArgParseSettings()

    @add_arg_table s begin
        "--tenxfile"
            help = "10XHDF5 file"
            arg_type = AbstractString
            required = true
        "--outdir", "-o"
            help = "The directory specified the directory you want to save the result."
            arg_type = AbstractString
            default = "."
            required = false
        "--scale"
            help = "{sqrt,log,raw}-scaling of the value."
            arg_type = AbstractString
            default = "sqrt"
        "--rowmeanlist", "-m"
            help = "The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--rowvarlist", "-v"
            help = "The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--colsumlist"
            help = "The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--dim", "-d"
            help = "The number of dimension of PCA."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--noversamples"
            help = "The number of over-sampling."
            arg_type = Union{Number,AbstractString}
            default = 5
        "--niter"
            help = "The number of power interation."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--chunksize"
            help = "The number of rows reading at once (e.g. 5000)."
            arg_type = Union{Number,AbstractString}
            default = 5000
        "--group"
            help = "The group name of 10XHDF5 (e.g. mm10)."
            arg_type = AbstractString
            required = true
        "--initW"
            help = "The CSV file saving the initial values of eigenvectors."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--initV"
            help = "The CSV file saving the initial values of loadings."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--logdir", "-l"
            help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--perm"
            help = "Whether the data matrix is shuffled at random"
            arg_type = Union{Bool,AbstractString}
            default = false
        "--cper"
            help = "Count per X (e.g. CPM: Count per million <1e+6>)"
            arg_type = Union{Number,AbstractString}
            default = 1f0
    end

    return parse_args(s)
end

# Parse command line options (ALGORITHM971, HALKO, SINGLEPASS, SINGLEPASS2)
function parse_commandline(pca::Union{ALGORITHM971,HALKO,SINGLEPASS,SINGLEPASS2})
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input", "-i"
            help = "Julia Binary file generated by `OnlinePCA.csv2bin` function."
            arg_type = AbstractString
            required = true
        "--outdir", "-o"
            help = "The directory specified the directory you want to save the result."
            arg_type = AbstractString
            default = "."
            required = false
        "--scale"
            help = "{log,ftt,raw}-scaling of the value."
            arg_type = AbstractString
            default = "ftt"
        "--pseudocount", "-p"
            help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
            arg_type = Union{Number,AbstractString}
            default = 1f0
        "--rowmeanlist", "-m"
            help = "The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--rowvarlist", "-v"
            help = "The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--colsumlist"
            help = "The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--dim", "-d"
            help = "The number of dimension of PCA."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--noversamples"
            help = "The number of over-sampling."
            arg_type = Union{Number,AbstractString}
            default = 5
        "--niter"
            help = "The number of power interation."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--initW"
            help = "The CSV file saving the initial values of eigenvectors."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--initV"
            help = "The CSV file saving the initial values of loadings."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--logdir", "-l"
            help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--perm"
            help = "Whether the data matrix is shuffled at random"
            arg_type = Union{Bool,AbstractString}
            default = false
        "--cper"
            help = "Count per X (e.g. CPM: Count per million <1e+6>)"
            arg_type = Union{Number,AbstractString}
            default = 1f0
    end

    return parse_args(s)
end

# Parse command line options (only ARNOLDI,LANCZOS)
function parse_commandline(pca::Union{ARNOLDI,LANCZOS})
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input", "-i"
            help = "Julia Binary file generated by `OnlinePCA.csv2bin` function."
            arg_type = AbstractString
            required = true
        "--outdir", "-o"
            help = "The directory specified the directory you want to save the result."
            arg_type = AbstractString
            default = "."
            required = false
        "--scale"
            help = "{log,ftt,raw}-scaling of the value."
            arg_type = AbstractString
            default = "ftt"
        "--pseudocount", "-p"
            help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
            arg_type = Union{Number,AbstractString}
            default = 1f0
        "--rowmeanlist", "-m"
            help = "The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--rowvarlist", "-v"
            help = "The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--colsumlist"
            help = "The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--dim", "-d"
            help = "The number of dimension of PCA."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--numepoch", "-e"
            help = "The number of epoch."
            arg_type = Union{Number,AbstractString}
            default = 5
        "--perm"
            help = "Whether the data matrix is shuffled at random"
            arg_type = Union{Bool,AbstractString}
            default = false
        "--expvar"
            help = "The calculation is determined as converged when captured variance is larger than this value (0 - 1)"
            arg_type = Union{Number,AbstractString}
            default = 0.1f0
        "--cper"
            help = "Count per X (e.g. CPM: Count per million <1e+6>)"
            arg_type = Union{Number,AbstractString}
            default = 1f0
    end

    return parse_args(s)
end

# Parse command line options (only ORTHITER, RBKITER)
function parse_commandline(pca::Union{ORTHITER,RBKITER})
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input", "-i"
            help = "Julia Binary file generated by `OnlinePCA.csv2bin` function."
            arg_type = AbstractString
            required = true
        "--outdir", "-o"
            help = "The directory specified the directory you want to save the result."
            arg_type = AbstractString
            default = "."
            required = false
        "--scale"
            help = "{log,ftt,raw}-scaling of the value."
            arg_type = AbstractString
            default = "ftt"
        "--pseudocount", "-p"
            help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
            arg_type = Union{Number,AbstractString}
            default = 1f0
        "--rowmeanlist", "-m"
            help = "The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--rowvarlist", "-v"
            help = "The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--colsumlist"
            help = "The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--dim", "-d"
            help = "The number of dimension of PCA."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--numepoch", "-e"
            help = "The number of epoch."
            arg_type = Union{Number,AbstractString}
            default = 5
        "--lower"
            help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 0
        "--upper"
            help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 1.0f+38
        "--initW"
            help = "The CSV file saving the initial values of eigenvectors."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--initV"
            help = "The CSV file saving the initial values of loadings."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--logdir", "-l"
            help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--perm"
            help = "Whether the data matrix is shuffled at random"
            arg_type = Union{Bool,AbstractString}
            default = false
        "--expvar"
            help = "The calculation is determined as converged when captured variance is larger than this value (0 - 1)"
            arg_type = Union{Number,AbstractString}
            default = 0.1f0
        "--cper"
            help = "Count per X (e.g. CPM: Count per million <1e+6>)"
            arg_type = Union{Number,AbstractString}
            default = 1f0
    end

    return parse_args(s)
end

# Parse command line options (only CCIPCA)
function parse_commandline(pca::CCIPCA)
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input", "-i"
            help = "Julia Binary file generated by `OnlinePCA.csv2bin` function."
            arg_type = AbstractString
            required = true
        "--outdir", "-o"
            help = "The directory specified the directory you want to save the result."
            arg_type = AbstractString
            default = "."
            required = false
        "--scale"
            help = "{log,ftt,raw}-scaling of the value."
            arg_type = AbstractString
            default = "ftt"
        "--pseudocount", "-p"
            help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
            arg_type = Union{Number,AbstractString}
            default = 1f0
        "--rowmeanlist", "-m"
            help = "The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--rowvarlist", "-v"
            help = "The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--colsumlist"
            help = "The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--dim", "-d"
            help = "The number of dimension of PCA."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--stepsize", "-s"
            help = "The parameter used in every iteration."
            arg_type = Union{Number,AbstractString}
            default = 1.0f3
        "--numepoch", "-e"
            help = "The number of epoch."
            arg_type = Union{Number,AbstractString}
            default = 5
        "--lower"
            help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 0
        "--upper"
            help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 1.0f+38
        "--evalfreq"
            help = "Evaluation Frequency of Reconstruction Error"
            arg_type = Union{Number,AbstractString}
            default = 5000
        "--offsetStoch"
            help = "Off set value for avoding overflow when calculating stochastic gradient"
            arg_type = Union{Number,AbstractString}
            default = 1f-15
        "--initW"
            help = "The CSV file saving the initial values of eigenvectors."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--initV"
            help = "The CSV file saving the initial values of loadings."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--logdir", "-l"
            help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--perm"
            help = "Whether the data matrix is shuffled at random"
            arg_type = Union{Bool,AbstractString}
            default = false
        "--expvar"
            help = "The calculation is determined as converged when captured variance is larger than this value (0 - 1)"
            arg_type = Union{Number,AbstractString}
            default = 0.1f0
        "--cper"
            help = "Count per X (e.g. CPM: Count per million <1e+6>)"
            arg_type = Union{Number,AbstractString}
            default = 1f0
    end

    return parse_args(s)
end

# Parse command line options (SGD、RSGD、SVRG、RSVRG)
function parse_commandline(pca::Union{SGD,RSGD,SVRG,RSVRG})
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input", "-i"
            help = "Julia Binary file generated by `OnlinePCA.csv2bin` function."
            arg_type = AbstractString
            required = true
        "--outdir", "-o"
            help = "The directory specified the directory you want to save the result."
            arg_type = AbstractString
            default = "."
            required = false
        "--scale"
            help = "{log,ftt,raw}-scaling of the value."
            arg_type = AbstractString
            default = "ftt"
        "--pseudocount", "-p"
            help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
            arg_type = Union{Number,AbstractString}
            default = 1f0
        "--rowmeanlist", "-m"
            help = "The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--rowvarlist", "-v"
            help = "The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--colsumlist"
            help = "The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--dim", "-d"
            help = "The number of dimension of PCA."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--stepsize", "-s"
            help = "The parameter used in every iteration."
            arg_type = Union{Number,AbstractString}
            default = 1.0f3
        "--numbatch", "-b"
            help = "The number of batch size."
            arg_type = Union{Number,AbstractString}
            default = 100
        "--numepoch", "-e"
            help = "The number of epoch."
            arg_type = Union{Number,AbstractString}
            default = 5
        "--scheduling"
            help = "Learning parameter scheduling. `robbins-monro`, `momentum`, `nag`, and `adagrad` are available."
            arg_type = AbstractString
            default = "robbins-monro"
        "-g"
            help = "The parameter that is used when scheduling is specified as nag."
            arg_type = Union{Number,AbstractString}
            default = 0.9f0
        "--epsilon"
            help = "The parameter that is used when scheduling is specified as adagrad."
            arg_type = Union{Number,AbstractString}
            default = 1.0f-8
        "--lower"
            help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 0
        "--upper"
            help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 1.0f+38
        "--evalfreq"
            help = "Evaluation Frequency of Reconstruction Error"
            arg_type = Union{Number,AbstractString}
            default = 5000
        "--offsetFull"
            help = "Off set value for avoding overflow when calculating full gradient"
            arg_type = Union{Number,AbstractString}
            default = 1.0f-20
        "--offsetStoch"
            help = "Off set value for avoding overflow when calculating stochastic gradient"
            arg_type = Union{Number,AbstractString}
            default = 1.0f-6
        "--initW"
            help = "The CSV file saving the initial values of eigenvectors."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--initV"
            help = "The CSV file saving the initial values of loadings."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--logdir", "-l"
            help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--perm"
            help = "Whether the data matrix is shuffled at random"
            arg_type = Union{Bool,AbstractString}
            default = false
        "--expvar"
            help = "The calculation is determined as converged when captured variance is larger than this value (0 - 1)"
            arg_type = Union{Number,AbstractString}
            default = 0.1f0
        "--cper"
            help = "Count per X (e.g. CPM: Count per million <1e+6>)"
            arg_type = Union{Number,AbstractString}
            default = 1f0
    end

    return parse_args(s)
end

# Parse command line options (other PCA)
function parse_commandline(pca::Union{OJA,GD})
    s = ArgParseSettings()

    @add_arg_table s begin
        "--input", "-i"
            help = "Julia Binary file generated by `OnlinePCA.csv2bin` function."
            arg_type = AbstractString
            required = true
        "--outdir", "-o"
            help = "The directory specified the directory you want to save the result."
            arg_type = AbstractString
            default = "."
            required = false
        "--scale"
            help = "{log,ftt,raw}-scaling of the value."
            arg_type = AbstractString
            default = "ftt"
        "--pseudocount", "-p"
            help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
            arg_type = Union{Number,AbstractString}
            default = 1f0
        "--rowmeanlist", "-m"
            help = "The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--rowvarlist", "-v"
            help = "The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--colsumlist"
            help = "The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions."
            arg_type = AbstractString
            default = ""
            required = false
        "--dim", "-d"
            help = "The number of dimension of PCA."
            arg_type = Union{Number,AbstractString}
            default = 3
        "--stepsize", "-s"
            help = "The parameter used in every iteration."
            arg_type = Union{Number,AbstractString}
            default = 1.0f3
        "--numepoch", "-e"
            help = "The number of epoch."
            arg_type = Union{Number,AbstractString}
            default = 5
        "--scheduling"
            help = "Learning parameter scheduling. `robbins-monro`, `momentum`, `nag`, and `adagrad` are available."
            arg_type = AbstractString
            default = "robbins-monro"
        "-g"
            help = "The parameter that is used when scheduling is specified as nag."
            arg_type = Union{Number,AbstractString}
            default = 0.9f0
        "--epsilon"
            help = "The parameter that is used when scheduling is specified as adagrad."
            arg_type = Union{Number,AbstractString}
            default = 1.0f-8
        "--lower"
            help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 0
        "--upper"
            help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
            arg_type = Union{Number,AbstractString}
            default = 1.0f+38
        "--evalfreq"
            help = "Evaluation Frequency of Reconstruction Error"
            arg_type = Union{Number,AbstractString}
            default = 5000
        "--offsetFull"
            help = "Off set value for avoding overflow when calculating full gradient"
            arg_type = Union{Number,AbstractString}
            default = 1.0f-20
        "--offsetStoch"
            help = "Off set value for avoding overflow when calculating stochastic gradient"
            arg_type = Union{Number,AbstractString}
            default = 1.0f-6
        "--initW"
            help = "The CSV file saving the initial values of eigenvectors."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--initV"
            help = "The CSV file saving the initial values of loadings."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--logdir", "-l"
            help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
            arg_type = Union{Nothing,AbstractString}
            default = nothing
        "--perm"
            help = "Whether the data matrix is shuffled at random"
            arg_type = Union{Bool,AbstractString}
            default = false
        "--expvar"
            help = "The calculation is determined as converged when captured variance is larger than this value (0 - 1)"
            arg_type = Union{Number,AbstractString}
            default = 0.1f0
        "--cper"
            help = "Count per X (e.g. CPM: Count per million <1e+6>)"
            arg_type = Union{Number,AbstractString}
            default = 1f0
    end

    return parse_args(s)
end

# Return N, M
function nm(input::AbstractString)
    N = zeros(UInt32, 1)
    M = zeros(UInt32, 1)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, N)
        read!(stream, M)
        close(stream)
    end
    return N[], M[]
end

# Initialization (ALGORITHM971, HALKO, SINGLEPASS, SINGLEPASS2)
function init(input::AbstractString, pseudocount::Number, dim::Number, rowmeanlist::AbstractString, rowvarlist::AbstractString, colsumlist::AbstractString, initW::Union{Nothing,AbstractString}, initV::Union{Nothing,AbstractString}, logdir::Union{Nothing,AbstractString}, pca::Union{ALGORITHM971,HALKO,SINGLEPASS,SINGLEPASS2}, cper::Number, scale::AbstractString="ftt")
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    pseudocount = Float32(pseudocount)
    # Eigen vectors
    # Eigen vectors
    if initW == nothing
        W = zeros(Float32, M, dim)
        for i=1:dim
            W[i, i] = 1
        end
    end
    if typeof(initW) == String
        if initV == nothing
            W = readcsv(initW, Float32)
        else
            error("initW and initV are not specified at once. You only have one choice.")
        end
    end
    if typeof(initV) == String
            V = readcsv(initV, Float32)
            V = V[:,1:dim]
    end
    D = Diagonal(reverse(1:dim)) # Diagonal Matrix
    x = zeros(UInt32, M)
    rowmeanvec = zeros(Float32, N, 1)
    rowvarvec = zeros(Float32, N, 1)
    colsumvec = zeros(Float32, M, 1)
    if rowmeanlist != ""
        rowmeanvec = readcsv(rowmeanlist, Float32)
    end
    if rowvarlist != ""
        rowvarvec = readcsv(rowvarlist, Float32)
    end
    if colsumlist != ""
        colsumvec = readcsv(colsumlist, Float32)
    end
    # N, M, All Variance
    TotalVar = 0
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            TotalVar = TotalVar + normx'normx
            # Creating W from V
            if typeof(initV) == String
                W = W .+ (V[n,:]*normx')'
            end
        end
        close(stream)
    end
    TotalVar = TotalVar / M
    # directory for log file
    if logdir isa String
        if(!isdir(logdir))
            mkdir(logdir)
        end
    end
    return pseudocount, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar
end

# Initialization (only ARNOLDI,LANCZOS)
function init(input::AbstractString, pseudocount::Number, dim::Number, rowmeanlist::AbstractString, rowvarlist::AbstractString, colsumlist::AbstractString, pca::Union{ARNOLDI,LANCZOS}, numepoch::Number, cper::Number, scale::AbstractString="ftt")
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    pseudocount = Float32(pseudocount)
    # Eigen vectors
    W = rand(Float32, M, numepoch+1)
    # Normalization
    for i=1:(numepoch+1)
        W[:, i] = W[:, i] / norm(W[:, i])
    end
    X = zeros(Float32, M, dim+1) # Temporal Vector (Same length
    D = Diagonal(reverse(1:dim)) # Diagonaml Matrix
    x = zeros(UInt32, M)
    rowmeanvec = zeros(Float32, N, 1)
    rowvarvec = zeros(Float32, N, 1)
    colsumvec = zeros(Float32, M, 1)
    if rowmeanlist != ""
        rowmeanvec = readcsv(rowmeanlist, Float32)
    end
    if rowvarlist != ""
        rowvarvec = readcsv(rowvarlist, Float32)
    end
    if colsumlist != ""
        colsumvec = readcsv(colsumlist, Float32)
    end
    # N, M, All Variance
    TotalVar = 0
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            TotalVar = TotalVar + normx'normx
        end
        close(stream)
    end
    TotalVar = TotalVar / M
    return pseudocount, W, X, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar
end

# Initialization (ORTHITER, RBKITER)
function init(input::AbstractString, pseudocount::Number, dim::Number, rowmeanlist::AbstractString, rowvarlist::AbstractString, colsumlist::AbstractString, initW::Union{Nothing,AbstractString}, initV::Union{Nothing,AbstractString}, logdir::Union{Nothing,AbstractString}, pca::Union{ORTHITER,RBKITER}, lower::Number, upper::Number, cper::Number, scale::AbstractString="ftt")
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    pseudocount = Float32(pseudocount)
    lower = Float32(lower)
    upper = Float32(upper)
    # Eigen vectors
    if initW == nothing
        W = rand(Float32, M, dim)
            # Normalization
            for i=1:dim
                W[:, i] = W[:, i] / norm(W[:, i])
            end
    end
    if typeof(initW) == String
        if initV == nothing
            W = readcsv(initW, Float32)
        else
            error("initW and initV are not specified at once. You only have one choice.")
        end
    end
    if typeof(initV) == String
            V = readcsv(initV, Float32)
            V = V[:,1:dim]
    end
    X = zeros(Float32, M, dim+1) # Temporal Vector (Same length
    D = Diagonal(reverse(1:dim)) # Diagonaml Matrix
    x = zeros(UInt32, M)
    rowmeanvec = zeros(Float32, N, 1)
    rowvarvec = zeros(Float32, N, 1)
    colsumvec = zeros(Float32, M, 1)
    if rowmeanlist != ""
        rowmeanvec = readcsv(rowmeanlist, Float32)
    end
    if rowvarlist != ""
        rowvarvec = readcsv(rowvarlist, Float32)
    end
    if colsumlist != ""
        colsumvec = readcsv(colsumlist, Float32)
    end
    # N, M, All Variance
    TotalVar = 0
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            TotalVar = TotalVar + normx'normx
            # Creating W from V
            if typeof(initV) == String
                W = W .+ (V[n,:]*normx')'
            end
        end
        close(stream)
    end
    TotalVar = TotalVar / M
    # directory for log file
    if logdir isa String
        if(!isdir(logdir))
            mkdir(logdir)
        end
    end
    return pseudocount, W, X, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper
end

# Initialization (only CCIPCA)
function init(input::AbstractString, pseudocount::Number, stepsize::Number, dim::Number, rowmeanlist::AbstractString, rowvarlist::AbstractString, colsumlist::AbstractString, initW::Union{Nothing,AbstractString}, initV::Union{Nothing,AbstractString}, logdir::Union{Nothing,AbstractString}, pca::CCIPCA, lower::Number, upper::Number, evalfreq::Number, offsetStoch::Number, cper::Number, scale::AbstractString="ftt")
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    pseudocount = Float32(pseudocount)
    stepsize = Float32(stepsize)
    lower = Float32(lower)
    upper = Float32(upper)
    evalfreq = Int64(evalfreq)
    offsetStoch = Float32(offsetStoch)
    # Eigen vectors
    if initW == nothing
        W = zeros(Float32, M, dim)
        for i=1:dim
            W[i, i] = 1
        end
    end
    if typeof(initW) == String
        if initV == nothing
            W = readcsv(initW, Float32)
        else
            error("initW and initV are not specified at once. You only have one choice.")
        end
    end
    if typeof(initV) == String
            V = readcsv(initV, Float32)
            V = V[:,1:dim]
    end
    X = zeros(Float32, M, dim+1) # Temporal Vector (Same length
    D = Diagonal(reverse(1:dim)) # Diagonaml Matrix
    x = zeros(UInt32, M)
    rowmeanvec = zeros(Float32, N, 1)
    rowvarvec = zeros(Float32, N, 1)
    colsumvec = zeros(Float32, M, 1)
    if rowmeanlist != ""
        rowmeanvec = readcsv(rowmeanlist, Float32)
    end
    if rowvarlist != ""
        rowvarvec = readcsv(rowvarlist, Float32)
    end
    if colsumlist != ""
        colsumvec = readcsv(colsumlist, Float32)
    end
    # N, M, All Variance
    TotalVar = 0
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            TotalVar = TotalVar + normx'normx
            # Creating W from V
            if typeof(initV) == String
                W = W .+ (V[n,:]*normx')'
            end
        end
        close(stream)
    end
    TotalVar = TotalVar / M
    # directory for log file
    if logdir isa String
        if(!isdir(logdir))
            mkdir(logdir)
        end
    end
    return pseudocount, stepsize, W, X, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetStoch
end

# Initialization (other PCA)
function init(input::AbstractString, pseudocount::Number, stepsize::Number, g::Number, epsilon::Number, dim::Number, rowmeanlist::AbstractString, rowvarlist::AbstractString, colsumlist::AbstractString, initW::Union{Nothing,AbstractString}, initV::Union{Nothing,AbstractString}, logdir::Union{Nothing,AbstractString}, pca::Union{OJA,SGD,GD,RSGD,SVRG,RSVRG}, lower::Number, upper::Number, evalfreq::Number, offsetFull::Number, offsetStoch::Number, cper::Number, scale::AbstractString="ftt")
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    pseudocount = Float32(pseudocount)
    stepsize = Float32(stepsize)
    g = Float32(g)
    epsilon = Float32(epsilon)
    lower = Float32(lower)
    upper = Float32(upper)
    evalfreq = Int64(evalfreq)
    offsetFull = Float32(offsetFull)
    offsetStoch = Float32(offsetStoch)
    # Eigen vectors
    if initW == nothing
        W = zeros(Float32, M, dim)
        for i=1:dim
            W[i, i] = 1
        end
    end
    if typeof(initW) == String
        if initV == nothing
            W = readcsv(initW, Float32)
        else
            error("initW and initV are not specified at once. You only have one choice.")
        end
    end
    if typeof(initV) == String
            V = readcsv(initV, Float32)
            V = V[:,1:dim]
    end
    v = zeros(Float32, M, dim) # Temporal Vector (Same length
    D = Diagonal(reverse(1:dim)) # Diagonaml Matrix
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    rowmeanvec = zeros(Float32, N, 1)
    rowvarvec = zeros(Float32, N, 1)
    colsumvec = zeros(Float32, M, 1)
    if rowmeanlist != ""
        rowmeanvec = readcsv(rowmeanlist, Float32)
    end
    if rowvarlist != ""
        rowvarvec = readcsv(rowvarlist, Float32)
    end
    if colsumlist != ""
        colsumvec = readcsv(colsumlist, Float32)
    end
    # N, M, All Variance
    TotalVar = 0
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            TotalVar = TotalVar + normx'normx
            # Creating W from V
            if typeof(initV) == String
                W = W .+ (V[n,:]*normx')'
            end
        end
        close(stream)
    end
    TotalVar = TotalVar / M
    # directory for log file
    if logdir isa String
        if !isdir(logdir)
            mkdir(logdir)
        end
    end
    return pseudocount, stepsize, g, epsilon, W, v, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, lower, upper, evalfreq, offsetFull, offsetStoch
end

# Eigen value, Loading, Scores
function WλV(W::AbstractArray, input::AbstractString, dim::Number, scale::AbstractString, pseudocount::Number, rowmeanlist::AbstractString, rowmeanvec::AbstractArray, rowvarlist::AbstractString, rowvarvec::AbstractArray, colsumlist::AbstractString, colsumvec::AbstractArray, TotalVar::Number, cper::Number)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    V = zeros(Float32, N, dim)
    Scores = zeros(Float32, M, dim)
    x = zeros(UInt32, M)
    normx = zeros(UInt32, M)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            V[n, :] = normx'W
        end
        close(stream)
    end
    # Eigen value
    σ = Float32[norm(V[:, x]) for x=1:dim]
    for n = 1:dim
        V[:, n] ./= σ[n]
    end
    λ = σ .* σ ./ M

    # Sort by Eigen value
    idx = sortperm(λ, rev=true)
    W .= W[:, idx]
    λ .= λ[idx]
    V .= V[:, idx]
    for n = 1:dim
        Scores[:, n] .= λ[n] .* W[:, n]
    end
    ExpVar = sum(λ) / TotalVar
    # Return
    return W, λ, V, Scores, ExpVar, TotalVar
end

# Output log file （only GD）
function outputlog(s::Number, input::AbstractString, dim::Number, logdir::AbstractString, W::AbstractArray, pca::GD, TotalVar::Number, scale::AbstractString, pseudocount::Number, rowmeanlist::AbstractString, rowmeanvec::AbstractArray, rowvarlist::AbstractString, rowvarvec::AbstractArray, colsumlist::AbstractString, colsumvec::AbstractArray, lower::Number, upper::Number, stop::Number, cper::Number)
    REs = RecError(W, input, TotalVar, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
    if s != 1
        old_E = readcsv(joinpath(logdir, "RecError_Epoch"*string(s-1)*".csv"))
        RelChange = abs(REs[1][2] - old_E[1,2]) / REs[1][2]
        REs = [REs[1], REs[2], REs[3], REs[4], REs[5], REs[6], "RelChange" => RelChange]
        if RelChange < lower
            println("Relative change of reconstruction error is below the lower value (no change)")
            stop = 1
        end
        if RelChange > upper
            println("Relative change of reconstruction error is above the upper value (unstable)")
            stop = 2
        end
    end
    writecsv(joinpath(logdir, "RecError_Epoch"*string(s)*".csv"), REs)
    writecsv(joinpath(logdir, "W_Epoch"*string(s)*".csv"), W)
    return stop
end

# Output log file (other PCA)
function outputlog(N::Number, s::Number, n::Number, input::AbstractString, dim::Number, logdir::AbstractString, W::AbstractArray, pca::Union{OJA,SGD,CCIPCA,RSGD,SVRG,RSVRG}, TotalVar::Number, scale::AbstractString, pseudocount::Number, rowmeanlist::AbstractString, rowmeanvec::AbstractArray, rowvarlist::AbstractString, rowvarvec::AbstractArray, colsumlist::AbstractString, colsumvec::AbstractArray, lower::Number, upper::Number, stop::Number, evalfreq::Number, cper::Number)
    if(mod((N*(s-1)+n), evalfreq) == 0)
        REs = RecError(W, input, TotalVar, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
        if n != evalfreq && (N*(s-1)+(n-evalfreq)) != 0
            old_E = readcsv(joinpath(logdir, "RecError_"*string((N*(s-1)+(n-evalfreq)))*".csv"))
            RelChange = abs(REs[1][2] - old_E[1,2]) / REs[1][2]
            REs = [REs[1], REs[2], REs[3], REs[4], REs[5], REs[6], "RelChange"=> RelChange]
            if RelChange < lower
                println("Relative change of reconstruction error is below the lower value (no change)")
                stop = 1
            end
            if RelChange > upper
                println("Relative change of reconstruction error is above the upper value (unstable)")
                stop = 2
            end
        end
        writecsv(joinpath(logdir, "W_"*string((N*(s-1)+n))*".csv"), W)
        writecsv(joinpath(logdir, "RecError_"*string((N*(s-1)+n))*".csv"), REs)
    end
    return stop
end

# Reconstuction Error
function RecError(W::AbstractArray, input::AbstractString, TotalVar::Number, scale::AbstractString, pseudocount::Number, rowmeanlist::AbstractString, rowmeanvec::AbstractArray, rowvarlist::AbstractString, rowvarvec::AbstractArray, colsumlist::AbstractString, colsumvec::AbstractArray, cper::Number)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    dim = size(W)[2]
    V = zeros(Float32, N, dim)
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    E = 0.0f0
    AE = 0.0f0
    RMSE = 0.0f0
    ARE = 0.0f0
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            pc = W'normx
            E = E + dot(normx, normx) - dot(pc, pc)
            V[n, :] = pc
        end
        close(stream)
    end
    # Eigen value
    σ = Float32[norm(V[:, x]) for x=1:dim]
    for n = 1:dim
        V[:, n] ./= σ[n]
    end
    λ = σ .* σ ./ M
    ExpVar = sum(λ) / TotalVar

    AE = E / M
    RMSE = sqrt(E / (N * M))
    ARE = sqrt(E / TotalVar)
    @assert E isa Float32
    @assert AE isa Float32
    @assert RMSE isa Float32
    @assert ARE isa Float32
    # Return
    return ["E"=>E, "AE"=>AE, "RMSE"=>RMSE, "ARE"=>ARE, "Explained Variance"=>ExpVar, "Total Variance"=>TotalVar]
end

# Row vector
function normalizex(x::Array{UInt32,1}, n::Number, stream, scale::AbstractString, pseudocount::Number, rowmeanlist::AbstractString, rowmeanvec::AbstractArray, rowvarlist::AbstractString, rowvarvec::AbstractArray, colsumlist::AbstractString, colsumvec::AbstractArray, cper::Number)
    # Arugment Check
    if rowmeanlist == ""
        error("No rowmeanlist is specified!")
    end
    if !(scale in ["log", "ftt", "raw"])
        error("scale must be specified as log, ftt, or raw")
    end

    # Input
    xx = convert(Vector{Float32}, x)

    # Normalization
    if colsumlist != ""
        @inbounds for i in 1:length(xx)
            xx[i] = cper * xx[i] / colsumvec[i, 1]
        end
    end

    # Variance Stabilizing Transformation
    if scale == "log"
        @inbounds for i in 1:length(xx)
            xx[i] = log10(xx[i] + pseudocount)
        end
    end
    if scale == "ftt"
        @inbounds for i in 1:length(xx)
            xx[i] = sqrt(xx[i]) + sqrt(xx[i] + 1.0f0)
        end
    end

    # Centering
    @inbounds for i in 1:length(xx)
        xx[i] = xx[i] - rowmeanvec[n, 1]
    end

    # Scaling
    if rowvarlist != ""
        @inbounds for i in 1:length(xx)
            xx[i] = xx[i] / rowvarvec[n, 1]
        end
    end

    # Return
    return xx
end

# Full Gradient
function ∇f(W::AbstractArray, input::AbstractString, D::AbstractArray, scale::AbstractString, pseudocount::Number, rowmeanlist::AbstractString, rowmeanvec::AbstractArray, rowvarlist::AbstractString, rowvarvec::AbstractArray, colsumlist::AbstractString, colsumvec::AbstractArray, stepsize::Number, offsetFull::Number, offsetStoch::Number, perm::Bool, cper::Number)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    tmpW = zeros(Float32, size(W)[1], size(W)[2])
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        for n = 1:N
            # Data Import
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            if perm
                normx .= normx[randperm(length(normx))]
            end
            # Full Gradient
            tmpW .+= offsetFull * ∇fn(W, normx, D, M, stepsize, offsetStoch)
        end
        close(stream)
    end
    return tmpW / offsetFull
end

# Stochastic Gradient
function ∇fn(W::AbstractArray, x::Array{Float32,1}, D::AbstractArray, M::Number, stepsize::Number, offsetStoch::Number)
    return 1/offsetStoch * stepsize * Float32(2 / M) * x * (offsetStoch * x'W * D)
end

# sym
function sym(Y::AbstractArray)
    return (Y + Y') / 2
end

# Riemannian Gradient
function Pw(Z::AbstractArray, W::AbstractArray)
    return Z - W * sym(W'Z)
end


text:
# __precompile__(false)
module HOQSTTutorials

using Weave, Pkg, InteractiveUtils, IJulia

repo_directory = joinpath(@__DIR__, "..")
cssfile = joinpath(@__DIR__, "..", "templates", "skeleton_css.css")
latexfile = joinpath(@__DIR__, "..", "templates", "julia_tex.tpl")

function weave_file(folder, file, build_list=(:script, :html, :pdf, :github, :notebook); kwargs...)
    tmp = joinpath(repo_directory, "tutorials", folder, file)
    Pkg.activate(dirname(tmp))
    Pkg.instantiate()
    args = Dict{Symbol,String}(:folder => folder, :file => file)
    if :script ∈ build_list
        println("Building Script")
        dir = joinpath(repo_directory, "script", folder)
        isdir(dir) || mkdir(dir)
        args[:doctype] = "script"
        tangle(tmp;out_path=dir)
    end
    if :html ∈ build_list
        println("Building HTML")
        dir = joinpath(repo_directory, "html", folder)
        isdir(dir) || mkdir(dir)
        args[:doctype] = "html"
        weave(tmp, doctype="md2html", out_path=dir, args=args; fig_ext=".svg", css=cssfile, kwargs...)
    end
    if :pdf ∈ build_list
        println("Building PDF")
        dir = joinpath(repo_directory, "pdf", folder)
        isdir(dir) || mkdir(dir)
        args[:doctype] = "pdf"
        try
            weave(tmp, doctype="md2pdf", out_path=dir, args=args; template=latexfile, kwargs...)
        catch ex
            @warn "PDF generation failed" exception = (ex, catch_backtrace())
        end
    end
    if :github ∈ build_list
        println("Building Github Markdown")
        dir = joinpath(repo_directory, "markdown", folder)
        isdir(dir) || mkdir(dir)
        args[:doctype] = "github"
        weave(tmp, doctype="github", out_path=dir, args=args; kwargs...)
    end
    if :notebook ∈ build_list
        println("Building Notebook")
        dir = joinpath(repo_directory, "notebook", folder)
        isdir(dir) || mkdir(dir)
        args[:doctype] = "notebook"
        Weave.convert_doc(tmp, joinpath(dir, file[1:end - 4] * ".ipynb"))
    end
end

function weave_all()
    for folder in readdir(joinpath(repo_directory, "tutorials"))
        folder == "test.jmd" && continue
        weave_folder(folder)
    end
end

function weave_folder(folder)
    for file in readdir(joinpath(repo_directory, "tutorials", folder))
        println("Building $(joinpath(folder, file)))")
        try
            weave_file(folder, file)
        catch
        end
    end
end

function tutorial_footer(folder=nothing, file=nothing; remove_homedir=true)
    display("text/markdown", """
    ## Appendix
     This tutorial is part of the HOQSTTutorials.jl repository, found at: <https://github.com/USCqserver/HOQSTTutorials.jl>.
    """)
    if folder !== nothing && file !== nothing
        display("text/markdown", """
        To locally run this tutorial, do the following commands:
        ```
        using HOQSTTutorials
        HOQSTTutorials.weave_file("$folder","$file")
        ```
        """)
    end
    display("text/markdown", "Computer Information:")
    vinfo = sprint(InteractiveUtils.versioninfo)
    display("text/markdown",  """
    ```
    $(vinfo)
    ```
    """)
    ctx = Pkg.API.Context()
    pkgs = Pkg.dependencies();
    projfile = ctx.env.project_file
    remove_homedir && (projfile = projfile[findlast("tutorials", projfile)[1]:end])

    display("text/markdown","""
    Package Information:
    """)

    md = ""
    md *= "```\nStatus `$(projfile)`\n"

    for (uuid, pkg) in pkgs
        if pkg.is_direct_dep == true
            md *= "[$(string(uuid))] $(string(pkg.name)) $(string(pkg.version))\n"
        end
    end
    md *= "```"
    display("text/markdown", md)
end

function open_notebooks()
    Base.eval(Main, Meta.parse("import IJulia"))
    path = joinpath(repo_directory, "notebook")
    IJulia.notebook(;dir=path)
end

end


text:
# deprecated functions

import Base: depwarn

function add!(a::Matrix, b::AbstractPDMat)
    depwarn("add! is deprecated in favor of pdadd!", :add!)
    pdadd!(a, b)
end

function add_scal!(a::Matrix, b::AbstractPDMat, c::Real)
    depwarn("add! is deprecated in favor of pdadd!", :add_scal!)
    pdadd!(a, b, c)
end

function add_scal(a::Matrix, b::AbstractPDMat, c::Real)
    depwarn("add_scal is deprecated in favor of pdadd", :add_scal)
    pdadd(a, b, c)
end


text:
##################################################

# type to throw on succesful convergence
mutable struct StateConverged
    x0::Number
end

# type to throw on failure
mutable struct ConvergenceFailed
    reason::AbstractString
end

##################################################
## Helpers for the various methods

_unitless(x) = x / oneunit(x)

# NaN of correct type
nan(::Type{Float16}) = NaN16
nan(::Type{Float32}) = NaN32
nan(::Type{Float64}) = NaN
nan(x::T) where {T <: Real} = NaN*one(T)
nan(::Any) = NaN

## issue with approx derivative
isissue(x) = iszero(x) || isnan(x) || isinf(x)

# of (a,fa), (b,fb) choose pair where |f| is smallest
@inline choose_smallest(a, b, fa, fb) = abs(fa) < abs(fb) ? (a,fa) : (b,fb)
@inline sort_smallest(a, b, fa, fb) = abs(fa) < abs(fb) ? (a, b, fa, fb) : (b, a, fb, fa)

# from an interable extract a bracketing interval and promote to floating point values
# used when an interval is desired
_extrema(x::Number) = throw(ArgumentError("Need extrema to return two distinct values"))
function _extrema(x)
    a, b = extrema(x)
    a == b && throw(ArgumentError("Need extrema to return two distinct values"))
    a, b
end
# fix type instability issues of tuples of mixed types
function _extrema(x::Tuple{<:Number,<:Number})
    a, b = x
    a == b && throw(ArgumentError("Need extrema to return two distinct values"))
    extrema(promote(a, b))
end

# used by secant. Get x₀, x₁ for x
function x₀x₁(x::Number)
    x₁ = float(x)
    promote(_default_secant_step(x₁), x₁)
end

function x₀x₁(x)
    x₀, x₁′ = x
    x₁ = first(x₁′)
    promote(float(x₀), float(x₁))
end

## find a default secant step
function _default_secant_step(x1)
    ϵ = eps(one(real(x1)))
    h = cbrt(ϵ)
    dx = h*oneunit(x1) + abs(x1)*h^2 # adjust for if eps(x1) > h
    x0 = x1 + dx
    x0
end

# a bit better than a - fa/f_ab
@inline secant_step(a, b, fa, fb) =  a - fa * (b - a) / (fb - fa)

function guarded_secant_step(alpha, beta, falpha, fbeta)

    fp = (fbeta - falpha) /  (beta - alpha)
    Δ = fbeta / fp
    ## odd, we get allocations if we define Delta, then beta - Delta
    ## Δ = beta - fbeta * (beta - alpha) / (fbeta - falpha)

    if isissue(Δ)
        Δ = oneunit(alpha)/1000
    elseif abs(Δ) >= 100 * abs(alpha - beta) # guard runaway
        Δ = sign(Δ) * 100 * min(oneunit(alpha), abs(alpha - beta))
    end

    if isissue(Δ)
        return (alpha + (beta - alpha)*(0.5), true) # midpoint
    else
        return (beta - Δ, false)
    end
end

"""
    steff_step(M, x, fx)

Return first Steffensen step x + fx (with proper units).
May be overriden to provide a guard when fx is too large.

"""
function steff_step(M::Any, x::T, fx::S) where {T, S}
    x + fx * oneunit(T) / oneunit(S)
end

# return vertex of parabola through (a,fa),(b,fb),(c,fc)
# first time trhough, we have picture of a > b > c; |fa|, |fc| > |fb|, all same sign
function quad_vertex(c,fc,b,fb,a,fa)
    fba = (fb-fa)/(b-a)
    fbc = (fb-fc)/(b-c)

    1/2 * ((a+b) - fba/(fbc - fba)*(c-a))

end

## inverse quadratic
function inverse_quadratic_step(a::T,b,c,fa,fb,fc) where {T}
    s = zero(T)
    s += a * fb * fc / (fa - fb) / (fa - fc) # quad step
    s += b * fa * fc / (fb - fa) / (fb - fc)
    s += c * fa * fb / (fc - fa) / (fc - fb)
    s
end

## Different functions for approximating f'(xn)
## return fpxn and whether it is an issue

## use f[a,b] to approximate f'(x)
function _fbracket(a, b, fa, fb)
    num, den = fb - fa, b - a
    iszero(num) && iszero(den) && return Inf, true
    out = num / den
    out, isissue(out)
end

## use f[y,z] - f[x,y] + f[x,z] to approximate
function _fbracket_diff(a,b,c, fa, fb, fc)
    x1, issue = _fbracket(b, c, fb,  fc)
    issue && return x1, issue
    x2, issue = _fbracket(a, b, fa,  fb)
    issue && return x2, issue
    x3, issue = _fbracket(a, c, fa,  fc)
    issue && return x3, issue

    out = x1 - x2 + x3
    out, isissue(out)
end

## use f[a,b] * f[a,c] / f[b,c]
function _fbracket_ratio(a, b, c, fa, fb, fc)
    x1, _ = _fbracket(a, b, fa, fb)
    x2, _ = _fbracket(a, c, fa, fc)
    x3, _ = _fbracket(b, c, fb, fc)
    out = (x1 * x2) / x3
    out, isissue(out)
end

## from https://core.ac.uk/download/pdf/82282744.pdf
## signum based iteration allows one to guarantee a valid starting point
## if N is big enough. (THough this depends on algorithm, a, b and function)
## N here would need to be tuned. But, as is, this may prove useful anyways.
function identify_starting_point(f, a, b, N)
    pts = range(a, stop=b, length=N+1)
    fxs = f.(pts)
    sfxs = sign.(f.(pts))
    identify_starting_point(a,b,sfxs)
end

function identify_starting_point(a, b, sfxs)
    N = length(sfxs) - 1
    p0 = a + (b-a)/2
    p1 = p0 + (b-a)/(2N) * sfxs[1] * sum(s for s in sfxs[2:end-1])
    p1
end

## ----

function unicode_subscript(io, j)
    a = ("⁻","","","₀","₁","₂","₃","₄","₅","₆","₇","₈","₉")
    for i in string(j)
        print(io, a[Int(i)-44])
    end
end


text:
"""
    DynamicAggregation

Functor object for using dynamic aggregation in an L-shaped algorithm. Create by supplying a [`DynamicAggregate`](@ref) object through `aggregate ` in the `LShapedSolver` factory function and then pass to a `StochasticPrograms.jl` model.

The following selection rules are available
- [`SelectUniform`](@ref)
- [`SelectDecaying`](@ref)
- [`SelectRandom`](@ref
- [`SelectClosest`](@ref)
- [`SelectClosestToReference`](@ref)

...
# Parameters
- `naggregates::Int`: Number of aggregates
- `rule::SelectionRule`: Rule that determines which aggregate an incoming cut should be placed in
- `lock_after::Function = (τ,n)->false`: Function that determines if the current aggregation scheme should be fixed, based on the current optimality gap `τ` and the number of iterations `n`
...
"""
struct DynamicAggregation{T <: AbstractFloat, S <: SelectionRule} <: AbstractAggregation
    aggregates::Vector{AggregatedOptimalityCut{T}}
    rule::S
    partitioning::Dict{Int,Int}
    lock::Function

    function DynamicAggregation(naggregates::Integer, rule::SelectionRule, lock::Function, ::Type{T}) where T <: AbstractFloat
        S = typeof(rule)
        aggregates = [zero(AggregatedOptimalityCut{T}) for _ = 1:naggregates]
        return new{T,S}(aggregates, rule, Dict{Int,Int}(), lock)
    end
end

function aggregate_cut!(lshaped::AbstractLShapedSolver, ::DynamicAggregation, cut::HyperPlane)
    return add_cut!(lshaped, cut)
end

function aggregate_cut!(lshaped::AbstractLShapedSolver, aggregation::DynamicAggregation{T}, cut::HyperPlane{OptimalityCut}) where T <: AbstractFloat
    if aggregation.lock(gap(lshaped),niterations(lshaped)) && haskey(aggregation.partitioning, cut.id)
        aggregation.aggregates[aggregation.partitioning[cut.id]] += cut
        return false
    end
    (idx, full) = select(aggregation.rule, aggregation.aggregates, cut)
    aggregation.aggregates[idx] += cut
    aggregation.partitioning[cut.id] = idx
    if full
        aggregate = aggregation.aggregates[idx]
        aggregation.aggregates[idx] = zero(AggregatedOptimalityCut{T})
        return add_cut!(lshaped, aggregate)
    end
    return false
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::DynamicAggregation, ::MetaData, t::Integer, cut::HyperPlane, x::AbstractArray)
    put!(cutqueue, (t, cut))
    return nothing
end

function aggregate_cut!(cutqueue::CutQueue, aggregation::DynamicAggregation{T}, metadata::MetaData, t::Integer, cut::HyperPlane{OptimalityCut}, x::AbstractArray) where T <: AbstractFloat
    gap = fetch(metadata, t, :gap)
    if aggregation.lock(gap, t) && haskey(aggregation.partitioning, cut.id)
        aggregation.aggregates[aggregation.partitioning[cut.id]] += cut
        return nothing
    end
    (idx, full) = select(aggregation.rule, aggregation.aggregates, cut)
    aggregation.partitioning[cut.id] = idx
    aggregation.aggregates[idx] += cut
    if full
        aggregate = aggregation.aggregates[idx]
        put!(cutqueue, (t, aggregate))
        aggregation.aggregates[idx] = zero(AggregatedOptimalityCut{T})
    end
    return nothing
end

function nthetas(nscenarios::Integer, ::DynamicAggregation)
    return nscenarios
end

function nthetas(nscenarios::Integer, ::DynamicAggregation, ::AbstractScenarioProblems)
    return nscenarios
end

function flush!(lshaped::AbstractLShapedSolver, aggregation::DynamicAggregation{T}) where T <: AbstractFloat
    added = false
    for (i,aggregate) in enumerate(aggregation.aggregates)
        if !iszero(aggregate)
            added |= add_cut!(lshaped, aggregate)
            aggregation.aggregates[i] = zero(AggregatedOptimalityCut{T})
        end
    end
    reset!(aggregation.rule)
    return added
end

function flush!(cutqueue::CutQueue, aggregation::DynamicAggregation{T}, ::MetaData, t::Integer, x::AbstractArray) where T <: AbstractFloat
    for (i,aggregate) in enumerate(aggregation.aggregates)
        if !iszero(aggregate)
            put!(cutqueue, (t, aggregate))
            aggregation.aggregates[i] = zero(AggregatedOptimalityCut{T})
        end
    end
    reset!(aggregation.rule)
    return nothing
end

# API
# ------------------------------------------------------------
"""
    DynamicAggregate(naggregates::Integer, rule::SelectionRule; lock_after::Function = (τ,n)->false)

Factory object for [`DynamicAggregation`](@ref). Pass to `aggregate ` in the `LShapedSolver` factory function. See ?DynamicAggregation for parameter descriptions.

"""
struct DynamicAggregate{S <: SelectionRule} <: AbstractAggregator
    naggregates::Int
    rule::S
    lock::Function

    function DynamicAggregate(naggregates::Integer, rule::SelectionRule; lock_after = (τ,n)->false)
        S = typeof(rule)
        return new{S}(naggregates, rule, lock_after)
    end
end

function (aggregator::DynamicAggregate)(nscenarios::Integer, T::Type{<:AbstractFloat})
    if aggregator.rule isa SelectRandom && aggregator.rule.max == 1
        return NoAggregation()
    end
    return DynamicAggregation(aggregator.naggregates, aggregator.rule, aggregator.lock, T)
end

function remote_aggregator(aggregation::DynamicAggregation, ::AbstractScenarioProblems, ::Integer)
    return DynamicAggregate(length(aggregation.aggregates), aggregation.rule; lock_after = aggregation.lock)
end

function str(aggregator::DynamicAggregate)
    return "dynamic aggregation ruled by $(str(aggregator.rule))"
end


text:
export Cuboid

# Below, the cuboid axes describe the directions of edges.  For example, if a₁, a₂, a₃ are
# the axes of a 3D cuboid, then the volume inside the cuboid is the collection of points
# x₁a₁ + x₂a₂ + x₃a₃, where -rₙ ≤ xₙ ≤ rₙ.  This means aₙ and aₘ with n ≠ m span a face of
# the cuboid.
#
# The cuboid projection matrix (Cuboid.p) is the matrix that produces the coordinates
# (x₁, x₂, x₃) of a given point.  In other words, for a point x inside the cuboid, the nth
# entry of Cuboid.p * x must have magnitude ≤ rₙ.  This means that Cuboid.p * aₙ = eₙ,
# because aₙ has only the aₙ-component and the component must be 1.  This also means that
# each row of Cuboid.p is orthogononal to two cuboid axes, and therefore normal to the face
# spanned by the two cuboid axes.  (Note that the rows of Cuboid.p are not unit normals to
# the faces, becuase they are not unit vectors.)
mutable struct Cuboid{N,N²} <: Shape{N,N²}
    c::SVector{N,Float64}  # center of cuboid
    r::SVector{N,Float64}  # "radii" (semi-axes) in axis directions
    p::SMatrix{N,N,Float64,N²}  # projection matrix to cuboid coordinates
    Cuboid{N,N²}(c,r,p) where {N,N²} = new(c,r,p)  # suppress default outer constructor
end

Cuboid(c::SVector{N,<:Real},
       s::SVector{N,<:Real},
       axes::SMatrix{N,N,<:Real}=SMatrix{N,N,Float64}(I)
       ) where {N} =
    Cuboid{N,N*N}(c, 0.5s, inv(axes ./ sqrt.(sum(abs2,axes,dims=Val(1)))))

Cuboid(c::AbstractVector{<:Real},  # center of cuboid
       s::AbstractVector{<:Real},  # size of cuboid in axis directions
       axes::AbstractMatrix{<:Real}=Matrix{Float64}(I,length(c),length(c))) =  # columns are axes vectors (each being parallel to two sets of faces in 3D)
    (N = length(c); Cuboid(SVector{N}(c), SVector{N}(s), SMatrix{N,N}(axes)))

Cuboid(d::NTuple{2,AbstractVector{<:Real}}) =  # end points of diagonal
    Cuboid((d[1]+d[2])/2, abs.(d[2]-d[1]))

Base.:(==)(b1::Cuboid, b2::Cuboid) = b1.c==b2.c && b1.r==b2.r && b1.p==b2.p
Base.isapprox(b1::Cuboid, b2::Cuboid) = b1.c≈b2.c && b1.r≈b2.r && b1.p≈b2.p
Base.hash(b::Cuboid, h::UInt) = hash(b.c, hash(b.r, hash(b.p, hash(:Cuboid, h))))

function level(x::SVector{N,<:Real}, b::Cuboid{N}) where {N}
    d = b.p * (x - b.c)

    return maximum(abs.(d) ./ b.r) - 1.0
end

function surfpt_nearby(x::SVector{N,<:Real}, b::Cuboid{N}) where {N}
    ax = inv(b.p)  # axes: columns are unit vectors

    # Below, the rows of n are the unit normals to the faces of the cuboid.  Appropriate
    # signs will be multiplied later.  (Note that the signs will be such that the rows of n
    # are always outward directions, even if x is inside the cuboid.)
    n = (b.p ./ sqrt.(sum(abs2,b.p,dims=Val(2))[:,1]))  # b.p normalized in row direction

    # Below, θ[i], the angle between ax[:,i] and n[i,:], is always acute (i.e, cosθ .≥ 0),
    # because the diagonal entries of b.p * ax are positive (= 1) and the diagonal entries
    # of n * ax are the scaled version of the diagonal entries of b.p * ax with positive
    # scale factors.
    cosθ = sum(ax.*n', dims=Val(1))[1,:]  # equivalent to diag(n*ax)
    # cosθ = diag(n*ax)  # faster than SVector(ntuple(i -> ax[:,i]⋅n[i,:], Val(N)))
    # @assert all(cosθ .≥ 0)

    d = b.p * (x - b.c)
    n = n .* copysign.(1.0,d)  # operation returns SMatrix (reason for leaving n untransposed)
    absd = abs.(d)
    onbnd = abs.(b.r.-absd) .≤ Base.rtoldefault(Float64) .* b.r  # basically b.r .≈ absd but faster
    isout = (b.r.<absd) .| onbnd
    ∆ = (b.r .- absd) .* cosθ  # entries can be negative
    if count(isout) == 0  # x strictly inside cuboid; ∆ all positive
        l∆x, i = findmin(∆)  # find closest face
        nout = n[i,:]
        ∆x = l∆x * nout
    else  # x outside cuboid or on boundary in one or multiple directions
        ∆x = n' * (∆ .* isout)  # project out .!isout directions

        # Below, (.!isout .| onbnd) tests if each entry of ∆ (not ∆x) is either projected
        # out or on the boundary.  If the ith entry of ∆ is projected out, ∆x (not ∆) does
        # not have any component in the n[i,:] direction.  Even if the ith entry of ∆ is not
        # projected out, if it is too small then ∆x barely has any component in the n[i,:]
        # direction.  If all the entries of ∆ satisfy one of the two conditions, ∆x ≈ 0 and
        # we cannot use -∆x for nout.  In that case, take n'*onbnd as nout.  When x is
        # outside the cuboid only in one dimension, n'*onbnd is basically the row of n along
        # that dimension, which reduces to normalize(-∆x) even if ∆x is very small.
        nout = all(.!isout .| onbnd) ? n'*onbnd : -∆x
        nout = normalize(nout)
    end

    return x+∆x, nout
end

signmatrix(b::Cuboid{1}) = SMatrix{1,1}(1)
signmatrix(b::Cuboid{2}) = SMatrix{2,2}(1,1, -1,1)
signmatrix(b::Cuboid{3}) = SMatrix{3,4}(1,1,1, -1,1,1, 1,-1,1, 1,1,-1)

function bounds(b::Cuboid)
    A = inv(b.p) .* b.r'
    m = maximum(abs.(A * signmatrix(b)), dims=Val(2))[:,1] # extrema of all 2^N corners of the cuboid
    return (b.c-m,b.c+m)
end


text:
@testset "Lattice" begin

    param = Dict{String,Any}("L" => 3)
    @testset "$lname" for (lname, latt) in [("dimer", dimer_lattice),
                                            ("chain", chain_lattice),
                                            ("square", square_lattice),
                                            ("triangular", triangular_lattice),
                                            ("cubic", cubic_lattice),
                                           ]
        lat = latt(param)
        for s in 1:numsites(lat)
            cs = sitecoordinate(lat,s)
            for (ns, nb) in neighbors(lat, s)
                if source(lat,nb) == s
                    @test target(lat,nb) == ns
                    src = s
                    tgt = ns
                else
                    @test source(lat,nb) == ns
                    @test target(lat,nb) == s
                    src = ns
                    tgt = s
                end
                cs = lattice_sitecoordinate(lat,src)
                ct = lattice_sitecoordinate(lat,tgt)
                cb  = lattice_bonddirection(lat,nb)
                ct2 = (cs .+ cb .+ size(lat)) .% size(lat)
                @test ct == ct2
            end
        end
    end
end


text:
struct TikzPlotFormula <: AbstractTikzPlot
	x::AbstractVector{Number}
	y::AbstractVector{Number}
	attributes::Dict{String, String}
end

# A fomula specifying only a function of the form y = f(x)
function formula(y::String, attributes::Dict{String,String}=Dict())

end

# A fomula with a range of x values and y specified as a fomula
function formula(x::AbstractVector{<:Real}, y::String, attributes::Dict{String,String}=Dict())
	return TikzPlotScatter(x, y, attributes)
end

function ToTikzString(plot::TikzPlotScatter)
	retString = ""
	retString *= "\\addplot[\n"
	for (k,v) in plot.attributes
		retString *= "$k = $v,\n"
	end
	retString *= "]\n"
	retString *= "coordinates\n{\n"
	for pt_i in range(1, length(plot.x))
		retString *= "($(plot.x[pt_i]), $(plot.y[pt_i]))\n"
	end
	retString *= "};\n"
	return retString
end


text:
include("_common.jl")

if nworkers() < 2
    addprocs(min(length(APPROACHES), CPU_CORES), exeflags="--project=..")
end

include_everywhere("_common.jl")
@everywhere baseline_path = joinpath(DATA_PATH, "baselines") * "/"
mkpath(baseline_path)

@everywhere calib_fig_path = joinpath(FIG_PATH, "calib_baseline") * "/"
mkpath(calib_fig_path)

@everywhere function run_calibration(metric)
    approach, objfunc = metric

    wid = Distributed.myid()
    @info "Worker $wid : Calibrating $approach"

    sn, n_id = setup_network("410730_IHACRES.yml")
    node = sn[n_id]

    # Calibrate network using the BlackBoxOptim package
    # keyword arguments will be passed to the `bboptimize()` function
    # This will set node parameters to the optimal values found
    calib_start = 366
    func = nothing
    if !occursin("RMSE", approach)
        func = (obs, sim) -> abs(1.0 - objfunc(obs[calib_start:end], sim[calib_start:end]))
    else
        func = (obs, sim) -> objfunc(obs[calib_start:end], sim[calib_start:end])
    end

    calibrate!(sn, CALIB_CLIMATE, CALIB_FLOW; metric=func,
               TraceInterval=300, MaxTime=60*60*6, PopulationSize=1000)

    # Save calibrated network spec to file
    spec_fn = "$(baseline_path)cotter_baseline_IHACRES_$(approach).yml"
    Streamfall.save_network_spec(sn, spec_fn)

    reset!(node)

    # Run calibrated node
    Streamfall.run_node!(sn, n_id, FULL_CLIMATE)
    sim_flow = node.outflow

    calib_data = FULL_DATASET[calib_start:CALIB_LN, "410730_Q"]
    calib_sim = sim_flow[calib_start:CALIB_LN]
    rmse = round(Streamfall.RMSE(calib_data, calib_sim), digits=2)
    nse = round(Streamfall.NSE(calib_data, calib_sim), digits=2)
    kge = round(Streamfall.KGE(calib_data, calib_sim), digits=2)
    plot(FULL_DATASET.Date[calib_start:CALIB_LN], calib_data,
         title="$approach\n(RMSE: $rmse; NSE: $nse; KGE: $kge)",
         legend=:topright,
         label="Historic",
         xlabel="Date",
         ylabel="Streamflow [ML]")
    plot!(CALIB_DATES[calib_start:end], calib_sim, label="Calibration", alpha=0.7)
    savefig(joinpath(calib_fig_path, "$(approach)_calibration.png"))

    valid_data = FULL_DATASET[CALIB_LN+1:end, "410730_Q"]
    valid_sim = sim_flow[CALIB_LN+1:end]
    rmse = round(Streamfall.RMSE(valid_data, valid_sim), digits=2)
    nse = round(Streamfall.NSE(valid_data, valid_sim), digits=2)
    kge = round(Streamfall.KGE(valid_data, valid_sim), digits=2)
    plot(FULL_DATASET.Date[CALIB_LN+1:end], valid_data,
         title="$approach\n(RMSE: $rmse; NSE: $nse; KGE: $kge)",
         legend=:topleft,
         label="Historic",
         xlabel="Date",
         ylabel="Streamflow [ML]")
    plot!(VALID_DATES, valid_sim, label="Validation", alpha=0.7)
    savefig(joinpath(calib_fig_path, "$(approach)_validation.png"))
end

pmap(run_calibration, zip(APPROACHES, OBJFUNCS))


text:
julia> "My String"[2:end] # without first character
"y String"

julia> "My String"[1:end-1] # without last character
"My Strin"

julia> "My String"[2:end-1] # without first and last characters
"y Strin"

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Utility functions for Probabilistic Programming - 1
Wouter Kouw
Last update: 09-01-2020
"""

import SpecialFunctions: gamma
import Distributions: Normal, Gamma
import DataFrames: DataFrame, select!, insertcols!
import MCMCChains: Chains
import Plots: contour, contour!, xlabel!, ylabel!

function prob_normalgamma(μ::Number, τ::Number; m=0, l=1, a=1, b=1)
    "Compute probability of sample under Normal-Gamma distribution."
    return b^a*sqrt(l)/(gamma(a)*sqrt(2*pi)) * τ^(a - 1/2) * exp(-b*τ) * exp( -(l*τ*(μ - m)^2)/2)
end

function prob_normal_gamma(μ::Number, τ::Number; m=0, v=1, a=1, b=1)
    "Compute probability of sample under a bivariate [Normal, Gamma] distribution."

    # Probablity under Normal distribution (Distributions.jl uses mean-standard-deviation parameterization)
    pμ = 1/(sqrt(2*pi*v))*exp(-(μ - m)^2/(2*v))

    # Probablity under Gamma distribution (Distributions.jl uses shape-scale parameterization)
    pτ = b^a / gamma(a)*τ^(a-1)*exp(-b*τ)

    # Return product of marginal distributions
    return pμ * pτ
end

function plot_normalgamma(m, l, a, b;
                          label=nothing,
                          x1=range(-5, stop=5),
                          x2=range(0.1, stop=5),
                          overlay=false)
    "Contour plot of Normal-Gamma distribution with given parameters."

    if overlay
        # Overlay the posterior probability
        contour!(x1, x2, (x, y) -> prob_normalgamma(x, y, m=m, l=l, a=a, b=b), label=label, cbar=nothing)
    else
        # Overlay the posterior probability
        contour(x1, x2, (x, y) -> prob_normalgamma(x, y, m=m, l=l, a=a, b=b), label=label, cbar=nothing)
    end

    # Axis properties
    xlabel!("μ")
    ylabel!("τ")
end

function plot_normal_gamma(m_μ, v_μ, a_τ, b_τ;
                           label=nothing,
                           color=nothing,
                           x1=range(-5, stop=5),
                           x2=range(0.1, stop=5),
                           overlay=false)
    "Contour plot of bivariate [Normal, Gamma] distribution with given parameters."

    if overlay
        # Overlay the posterior probability
        contour!(x1, x2,
                 (x, y) -> prob_normal_gamma(x, y, m=m_μ, v=v_μ, a=a_τ, b=b_τ),
                 label=label,
                 color=color,
                 cbar=nothing)
    else
        # Overlay the posterior probability
        contour(x1, x2,
                (x, y) -> prob_normal_gamma(x, y, m=m_μ, v=v_μ, a=a_τ, b=b_τ),
                label=label,
                color=color,
                cbar=nothing)
    end

    # Axis properties
    xlabel!("μ")
    ylabel!("τ")
end

function plot_decisionboundary(weights;
                               label=nothing,
                               color="black",
                               x1=range(-4, step=0.01, stop=4),
                               x2=range(-4, step=0.01, stop=4),
                               overlay=true)
    "Contour plot of linear function"

    # Logistic function
    ϕ(x) = 1/(1 + exp(-x))

    if overlay
        contour!(x1, x2,
                 (x,y) -> ϕ(weights[1]*x + weights[2]*y + weights[3]) > 0.5,
                 label=label,
                 color=color,
                 cbar=nothing)
    else
        contour(x1, x2,
                (x,y) -> ϕ(weights[1]*x + weights[2]*y + weights[3]) > 0.5,
                label=label,
                color=color,
                cbar=nothing)
    end
end

function dfcol_str2int(df::DataFrame, colname::Symbol)
    "Map column of dataframe from type string to type integer."

    # Create dummy array
    dummy = zeros(size(df,1),)

    # Map string elements to integer
    for (i, element) in enumerate(unique(df[!, colname]))
        dummy[df[!,colname] .== element] .= i-1
    end

    # Remove column
    select!(df, Not(colname))

    # Insert dummy var as column
    insertcols!(df, size(df,2), (colname => dummy))

    return df
end

function mean_chain(chain::Chains, param::Symbol)
    "Quick function to extract sample average from MCMCChains object"

    # Mean of chain
    mchain = mean(chain[param].value.data, dims=[1, 3])

    # Return as matrix
    return vec(convert(Array{Float64,2}, mchain[:,:,1]))
end

max_stars_repo_name: lisakoolen/BMLIP

max_issues_repo_name: lisakoolen/BMLIP

max_forks_repo_name: lisakoolen/BMLIP


text:
using CompScienceMeshes, BEAST

width, height, h = 1.0, 0.5, 0.05
G1 = meshrectangle(width, height, h)
G2 = CompScienceMeshes.rotate(G1, 1.0π * x̂)
G3 = CompScienceMeshes.rotate(G1, 1.5π * x̂)

G12 = weld(G1,-G2)
G23 = weld(G2,-G3)
G31 = weld(G3,-G1)

G21 = weld(G2,-G1)

X12 = raviartthomas(G12, sort=:none)
X23 = raviartthomas(G23, sort=:none)
X31 = raviartthomas(G31, sort=:none)

X21 = raviartthomas(G21, sort=:none)

Y12 = buffachristiansen(G12, sort=:none)
Y23 = buffachristiansen(G23, sort=:none)
Y31 = buffachristiansen(G31, sort=:none)

Y21 = buffachristiansen(G21)

κ = 1.0
SL = Maxwell3D.singlelayer(wavenumber=κ)
N = NCross()

X = X12 × X23 × X31
Y = Y12 × Y23 × Y31

# X = X12 × X21
# Y = Y12 × Y21

E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
e = (n × E) × n

@hilbertspace j
@hilbertspace k
mtefie = @discretise(
    SL[k,j] == e[k],
    j ∈ X, k ∈ X)

# u = gmres(mtefie)
#
# offset = 1; u12 = u[offset:offset+numfunctions(X12)-1]
# offset += numfunctions(X12); u23 = u[offset:offset+numfunctions(X23)-1]
# offset += numfunctions(X23); u31 = u[offset:offset+numfunctions(X31)-1]
#
# fcr12, _ = facecurrents(u12, X12)
# fcr23, _ = facecurrents(u23, X23)
# fcr31, _ = facecurrents(u31, X31)

import Plotly
using LinearAlgebra

function blkdiagm(blks...)
    m = sum(size(b,1) for b in blks)
    n = sum(size(b,2) for b in blks)
    A = zeros(eltype(first(blks)), m, n)
    o1 = 1
    o2 = 1
    for b in blks
        m1 = size(b,1)
        n1 = size(b,2)
        A[o1:o1+m1-1,o2:o2+n1-1] .= b
        o1 += m1
        o2 += n1
    end
    A
end

Sxx = BEAST.sysmatrix(mtefie)
N1 = assemble(N, X12, Y12)
N2 = assemble(N, X23, Y23)
N3 = assemble(N, X31, Y31)
Nxy = blkdiagm(N1,N2,N3)
Syy = assemble(SL, Y, Y)
iNxy = blkdiagm(inv(Matrix(N1)), inv(Matrix(N2)), inv(Matrix(N3)))

# N21 = assemble(N, X21, Y21)
# Nxy = blkdiagm(N1, N21)
# iNxy = blkdiagm(inv(Matrix(N1)), inv(Matrix(N21)))

ex = assemble(e,X)
P = transpose(iNxy) * Syy * iNxy
Q = P * Sxx;
R = P * ex;

u1, ch1 = solve(BEAST.GMRESSolver(Sxx), ex)
u2, ch2 = solve(BEAST.GMRESSolver(Q), R)

@show ch1.iters
@show ch2.iters

G123 = weld(G1,G2,G3)
X123 = raviartthomas(G123)

stefie = @discretise(
    SL[k,j] == e[k],
    j ∈ X123, k ∈ X123)

Sst = assemble(SL, X123, X123)
bst = assemble(e, X123)
ust, chst = solve(BEAST.GMRESSolver(Sst),bst)

fcrst, geost = facecurrents(ust, X123)
Plotly.plot(patch(geost, norm.(fcrst)))

Φ, Θ = [0.0], range(0,stop=π,length=100)
pts = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for ϕ in Φ for θ in Θ]

ffd_mt    = potential(MWFarField3D(wavenumber=κ), pts, u1, X)
ffd_mt_pc = potential(MWFarField3D(wavenumber=κ), pts, u2, X)
ffd_err   = potential(MWFarField3D(wavenumber=κ), pts, u1-u2, X)
ffd_st    = potential(MWFarField3D(wavenumber=κ), pts, ust, X123)

using Plots
plot(norm.(ffd_mt))
scatter!(norm.(ffd_mt_pc))
scatter!(norm.(ffd_st))
scatter!(norm.(ffd_err))


text: # Most functions are moved to src/utils.jl


text:
using DifferentialEquations
using ProgressMeter
using PyCall
using MATLAB
using Printf
using BSON: @save, @load
# pygui(:qt)
using PyPlot
include("PSS_Definitions_Dimer_X.jl")
# pygui(true);

H=range(.08,stop=0.15,length=25)
H=.1
@showprogress for Energy in H
location_fig="/media/brandon_behring/Extra_Space/MATLAB_FIGURES/"
location_BSON="/media/brandon_behring/Extra_Space/BSON_FILES/"
h=replace(@sprintf("%.15f",Energy),"."=>"_")

t_end=1e5;
N_iter_Q=25
Q_start=0.0
Q_end=.9
N_iter_P=25
P_start=0.0
P_end=.45

# t_end=1e5;
# N_iter_Q=50
# Q_start=0.1
# Q_end=.3
# N_iter_P=50
# P_start=0.0
# P_end=.3

nQ=@sprintf("_%d",N_iter_Q)
nP=@sprintf("_%d",N_iter_P)
h_title=@sprintf("%.8f",Energy)
println(h_title)
file_name=location_fig*"DIMER_PSS_X"*h*nQ*nP*".fig"
h_BSON=location_BSON*"DIMER_PSS_X"*h*nQ*nP*".bson"

ArrP=range(P_start,stop=P_end,length=N_iter_P)
ArrQ=range(Q_start,stop=Q_end,length=N_iter_Q)
Brandons_Colors=["#393b79" ,"#5254a3","#6b6ecf","#9c9ede" ,"#637939","#8ca252" ,"#b5cf6b" ,"#cedb9c" ,"#8c6d31","#bd9e39" ,"#e7ba52","#e7cb94","#843c39","#ad494a" ,"#d6616b","#e7969c" ,"#7b4173" ,"#a55194","#ce6dbd" ,"#de9ed6"];
SAVE_DATA=Dict()
# figure()
@showprogress for k=1:N_iter_Q
    for j=1:N_iter_P
        Q_PSS,P_PSS=PSS_function(ArrQ[k], ArrP[j], Energy, t_end);
        SAVE_DATA[(ArrQ[k],ArrP[j])]=Q_PSS, P_PSS
        # current_color=Brandons_Colors[mod(k,length(Brandons_Colors))+1]
    end
end
# @showprogress for k=1:N_iter_P
#     Q_PSS,P_PSS=PSS_function(0,ArrP[k], Energy, t_end);
#     current_color=Brandons_Colors[mod(k,length(Brandons_Colors))+1]
#     SAVE_DATA[k+N_iter_Q]=Q_PSS, P_PSS
#     if Q_PSS!=0
#         # mat"plot($Q_PSS,$P_PSS,'.','MarkerSize',3,'color',$current_color); hold on;"
#         # mat"plot(-$Q_PSS,-$P_PSS,'.','MarkerSize',3,'color',$current_color); hold on;"
#     end
# end

@save h_BSON SAVE_DATA
Yfind(h)=sqrt(h/(2h+1))
P=Yfind(Energy)
mat"figure();set(gcf, 'Position',  [0, 0, 2500, 2000]); hold on;"
mat"xlabel('Q')"
mat"ylabel('P')"
mat"title($h_title)"
mat"axis([-.9, .9, -.45, .45])"
# mat"axis([-1, 1, -.5, .5])"
@showprogress for k=1:N_iter_Q
    for j=1:N_iter_P
        Q_PSS,P_PSS= SAVE_DATA[(ArrQ[k],ArrP[j])]
        current_color=Brandons_Colors[mod(k*j,length(Brandons_Colors))+1]
        mat"plot($Q_PSS,$P_PSS,'.','MarkerSize',2,'color',$current_color); hold on;"
        # mat"plot(-$Q_PSS,$P_PSS,'.','MarkerSize',2,'color',$current_color); hold on;"
        # mat"plot($Q_PSS,-$P_PSS,'.','MarkerSize',2,'color',$current_color); hold on;"
        mat"plot(-$Q_PSS,-$P_PSS,'.','MarkerSize',2,'color',$current_color); hold on;"
    end
end
Q_fix=sqrt(6)/3
mat"plot($Q_fix,0,'b.','MarkerSize',30)"
mat"plot(-$Q_fix,0,'b.','MarkerSize',30)"
mat"plot(0,$P,'r.','MarkerSize',30)"
mat"plot(0,-$P,'r.','MarkerSize',30)"
mat"savefig($file_name)"
end


text:
using DelimitedFiles

# returns solution of part 1
function x()
	t = vec(readdlm("10.txt", Int64))
	sort!(t)

	a=b=0
	prev = 0
	for ch in t
		if ch - prev == 1
			a += 1
		elseif ch - prev == 3
			b += 1
		end
		prev = ch
	end
	b += 1
	a*b
end

# returns solution of part 2
function y()
	t = sort(vec(readdlm("10.txt", Int64)))

	nums = Array{Int64,1}(zeros(maximum(t)+3))

	for ch in t
		s = ch <= 3 ? 1 : 0
		for i = max(1, ch-3):ch-1
			s += nums[i]
		end
		nums[ch] = s
	end
	nums[maximum(t)]
end