diff options
author | wreise <wojciech.reise@epfl.ch> | 2022-06-02 15:38:38 +0200 |
---|---|---|
committer | wreise <wojciech.reise@epfl.ch> | 2022-06-02 15:38:38 +0200 |
commit | a70ad064ac3e0aee5c9b8084d35b7ce329c6bddc (patch) | |
tree | 1b2bba69bea5be293bdc0275bf51f49f555361f6 /src/python/gudhi/representations | |
parent | 775cd158c869f75a4925ac6f7a94a1847f7d7788 (diff) |
Start landscape optimization
Diffstat (limited to 'src/python/gudhi/representations')
-rw-r--r-- | src/python/gudhi/representations/vector_methods.py | 64 |
1 files changed, 22 insertions, 42 deletions
diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py index c250c98c..acc62943 100644 --- a/src/python/gudhi/representations/vector_methods.py +++ b/src/python/gudhi/representations/vector_methods.py @@ -141,6 +141,16 @@ class Landscape(BaseEstimator, TransformerMixin): self.nan_in_range = np.isnan(np.array(self.sample_range)) self.new_resolution = self.resolution + self.nan_in_range.sum() + landscape_formula = "(-1)*ReLU(heights - Abs(x_values - midpoints))" + variables = [ + "heights = Vi(1)", + "midpoints = Vi(1)", + "x_values = Vj(1)", + ] + from pykeops.numpy import Genred + self.landscape = Genred(landscape_formula, variables, reduction_op="KMin", + axis=0, opt_arg=self.num_landscapes) + def fit(self, X, y=None): """ Fit the Landscape class on a list of persistence diagrams: if any of the values in **sample_range** is numpy.nan, replace it with the corresponding value computed on the given list of persistence diagrams. @@ -162,53 +172,23 @@ class Landscape(BaseEstimator, TransformerMixin): Returns: numpy array with shape (number of diagrams) x (number of samples = **num_landscapes** x **resolution**): output persistence landscapes. """ - num_diag, Xfit = len(X), [] - x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution) - step_x = x_values[1] - x_values[0] - - for i in range(num_diag): - - diagram, num_pts_in_diag = X[i], X[i].shape[0] - ls = np.zeros([self.num_landscapes, self.new_resolution]) - - events = [] - for j in range(self.new_resolution): - events.append([]) - - for j in range(num_pts_in_diag): - [px,py] = diagram[j,:2] - min_idx = np.clip(np.ceil((px - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - mid_idx = np.clip(np.ceil((0.5*(py+px) - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - max_idx = np.clip(np.ceil((py - self.sample_range[0]) / step_x).astype(int), 0, self.new_resolution) - - if min_idx < self.new_resolution and max_idx > 0: - - landscape_value = self.sample_range[0] + min_idx * step_x - px - for k in range(min_idx, mid_idx): - events[k].append(landscape_value) - landscape_value += step_x - - landscape_value = py - self.sample_range[0] - mid_idx * step_x - for k in range(mid_idx, max_idx): - events[k].append(landscape_value) - landscape_value -= step_x - - for j in range(self.new_resolution): - events[j].sort(reverse=True) - for k in range( min(self.num_landscapes, len(events[j])) ): - ls[k,j] = events[j][k] + Xfit = [] + x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.new_resolution) + for i, diag in enumerate(X): + midpoints, heights = (diag[:, 0] + diag[:, 1]) / 2., (diag[:, 1] - diag[:, 0]) / 2. + tent_functions = np.maximum(heights[None, :] - np.abs(x_values[:, None] - midpoints[None, :]), 0) + tent_functions.partition(diag.shape[0] - self.num_landscapes, axis=1) + landscapes = np.sort(tent_functions[-self.num_landscapes:, :])[::-1].T if self.nan_in_range[0]: - ls = ls[:,1:] + landscapes = landscapes[:,1:] if self.nan_in_range[1]: - ls = ls[:,:-1] - ls = np.sqrt(2)*np.reshape(ls,[1,-1]) - Xfit.append(ls) - - Xfit = np.concatenate(Xfit,0) + landscapes = landscapes[:,:-1] + landscapes = np.sqrt(2)*np.reshape(landscapes, [1, -1]) + Xfit.append(landscapes) - return Xfit + return np.stack(Xfit, axis=0) def __call__(self, diag): """ |