Merge pull request #57 from unitaryfund/md/density_canvas

added triangular and hexagonal lattices

aquapointer/density_canvas/DensityCanvas.py

@@ -336,6 +336,14 @@ def set_rectangular_lattice(self, num_x, num_y, spacing, rotation = None):
         rotation = self.set_lattice_rotation(rotation)
         lattice = Lattice.rectangular(num_x=num_x, num_y=num_y, spacing=spacing, rotation=rotation)
         self.set_lattice(lattice, centering=True)
+
+    def set_triangular_lattice(self, nrows, ncols, spacing):
+        lattice = Lattice.triangular(nrows=nrows, ncols=ncols, spacing=spacing)
+        self.set_lattice(lattice, centering=True)
+
+    def set_hexagonal_lattice(self, nrows, ncols, spacing):
+        lattice = Lattice.hexagonal(nrows=nrows, ncols=ncols, spacing=spacing)
+        self.set_lattice(lattice, centering=True)
 
     def set_poisson_disk_lattice(self, spacing: tuple, rotation = None):
         rotation = self.set_lattice_rotation(rotation)
@@ -414,10 +422,9 @@ def lattice_dynamics(
         except AttributeError:
             pass
 
-    def calculate_pubo_coefficients(
-        self, p: int, params: ArrayLike, high=None, low=None
-    ):
-        """Calcuates the coefficients of the cost function.
+
+    def calculate_pubo_coefficients(self, p: int, params: ArrayLike, high=None, low=None, efficient_qubo=True):
+        """ Calcuates the coefficients of the cost function.
         The coefficients are stored in a dictionary {1:{}, 2:{}, ..., p:{}} where the key
         represents the interaction order and the values are dictionaries.
         The dictionary associated to a key is of the type {(0,1): val, (0,2): val, ...}
@@ -465,9 +472,9 @@ def _component(i: int, pms: ArrayLike) -> Self:
 
         # calculate using formula
         self._pubo = {
-            "coeffs": lpn.Lp_coefficients(
-                len(lattice._coords), p, _base, _component, params, high, low
-            ),
+
+            "coeffs": lpn.Lp_coefficients(lattice._coords, p, _base, _component, params, high, low, efficient_qubo),
+
             "p": p,
             "params": params,
             "high": high,

aquapointer/density_canvas/Lattice.py

@@ -69,6 +69,29 @@ def rectangular(cls, num_x: int, num_y: int, spacing: numbers.Number):
 
         return cls(np.array(coords, dtype=float), min_spacing=spacing, length_x=length_x, length_y=length_y, center=center, type="rectangular")
 
+    @classmethod
+    def triangular(cls, nrows: int, ncols: int, spacing: numbers.Number):
+        coords = np.mgrid[:ncols, :nrows].transpose().reshape(-1, 2).astype(float)
+        coords[:, 0] += 0.5 * np.mod(coords[:, 1], 2)
+        coords[:, 1] *= np.sqrt(3) / 2
+        coords *= spacing
+        center = np.array([np.mean(coords[:,0]), np.mean(coords[:,1])])
+
+        return cls(np.array(coords, dtype=float), min_spacing=spacing, center=center, type="triangular")
+
+    @classmethod
+    def hexagonal(cls, nrows: int, ncols: int, spacing: numbers.Number):
+        rows = range(nrows)
+        cols = range(ncols)
+        xx = (0.5 + i + i // 2 + (j % 2) * ((i % 2) - 0.5) for i in cols for j in rows)
+        h = np.sqrt(3) / 2
+        yy = (h * j for i in cols for j in rows)
+
+        coords = spacing*np.array([(x, y) for x, y in zip(xx, yy)])
+        center = np.array([np.mean(coords[:,0]), np.mean(coords[:,1])])
+
+        return cls(np.array(coords, dtype=float), min_spacing=spacing, center=center, type="hexagonal")
+
     @classmethod
     def poisson_disk(cls, density: ArrayLike, length: tuple, spacing: tuple, max_num: int = 8000):
         """

aquapointer/density_canvas/Lp_norm.py

@@ -161,11 +161,13 @@ def Lp_ising_energy_first_order_truncation(x, p, base, component, params):
 
     return res
 
-def Lp_coefficients(M, p, base, component, params, high=None, low=None):
+def Lp_coefficients(coords, p, base, component, params, high=None, low=None, efficient_qubo=True):
     """Calculates the Lp coefficients and stores them in a dictionary
     high and low are optional parameters that can truncate the computation at some
     order, from above or from below"""
 
+    M = len(coords)
+
     if high is None:
         high = p
     if low is None:
@@ -186,6 +188,7 @@ def Lp_coefficients(M, p, base, component, params, high=None, low=None):
         for idx in idxs:
             coefficients[i][idx] = 0
 
+
     # calculate coefficients
     for k in range(1, p+1):
         fact_k = np.math.factorial(k)
@@ -201,17 +204,26 @@ def Lp_coefficients(M, p, base, component, params, high=None, low=None):
             # cycle over all permutations
             index_permutations = it.permutations(range(M), r=len(r))
             for indices in index_permutations:
-                prefactor = (-1)**k * binom(p, k) * multiplicity 
-                value = integral(p, k, list(indices), r, base, component, params) 
+                prefactor = (-1)**k * binom(p, k) * multiplicity
+                if p==2 and k==2 and efficient_qubo:
+                    if len(r)==1:
+                        d = np.linalg.norm(coords[indices[0]] - coords[indices[0]])
+                    else:
+                        d = np.linalg.norm(coords[indices[0]] - coords[indices[1]])
+                    value = params[0]*params[0]*np.exp(-d**2/(2*2*params[1]))/(2*np.pi*2*params[1])
+                else:
+                    value = integral(p, k, list(indices), r, base, component, params)
                 coefficients[len(r)][tuple(sorted(indices))] += prefactor*value
     return coefficients
 
-def Lp_coefficients_first_order(M, p, base, component, params):
+def Lp_coefficients_first_order(coords, p, base, component, params):
     """Calculates the first order Lp coefficients"""
-    
+
     if p//2 != p/2:
         raise ValueError("power p must be even")
-
+
+    M = len(coords)
+
     # calculate first order
     coefficients = []
     for i1 in range(M):