# Path: fdasrsf/fPLS.py
# 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]))
# 
#     wqf = wf
#     wqg = wg
# 
#     N = qf.shape[1]
#     rfi = np.zeros(N)
#     rgi = np.zeros(N)
# 
#     for l in range(0, N):
#         rfi[l] = innerprod_q(time, qf[:, l], wqf[:, 0])
#         rgi[l] = innerprod_q(time, qg[:, l], wqg[:, 0])
# 
#     cost = np.cov(rfi, rgi)[1, 0]
# 
#     return wqf, wqg, alpha, values, cost


import numpy as np
import matplotlib.pyplot as plt
import fdasrsf.utility_functions as uf
import fdasrsf.plot_style as plot
import fpls_warp as fpls
import collections
from scipy.integrate import trapz, cumtrapz
from scipy.linalg import svd
from numpy.linalg import norm
from joblib import Parallel, delayed
from fdasrsf.fPLS import pls_svd
    """
    print ("Initializing...")
    binsize = np.diff(time)
    binsize = binsize.mean()
    eps = np.finfo(np.double).eps
    M = f.shape[0]
    N = f.shape[1]
    f0 = f
    g0 = g

    if showplot:
        plot.f_plot(time, f, title="f Original Data")
        plot.f_plot(time, g, title="g Original Data")

    # Compute q-function of f and g
    f, g1, g2 = uf.gradient_spline(time, f, smoothdata)
    qf = g1 / np.sqrt(abs(g1) + eps)
    g, g1, g2 = uf.gradient_spline(time, g, smoothdata)
    qg = g1 / np.sqrt(abs(g1) + eps)

    print("Calculating fPLS weight functions for %d Warped Functions..." % N)
    itr = 0
    fi = np.zeros((M, N, max_itr + 1))
    fi[:, :, itr] = f
    gi = np.zeros((M, N, max_itr + 1))
    gi[:, :, itr] = g
    qfi = np.zeros((M, N, max_itr + 1))
    qfi[:, :, itr] = qf
    qgi = np.zeros((M, N, max_itr + 1))
    qgi[:, :, itr] = qg

# Path: fdasrsf/fPLS.py
# 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]))
# 
#     wqf = wf
#     wqg = wg
# 
#     N = qf.shape[1]
#     rfi = np.zeros(N)
#     rgi = np.zeros(N)
# 
#     for l in range(0, N):
#         rfi[l] = innerprod_q(time, qf[:, l], wqf[:, 0])
#         rgi[l] = innerprod_q(time, qg[:, l], wqg[:, 0])
# 
#     cost = np.cov(rfi, rgi)[1, 0]
# 
#     return wqf, wqg, alpha, values, cost


# Path: fdasrsf/time_warping.py
import numpy as np
import matplotlib.pyplot as plt
import fdasrsf.utility_functions as uf
import fdasrsf.plot_style as plot
import fpls_warp as fpls
import collections
from scipy.integrate import trapz, cumtrapz
from scipy.linalg import svd
from numpy.linalg import norm
from joblib import Parallel, delayed
from fdasrsf.fPLS import pls_svd

    """
    print ("Initializing...")
    binsize = np.diff(time)
    binsize = binsize.mean()
    eps = np.finfo(np.double).eps
    M = f.shape[0]
    N = f.shape[1]
    f0 = f
    g0 = g

    if showplot:
        plot.f_plot(time, f, title="f Original Data")
        plot.f_plot(time, g, title="g Original Data")

    # Compute q-function of f and g
    f, g1, g2 = uf.gradient_spline(time, f, smoothdata)
    qf = g1 / np.sqrt(abs(g1) + eps)
    g, g1, g2 = uf.gradient_spline(time, g, smoothdata)
    qg = g1 / np.sqrt(abs(g1) + eps)

    print("Calculating fPLS weight functions for %d Warped Functions..." % N)
    itr = 0
    fi = np.zeros((M, N, max_itr + 1))
    fi[:, :, itr] = f
    gi = np.zeros((M, N, max_itr + 1))
    gi[:, :, itr] = g
    qfi = np.zeros((M, N, max_itr + 1))
    qfi[:, :, itr] = qf
    qgi = np.zeros((M, N, max_itr + 1))
    qgi[:, :, itr] = qg

    wqf1, wqg1, alpha, values, costmp = pls_svd(time, qfi[:, :, itr],