"""
    nipals(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing, AbstractString,AbstractArray}="", dim::Number=3, numepoch::Number=3, lower::Number=0, upper::Number=1.0f+38, offsetStoch::Number=1f-20, logdir::Union{Nothing,AbstractString}=nothing)
Online CCA solved by NIPALS.
Input Arguments
---------
- `input` : Julia Binary file generated by `OnlineCCA.csv2bin` function. Only two files can be specified (e.g. X.zst and Y.zst)
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `colmeanlist` : The mean of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `colvarlist` : The variance of each column of matrix. The CSV file is generated by `OnlineCCA.sumr` functions.
- `dim` : The number of dimension of CCA.
- `numepoch` : The number of epoch.
- `lower` : Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)
- `upper` : Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)
- `evalfreq` : Evaluation Frequency of Reconstruction Error
- `offsetStoch` : Off set value for avoding overflow when calculating stochastic gradient
---------
- `T` : = (t1,t2,...,tK) (No. columns of X × dim), the score vectors of X
- `P` : = (p1,p2,...,pK) (No. rows of X × dim), the loading of X
- `U` : = (u1,u2,...,uK) (No. columns of Y × dim), the score vectors of Y
- `Q` : = (q1,q2,...,qK) (No. rows of Y × dim), the loading of Y
- `β` : The Regression of Coefficient (dim × dim), Y is predicted by the new X by Y_new = X_new P β Q'
"""
function nipals(;input::AbstractArray=[], outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractArray=[], pseudocount::AbstractArray=[], colmeanlist::AbstractArray="", colvarlist::Union{Nothing,AbstractString,AbstractArray}="", dim::Number=3, numepoch::Number=3, lower::Number=0, upper::Number=1.0f+38, offsetStoch::Number=1f-20, logdir::Union{Nothing,AbstractString}=nothing)
    # Initial Setting
    cca = NIPALS()
    pseudocount, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, lower, upper, offsetStoch = init(input, pseudocount, dim, colmeanlist, colvarlist, cca, lower, upper, offsetStoch, scale)
    # Perform CCA
    out = nipals(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, numepoch, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, lower, upper, offsetStoch)
    if outdir isa String
        output(outdir, out, cca)
    end
    return out
end

function nipals(input, outdir, scale, pseudocount, colmeanlist, colvarlist, dim, numepoch, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, lower, upper, offsetStoch)
    @assert length(input) == 2
    l = 2
    N, M = nm(input)
    for i=1:l
        @assert 0 < dim ≤ min(N[i], M[i]) - 1
    end
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    X = input[1]
    Y = input[2]
    x = zeros(UInt32, M[1])
    normx = zeros(UInt32, M[1])
    defx = zeros(UInt32, M[1])
    y = zeros(UInt32, M[2])
    normy = zeros(UInt32, M[2])
    defy = zeros(UInt32, M[2])
    T = zeros(Float32, N[1], dim)
    P = zeros(Float32, M[1], dim)
    U = zeros(Float32, N[2], dim)
    Q = zeros(Float32, M[2], dim)
    tmp = tempdir()
    old_t = zeros(Float32, N[1])
    # NIPALS Iteration
    for k=1:dim
        println("#######################################")
        println(string(k)*"-th u, t, p, and q are calculated...")
        println("#######################################")
        outsvd = oocmcca(input, outdir, scale, pseudocount, colmeanlist, colvarlist, 1, cca, W, v, D, colmeanvec, colvarvec, N, M, TotalCorVar, offsetStoch, false)
        u = outsvd[3][1][:,1]
        t = outsvd[3][2][:,1]
        p = zeros(Float32, M[1])
        q = zeros(Float32, M[2])
        iter = 1
        relchange = Any[]
        push!(relchange, 10^20)
        while (relchange[iter] > lower) && (iter <= numepoch)
            println(string(iter)*"-th iteration : ||t_old - t_new|| = "*string(relchange[iter]))
            old_t .= t
            ####################
            ##### p = X'u ######
            ####################
            open(input[1]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[1]
                    # Row vector of data matrix
                    read!(stream, x)
                    if k == 1
                        normx = normalizex(x, n, stream, scale[1], pseudocount[1], colmeanlist[1], colmeanvec[1], colvarlist[1], colvarvec[1])
                    else
                        normx = x
                    end
                    # p = X'u
                    p .+= normx[:,1] .* u[n]
                end
                close(stream)
            end
            p = p ./ norm(p)

            ####################
            ###### t = Xp ######
            ####################
            open(input[1]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[1]
                    # Row vector of data matrix
                    read!(stream, x)
                    if k == 1
                        normx = normalizex(x, n, stream, scale[1], pseudocount[1], colmeanlist[1], colmeanvec[1], colvarlist[1], colvarvec[1])
                    else
                        normx = x
                    end
                    # t = Xp
                    t[n] = normx[:,1]'p
                end
                close(stream)
            end

            ####################
            ###### q = Y't #####
            ####################
            open(input[2]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[2]
                    # Row vector of data matrix
                    read!(stream, y)
                    if k == 1
                        normy = normalizex(y, n, stream, scale[2], pseudocount[2], colmeanlist[2], colmeanvec[2], colvarlist[2], colvarvec[2])
                    else
                        normy = y
                    end
                    # q = Y't
                    q .+= normy[:,1] .* t[n]
                end
                close(stream)
            end
            q = q ./ norm(q)

            ####################
            ###### u = Yq ######
            ####################
            open(input[2]) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                for n=1:N[2]
                    # Row vector of data matrix
                    read!(stream, y)
                    if k == 1
                        normy = normalizex(y, n, stream, scale[2], pseudocount[2], colmeanlist[2], colmeanvec[2], colvarlist[2], colvarvec[2])
                    else
                        normy = y
                    end
                    # u = Yq
                    u[n] = normy[:,1]'q
                end
                close(stream)
            end
            # Relative Change
            iter = iter + 1
            push!(relchange, abs(norm(old_t - t)))
        end
        # Update
        T[:,k] = t
        P[:,k] = p
        U[:,k] = u
        Q[:,k] = q
        ####################
        #### Deflation #####
        ####################
        if k < dim
            defXfile = tmp*"/defX"*string(k)*".zst"
            defYfile = tmp*"/defY"*string(k)*".zst"
            open(defXfile, "w") do deffile
                defstream = ZstdCompressorStream(deffile)
                open(input[1]) do file
                    stream = ZstdDecompressorStream(file)
                    read!(stream, tmpN)
                    read!(stream, tmpM)
                    write(defstream, tmpN)
                    write(defstream, tmpM)
                    # Each step n
                    for n=1:N[1]
                        # Row vector of data matrix
                        read!(stream, x)
                        normx = normalizex(x, n, stream, scale[1], pseudocount[1], colmeanlist[1], colmeanvec[1], colvarlist[1], colvarvec[1])
                        # X = X - tp'
                        write(defstream, normx[:,1]' - t[n]*p')
                    end
                    close(stream)
                end
                close(defstream)
            end

            open(defYfile, "w") do deffile
                defstream = ZstdCompressorStream(deffile)
                open(input[2]) do file
                    stream = ZstdDecompressorStream(file)
                    read!(stream, tmpN)
                    read!(stream, tmpM)
                    write(defstream, tmpN)
                    write(defstream, tmpM)
                    # Each step n
                    for n=1:N[2]
                        # Row vector of data matrix
                        read!(stream, y)
                        normy = normalizex(y, n, stream, scale[2], pseudocount[2], colmeanlist[2], colmeanvec[2], colvarlist[2], colvarvec[2])
                        # Y = Y - uq'
                        write(defstream, normy[:,1]' - u[n]*q')
                    end
                    close(stream)
                end
                close(defstream)
            end
            input[1] = defXfile
            input[2] = defYfile
        end
    end
    # Regression Coefficient
    W,Σ,V = svd(T'T)
    β = V*Diagonal(1/Σ)*W'T'U
    return (T, P, U, Q, β)
end
