Merge pull request #5 from arnaudon/weighted

correct weights

examples/weighted_complex/plot.py

@@ -0,0 +1,9 @@
+from simplicial_kuramoto.frustration_scan import plot_projections, plot_order
+import matplotlib.pyplot as plt
+
+if __name__ == "__main__":
+    for i in range(50):
+        path = f"results/square_{i}.pkl"
+        plot_projections(path, f"projections_{i:02d}.png")
+        plot_order(path, f"order_{i:02d}.png")
+        plt.close('all')

examples/weighted_complex/scan.py

@@ -0,0 +1,39 @@
+import numpy as np
+
+from simplicial_kuramoto import SimplicialComplex
+from simplicial_kuramoto.frustration_scan import scan_frustration_parameters
+from simplicial_kuramoto.graph_generator import delaunay_with_holes
+from simplicial_kuramoto.plotting import draw_simplicial_complex
+
+if __name__ == "__main__":
+    points = [[-1, -1], [-1, 1], [1, 1], [1, -1]]
+    w = 0.1
+    ws = np.linspace(0.1, 1.5, 50)
+    for i, w in enumerate(ws):
+        graph, points = delaunay_with_holes(points=points)
+        for _i, (u, v) in enumerate(graph.edges):
+            if _i == 0:
+                graph[u][v]["weight"] = w
+            else:
+                graph[u][v]["weight"] = 1.0
+
+        Gsc = SimplicialComplex(graph=graph, face_weights=[w, w])
+        draw_simplicial_complex(Gsc, "complex.pdf")
+
+        alpha1 = np.linspace(0, 2.5, 80)
+        alpha2 = np.linspace(0, np.pi / 2.0, 120)
+        n_repeats = 1
+        t_max = 1000
+        n_t = 500
+        n_workers = 80
+
+        scan_frustration_parameters(
+            Gsc,
+            filename=f"square_{i}.pkl",
+            alpha1=alpha1,
+            alpha2=alpha2,
+            repeats=n_repeats,
+            n_workers=n_workers,
+            t_max=t_max,
+            n_t=n_t,
+        )

simplicial_kuramoto/frustration_scan.py

@@ -412,7 +412,7 @@ def _get_scan_boundary(vec, axis=0):
         vec = np.diff(vec, axis=axis) > 0
         return a2[vec.T] - step1 / 2.0, a1[vec.T]
 
-    plt.figure()  # figsize=(5, 4))
+    plt.figure(figsize=(5, 4))
     plt.imshow(harm_order, origin="lower", extent=extent, aspect="auto")
     if with_proj:
         plt.plot(*_get_scan_boundary(grad), c="k", lw=1)

simplicial_kuramoto/simplicial_complex.py

@@ -110,30 +110,31 @@ def flip_edge_orientation(self, edge_indices):
     def W0(self):
         """Create node weight matrix."""
         if self._W0 is None:
-            if "weight" in self.graph.nodes[list(self.graph)[0]]:
-                node_weights = [self.graph.nodes[u]["weight"] for u in self.graph]
-            else:
-                node_weights = np.ones(self.n_nodes)
-            self._W0 = sc.sparse.spdiags(node_weights, 0, self.n_nodes, self.n_nodes)
+            self._W0 = sc.sparse.spdiags(
+                [self.graph.nodes[u].get("weight", 1.0) for u in self.graph],
+                0,
+                self.n_nodes,
+                self.n_nodes,
+            )
         return self._W0
 
     @property
     def W1(self):
         """Create edge weight matrix."""
         if self._W1 is None:
-            if "weight" in self.graph.nodes[list(self.graph)[0]]:
-                edge_weights = [self.graph[u][v]["weight"] for u, v in self.graph.edges]
-            else:
-                edge_weights = np.ones(self.n_edges)
-            self._W1 = sc.sparse.spdiags(edge_weights, 0, self.n_edges, self.n_edges)
+            self._W1 = sc.sparse.spdiags(
+                [self.graph[u][v].get("weight", 1.0) for u, v in self.graph.edges],
+                0,
+                self.n_edges,
+                self.n_edges,
+            )
         return self._W1
 
     @property
     def W2(self):
         """Create face weight matrix."""
-        if self._W2 is None:
-            if self.faces is not None:
-                self._W2 = sc.sparse.spdiags(self.face_weights, 0, self.n_faces, self.n_faces)
+        if self._W2 is None and self.faces is not None:
+            self._W2 = sc.sparse.spdiags(self.face_weights, 0, self.n_faces, self.n_faces)
         return self._W2
 
     @property
@@ -162,21 +163,20 @@ def N0s(self):
     @property
     def B1(self):
         """Create edge incidence matrix."""
-        if self._B1 is None:
-            if self.faces is not None:
-                self._B1 = sc.sparse.lil_matrix((self.n_faces, self.n_edges))
-                for face_index, face in enumerate(self.faces):
-                    for i in range(3):
-                        edge = tuple(np.roll(face, i)[:2])
-                        edge_rev = tuple(np.roll(face, i)[1::-1])
-                        if edge in self.edgelist:
-                            edge_index = self.edgelist.index(edge)
-                            self._B1[face_index, edge_index] = 1.0
-                        elif edge_rev in self.edgelist:
-                            edge_index = self.edgelist.index(edge_rev)
-                            self._B1[face_index, edge_index] = -1.0
-                        else:
-                            raise Exception("The face is not a triangle in the graph")
+        if self._B1 is None and self.faces is not None:
+            self._B1 = sc.sparse.lil_matrix((self.n_faces, self.n_edges))
+            for face_index, face in enumerate(self.faces):
+                for i in range(3):
+                    edge = tuple(np.roll(face, i)[:2])
+                    edge_rev = tuple(np.roll(face, i)[1::-1])
+                    if edge in self.edgelist:
+                        edge_index = self.edgelist.index(edge)
+                        self._B1[face_index, edge_index] = 1.0
+                    elif edge_rev in self.edgelist:
+                        edge_index = self.edgelist.index(edge_rev)
+                        self._B1[face_index, edge_index] = -1.0
+                    else:
+                        raise Exception("The face is not a triangle in the graph")
         return self._B1
 
     @property
@@ -221,7 +221,6 @@ def V1(self):
             )
         return self._V1
 
-
     @property
     def V2(self):
         """Lift operator on faces."""