# Path: fdasrsf/utility_functions.py
# def diffop(n, binsize=1):
#     """
#     Creates a second order differential operator
# 
#     :param n: dimension
#     :param binsize: dx (default = 1)
# 
#     :rtype: numpy ndarray
#     :return m: matrix describing differential operator
# 
#     """
#     m = diagflat(ones(n - 1), k=1) + diagflat(ones(n - 1), k=-1) + diagflat(2 * ones(n))
#     m = inner(m.transpose(), m)
#     m[0, 0] = 6
#     m[-1, -1] = 6
#     m /= (binsize ** 4.)
# 
#     return m
# 
# def geigen(Amat, Bmat, Cmat):
#     """
#     generalized eigenvalue problem of the form
# 
#     max tr L'AM / sqrt(tr L'BL tr M'CM) w.r.t. L and M
# 
#     :param Amat numpy ndarray of shape (M,N)
#     :param Bmat numpy ndarray of shape (M,N)
#     :param Bmat numpy ndarray of shape (M,N)
# 
#     :rtype: numpy ndarray
#     :return values: eigenvalues
#     :return Lmat: left eigenvectors
#     :return Mmat: right eigenvectors
# 
#     """
#     if Bmat.shape[0] != Bmat.shape[1]:
#         print("BMAT is not square.\n")
#         sys.exit(1)
# 
#     if Cmat.shape[0] != Cmat.shape[1]:
#         print("CMAT is not square.\n")
#         sys.exit(1)
# 
#     p = Bmat.shape[0]
#     q = Cmat.shape[0]
# 
#     s = min(p, q)
#     tmp = fabs(Bmat - Bmat.transpose())
#     tmp1 = fabs(Bmat)
#     if tmp.max() / tmp1.max() > 1e-10:
#         print("BMAT not symmetric..\n")
#         sys.exit(1)
# 
#     tmp = fabs(Cmat - Cmat.transpose())
#     tmp1 = fabs(Cmat)
#     if tmp.max() / tmp1.max() > 1e-10:
#         print("CMAT not symmetric..\n")
#         sys.exit(1)
# 
#     Bmat = (Bmat + Bmat.transpose()) / 2.
#     Cmat = (Cmat + Cmat.transpose()) / 2.
#     Bfac = cholesky(Bmat)
#     Cfac = cholesky(Cmat)
#     Bfacinv = inv(Bfac)
#     Bfacinvt = Bfacinv.transpose()
#     Cfacinv = inv(Cfac)
#     Dmat = Bfacinvt.dot(Amat).dot(Cfacinv)
#     if p >= q:
#         u, d, v = svd(Dmat)
#         values = d
#         Lmat = Bfacinv.dot(u)
#         Mmat = Cfacinv.dot(v.transpose())
#     else:
#         u, d, v = svd(Dmat.transpose())
#         values = d
#         Lmat = Bfacinv.dot(u)
#         Mmat = Cfacinv.dot(v.transpose())
# 
#     return values, Lmat, Mmat
# 
# def innerprod_q(time, q1, q2):
#     """
#     calculates the innerproduct between two srsfs
# 
#     :param time vector descrbing time samples
#     :param q1 vector of srsf 1
#     :param q2 vector of srsf 2
# 
#     :rtype: scalar
#     :return val: inner product value
# 
#     """
#     val = trapz(q1 * q2, time)
#     return val


import numpy as np
from fdasrsf.utility_functions import diffop, geigen, innerprod_q
def pls_svd(time, qf, qg, no, alpha=0.0):
    """
    This function computes the partial least squares using SVD

    :param time: vector describing time samples
    :param qf: numpy ndarray of shape (M,N) of N functions with M samples
    :param qg: numpy ndarray of shape (M,N) of N functions with M samples
    :param no: number of components
    :param alpha: amount of smoothing (Default = 0.0 i.e., none)

    :rtype: numpy ndarray
    :return wqf: f weight function
    :return wqg: g weight function
    :return alpha: smoothing value
    :return values: singular values

    """
    binsize = np.diff(time)
    binsize = binsize.mean()

    Kfg = np.cov(qf, qg)
    Kfg = Kfg[0:qf.shape[0], qf.shape[0]:Kfg.shape[0]]
    nx = Kfg.shape[0]
    D4x = diffop(nx, binsize)
    values, Lmat, Mmat = geigen(Kfg, np.eye(nx) + alpha * D4x, np.eye(nx) + alpha * D4x)
    wf = Lmat[:, 0:no]
    wg = Mmat[:, 0:no]

    for ii in range(0, no):
        wf[:, ii] = wf[:, ii] / np.sqrt(innerprod_q(time, wf[:, ii], wf[:, ii]))

# Path: fdasrsf/utility_functions.py
# def diffop(n, binsize=1):
#     """
#     Creates a second order differential operator
# 
#     :param n: dimension
#     :param binsize: dx (default = 1)
# 
#     :rtype: numpy ndarray
#     :return m: matrix describing differential operator
# 
#     """
#     m = diagflat(ones(n - 1), k=1) + diagflat(ones(n - 1), k=-1) + diagflat(2 * ones(n))
#     m = inner(m.transpose(), m)
#     m[0, 0] = 6
#     m[-1, -1] = 6
#     m /= (binsize ** 4.)
# 
#     return m
# 
# def geigen(Amat, Bmat, Cmat):
#     """
#     generalized eigenvalue problem of the form
# 
#     max tr L'AM / sqrt(tr L'BL tr M'CM) w.r.t. L and M
# 
#     :param Amat numpy ndarray of shape (M,N)
#     :param Bmat numpy ndarray of shape (M,N)
#     :param Bmat numpy ndarray of shape (M,N)
# 
#     :rtype: numpy ndarray
#     :return values: eigenvalues
#     :return Lmat: left eigenvectors
#     :return Mmat: right eigenvectors
# 
#     """
#     if Bmat.shape[0] != Bmat.shape[1]:
#         print("BMAT is not square.\n")
#         sys.exit(1)
# 
#     if Cmat.shape[0] != Cmat.shape[1]:
#         print("CMAT is not square.\n")
#         sys.exit(1)
# 
#     p = Bmat.shape[0]
#     q = Cmat.shape[0]
# 
#     s = min(p, q)
#     tmp = fabs(Bmat - Bmat.transpose())
#     tmp1 = fabs(Bmat)
#     if tmp.max() / tmp1.max() > 1e-10:
#         print("BMAT not symmetric..\n")
#         sys.exit(1)
# 
#     tmp = fabs(Cmat - Cmat.transpose())
#     tmp1 = fabs(Cmat)
#     if tmp.max() / tmp1.max() > 1e-10:
#         print("CMAT not symmetric..\n")
#         sys.exit(1)
# 
#     Bmat = (Bmat + Bmat.transpose()) / 2.
#     Cmat = (Cmat + Cmat.transpose()) / 2.
#     Bfac = cholesky(Bmat)
#     Cfac = cholesky(Cmat)
#     Bfacinv = inv(Bfac)
#     Bfacinvt = Bfacinv.transpose()
#     Cfacinv = inv(Cfac)
#     Dmat = Bfacinvt.dot(Amat).dot(Cfacinv)
#     if p >= q:
#         u, d, v = svd(Dmat)
#         values = d
#         Lmat = Bfacinv.dot(u)
#         Mmat = Cfacinv.dot(v.transpose())
#     else:
#         u, d, v = svd(Dmat.transpose())
#         values = d
#         Lmat = Bfacinv.dot(u)
#         Mmat = Cfacinv.dot(v.transpose())
# 
#     return values, Lmat, Mmat
# 
# def innerprod_q(time, q1, q2):
#     """
#     calculates the innerproduct between two srsfs
# 
#     :param time vector descrbing time samples
#     :param q1 vector of srsf 1
#     :param q2 vector of srsf 2
# 
#     :rtype: scalar
#     :return val: inner product value
# 
#     """
#     val = trapz(q1 * q2, time)
#     return val


# Path: fdasrsf/fPLS.py
import numpy as np
from fdasrsf.utility_functions import diffop, geigen, innerprod_q

def pls_svd(time, qf, qg, no, alpha=0.0):
    """
    This function computes the partial least squares using SVD

    :param time: vector describing time samples
    :param qf: numpy ndarray of shape (M,N) of N functions with M samples
    :param qg: numpy ndarray of shape (M,N) of N functions with M samples
    :param no: number of components
    :param alpha: amount of smoothing (Default = 0.0 i.e., none)

    :rtype: numpy ndarray
    :return wqf: f weight function
    :return wqg: g weight function
    :return alpha: smoothing value
    :return values: singular values

    """
    binsize = np.diff(time)
    binsize = binsize.mean()

    Kfg = np.cov(qf, qg)
    Kfg = Kfg[0:qf.shape[0], qf.shape[0]:Kfg.shape[0]]
    nx = Kfg.shape[0]
    D4x = diffop(nx, binsize)
    values, Lmat, Mmat = geigen(Kfg, np.eye(nx) + alpha * D4x, np.eye(nx) + alpha * D4x)
    wf = Lmat[:, 0:no]
    wg = Mmat[:, 0:no]

    for ii in range(0, no):
        wf[:, ii] = wf[:, ii] / np.sqrt(innerprod_q(time, wf[:, ii], wf[:, ii]))

        wg[:, ii] = wg[:, ii] / np.sqrt(innerprod_q(time, wg[:, ii], wg[:, ii]))