From 767f8ae8bcdc6bd0bdff6ad677d60a36cc265cd9 Mon Sep 17 00:00:00 2001 From: "J. Derek Tucker" Date: Thu, 30 Nov 2023 20:22:56 -0700 Subject: [PATCH] add prediction interval --- fdasrsf/pcr_regression.py | 643 +++++++++++++++++++++----------------- 1 file changed, 353 insertions(+), 290 deletions(-) diff --git a/fdasrsf/pcr_regression.py b/fdasrsf/pcr_regression.py index 492de514..acbf92d0 100644 --- a/fdasrsf/pcr_regression.py +++ b/fdasrsf/pcr_regression.py @@ -11,17 +11,18 @@ import fdasrsf.fPCA as fpca import fdasrsf.regression as rg import fdasrsf.geometry as geo -from scipy.linalg import inv, norm -from scipy.integrate import trapz, cumtrapz +from scipy.linalg import inv +from scipy.stats import t from scipy.optimize import fmin_l_bfgs_b + class elastic_pcr_regression: """ This class provides elastic pcr regression for functional data using the SRVF framework accounting for warping - + Usage: obj = elastic_pcr_regression(f,y,time) - + :param f: (M,N) % matrix defining N functions of M samples :param y: response vector of length N :param warp_data: fdawarp object of alignment @@ -44,15 +45,21 @@ def __init__(self, f, y, time): a = time.shape[0] if f.shape[0] != a: - raise Exception('Columns of f and time must be equal') + raise Exception("Columns of f and time must be equal") self.f = f self.y = y self.time = time - def calc_model(self, pca_method="combined", no=5, - smooth_data=False, sparam=25, parallel=False, - C=None): + def calc_model( + self, + pca_method="combined", + no=5, + smooth_data=False, + sparam=25, + parallel=False, + C=None, + ): """ This function identifies a regression model with phase-variability using elastic pca @@ -67,59 +74,65 @@ def calc_model(self, pca_method="combined", no=5, """ if smooth_data: - self.f = fs.smooth_data(self.f,sparam) - + self.f = fs.smooth_data(self.f, sparam) + N1 = self.f.shape[1] # Align Data - self.warp_data = fs.fdawarp(self.f,self.time) + self.warp_data = fs.fdawarp(self.f, self.time) self.warp_data.srsf_align(parallel=parallel) # Calculate PCA - if pca_method=='combined': + if pca_method == "combined": self.pca = fpca.fdajpca(self.warp_data) - elif pca_method=='vert': + elif pca_method == "vert": self.pca = fpca.fdavpca(self.warp_data) - elif pca_method=='horiz': + elif pca_method == "horiz": self.pca = fpca.fdahpca(self.warp_data) else: - raise Exception('Invalid fPCA Method') + raise Exception("Invalid fPCA Method") self.pca.calc_fpca(no) - + # OLS using PCA basis lam = 0 R = 0 - Phi = np.ones((N1, no+1)) - Phi[:,1:(no+1)] = self.pca.coef + Phi = np.ones((N1, no + 1)) + Phi[:, 1: (no + 1)] = self.pca.coef xx = np.dot(Phi.T, Phi) inv_xx = inv(xx + lam * R) xy = np.dot(Phi.T, self.y) b = np.dot(inv_xx, xy) alpha = b[0] - b = b[1:no+1] + b = b[1: no + 1] # compute the SSE int_X = np.zeros(N1) - for ii in range(0,N1): - int_X[ii] = np.sum(self.pca.coef[ii,:]*b) - - SSE = np.sum((self.y-alpha-int_X)**2) + for ii in range(0, N1): + int_X[ii] = np.sum(self.pca.coef[ii, :] * b) + + y_pred = alpha + int_X + SSE = np.sum((self.y - y_pred) ** 2) + + self.Ntr = len(self.y) + stdev = np.sqrt(sum((y_pred - self.y)**2) / (self.Ntr - 2)) self.alpha = alpha self.b = b self.SSE = SSE + self.stdev = stdev self.pca_method = pca_method return - - def predict(self, newdata=None): + def predict(self, newdata=None, alpha=.05): """ - This function performs prediction on regression model on new data if available or current stored data in object + This function performs prediction on regression model on new + data if available or current stored data in object Usage: obj.predict() obj.predict(newdata) - :param newdata: dict containing new data for prediction (needs the keys below, if None predicts on training data) + :param newdata: dict containing new data for prediction (needs + the keys below, if None predicts on training data) :type newdata: dict :param f: (M,N) matrix of functions :param time: vector of time points @@ -132,102 +145,119 @@ def predict(self, newdata=None): lam = self.warp_data.lam M = self.time.shape[0] - if newdata != None: - f = newdata['f'] - time = newdata['time'] - y = newdata['y'] - if newdata['smooth']: - sparam = newdata['sparam'] - f = fs.smooth_data(f,sparam) - - q1 = fs.f_to_srsf(f,time) + if newdata is not None: + f = newdata["f"] + time = newdata["time"] + y = newdata["y"] + if newdata["smooth"]: + sparam = newdata["sparam"] + f = fs.smooth_data(f, sparam) + + q1 = fs.f_to_srsf(f, time) n = q1.shape[1] self.y_pred = np.zeros(n) mq = self.warp_data.mqn - fn = np.zeros((M,n)) - qn = np.zeros((M,n)) - gam = np.zeros((M,n)) - for ii in range(0,n): - gam[:,ii] = uf.optimum_reparam(mq,time,q1[:,ii],omethod,lam) - fn[:,ii] = uf.warp_f_gamma(time,f[:,ii],gam[:,ii]) - qn[:,ii] = uf.f_to_srsf(fn[:,ii],time) - + fn = np.zeros((M, n)) + qn = np.zeros((M, n)) + gam = np.zeros((M, n)) + for ii in range(0, n): + gam[:, ii] = uf.optimum_reparam(mq, time, q1[:, ii], omethod, lam) + fn[:, ii] = uf.warp_f_gamma(time, f[:, ii], gam[:, ii]) + qn[:, ii] = uf.f_to_srsf(fn[:, ii], time) + U = self.pca.U no = U.shape[1] - if self.pca.__class__.__name__ == 'fdajpca': - m_new = np.sign(fn[self.pca.id,:])*np.sqrt(np.abs(fn[self.pca.id,:])) + if self.pca.__class__.__name__ == "fdajpca": + m_new = np.sign(fn[self.pca.id, :]) * np.sqrt( + np.abs(fn[self.pca.id, :]) + ) qn1 = np.vstack((qn, m_new)) C = self.pca.C TT = self.time.shape[0] mu_g = self.pca.mu_g mu_psi = self.pca.mu_psi - vec = np.zeros((M,n)) - psi = np.zeros((TT,n)) + vec = np.zeros((M, n)) + psi = np.zeros((TT, n)) binsize = np.mean(np.diff(self.time)) - for i in range(0,n): - psi[:,i] = np.sqrt(np.gradient(gam[:,i],binsize)) - out, theta = geo.inv_exp_map(mu_psi, psi[:,i]) - vec[:,i] = out - - g = np.vstack((qn1, C*vec)) - a = np.zeros((n,no)) - for i in range(0,n): - for j in range(0,no): - tmp = (g[:,i]-mu_g) - a[i,j] = np.dot(tmp.T, U[:,j]) - - elif self.pca.__class__.__name__ == 'fdavpca': - m_new = np.sign(fn[self.pca.id,:])*np.sqrt(np.abs(fn[self.pca.id,:])) + for i in range(0, n): + psi[:, i] = np.sqrt(np.gradient(gam[:, i], binsize)) + out, theta = geo.inv_exp_map(mu_psi, psi[:, i]) + vec[:, i] = out + + g = np.vstack((qn1, C * vec)) + a = np.zeros((n, no)) + for i in range(0, n): + for j in range(0, no): + tmp = g[:, i] - mu_g + a[i, j] = np.dot(tmp.T, U[:, j]) + + elif self.pca.__class__.__name__ == "fdavpca": + m_new = np.sign(fn[self.pca.id, :]) * np.sqrt( + np.abs(fn[self.pca.id, :]) + ) qn1 = np.vstack((qn, m_new)) - a = np.zeros((n,no)) - for i in range(0,n): - for j in range(0,no): - tmp = (qn1[:,i]-self.pca.mqn) - a[i,j] = np.dot(tmp.T, U[:,j]) - - elif self.pca.__class__.__name__ == 'fdahpca': - a = np.zeros((n,no)) + a = np.zeros((n, no)) + for i in range(0, n): + for j in range(0, no): + tmp = qn1[:, i] - self.pca.mqn + a[i, j] = np.dot(tmp.T, U[:, j]) + + elif self.pca.__class__.__name__ == "fdahpca": + a = np.zeros((n, no)) mu_psi = self.pca.psi_mu - vec = np.zeros((M,n)) + vec = np.zeros((M, n)) TT = self.time.shape[0] - psi = np.zeros((TT,n)) + psi = np.zeros((TT, n)) binsize = np.mean(np.diff(self.time)) - for i in range(0,n): - psi[:,i] = np.sqrt(np.gradient(gam[:,i],binsize)) - out, theta = geo.inv_exp_map(mu_psi, psi[:,i]) - vec[:,i] = out - + for i in range(0, n): + psi[:, i] = np.sqrt(np.gradient(gam[:, i], binsize)) + out, theta = geo.inv_exp_map(mu_psi, psi[:, i]) + vec[:, i] = out + vm = self.pca.vec.mean(axis=1) - for i in range(0,n): - for j in range(0,no): - a[i,j] = np.sum(np.dot(vec[:,i]-vm,U[:,j])) - else: - raise Exception('Invalid fPCA Method') + for i in range(0, n): + for j in range(0, no): + a[i, j] = np.sum(np.dot(vec[:, i] - vm, U[:, j])) + else: + raise Exception("Invalid fPCA Method") + + self.y_int = np.zeros((n, 2)) + X = (self.pca.coef.T@self.pca.coef) + df = self.Ntr - a.shape[1] - 1 + for ii in range(0, n): + self.y_pred[ii] = self.alpha + np.dot(a[ii, :], self.b) + interval = t.ppf(alpha/2, df)*self.stdev*np.sqrt(1+a[ii, :]@X@a[ii, :]) + self.y_int[ii, 0] = -interval + self.y_int[ii, 1] = interval - for ii in range(0,n): - self.y_pred[ii] = self.alpha + np.dot(a[ii,:],self.b) - if y is None: self.SSE = np.nan else: - self.SSE = np.sum((y-self.y_pred)**2) + self.SSE = np.sum((y - self.y_pred) ** 2) else: n = self.pca.coef.shape[0] self.y_pred = np.zeros(n) - for ii in range(0,n): - self.y_pred[ii] = self.alpha + np.dot(self.pca.coef[ii,:],self.b) - - self.SSE = np.sum((self.y-self.y_pred)**2) + self.y_int = np.zeros((n, 2)) + X = (self.pca.coef.T@self.pca.coef) + df = self.Ntr - a.shape[1] - 1 + for ii in range(0, n): + self.y_pred[ii] = self.alpha + np.dot(self.pca.coef[ii, :], self.b) + interval = t.ppf(alpha/2, df)*self.stdev*np.sqrt(1+self.pca.coef[ii, :]@X@self.pca.coef[ii, :]) + self.y_int[ii, 0] = -interval + self.y_int[ii, 1] = interval + + self.SSE = np.sum((self.y - self.y_pred) ** 2) return + class elastic_lpcr_regression: """ - This class provides elastic logistic pcr regression for functional + This class provides elastic logistic pcr regression for functional data using the SRVF framework accounting for warping - + Usage: obj = elastic_lpcr_regression(f,y,time) :param f: (M,N) % matrix defining N functions of M samples @@ -240,7 +270,7 @@ class elastic_lpcr_regression: :param Loss: logistic loss :param PC: probability of classification :param ylabels: predicted labels - + Author : J. D. Tucker (JDT) Date : 18-Mar-2018 """ @@ -255,14 +285,15 @@ def __init__(self, f, y, time): a = time.shape[0] if f.shape[0] != a: - raise Exception('Columns of f and time must be equal') + raise Exception("Columns of f and time must be equal") self.f = f self.y = y self.time = time - def calc_model(self, pca_method="combined", no=5, - smooth_data=False, sparam=25, parallel=False): + def calc_model( + self, pca_method="combined", no=5, smooth_data=False, sparam=25, parallel=False + ): """ This function identifies a logistic regression model with phase-variability using elastic pca @@ -278,43 +309,50 @@ def calc_model(self, pca_method="combined", no=5, """ if smooth_data: - self.f = fs.smooth_data(self.f,sparam) - + self.f = fs.smooth_data(self.f, sparam) + N1 = self.f.shape[1] # Align Data - self.warp_data = fs.fdawarp(self.f,self.time) + self.warp_data = fs.fdawarp(self.f, self.time) self.warp_data.srsf_align(parallel=parallel) # Calculate PCA - if pca_method=='combined': + if pca_method == "combined": self.pca = fpca.fdajpca(self.warp_data) - elif pca_method=='vert': + elif pca_method == "vert": self.pca = fpca.fdavpca(self.warp_data) - elif pca_method=='horiz': + elif pca_method == "horiz": self.pca = fpca.fdahpca(self.warp_data) else: - raise Exception('Invalid fPCA Method') + raise Exception("Invalid fPCA Method") self.pca.calc_fpca(no) - + # OLS using PCA basis lam = 0 R = 0 - Phi = np.ones((N1, no+1)) - Phi[:,1:(no+1)] = self.pca.coef + Phi = np.ones((N1, no + 1)) + Phi[:, 1 : (no + 1)] = self.pca.coef # Find alpha and beta using l_bfgs - b0 = np.zeros(no+1) - out = fmin_l_bfgs_b(rg.logit_loss, b0, fprime=rg.logit_gradient, - args=(Phi, self.y), pgtol=1e-10, maxiter=200, - maxfun=250, factr=1e-30) + b0 = np.zeros(no + 1) + out = fmin_l_bfgs_b( + rg.logit_loss, + b0, + fprime=rg.logit_gradient, + args=(Phi, self.y), + pgtol=1e-10, + maxiter=200, + maxfun=250, + factr=1e-30, + ) b = out[0] alpha = b[0] # compute the Loss - LL = rg.logit_loss(b,Phi,self.y) + LL = rg.logit_loss(b, Phi, self.y) - b = b[1:no+1] + b = b[1 : no + 1] self.alpha = alpha self.b = b @@ -343,82 +381,86 @@ def predict(self, newdata=None): M = self.time.shape[0] if newdata != None: - f = newdata['f'] - time = newdata['time'] - y = newdata['y'] - if newdata['smooth']: - sparam = newdata['sparam'] - f = fs.smooth_data(f,sparam) - - q1 = fs.f_to_srsf(f,time) + f = newdata["f"] + time = newdata["time"] + y = newdata["y"] + if newdata["smooth"]: + sparam = newdata["sparam"] + f = fs.smooth_data(f, sparam) + + q1 = fs.f_to_srsf(f, time) n = q1.shape[1] self.y_pred = np.zeros(n) mq = self.warp_data.mqn - fn = np.zeros((M,n)) - qn = np.zeros((M,n)) - gam = np.zeros((M,n)) - for ii in range(0,n): - gam[:,ii] = uf.optimum_reparam(mq,time,q1[:,ii],omethod) - fn[:,ii] = uf.warp_f_gamma(time,f[:,ii],gam[:,ii]) - qn[:,ii] = uf.f_to_srsf(fn[:,ii],time) - + fn = np.zeros((M, n)) + qn = np.zeros((M, n)) + gam = np.zeros((M, n)) + for ii in range(0, n): + gam[:, ii] = uf.optimum_reparam(mq, time, q1[:, ii], omethod) + fn[:, ii] = uf.warp_f_gamma(time, f[:, ii], gam[:, ii]) + qn[:, ii] = uf.f_to_srsf(fn[:, ii], time) + U = self.pca.U no = U.shape[1] - if self.pca.__class__.__name__ == 'fdajpca': - m_new = np.sign(fn[self.pca.id,:])*np.sqrt(np.abs(fn[self.pca.id,:])) + if self.pca.__class__.__name__ == "fdajpca": + m_new = np.sign(fn[self.pca.id, :]) * np.sqrt( + np.abs(fn[self.pca.id, :]) + ) qn1 = np.vstack((qn, m_new)) C = self.pca.C TT = self.time.shape[0] mu_g = self.pca.mu_g mu_psi = self.pca.mu_psi - vec = np.zeros((M,n)) - psi = np.zeros((TT,n)) + vec = np.zeros((M, n)) + psi = np.zeros((TT, n)) binsize = np.mean(np.diff(self.time)) - for i in range(0,n): - psi[:,i] = np.sqrt(np.gradient(gam[:,i],binsize)) - out, theta = geo.inv_exp_map(mu_psi, psi[:,i]) - vec[:,i] = out - - g = np.vstack((qn1, C*vec)) - a = np.zeros((n,no)) - for i in range(0,n): - for j in range(0,no): - tmp = (g[:,i]-mu_g) - a[i,j] = np.dot(tmp.T, U[:,j]) - - elif self.pca.__class__.__name__ == 'fdavpca': - m_new = np.sign(fn[self.pca.id,:])*np.sqrt(np.abs(fn[self.pca.id,:])) + for i in range(0, n): + psi[:, i] = np.sqrt(np.gradient(gam[:, i], binsize)) + out, theta = geo.inv_exp_map(mu_psi, psi[:, i]) + vec[:, i] = out + + g = np.vstack((qn1, C * vec)) + a = np.zeros((n, no)) + for i in range(0, n): + for j in range(0, no): + tmp = g[:, i] - mu_g + a[i, j] = np.dot(tmp.T, U[:, j]) + + elif self.pca.__class__.__name__ == "fdavpca": + m_new = np.sign(fn[self.pca.id, :]) * np.sqrt( + np.abs(fn[self.pca.id, :]) + ) qn1 = np.vstack((qn, m_new)) - a = np.zeros((n,no)) - for i in range(0,n): - for j in range(0,no): - tmp = (qn1[:,i]-self.pca.mqn) - a[i,j] = np.dot(tmp.T, U[:,j]) - - elif self.pca.__class__.__name__ == 'fdahpca': - a = np.zeros((n,no)) + a = np.zeros((n, no)) + for i in range(0, n): + for j in range(0, no): + tmp = qn1[:, i] - self.pca.mqn + a[i, j] = np.dot(tmp.T, U[:, j]) + + elif self.pca.__class__.__name__ == "fdahpca": + a = np.zeros((n, no)) mu_psi = self.pca.psi_mu - vec = np.zeros((M,n)) + vec = np.zeros((M, n)) TT = self.time.shape[0] - psi = np.zeros((TT,n)) + psi = np.zeros((TT, n)) binsize = np.mean(np.diff(self.time)) - for i in range(0,n): - psi[:,i] = np.sqrt(np.gradient(gam[:,i],binsize)) - out, theta = geo.inv_exp_map(mu_psi, psi[:,i]) - vec[:,i] = out - + for i in range(0, n): + psi[:, i] = np.sqrt(np.gradient(gam[:, i], binsize)) + out, theta = geo.inv_exp_map(mu_psi, psi[:, i]) + vec[:, i] = out + vm = self.pca.vec.mean(axis=1) - for i in range(0,n): - for j in range(0,no): - a[i,j] = np.sum(np.dot(vec[:,i]-vm,U[:,j])) - else: - raise Exception('Invalid fPCA Method') + for i in range(0, n): + for j in range(0, no): + a[i, j] = np.sum(np.dot(vec[:, i] - vm, U[:, j])) + else: + raise Exception("Invalid fPCA Method") + + for ii in range(0, n): + self.y_pred[ii] = self.alpha + np.sum(a[ii, :] * self.b) - for ii in range(0,n): - self.y_pred[ii] = self.alpha + np.sum(a[ii,:]*self.b) - if y is None: self.y_pred = rg.phi(self.y_pred) self.y_labels = np.ones(n) @@ -432,13 +474,13 @@ def predict(self, newdata=None): FP = np.sum(y[self.y_labels == -1] == 1) TN = np.sum(y[self.y_labels == -1] == -1) FN = np.sum(y[self.y_labels == 1] == -1) - self.PC = (TP+TN)/(TP+FP+FN+TN) + self.PC = (TP + TN) / (TP + FP + FN + TN) else: n = self.pca.coef.shape[0] self.y_pred = np.zeros(n) - for ii in range(0,n): - self.y_pred[ii] = self.alpha + np.dot(self.pca.coef[ii,:],self.b) - + for ii in range(0, n): + self.y_pred[ii] = self.alpha + np.dot(self.pca.coef[ii, :], self.b) + self.y_pred = rg.phi(self.y_pred) self.y_labels = np.ones(n) self.y_labels[self.y_pred < 0.5] = -1 @@ -446,16 +488,16 @@ def predict(self, newdata=None): FP = np.sum(self.y[self.y_labels == -1] == 1) TN = np.sum(self.y[self.y_labels == -1] == -1) FN = np.sum(self.y[self.y_labels == 1] == -1) - self.PC = (TP+TN)/(TP+FP+FN+TN) + self.PC = (TP + TN) / (TP + FP + FN + TN) return - + class elastic_mlpcr_regression: """ This class provides elastic multinomial logistic pcr regression for functional data using the SRVF framework accounting for warping - + Usage: obj = elastic_mlpcr_regression(f,y,time) :param f: (M,N) % matrix defining N functions of M samples @@ -484,7 +526,7 @@ def __init__(self, f, y, time): a = time.shape[0] if f.shape[0] != a: - raise Exception('Columns of f and time must be equal') + raise Exception("Columns of f and time must be equal") self.f = f self.y = y @@ -496,10 +538,11 @@ def __init__(self, f, y, time): self.n_classes = m self.Y = np.zeros((N1, m), dtype=int) for ii in range(0, N1): - self.Y[ii, y[ii]-1] = 1 + self.Y[ii, y[ii] - 1] = 1 - def calc_model(self, pca_method="combined", no=5, - smooth_data=False, sparam=25, parallel=False): + def calc_model( + self, pca_method="combined", no=5, smooth_data=False, sparam=25, parallel=False + ): """ This function identifies a logistic regression model with phase-variability using elastic pca @@ -518,44 +561,51 @@ def calc_model(self, pca_method="combined", no=5, """ if smooth_data: - self.f = fs.smooth_data(self.f,sparam) - + self.f = fs.smooth_data(self.f, sparam) + N1 = self.f.shape[1] # Align Data - self.warp_data = fs.fdawarp(self.f,self.time) + self.warp_data = fs.fdawarp(self.f, self.time) self.warp_data.srsf_align(parallel=parallel) # Calculate PCA - if pca_method=='combined': + if pca_method == "combined": self.pca = fpca.fdajpca(self.warp_data) - elif pca_method=='vert': + elif pca_method == "vert": self.pca = fpca.fdavpca(self.warp_data) - elif pca_method=='horiz': + elif pca_method == "horiz": self.pca = fpca.fdahpca(self.warp_data) else: - raise Exception('Invalid fPCA Method') + raise Exception("Invalid fPCA Method") self.pca.calc_fpca(no) - + # OLS using PCA basis lam = 0 R = 0 - Phi = np.ones((N1, no+1)) - Phi[:,1:(no+1)] = self.pca.coef + Phi = np.ones((N1, no + 1)) + Phi[:, 1 : (no + 1)] = self.pca.coef # Find alpha and beta using l_bfgs - b0 = np.zeros(self.n_classes*(no+1)) - out = fmin_l_bfgs_b(rg.mlogit_loss, b0, fprime=rg.mlogit_gradient, - args=(Phi, self.Y), pgtol=1e-10, maxiter=200, - maxfun=250, factr=1e-30) + b0 = np.zeros(self.n_classes * (no + 1)) + out = fmin_l_bfgs_b( + rg.mlogit_loss, + b0, + fprime=rg.mlogit_gradient, + args=(Phi, self.Y), + pgtol=1e-10, + maxiter=200, + maxfun=250, + factr=1e-30, + ) b = out[0] - B0 = b.reshape(no+1, self.n_classes) + B0 = b.reshape(no + 1, self.n_classes) alpha = B0[0, :] # compute the Loss - LL = rg.mlogit_loss(b,Phi,self.y) + LL = rg.mlogit_loss(b, Phi, self.y) - b = B0[1:no+1,:] + b = B0[1 : no + 1, :] self.alpha = alpha self.b = b @@ -585,124 +635,137 @@ def predict(self, newdata=None): M = self.time.shape[0] if newdata != None: - f = newdata['f'] - time = newdata['time'] - y = newdata['y'] - if newdata['smooth']: - sparam = newdata['sparam'] - f = fs.smooth_data(f,sparam) - - q1 = fs.f_to_srsf(f,time) + f = newdata["f"] + time = newdata["time"] + y = newdata["y"] + if newdata["smooth"]: + sparam = newdata["sparam"] + f = fs.smooth_data(f, sparam) + + q1 = fs.f_to_srsf(f, time) n = q1.shape[1] - self.y_pred = np.zeros((n,m)) + self.y_pred = np.zeros((n, m)) mq = self.warp_data.mqn - fn = np.zeros((M,n)) - qn = np.zeros((M,n)) - gam = np.zeros((M,n)) - for ii in range(0,n): - gam[:,ii] = uf.optimum_reparam(mq,time,q1[:,ii],omethod) - fn[:,ii] = uf.warp_f_gamma(time,f[:,ii],gam[:,ii]) - qn[:,ii] = uf.f_to_srsf(fn[:,ii],time) - + fn = np.zeros((M, n)) + qn = np.zeros((M, n)) + gam = np.zeros((M, n)) + for ii in range(0, n): + gam[:, ii] = uf.optimum_reparam(mq, time, q1[:, ii], omethod) + fn[:, ii] = uf.warp_f_gamma(time, f[:, ii], gam[:, ii]) + qn[:, ii] = uf.f_to_srsf(fn[:, ii], time) + U = self.pca.U no = U.shape[1] - if self.pca.__class__.__name__ == 'fdajpca': - m_new = np.sign(fn[self.pca.id,:])*np.sqrt(np.abs(fn[self.pca.id,:])) + if self.pca.__class__.__name__ == "fdajpca": + m_new = np.sign(fn[self.pca.id, :]) * np.sqrt( + np.abs(fn[self.pca.id, :]) + ) qn1 = np.vstack((qn, m_new)) C = self.pca.C TT = self.time.shape[0] mu_g = self.pca.mu_g mu_psi = self.pca.mu_psi - vec = np.zeros((M,n)) - psi = np.zeros((TT,n)) + vec = np.zeros((M, n)) + psi = np.zeros((TT, n)) binsize = np.mean(np.diff(self.time)) - for i in range(0,n): - psi[:,i] = np.sqrt(np.gradient(gam[:,i],binsize)) - out, theta = geo.inv_exp_map(mu_psi, psi[:,i]) - vec[:,i] = out - - g = np.vstack((qn1, C*vec)) - a = np.zeros((n,no)) - for i in range(0,n): - for j in range(0,no): - tmp = (g[:,i]-mu_g) - a[i,j] = np.dot(tmp.T, U[:,j]) - - elif self.pca.__class__.__name__ == 'fdavpca': - m_new = np.sign(fn[self.pca.id,:])*np.sqrt(np.abs(fn[self.pca.id,:])) + for i in range(0, n): + psi[:, i] = np.sqrt(np.gradient(gam[:, i], binsize)) + out, theta = geo.inv_exp_map(mu_psi, psi[:, i]) + vec[:, i] = out + + g = np.vstack((qn1, C * vec)) + a = np.zeros((n, no)) + for i in range(0, n): + for j in range(0, no): + tmp = g[:, i] - mu_g + a[i, j] = np.dot(tmp.T, U[:, j]) + + elif self.pca.__class__.__name__ == "fdavpca": + m_new = np.sign(fn[self.pca.id, :]) * np.sqrt( + np.abs(fn[self.pca.id, :]) + ) qn1 = np.vstack((qn, m_new)) - a = np.zeros((n,no)) - for i in range(0,n): - for j in range(0,no): - tmp = (qn1[:,i]-self.pca.mqn) - a[i,j] = np.dot(tmp.T, U[:,j]) - - elif self.pca.__class__.__name__ == 'fdahpca': - a = np.zeros((n,no)) + a = np.zeros((n, no)) + for i in range(0, n): + for j in range(0, no): + tmp = qn1[:, i] - self.pca.mqn + a[i, j] = np.dot(tmp.T, U[:, j]) + + elif self.pca.__class__.__name__ == "fdahpca": + a = np.zeros((n, no)) mu_psi = self.pca.psi_mu - vec = np.zeros((M,n)) + vec = np.zeros((M, n)) TT = self.time.shape[0] - psi = np.zeros((TT,n)) + psi = np.zeros((TT, n)) binsize = np.mean(np.diff(self.time)) - for i in range(0,n): - psi[:,i] = np.sqrt(np.gradient(gam[:,i],binsize)) - out, theta = geo.inv_exp_map(mu_psi, psi[:,i]) - vec[:,i] = out - + for i in range(0, n): + psi[:, i] = np.sqrt(np.gradient(gam[:, i], binsize)) + out, theta = geo.inv_exp_map(mu_psi, psi[:, i]) + vec[:, i] = out + vm = self.pca.vec.mean(axis=1) - for i in range(0,n): - for j in range(0,no): - a[i,j] = np.sum(np.dot(vec[:,i]-vm,U[:,j])) - else: - raise Exception('Invalid fPCA Method') - - for ii in range(0,n): - for jj in range(0,m): - self.y_pred[ii,jj] = self.alpha[jj] + np.sum(a[ii,:]*self.b[:,jj]) - - + for i in range(0, n): + for j in range(0, no): + a[i, j] = np.sum(np.dot(vec[:, i] - vm, U[:, j])) + else: + raise Exception("Invalid fPCA Method") + + for ii in range(0, n): + for jj in range(0, m): + self.y_pred[ii, jj] = self.alpha[jj] + np.sum( + a[ii, :] * self.b[:, jj] + ) + if y is None: - self.y_pred = rg.phi(self.y_pred.reshape((1,n*m))) - self.y_pred = self.y_pred.reshape((n,m)) - self.y_labels = np.argmax(self.y_pred,axis=1) + self.y_pred = rg.phi(self.y_pred.reshape((1, n * m))) + self.y_pred = self.y_pred.reshape((n, m)) + self.y_labels = np.argmax(self.y_pred, axis=1) self.PC = np.nan else: - self.y_pred = rg.phi(self.y_pred.reshape((1,n*m))) - self.y_pred = self.y_pred.reshape((n,m)) - self.y_labels = np.argmax(self.y_pred,axis=1) + self.y_pred = rg.phi(self.y_pred.reshape((1, n * m))) + self.y_pred = self.y_pred.reshape((n, m)) + self.y_labels = np.argmax(self.y_pred, axis=1) self.PC = np.zeros(m) - cls_set = np.arange(0,m) - for ii in range(0,m): - cls_sub = np.setdiff1d(cls_set,ii) + cls_set = np.arange(0, m) + for ii in range(0, m): + cls_sub = np.setdiff1d(cls_set, ii) TP = np.sum(y[self.y_labels == ii] == ii) - FP = np.sum(y[np.in1d(self.y_labels,cls_sub)] == ii) - TN = np.sum(y[np.in1d(self.y_labels,cls_sub)] == self.y_labels[np.in1d(self.y_labels,cls_sub)]) - FN = np.sum(np.in1d(y[self.y_labels==ii],cls_sub)) - self.PC[ii] = (TP+TN)/(TP+FP+FN+TN) - - self.PCo = np.sum(y == self.y_labels)/self.y_labels.shape[0] + FP = np.sum(y[np.in1d(self.y_labels, cls_sub)] == ii) + TN = np.sum( + y[np.in1d(self.y_labels, cls_sub)] + == self.y_labels[np.in1d(self.y_labels, cls_sub)] + ) + FN = np.sum(np.in1d(y[self.y_labels == ii], cls_sub)) + self.PC[ii] = (TP + TN) / (TP + FP + FN + TN) + + self.PCo = np.sum(y == self.y_labels) / self.y_labels.shape[0] else: n = self.pca.coef.shape[0] - self.y_pred = np.zeros((n,m)) - for ii in range(0,n): - for jj in range(0,m): - self.y_pred[ii,jj] = self.alpha[jj] + np.sum(self.pca.coef[ii,:]*self.b[:,jj]) - - self.y_pred = rg.phi(self.y_pred.reshape((1,n*m))) - self.y_pred = self.y_pred.reshape((n,m)) - self.y_labels = np.argmax(self.y_pred,axis=1) + self.y_pred = np.zeros((n, m)) + for ii in range(0, n): + for jj in range(0, m): + self.y_pred[ii, jj] = self.alpha[jj] + np.sum( + self.pca.coef[ii, :] * self.b[:, jj] + ) + + self.y_pred = rg.phi(self.y_pred.reshape((1, n * m))) + self.y_pred = self.y_pred.reshape((n, m)) + self.y_labels = np.argmax(self.y_pred, axis=1) self.PC = np.zeros(m) - cls_set = np.arange(0,m) - for ii in range(0,m): - cls_sub = np.setdiff1d(cls_set,ii) + cls_set = np.arange(0, m) + for ii in range(0, m): + cls_sub = np.setdiff1d(cls_set, ii) TP = np.sum(self.y[self.y_labels == ii] == ii) - FP = np.sum(self.y[np.in1d(self.y_labels,cls_sub)] == ii) - TN = np.sum(self.y[np.in1d(self.y_labels,cls_sub)] == self.y_labels[np.in1d(self.y_labels,cls_sub)]) - FN = np.sum(np.in1d(y[self.y_labels==ii],cls_sub)) - self.PC[ii] = (TP+TN)/(TP+FP+FN+TN) - - self.PCo = np.sum(y == self.y_labels)/self.y_labels.shape[0] + FP = np.sum(self.y[np.in1d(self.y_labels, cls_sub)] == ii) + TN = np.sum( + self.y[np.in1d(self.y_labels, cls_sub)] + == self.y_labels[np.in1d(self.y_labels, cls_sub)] + ) + FN = np.sum(np.in1d(y[self.y_labels == ii], cls_sub)) + self.PC[ii] = (TP + TN) / (TP + FP + FN + TN) + + self.PCo = np.sum(y == self.y_labels) / self.y_labels.shape[0] return