Added a few new functions, fixed various bugs, and bumped to v0.4.0

Makefile

@@ -21,31 +21,24 @@ endif
 USERID=$(shell id -u)
 USERIDN=$(shell id -u -n)
 
-.PHONY: all build rebuild install uninstall run pull update test build-with-root-user
+.PHONY: all build install uninstall run test build-with-root-user
 
 all:
- @echo "Please specify an instruction (e.g make build)"
+ @ echo "Please specify an instruction (e.g make install)"
 
 build:
- @if [ "$(USERID)" = "0" ]; then \
+ @ if [ "$(USERID)" = "0" ]; then \
  echo "Please run make as a non-root user and without sudo!"; \
  false; \
  fi
+ @ echo "Deleting old CYTools image..."
+ sudo docker rmi cytools:uid-$(USERID) | echo "Old CYTools image does not exist or cannot be deleted"
  @ echo "Building CYTools image for user $(USERIDN)..."
- sudo docker build -t cytools:uid-$(USERID) --build-arg USERNAME=cytools\
+ sudo docker build --force-rm -t cytools:uid-$(USERID) --build-arg USERNAME=cytools\
  --build-arg USERID=$(USERID) --build-arg ARCH=$(arch) --build-arg AARCH=$(aarch)\
  --build-arg VIRTUAL_ENV=/home/cytools/cytools-venv/ --build-arg ALLOW_ROOT_ARG=" " .
  @ echo "Successfully built CYTools image for user $(USERIDN)"
 
-rebuild:
- @if [ "$(USERID)" = "0" ]; then \
- echo "Please run make as a non-root user and without sudo!"; \
- false; \
- fi
- sudo docker build --no-cache -t cytools:uid-$(USERID) --build-arg USERNAME=cytools\
- --build-arg USERID=$(USERID) --build-arg ARCH=$(arch) --build-arg AARCH=$(aarch)\
- --build-arg VIRTUAL_ENV=/home/cytools/cytools-venv/ --build-arg ALLOW_ROOT_ARG=" " .
-
 install: build
  @if [ "$(USERID)" = "0" ]; then \
  echo "Please run make as a non-root user and without sudo!"; \
@@ -97,11 +90,6 @@ run:
  bash scripts/linux/cytools; \
  fi
 
-pull:
- git pull
-
-update: pull install
-
 unittests:
  @if [ "$(USERID)" = "0" ]; then \
  echo "Please run make as a non-root user and without sudo!"; \
@@ -127,6 +115,9 @@ build-with-root-user:
  @ echo "********************************************************************"
  @ echo " "
  @ read -p "Press enter to continue or ctrl+c to cancel"
+ @ echo "Deleting old CYTools image..."
+ sudo docker rmi cytools:root | echo "Old CYTools image does not exist or cannot be deleted"
+ @ echo "Building CYTools image for root user..."
  sudo docker build -t cytools:root --build-arg USERNAME=root\
  --build-arg USERID=0 --build-arg ARCH=$(arch) --build-arg AARCH=$(aarch)\
  --build-arg VIRTUAL_ENV=/opt/cytools/cytools-venv/ --build-arg ALLOW_ROOT_ARG="--allow-root" .

cytools/__init__.py

@@ -19,7 +19,7 @@
 from cytools.utils import read_polytopes, fetch_polytopes
 
 # Latest version
-version = "0.3.5"
+version = "0.4.0"
 versions_with_serious_bugs = []
 
 # Check for more recent versions of CYTools
@@ -59,7 +59,8 @@ def check_for_updates():
  ver = tuple(int(c) for c in version.split("."))
  if latest_ver <= ver:
  continue
- print("\nInfo: A more recent version of CYTools is available. "
+ print("\nInfo: A more recent version of CYTools is available: "
+ f"{ver} -> {latest_ver}.\n"
  "We recommend upgrading before continuing.\n")
  elif not checked_bugs and "versions_with_serious_bugs =" in l:
  checked_bugs = True

cytools/calabiyau.py

@@ -1406,7 +1406,7 @@ def prime_toric_divisors(self):
  self.intersection_numbers()
  return self._prime_divs
 
- def second_chern_class(self, in_basis=False):
+ def second_chern_class(self, in_basis=False, include_origin=True):
  """
  **Description:**
  Computes the second Chern class of the CY hypersurface. Returns the
@@ -1419,6 +1419,8 @@ def second_chern_class(self, in_basis=False):
  **Arguments:**
  - `in_basis` *(bool, optional, default=False)*: Only return the
  integrals over a basis of divisors.
+ - `include_origin` *(bool, optional, default=True)*: Include the origin
+ in the vector, which corresponds to the canonical divisor.
 
  **Returns:**
  *(numpy.ndarray)* A vector containing the integrals.
@@ -1431,36 +1433,37 @@ def second_chern_class(self, in_basis=False):
  t = p.triangulate()
  cy = t.get_cy()
  cy.second_chern_class()
- # array([ 36, 306, 204, 36, 36, -6])
+ # array([-612, 36,  306,  204,  36,  36,  -6])
  ```
  """
  if self.dim() != 3:
  raise Exception("This function currently only supports 3-folds.")
  if self._second_chern_class is None:
- c2 = np.zeros(self.h11()+self.dim()+1, dtype=int)
+ c2 = np.zeros(len(self.prime_toric_divisors())+1, dtype=int)
  intnums = self.intersection_numbers(in_basis=False)
  for ii in intnums:
  if ii[0] == 0:
  continue
  if ii[0] == ii[1] == ii[2]:
  continue
  elif ii[0] == ii[1]:
- c2[ii[0]-1] += intnums[ii]
+ c2[ii[0]] += intnums[ii]
  elif ii[0] == ii[2]:
- c2[ii[0]-1] += intnums[ii]
+ c2[ii[0]] += intnums[ii]
  elif ii[1] == ii[2]:
- c2[ii[1]-1] += intnums[ii]
+ c2[ii[1]] += intnums[ii]
  else:
- c2[ii[0]-1] += intnums[ii]
- c2[ii[1]-1] += intnums[ii]
- c2[ii[2]-1] += intnums[ii]
+ c2[ii[0]] += intnums[ii]
+ c2[ii[1]] += intnums[ii]
+ c2[ii[2]] += intnums[ii]
+ c2[0] = -np.sum(c2)
  self._second_chern_class = c2
  if in_basis:
- basis = self.divisor_basis(include_origin=False)
+ basis = self.divisor_basis()
  if len(basis.shape) == 2: # If basis is matrix
  return self._second_chern_class.dot(basis.T)
  return np.array(self._second_chern_class[basis])
- return np.array(self._second_chern_class)
+ return np.array(self._second_chern_class)[(0 if include_origin else 1):]
 
  def is_smooth(self):
  """
@@ -1564,7 +1567,7 @@ def toric_kahler_cone(self):
 
  **Example:**
  We construct a Calabi-Yau hypersurface and find its Kähler cone.
- ```python {6}
+ ```python {4}
  p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
  t = p.triangulate()
  cy = t.get_cy()
@@ -1616,7 +1619,7 @@ def compute_cy_volume(self, tloc):
  **Example:**
  We construct a Calabi-Yau hypersurface and find its volume at the tip
  of the stretched Kähler cone.
- ```python {6}
+ ```python {5}
  p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
  t = p.triangulate()
  cy = t.get_cy()
@@ -1641,7 +1644,7 @@ def compute_cy_volume(self, tloc):
  xvol += intnums[ii]*np.prod([tloc[int(j)] for j in ii])/mult
  return xvol
 
- def compute_divisor_volumes(self, tloc):
+ def compute_divisor_volumes(self, tloc, in_basis=False):
  """
  **Description:**
  Computes the volume of the basis divisors at a location in the Kähler
@@ -1650,25 +1653,41 @@ def compute_divisor_volumes(self, tloc):
  **Arguments:**
  - `tloc` *(array_like)*: A vector specifying a location in the
  Kähler cone.
+ - `in_basis` *(bool, optional, default=False)*: When set to True, the
+ volumes of the current basis of divisors are computed. Otherwise, the
+ volumes of all prime toric divisors are computed.
 
  **Returns:**
- *(numpy.ndarray)* The list of volumes of the basis divisors at the
- specified location.
+ *(numpy.ndarray)* The list of volumes of the prime toric divisors or of
+ the basis divisors at the specified location.
 
  **Example:**
  We construct a Calabi-Yau hypersurface and find the volumes of the
- basis divisors at the tip of the stretched Kähler cone.
- ```python {6}
+ prime toric divisors at the tip of the stretched Kähler cone.
+ ```python {5}
  p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
  t = p.triangulate()
  cy = t.get_cy()
  tip = cy.toric_kahler_cone().tip_of_stretched_cone(1)
  cy.compute_divisor_volumes(tip)
- # array([2.5, 0.5])
+ # array([ 2.5 , 23.99999999, 16. , 2.5 , 2.5 ,
+ # 0.5 ])
  ```
  """
- intnums = self.intersection_numbers(in_basis=True,
- exact_arithmetic=False)
+ if not in_basis:
+ tloc_new = np.array(tloc).dot(self.divisor_basis(as_matrix=True, include_origin=False))
+ intnums = self.intersection_numbers(in_basis=False, exact_arithmetic=False)
+ tau = np.zeros(len(self.prime_toric_divisors()), dtype=float)
+ for ii in intnums:
+ if 0 in ii:
+ continue
+ c = Counter(ii)
+ for j in c.keys():
+ tau[j-1] += intnums[ii] * np.prod(
+ [tloc_new[k-1]**(c[k]-(j==k))/factorial(c[k]-(j==k))
+ for k in c.keys()])
+ return np.array(tau)
+ intnums = self.intersection_numbers(in_basis=True, exact_arithmetic=False)
  basis = self.divisor_basis()
  if len(basis.shape) == 2: # If basis is matrix
  tmp = np.array(intnums)
@@ -1685,6 +1704,43 @@ def compute_divisor_volumes(self, tloc):
  for k in c.keys()])
  return np.array(tau)
 
+ def compute_curve_volumes(self, tloc, only_extremal=False):
+ """
+ **Description:**
+ Computes the volume of the curves corresponding to (not necessarily
+ minimal) generators of the Mori cone inferred from toric geometry (i.e.
+ the cone obtained with the [`toric_mori_cone`](#toric_mori_cone)
+ function).
+
+ **Arguments:**
+ - `tloc` *(array_like)*: A vector specifying a location in the
+ Kähler cone.
+ - `only_extremal` *(bool, optional, default=False)*: Use only the
+ extremal rays of the Mori cone.
+
+ **Returns:**
+ *(numpy.ndarray)* The list of volumes of the curves.
+
+ **Example:**
+ We construct a Calabi-Yau hypersurface and find the volumes of the
+ generators of the Mori cone at the tip of the stretched Kähler cone.
+ ```python {5}
+ p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
+ t = p.triangulate()
+ cy = t.get_cy()
+ tip = cy.toric_kahler_cone().tip_of_stretched_cone(1)
+ cy.compute_curve_volumes(tip)
+ # array([0.99997511, 3.99992091, 0.99998193])
+ ```
+ As expected, all generators of the Mori cone have volumes greater than
+ or equal to 1 (up to rounding errors) at the tip of the stretched
+ Kähler cone.
+ """
+ c = self.toric_mori_cone(in_basis=True)
+ if only_extremal:
+ return c.extremal_rays().dot(tloc)
+ return c.rays().dot(tloc)
+
  def compute_AA(self, tloc):
  """
  **Description:**
@@ -1705,7 +1761,7 @@ def compute_AA(self, tloc):
  **Example:**
  We construct a Calabi-Yau hypersurface and compute this matrix at the
  tip of the stretched Kähler cone.
- ```python {6}
+ ```python {5}
  p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
  t = p.triangulate()
  cy = t.get_cy()
@@ -1761,7 +1817,7 @@ def compute_Kinv(self, tloc):
  **Example:**
  We construct a Calabi-Yau hypersurface and compute the inverse Kähler
  metric at the tip of the stretched Kähler cone.
- ```python {6}
+ ```python {5}
  p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
  t = p.triangulate()
  cy = t.get_cy()

cytools/polytope.py

@@ -1811,6 +1811,8 @@ def glsm_charge_matrix(self, include_origin=True,
  indices[:] = np.array([0] + list(range(1,linrel.shape[1]))[::-1])
  indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
  elif n_try > 3:
+ if n_try == 4:
+ np.random.seed(1337)
  np.random.shuffle(indices[1:])
  indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
  for ctr in range(np.prod(linrel.shape)+1):
@@ -1947,10 +1949,10 @@ def glsm_linear_relations(self, include_origin=True,
  ```
  """
  if points is not None:
- pts_ind = tuple(set(points))
+ pts_ind = tuple(set(list(points)+[0]))
  if min(pts_ind) < 0 or max(pts_ind) > self.points().shape[0]:
  raise Exception("An index is out of the allowed range.")
- include_origin = 0 in pts_ind
+ include_origin = 0 in points
  elif include_points_interior_to_facets:
  pts_ind = tuple(range(self.points().shape[0]))
  else:
@@ -2011,10 +2013,10 @@ def glsm_basis(self, include_origin=True,
  ```
  """
  if points is not None:
- pts_ind = tuple(set(points))
+ pts_ind = tuple(set(list(points)+[0]))
  if min(pts_ind) < 0 or max(pts_ind) > self.points().shape[0]:
  raise Exception("An index is out of the allowed range.")
- include_origin = 0 in pts_ind
+ include_origin = 0 in points
  elif include_points_interior_to_facets:
  pts_ind = tuple(range(self.points().shape[0]))
  else:
@@ -2486,6 +2488,65 @@ def automorphisms(self, square_to_one=False):
  self._autos = [np.array(autos), np.array(autos2)]
  return np.array(self._autos[args_id])
 
+ def find_2d_reflexive_subpolytopes(self):
+ """
+ **Description:**
+ Use the algorithm by Huang and Taylor described in
+ [1907.09482](https://arxiv.org/abs/1907.09482) to find 2D reflexive
+ subpolytopes in 4D polytopes.
+
+ **Arguments:**
+ None.
+
+ **Returns:**
+ *(list)* The list of 2D reflexive subpolytopes.
+
+ **Example:**
+ We construct a polytope and find its 2D reflexive subpolytopes.
+ ```python {2}
+ p = Polytope([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1],[-1,-1,-6,-9]])
+ p.find_2d_reflexive_subpolytopes()
+ # [A 2-dimensional lattice polytope in ZZ^4]
+ ```
+ """
+ if not self.is_reflexive() or self.dim() != 4:
+ raise Exception("Only 4D reflexive polytopes are supported.")
+ pts = self.points()
+ dual_vert = self.dual().vertices()
+ # Construct the sets S_i by finding the maximum dot product with dual vertices
+ S_i = [[]]*3
+ for p in pts:
+ m = max(p.dot(v) for v in dual_vert)
+ if m in (1,2,3):
+ S_i[m-1].append(tuple(p))
+ # Check each of the three conditions
+ gen_pts = []
+ for i in range(len(S_i[0])):
+ if tuple(-np.array(S_i[0][i])) in S_i[0]:
+ for j in range(i+1,len(S_i[0])):
+ if (tuple(-np.array(S_i[0][j])) in S_i[0]
+ and tuple(-np.array(S_i[0][i]))!=S_i[0][j]):
+ gen_pts.append((S_i[0][i],S_i[0][j]))
+ for i in range(len(S_i[1])):
+ for j in range(i+1,len(S_i[1])):
+ p = tuple(-np.array(S_i[1][i])-np.array(S_i[1][j]))
+ if p in S_i[0] or p in S_i[1]:
+ gen_pts.append((S_i[1][i],S_i[1][j]))
+ for i in range(len(S_i[2])):
+ for j in range(i+1,len(S_i[2])):
+ p = -np.array(S_i[2][i])-np.array(S_i[2][j])
+ if all(c%2 == 0 for c in p) and tuple(p//2) in S_i[0]:
+ gen_pts.append((S_i[2][i],S_i[2][j]))
+ polys_2d = set()
+ for p1,p2 in gen_pts:
+ pts_2d = set()
+ for p in pts:
+ if np.linalg.matrix_rank((p1,p2,p)) == 2:
+ pts_2d.add(tuple(p))
+ if np.linalg.matrix_rank(list(pts_2d)) == 2:
+ polys_2d.add(tuple(sorted(pts_2d)))
+ return [Polytope(pp) for pp in polys_2d]
+
  def minkowski_sum(self, other):
  """
  **Description:**