Start landscape optimization

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):
         """