From f7d518ca120352f2e33ad627274de8e8666dd40e Mon Sep 17 00:00:00 2001
From: Kiran Vaddi <kiranvad@cpn-f07-09.cbls.ccr.buffalo.edu>
Date: Sat, 25 Jan 2020 16:58:08 -0500
Subject: [PATCH] add weighted delaunay triangulation

---
 cechmate/filtrations/alpha.py   | 64 +++++++++++++++++++++++++++++----
 docs/notebooks/BasicUsage.ipynb | 23 +++++++++++-
 2 files changed, 80 insertions(+), 7 deletions(-)

diff --git a/cechmate/filtrations/alpha.py b/cechmate/filtrations/alpha.py
index 8b33fc0..4cf95d9 100644
--- a/cechmate/filtrations/alpha.py
+++ b/cechmate/filtrations/alpha.py
@@ -4,6 +4,7 @@
 import numpy as np
 import numpy.linalg as linalg
 from scipy import spatial
+from scipy.spatial import ConvexHull
 
 from .base import BaseFiltration
 
@@ -26,7 +27,7 @@ class Alpha(BaseFiltration):
 
     """
 
-    def build(self, X):
+    def build(self, X, weights=[]):
         """
         Do the Alpha filtration of a Euclidean point set (requires scipy)
         
@@ -34,8 +35,13 @@ def build(self, X):
         ===========
         X: Nxd array
             Array of N Euclidean vectors in d dimensions
+        weights : Nx1 array
+            Array of weights (default is non-weighted)
+                  
         """
-
+        self.weights = weights
+        
+        
         if X.shape[0] < X.shape[1]:
             warnings.warn(
                 "The input point cloud has more columns than rows; "
@@ -47,14 +53,17 @@ def build(self, X):
 
         ## Step 1: Figure out the filtration
         if self.verbose:
-            print("Doing spatial.Delaunay triangulation...")
+            print("Doing Delaunay triangulation...")
             tic = time.time()
-
-        delaunay_faces = spatial.Delaunay(X).simplices
+        
+        if not np.shape(self.weights)[0]==0:
+            delaunay_faces = self._weighted_delaunay(X)
+        else:
+            delaunay_faces = spatial.Delaunay(X).simplices
         
         if self.verbose:
             print(
-                "Finished spatial.Delaunay triangulation (Elapsed Time %.3g)"
+                "Finished Delaunay triangulation (Elapsed Time %.3g)"
                 % (time.time() - tic)
             )
             print("Building alpha filtration...")
@@ -169,3 +178,46 @@ def _get_circumcenter(self, X):
                 x = x.dot(V.T) + muV
             return (x, rSqr)
         return (np.inf, np.inf)  # SC2 (Points not in general position)
+    
+    def _weighted_delaunay(self, X):
+        """
+        Computes a weighted delaunay triangulation using the duality of delaunay triangulation with a convex hull of a dimension higher
+        This function is triggered if a weight function is input to build the Alpha filtration 
+        and replaces the delaunay triangulation of scipy.
+        
+        Returns
+        -------
+        Triangular faces from a weighted delaunay triangulation in the form of a numpy array. 
+        Similar to the output from scipy.spatial.Delaunay
+
+        """
+        num, dim = np.shape(X)
+
+        lifted = np.zeros((num,dim+1))
+
+        for i in range(num):
+            p = X[i,:]
+            lifted[i,:] = np.append(p,np.sum(p**2) - self.weights[i]**2)
+        pinf = np.append(np.zeros((1,dim)),1e10);    
+        lifted = np.vstack((lifted, pinf))
+        hull = ConvexHull(lifted)
+        delaunay = []
+        for simplex in hull.simplices:
+            if num not in simplex:
+                delaunay.append(simplex.tolist())
+
+        return np.asarray(delaunay)
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        
+        
diff --git a/docs/notebooks/BasicUsage.ipynb b/docs/notebooks/BasicUsage.ipynb
index b1d9af4..ad46d7e 100644
--- a/docs/notebooks/BasicUsage.ipynb
+++ b/docs/notebooks/BasicUsage.ipynb
@@ -331,6 +331,27 @@
     "print(\"H1:\\n\", dgms[1])"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Weighted Alpha filtrations\n",
+    "He we showcase an usage of Alpha filtrations when the points have some weights. \n",
+    "This is done by computing a weighted delaunay triangulation and changing the circum radius accordingly for Alpha filtratiion computation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_points = 10\n",
+    "points = np.random.rand(num_points, 2)\n",
+    "weights = np.random.rand(num_points,1)\n",
+    "filtration = alpha.build(points, weights=weights)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -355,7 +376,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.1"
+   "version": "3.7.3"
   }
  },
  "nbformat": 4,