From 68f7790792454bb638816d3e03cfcf6f3216716b Mon Sep 17 00:00:00 2001 From: "J. Derek Tucker" Date: Mon, 29 Sep 2025 12:09:30 -0600 Subject: [PATCH 1/4] add PPD --- fdasrsf/PPD.py | 298 ++++++++++++++++++++++++++++++++++++++++ fdasrsf/time_warping.py | 82 +++++++++++ 2 files changed, 380 insertions(+) create mode 100644 fdasrsf/PPD.py diff --git a/fdasrsf/PPD.py b/fdasrsf/PPD.py new file mode 100644 index 0000000..698c098 --- /dev/null +++ b/fdasrsf/PPD.py @@ -0,0 +1,298 @@ +import numpy as np +from scipy.signal import find_peaks +from scipy.spatial.distance import pdist +from scipy.cluster.hierarchy import linkage, fcluster +import matplotlib.pyplot as plt +from matplotlib.patches import Rectangle +from mpl_toolkits.mplot3d import Axes3D +from matplotlib.colors import LinearSegmentedColormap + + +def getPPDinfo(t, Fa, lam, th): + n_lams = lam.shape[0] + + # Compute the mean of each function in FN_temp over its rows + FNm = np.zeros(t.shape[0], len(Fa)) + FNm[:, 0] = Fa[0].mean(axis=1) + + # Find indices of local maxima in the first function's mean + idxMaxFirst = find_peaks(FNm[:, 0]) + + # Initalize labels and Locations for the first function + Labels = [] + Locs = [] + Labels.append(np.arange(0, idxMaxFirst.shape[0])) + Locs.append(idxMaxFirst) + + # Initialize the maximum label number + labelMax = Labels[0].max() + + # Process each function to assign labels and locate peaks + for i in range(n_lams - 1): + currentLabel = Labels[i] + Labelst, labelMax = peak_successor(Fa[i], Fa[i+1], currentLabel, labelMax, 1) + Labels.append(Labelst) + + # Find peak locations in the next function's mean + FnmNextMean = Fa[i+1].mean(axis=1) + idxMaxNext = find_peaks(FnmNextMean) + Locs.append(idxMaxNext) + + # Update the mean function + FNm[:, i+1] = FnmNextMean + + IndicatorMatrix, Heights, Heights2 = PreprocessingForPPD(t, lam, Labels, Locs, labelMax, FNm, th) + + if np.all(np.isnan(Heights2)): + NameError("All peaks are ignored. A smaller threshold is required.") + + return (IndicatorMatrix, Heights, Locs, Labels, FNm) + + +def peak_successor(f1, f2, labels1, labelMax, smooth_parameter): + # Combine f1 and f2 into a 3D array and compute the mean across + # the second dimension + F = np.concatenate((f1[:,:,np.newaxis], f2[:,:,np.newaxis]), axis=2) + fm = F.mean(axis=1) + + # compute peak ranges and labes for fm[:,0] + ranges, idx_max1 = computePeakRanges(fm[:, 0]) + + # compute indices of local maxima in fm[:,1] + idx_max2 = find_peaks(fm[:,1]) + + if idx_max1.size==0: + labels2 = labelMax + np.arange(1, idx_max2.shape[0]+1) + else: + # assign labels to peaks in fm[:,1] based on matching ranges in fm[:,0] + labels2 = assignLabelsToPeaks(idx_max2, ranges, idx_max1, labels1) + + # ensure no overlapping labels in lables2 + labels2 = resolveOverlappingLabels(labels2, idx_max1, idx_max2, labels1) + + # assign new labels to unmatched peaks in fm[:,1] + unmatched = labels2 == 0 + labels2[unmatched] = labelMax + np.arange(1,np.sum(unmatched)+1) + labelMax = labelMax + sum(unmatched) + + return(labels2, labelMax) + + +def computePeakRanges(data): + # computes peak ranges defined by adjacent minima in the data + idx_max = find_peaks(data) + if idx_max.size == 0: + idx_max = [] + ranges = [] + + idx_min = find_peaks(-1*data) + idx_min = np.unique(np.concatenate((0, data.shape[0], idx_min))) + ranges = np.zeros((idx_max.shape[0], 2)) + for i in range(idx_max.shape[0]): + ranges[i,0] = np.max(idx_min[idx_minidx_max[i]]) + # remove degenerate ranges + idx = ranges[:,0] == ranges[:,1] + ranges = np.delete(ranges, idx, 0) + return(ranges, idx_max) + + +def assignLabelsToPeaks(idx_max2, ranges, idx_max1, labels1): + # assigns labels to peaks in idx_max2 based on matching ranges in idx_max1 + labels = np.zeros(idx_max2.shape[0]) + for i in range(idx_max2.shape): + # find the range in fm[:,0] that contains the current peak in fm[:,1] + in_range = np.logical_and(idx_max2[i] >= ranges[:,0], idx_max2[i] <= ranges[:,1]) + matching_ranges = np.argwhere(in_range) + + if matching_ranges.size != 0: + # choose the closest peak if multiple ranges match + if matching_ranges.size > 1: + tmp = np.abs(idx_max1[matching_ranges] - idx_max2[i]) + closest_idx = tmp.argmin() + else: + matching_range = matching_ranges.copy() + + labels[i] = labels1[matching_range] + + return labels + + +def resolveOverlappingLabels(labels, idx_max1, idx_max2, labels1): + # ensures no overlapping labels in the assigned labels + unique_labels = np.unique(labels[labels > 0]) + for label in unique_labels: + duplicates = np.argwhere(labels==label) + if duplicates.size > 1: + # keep the closest peak and reset others + distances = np.abs(idx_max1[labels1==label] - idx_max2[duplicates]) + min_idx = distances.arg_min() + duplicates = np.delete(duplicates, min_idx, 0) + labels[duplicates] = 0 + return labels + + +def PreprocessingForPPD(t, lam, Labels, Locs, labelMax, FNm, th): + K = lam.shape[0] + # initialize output matrices + IndicatorMatrix = np.full((K, labelMax), np.nan) + curvatures = np.zeros((K, labelMax)) + Heights = np.full((K, labelMax), np.nan) + + # assume t is uniformly spaced; compute time step + dx = t[1] - t[0] + + # process each function + for i in range(K): + # extract the function values for the current parameter lambda + fnm = FNm[:, i] + + # compute negative curvature (sec) + negCurvature = -1*np.gradient(np.gradient(fnm,dx), dx) + + # Ensure non-negative curvature values + negCurvature[negCurvature < 0] = 0 + + # normalize negative curvature to [0,1] if possible + maxNegCurvature = negCurvature.max() + if maxNegCurvature > 0: + negCurvature = negCurvature / negCurvature + + # retrieve peak locations and labels for the current function + locsCurrent = Locs[i] + labelsCurrent = Labels[i] + + # select negative curvature values at specified peak locations + negCurvSelected = negCurvature[locsCurrent] + + # update curvatures and heights matrices at the appropriate labels + curvatures[i, labelsCurrent] = negCurvSelected + Heights[i, labelsCurrent] = fnm[locsCurrent] + + # apply threshold to select significant peaks based on curvature + significantLabels = labelsCurrent[negCurvSelected >= th] + + # update the indicator matrix for signficiant peaks + IndicatorMatrix[i, significantLabels] = 1 + + # compute heighs2 by multiply heighs with the indicator matrix + Heights2 = IndicatorMatrix * Heights + + return(IndicatorMatrix, Heights, Heights2) + + +def getPersistentPeaks(IndicatorMatrix): + + if IndicatorMatrix.size == 0: + Clt2 = [] + return Clt2 + + # count the number of ones (occurences) for each peak (ignore nans) + occurrenceCounts = np.nansum(IndicatorMatrix, 1) + + data = np.append(occurrenceCounts, 0) + if np.all(data == 0) or np.unique(occurrenceCounts).size == 1: + Clt2 = [] + return Clt2 + + # compute pairwise distances between observations + pairwiseDistances = pdist(data, metric='euclidean') + Y = linkage(pairwiseDistances, 'ward') + + # Cluster the data into the specified number of clusters + clusterAssignments = fcluster(Y, 2, 'maxclust') + referenceCluster = clusterAssignments[-1] + clusterAssignments = clusterAssignments[:-1] + + # identify indices wehre cluster assignmetns differ from the reference + Clt2 = np.argwhere(clusterAssignments != referenceCluster) + + return Clt2 + + +def drawPPDBarChart(IndicatorMatrix, Heights, lam, idx_opt): + lam_diff = lam[1] - lam[0] + len_lam, labelMax = IndicatorMatrix.shape + + fig, ax = plt.subplots() + + for i in range(len_lam): + label_all_peaks = np.argwhere(np.logical_not(np.isnan(Heights[i,:]))) + label_persistent_peaks = np.argwhere(np.logical_not(np.isnan(IndicatorMatrix[i,:]))) + + if label_all_peaks.sum() == 0: + continue + + for j in range(label_all_peaks.shape[0]): + x = lam[i] + y = label_all_peaks[j] + + if y in label_persistent_peaks: + rect = Rectangle((x, y), lam_diff, 1, facecolor='black', edgecolor='none') + else: + rect = Rectangle((x, y), lam_diff, 0.5, 1, facecolor='grey', edgecolor='none') + + ax.add_patch(rect) + + for j in range(labelMax): + plt.hlines(y=j+0.5, xmin=plt.xlim()[0], xmax=plt.xlim()[1]) + + plt.axvline(x=lam[idx_opt], color="red", linestyles='dashed', linewdith=2) + + plt.xlim((lam[0], lam[-1])) + plt.ylim((0.5, labelMax+0.5)) + + plt.yticks(np.arange(1,labelMax+1)) + + plt.xlabel('$\lambda$') + plt.ylabel('Peak Index') + + +def drawPPDSurface(t,lam,FNm,Heights,Locs,IndicatorMatrix,Labels,idx_opt): + parula_colors = [ + (0.24, 0.15, 0.66), # Dark blue + (0.26, 0.44, 0.82), # Medium blue + (0.28, 0.70, 0.81), # Cyan + (0.60, 0.80, 0.44), # Greenish-yellow + (0.99, 0.81, 0.19) # Yellow + ] + + # Create the custom colormap + parula_cmap = LinearSegmentedColormap.from_list('parula_custom', parula_colors) + + + n_lams = lam.shape[0] + labelMax = IndicatorMatrix.shape[1] + + LocationMatrix_full = np.full((n_lams, labelMax), np.nan) + for i in range(n_lams): + LocationMatrix_full[i, Labels[i]] = Locs[i] + + LocationMatrix_sig = LocationMatrix_full * IndicatorMatrix + HeighMatrix_full = Heights.copy() + HeightMatrix_sig = HeighMatrix_full * IndicatorMatrix + + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + X, Y = np.meshgrid(t, lam) + ax.plot_surface(X, Y, FNm.T, cmap=parula_cmap, linewidth=0.5) + + plt.xlim(t[0], t[-1]) + + ax.plot(t, lam[idx_opt]*np.ones(t.shape[0]), FNm[:, idx_opt], linestyle='dashed', linewidth=2, color='magenta') + # loop through each label + for j in range(labelMax): + + # find non-nan indices for full location matrix and plot + idx_full = np.argwhere(np.logical_not(np.isnan(LocationMatrix_full[:, j]))) + ax.plot(t[LocationMatrix_full[idx_full, j]], lam[idx_full], HeighMatrix_full[idx_full, j], marker='o', color='black', linewidth=1.5, markersize=5, mfc='black') + + # find non-nan indices for significant location matrix and plot + idx_sig = np.argwhere(np.logical_not(np.isnan(LocationMatrix_sig[:, j]))) + ax.plot(t[LocationMatrix_sig[:, j]], lam[idx_sig], HeightMatrix_sig[idx_sig, j], linestyle='dashed', linewidth=2, color='magenta', markersize=6) + + plt.xlabel('$t$') + plt.ylabel('$\lambda$') + plt.zlabel('$\hat{g}_\lambda(t)$') + + plt.grid(True) diff --git a/fdasrsf/time_warping.py b/fdasrsf/time_warping.py index 2272ef1..d631b3c 100644 --- a/fdasrsf/time_warping.py +++ b/fdasrsf/time_warping.py @@ -11,11 +11,13 @@ import fdasrsf.fPCA as fpca import fdasrsf.geometry as geo from fdasrsf.gp import gp_posterior +from fdasrsf.PPD import getPPDinfo, getPersistentPeaks, drawPPDBarChart, drawPPDSurface from scipy.integrate import trapezoid, cumulative_trapezoid from scipy.interpolate import interp1d from scipy.linalg import svd, cholesky from scipy.cluster.hierarchy import linkage, fcluster from scipy.spatial.distance import squareform, pdist +from scipy.signal import find_peaks from numpy.linalg import norm, inv from numpy.random import rand, normal from joblib import Parallel, delayed @@ -358,6 +360,86 @@ def srsf_align( self.qun = qun[0:r] return + + def ppd(self, + max_lam = 2, + num_lam=10, + pt=0.15, + method="mean", + omethod="DP2", + center=True, + smoothdata=False, + MaxItr=20, + parallel=True, + pen="roughness", + cores=-1, + grid_dim=7, + verbose=True): + """ + This computes the peak persistance diagram over a range of + lambda. This can help determine the proper elasticity (penalty). + This can be slow and recommended to run in parallel. + + + :param max_lam: max regularization parameter (default 2) + :param num_lam: number of steps (default 10) + :param pt: the percentile of negative curvature of raw data (default .15) + :param method: (string) warp calculate Karcher Mean or Median + (options = "mean" or "median") (default="mean") + :param omethod: optimization method (DP, DP2, RBFGS, cRBFGS) + (default = DP2) + :param center: center warping functions (default = T) + :param smoothdata: Smooth the data using a box filter (default = F) + :param MaxItr: Maximum number of iterations (default = 20) + :param parallel: run in parallel (default = F) + :param penalty: penalty type (default="roughness") options are "roughness", + "l2gam", "l2psi", "geodesic". Only roughness implemented + in all methods. To use others method needs to be "RBFGS" + or "cRBFGS" + :param cores: number of cores for parallel (default = -1 (all)) + :param grid_dim: size of the grid, for the DP2 method only + (default = 7) + :param verbose: print status output (default = T) + :param thresh: cost function threshold (default = .01) + :type lam: double + :type smoothdata: bool + + """ + + # Obtain a lambda candidate set + lam_vec = np.linspace(0, max_lam, num_lam) + fns = [] + for i in range(num_lam): + self.srsf_align(method, + omethod, + center, + smoothdata, + MaxItr, + parallel, + lam_vec[i], + pen, + cores, + grid_dim, + verbose) + fns.append(self.fn) + + # Peak Persistent Diagrams + # Get the threshold for significant peak + diff_t = np.mean(np.diff(self.time)) + taus = [] + for i in range(self.f.shape[1]): + idx = find_peaks(self.f[:,i]) + df2 = np.gradient(np.gradient(self[:,i], diff_t), diff_t) + tau = -df2 / np.maximum(-df2) + tau[tau < 0] = 0 + taus.append(tau[idx]) + + th = np.percentile(taus, pt*100) + + IndicatorMatrix, Heights, Locs, Labels, FNm = getPPDinfo(self.time, fns, lam_vec, th) + + + def plot(self): """ From d1a44a813748f957eb91d4a534196b5e1ef2453a Mon Sep 17 00:00:00 2001 From: "J. Derek Tucker" Date: Mon, 29 Sep 2025 12:40:36 -0600 Subject: [PATCH 2/4] more updates to PPD code --- fdasrsf/PPD.py | 6 +++--- fdasrsf/time_warping.py | 36 ++++++++++++++++++++++++++++++++---- 2 files changed, 35 insertions(+), 7 deletions(-) diff --git a/fdasrsf/PPD.py b/fdasrsf/PPD.py index 698c098..ffd2026 100644 --- a/fdasrsf/PPD.py +++ b/fdasrsf/PPD.py @@ -244,7 +244,7 @@ def drawPPDBarChart(IndicatorMatrix, Heights, lam, idx_opt): plt.yticks(np.arange(1,labelMax+1)) - plt.xlabel('$\lambda$') + plt.xlabel('$\\lambda$') plt.ylabel('Peak Index') @@ -292,7 +292,7 @@ def drawPPDSurface(t,lam,FNm,Heights,Locs,IndicatorMatrix,Labels,idx_opt): ax.plot(t[LocationMatrix_sig[:, j]], lam[idx_sig], HeightMatrix_sig[idx_sig, j], linestyle='dashed', linewidth=2, color='magenta', markersize=6) plt.xlabel('$t$') - plt.ylabel('$\lambda$') - plt.zlabel('$\hat{g}_\lambda(t)$') + plt.ylabel('$\\lambda$') + plt.zlabel('$\\hat{g}_\\lambda(t)$') plt.grid(True) diff --git a/fdasrsf/time_warping.py b/fdasrsf/time_warping.py index d631b3c..ff1893d 100644 --- a/fdasrsf/time_warping.py +++ b/fdasrsf/time_warping.py @@ -428,18 +428,46 @@ def ppd(self, diff_t = np.mean(np.diff(self.time)) taus = [] for i in range(self.f.shape[1]): - idx = find_peaks(self.f[:,i]) - df2 = np.gradient(np.gradient(self[:,i], diff_t), diff_t) - tau = -df2 / np.maximum(-df2) + idx, _ = find_peaks(self.f[:,i]) + df2 = np.gradient(np.gradient(self.f[:,i], diff_t), diff_t) + tau = -df2 / np.max(-df2) tau[tau < 0] = 0 taus.append(tau[idx]) th = np.percentile(taus, pt*100) IndicatorMatrix, Heights, Locs, Labels, FNm = getPPDinfo(self.time, fns, lam_vec, th) - + persistent_peak_labels = getPersistentPeaks(IndicatorMatrix.T) + + # choose optimal lambda + ref_row = np.zeros(IndicatorMatrix.shape[1]) + ref_row[persistent_peak_labels] = 1 + + n_lams = IndicatorMatrix.shape[0] + hamming_distances = np.zeros(n_lams) + exact_match_indices = [] + for i in range(n_lams): + comp = np.logical_not(np.isnan(IndicatorMatrix[i, :])) + + # check for exact match + if comp == ref_row: + exact_match_indices.append(i) + + # calculate hamming distance + hamming_distances[i] = np.sum(comp != ref_row) + + if len(exact_match_indices) != 0: + idx_opt = min(exact_match_indices) + else: + min_distance_indices = hamming_distances.argmin() + idx_opt = min_distance_indices.min() + # draw PPD bar chart + drawPPDBarChart(IndicatorMatrix,Heights,lam_vec,idx_opt) + #draw PPD surface + drawPPDSurface(obj.time,lam_vec,FNm,Heights,Locs,IndicatorMatrix,Labels,idx_opt) + self.lam_opt = lam_vec[idx_opt] def plot(self): """ From 57ffa0fb4ac4f1b87fa6fc04a756450ae3e7b49e Mon Sep 17 00:00:00 2001 From: "J. Derek Tucker" Date: Mon, 29 Sep 2025 15:33:50 -0600 Subject: [PATCH 3/4] finished debugging --- fdasrsf/PPD.py | 47 ++++++++++++++++++++++------------------- fdasrsf/time_warping.py | 4 ++-- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/fdasrsf/PPD.py b/fdasrsf/PPD.py index ffd2026..0e0121e 100644 --- a/fdasrsf/PPD.py +++ b/fdasrsf/PPD.py @@ -12,16 +12,16 @@ def getPPDinfo(t, Fa, lam, th): n_lams = lam.shape[0] # Compute the mean of each function in FN_temp over its rows - FNm = np.zeros(t.shape[0], len(Fa)) + FNm = np.zeros((t.shape[0], len(Fa))) FNm[:, 0] = Fa[0].mean(axis=1) # Find indices of local maxima in the first function's mean - idxMaxFirst = find_peaks(FNm[:, 0]) + idxMaxFirst, _ = find_peaks(FNm[:, 0]) # Initalize labels and Locations for the first function Labels = [] Locs = [] - Labels.append(np.arange(0, idxMaxFirst.shape[0])) + Labels.append(np.arange(1, idxMaxFirst.shape[0]+1, dtype=int)) Locs.append(idxMaxFirst) # Initialize the maximum label number @@ -35,7 +35,7 @@ def getPPDinfo(t, Fa, lam, th): # Find peak locations in the next function's mean FnmNextMean = Fa[i+1].mean(axis=1) - idxMaxNext = find_peaks(FnmNextMean) + idxMaxNext, _ = find_peaks(FnmNextMean) Locs.append(idxMaxNext) # Update the mean function @@ -59,10 +59,10 @@ def peak_successor(f1, f2, labels1, labelMax, smooth_parameter): ranges, idx_max1 = computePeakRanges(fm[:, 0]) # compute indices of local maxima in fm[:,1] - idx_max2 = find_peaks(fm[:,1]) + idx_max2, _ = find_peaks(fm[:,1]) if idx_max1.size==0: - labels2 = labelMax + np.arange(1, idx_max2.shape[0]+1) + labels2 = labelMax + np.arange(1, idx_max2.shape[0]+1, dtype=int) else: # assign labels to peaks in fm[:,1] based on matching ranges in fm[:,0] labels2 = assignLabelsToPeaks(idx_max2, ranges, idx_max1, labels1) @@ -72,7 +72,7 @@ def peak_successor(f1, f2, labels1, labelMax, smooth_parameter): # assign new labels to unmatched peaks in fm[:,1] unmatched = labels2 == 0 - labels2[unmatched] = labelMax + np.arange(1,np.sum(unmatched)+1) + labels2[unmatched] = labelMax + np.arange(1,np.sum(unmatched)+1, dtype=int) labelMax = labelMax + sum(unmatched) return(labels2, labelMax) @@ -80,13 +80,15 @@ def peak_successor(f1, f2, labels1, labelMax, smooth_parameter): def computePeakRanges(data): # computes peak ranges defined by adjacent minima in the data - idx_max = find_peaks(data) + idx_max, _ = find_peaks(data) if idx_max.size == 0: idx_max = [] ranges = [] - idx_min = find_peaks(-1*data) - idx_min = np.unique(np.concatenate((0, data.shape[0], idx_min))) + idx_min, _ = find_peaks(-1*data) + idx_min = np.insert(idx_min, 0, 0) + idx_min = np.insert(idx_min, 0, data.shape[0]) + idx_min = np.unique(idx_min) ranges = np.zeros((idx_max.shape[0], 2)) for i in range(idx_max.shape[0]): ranges[i,0] = np.max(idx_min[idx_min= ranges[:,0], idx_max2[i] <= ranges[:,1]) matching_ranges = np.argwhere(in_range) @@ -110,6 +112,7 @@ def assignLabelsToPeaks(idx_max2, ranges, idx_max1, labels1): if matching_ranges.size > 1: tmp = np.abs(idx_max1[matching_ranges] - idx_max2[i]) closest_idx = tmp.argmin() + matching_range = matching_ranges[closest_idx] else: matching_range = matching_ranges.copy() @@ -166,14 +169,14 @@ def PreprocessingForPPD(t, lam, Labels, Locs, labelMax, FNm, th): negCurvSelected = negCurvature[locsCurrent] # update curvatures and heights matrices at the appropriate labels - curvatures[i, labelsCurrent] = negCurvSelected - Heights[i, labelsCurrent] = fnm[locsCurrent] + curvatures[i, labelsCurrent-1] = negCurvSelected + Heights[i, labelsCurrent-1] = fnm[locsCurrent] # apply threshold to select significant peaks based on curvature significantLabels = labelsCurrent[negCurvSelected >= th] # update the indicator matrix for signficiant peaks - IndicatorMatrix[i, significantLabels] = 1 + IndicatorMatrix[i, significantLabels-1] = 1 # compute heighs2 by multiply heighs with the indicator matrix Heights2 = IndicatorMatrix * Heights @@ -225,7 +228,7 @@ def drawPPDBarChart(IndicatorMatrix, Heights, lam, idx_opt): for j in range(label_all_peaks.shape[0]): x = lam[i] - y = label_all_peaks[j] + y = label_all_peaks[j,0] if y in label_persistent_peaks: rect = Rectangle((x, y), lam_diff, 1, facecolor='black', edgecolor='none') @@ -237,10 +240,10 @@ def drawPPDBarChart(IndicatorMatrix, Heights, lam, idx_opt): for j in range(labelMax): plt.hlines(y=j+0.5, xmin=plt.xlim()[0], xmax=plt.xlim()[1]) - plt.axvline(x=lam[idx_opt], color="red", linestyles='dashed', linewdith=2) + plt.axvline(x=lam[idx_opt], color="red", linestyle='dashed', linewidth=2) plt.xlim((lam[0], lam[-1])) - plt.ylim((0.5, labelMax+0.5)) + plt.ylim((0.5, labelMax)) plt.yticks(np.arange(1,labelMax+1)) @@ -259,16 +262,16 @@ def drawPPDSurface(t,lam,FNm,Heights,Locs,IndicatorMatrix,Labels,idx_opt): # Create the custom colormap parula_cmap = LinearSegmentedColormap.from_list('parula_custom', parula_colors) - n_lams = lam.shape[0] labelMax = IndicatorMatrix.shape[1] - LocationMatrix_full = np.full((n_lams, labelMax), np.nan) + LocationMatrix_full = np.full((n_lams, labelMax), np.nan, dtype=int) for i in range(n_lams): - LocationMatrix_full[i, Labels[i]] = Locs[i] + LocationMatrix_full[i, Labels[i]-1] = Locs[i] LocationMatrix_sig = LocationMatrix_full * IndicatorMatrix + LocationMatrix_sig = LocationMatrix_sig.astype(int) HeighMatrix_full = Heights.copy() HeightMatrix_sig = HeighMatrix_full * IndicatorMatrix @@ -293,6 +296,6 @@ def drawPPDSurface(t,lam,FNm,Heights,Locs,IndicatorMatrix,Labels,idx_opt): plt.xlabel('$t$') plt.ylabel('$\\lambda$') - plt.zlabel('$\\hat{g}_\\lambda(t)$') + ax.set_zlabel('$\\hat{g}_\\lambda(t)$') plt.grid(True) diff --git a/fdasrsf/time_warping.py b/fdasrsf/time_warping.py index ff1893d..fc5e509 100644 --- a/fdasrsf/time_warping.py +++ b/fdasrsf/time_warping.py @@ -451,7 +451,7 @@ def ppd(self, comp = np.logical_not(np.isnan(IndicatorMatrix[i, :])) # check for exact match - if comp == ref_row: + if np.all(comp == ref_row): exact_match_indices.append(i) # calculate hamming distance @@ -466,7 +466,7 @@ def ppd(self, # draw PPD bar chart drawPPDBarChart(IndicatorMatrix,Heights,lam_vec,idx_opt) #draw PPD surface - drawPPDSurface(obj.time,lam_vec,FNm,Heights,Locs,IndicatorMatrix,Labels,idx_opt) + drawPPDSurface(self.time,lam_vec,FNm,Heights,Locs,IndicatorMatrix,Labels,idx_opt) self.lam_opt = lam_vec[idx_opt] def plot(self): From 1094ec886fecbb6f2961986eaeee9d2b5c189cbc Mon Sep 17 00:00:00 2001 From: "J. Derek Tucker" Date: Mon, 29 Sep 2025 15:34:21 -0600 Subject: [PATCH 4/4] update changes --- CHANGES.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGES.txt b/CHANGES.txt index 2022d17..6df18f0 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -1,3 +1,4 @@ +v2.6.X, 2025-09-30 -- bugfixes, add PPD v2.6.4, 2025-07-03 -- bugfixes v2.6.3, 2025-06-27 -- bugfixes v2.6.2, 2025-04-03 -- bugfixes and add pns