Refactor p_norm throughout

The former approximate p_norm method used vector norm from numpy, giving incorrect results. Refactor exact p_norm and move to auxiliary as _p_norm. Now both exact and approximate p_norm methods call _p_norm

Author: catanzaromj

persim/landscapes/approximate.py

@@ -5,7 +5,7 @@
 import numpy as np
 from operator import itemgetter
 from .base import PersLandscape
-from .auxiliary import union_vals, ndsnap_regular
+from .auxiliary import union_vals, ndsnap_regular, _p_norm
 
 __all__ = ["PersLandscapeApprox"]
 
@@ -364,7 +364,8 @@ def p_norm(self, p: int = 2) -> float:
         p: float, default 2
             value p of the L_{`p`} norm
         """
-        return np.sum([np.linalg.norm(depth, p) for depth in self.values])
+        super().p_norm(p=p)
+        return _p_norm(p=p, critical_pairs=self.values_to_pairs())
 
     def sup_norm(self) -> float:
         """

persim/landscapes/auxiliary.py

@@ -91,10 +91,7 @@ def union_crit_pairs(A, B):
         else:
             result_pairs.append(
                 slope_to_pos_interp(
-                    sum_slopes(
-                        pos_to_slope_interp(a),
-                        pos_to_slope_interp(b),
-                    )
+                    sum_slopes(pos_to_slope_interp(a), pos_to_slope_interp(b),)
                 )
             )
     return result_pairs
@@ -187,3 +184,33 @@ def ndsnap_regular(points, *grid_axes):
         best = np.argmin(np.abs(diff), axis=0)
         snapped.append(ax[best])
     return np.array(snapped).T
+
+
+def _p_norm(p: float, critical_pairs: list = []):
+    """
+    Compute `p` norm of interpolated piecewise linear function defined from list of 
+    critical pairs.
+    """
+    result = 0.0
+    for l in critical_pairs:
+        for [[x0, y0], [x1, y1]] in zip(l, l[1:]):
+            if y0 == y1:
+                # horizontal line segment
+                result += (np.abs(y0) ** p) * (x1 - x0)
+                continue
+            # slope is well-defined
+            slope = (y1 - y0) / (x1 - x0)
+            b = y0 - slope * x0
+            # segment crosses the x-axis
+            if (y0 < 0 and y1 > 0) or (y0 > 0 and y1 < 0):
+                z = -b / slope
+                ev_x1 = (slope * x1 + b) ** (p + 1) / (slope * (p + 1))
+                ev_x0 = (slope * x0 + b) ** (p + 1) / (slope * (p + 1))
+                ev_z = (slope * z + +b) ** (p + 1) / (slope * (p + 1))
+                result += np.abs(ev_x1 + ev_x0 - 2 * ev_z)
+            # segment does not cross the x-axis
+            else:
+                ev_x1 = (slope * x1 + b) ** (p + 1) / (slope * (p + 1))
+                ev_x0 = (slope * x0 + b) ** (p + 1) / (slope * (p + 1))
+                result += np.abs(ev_x1 - ev_x0)
+    return (result) ** (1.0 / p)

persim/landscapes/base.py

@@ -38,7 +38,11 @@ def __init__(self, dgms: list = [], hom_deg: int = 0) -> None:
 
     @abstractmethod
     def p_norm(self, p: int = 2) -> float:
-        pass
+        if p < -1 or -1 < p < 0:
+            raise ValueError(f"p can't be negative, but {p} was passed")
+        self.compute_landscape()
+        if p == -1:
+            return self.sup_norm()
 
     @abstractmethod
     def sup_norm(self) -> float:

persim/landscapes/exact.py

@@ -7,7 +7,7 @@
 from operator import itemgetter
 
 from .approximate import PersLandscapeApprox
-from .auxiliary import union_crit_pairs
+from .auxiliary import union_crit_pairs, _p_norm
 from .base import PersLandscape
 
 __all__ = ["PersLandscapeExact"]
@@ -402,34 +402,8 @@ def p_norm(self, p: int = 2) -> float:
         p: float, default 2
             value p of the L_{`p`} norm
         """
-        if p == -1:
-            return self.infinity_norm()
-        if p < -1 or -1 < p < 0:
-            raise ValueError(f"p can't be negative, but {p} was passed")
-        self.compute_landscape()
-        result = 0.0
-        for l in self.critical_pairs:
-            for [[x0, y0], [x1, y1]] in zip(l, l[1:]):
-                if y0 == y1:
-                    # horizontal line segment
-                    result += (np.abs(y0) ** p) * (x1 - x0)
-                    continue
-                # slope is well-defined
-                slope = (y1 - y0) / (x1 - x0)
-                b = y0 - slope * x0
-                # segment crosses the x-axis
-                if (y0 < 0 and y1 > 0) or (y0 > 0 and y1 < 0):
-                    z = -b / slope
-                    ev_x1 = (slope * x1 + b) ** (p + 1) / (slope * (p + 1))
-                    ev_x0 = (slope * x0 + b) ** (p + 1) / (slope * (p + 1))
-                    ev_z = (slope * z + +b) ** (p + 1) / (slope * (p + 1))
-                    result += np.abs(ev_x1 + ev_x0 - 2 * ev_z)
-                # segment does not cross the x-axis
-                else:
-                    ev_x1 = (slope * x1 + b) ** (p + 1) / (slope * (p + 1))
-                    ev_x0 = (slope * x0 + b) ** (p + 1) / (slope * (p + 1))
-                    result += np.abs(ev_x1 - ev_x0)
-        return (result) ** (1.0 / p)
+        super().p_norm(p=p)
+        return _p_norm(p=p, critical_pairs=self.critical_pairs)
 
     def sup_norm(self) -> float:
         """