Fixed bug when finding a basis for a polytope whose points do not spa…

…n the lattice

cytools/__init__.py

@@ -19,7 +19,7 @@
 from cytools.utils import read_polytopes, fetch_polytopes
 
 # Latest version
-version = "0.3.3"
+version = "0.3.4"
 versions_with_serious_bugs = []
 
 # Check for more recent versions of CYTools

cytools/calabiyau.py

@@ -918,8 +918,7 @@ def glsm_charge_matrix(self, include_origin=True):
  """
  if self._glsm_charge_matrix is not None:
  return np.array(self._glsm_charge_matrix)[:,(0 if include_origin else 1):]
- toric_divs = [0]+list(self.prime_toric_divisors())
- pts = self.polytope().points_to_indices(self.polytope().points()[toric_divs])
+ pts = [0]+list(self.prime_toric_divisors())
  self._glsm_charge_matrix = self.polytope().glsm_charge_matrix(
  include_origin=True,
  points=pts)
@@ -956,8 +955,7 @@ def glsm_linear_relations(self, include_origin=True):
  """
  if self._glsm_linrels is not None:
  return np.array(self._glsm_linrels)[(0 if include_origin else 1):,(0 if include_origin else 1):]
- toric_divs = [0]+list(self.prime_toric_divisors())
- pts = self.polytope().points_to_indices(self.polytope().points()[toric_divs])
+ pts = [0]+list(self.prime_toric_divisors())
  self._glsm_linrels = self.polytope().glsm_linear_relations(
  include_origin=True,
  points=pts)
@@ -1004,12 +1002,11 @@ def divisor_basis(self, include_origin=True, as_matrix=False):
  ```
  """
  if self._divisor_basis is None:
+ pts = [0]+list(self.prime_toric_divisors())
  self.set_divisor_basis(self.polytope().glsm_basis(
  integral=True,
  include_origin=True,
- points=self.polytope().points_to_indices(
- self.polytope().points()[[0]+list(self.prime_toric_divisors())]))
- )
+ points=pts))
  if len(self._divisor_basis.shape) == 1:
  if 0 in self._divisor_basis and not include_origin:
  raise Exception("The basis was requested not including the "

cytools/polytope.py

@@ -1517,13 +1517,14 @@ def facets(self):
  """
  return self.faces(self._dim-1)
 
- def vertices(self):
+ def vertices(self, as_indices=False):
  """
  **Description:**
  Returns the vertices of the polytope.
 
  **Arguments:**
- None.
+ - `as_indices` *(bool)*: Return the points as indices of the full
+ list of points of the polytope.
 
  **Returns:**
  *(numpy.ndarray)* The list of vertices of the polytope.
@@ -1542,6 +1543,8 @@ def vertices(self):
  ```
  """
  if self._vertices is not None:
+ if as_indices:
+ return self.points_to_indices(self._vertices)
  return np.array(self._vertices)
  if self._dim == 0:
  self._vertices = np.array([self._input_pts[0]])
@@ -1607,6 +1610,8 @@ def vertices(self):
  if pt_tup not in input_pts:
  input_pts.append(pt_tup)
  self._vertices = np.array([list(pt) for pt in input_pts if pt in tmp_vert])
+ if as_indices:
+ return self.points_to_indices(self._vertices)
  return np.array(self._vertices)
 
  def dual(self):
@@ -1787,24 +1792,26 @@ def glsm_charge_matrix(self, include_origin=True,
  # If the result is not cached we do the computation
  # We start by finding a basis of columns
  linrel = self.points()[list(pts_ind)].T
- if int(round(np.linalg.det(np.array(fmpz_mat(linrel.tolist()).snf().tolist(), dtype=int)[:,:4]))) != 1:
- raise Exception("The points do not generate the lattice.")
+ sublat_ind = int(round(np.linalg.det(np.array(fmpz_mat(linrel.tolist()).snf().tolist(), dtype=int)[:,:linrel.shape[0]])))
  norms = [np.linalg.norm(p,1) for p in linrel.T]
  linrel = np.insert(linrel, 0, np.ones(linrel.shape[1], dtype=int), axis=0)
  good_exclusions = 0
  basis_exc = []
  indices = np.argsort(norms)
+ indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
  for n_try in range(4):
  if n_try == 1:
  indices[:] = np.array(range(linrel.shape[1]))
  elif n_try > 1:
- np.shuffle(indices[1:])
+ np.random.shuffle(indices[1:])
+ indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
  for ctr in range(np.prod(linrel.shape)+1):
  found_good_basis=True
  ctr += 1
  if ctr > 0:
  st = max([good_exclusions,1])
  indices[st:] = np.roll(indices[st:], -1)
+ indices[:linrel.shape[0]] = np.sort(indices[:linrel.shape[0]])
  linrel_rand = np.array(linrel[:,indices])
  try:
  linrel_hnf = fmpz_mat(linrel_rand.tolist()).hnf()
@@ -1813,12 +1820,14 @@ def glsm_charge_matrix(self, include_origin=True,
  linrel_rand = np.array(linrel_hnf.tolist(), dtype=int)
  good_exclusions = 0
  basis_exc = []
+ tmp_sublat_ind = 1
  for v in linrel_rand:
  for i,ii in enumerate(v):
  if ii != 0:
  if integral:
- if abs(ii) == 1:
- v *= ii
+ tmp_sublat_ind *= abs(ii)
+ if sublat_ind % tmp_sublat_ind == 0:
+ v *= ii//abs(ii)
  good_exclusions += 1
  else:
  found_good_basis = False
@@ -1829,7 +1838,7 @@ def glsm_charge_matrix(self, include_origin=True,
  if found_good_basis:
  break
  if found_good_basis:
-  break
+ break
  if not found_good_basis:
  raise Exception("failed")
  if ctr == np.prod(linrel.shape):
@@ -1845,16 +1854,16 @@ def glsm_charge_matrix(self, include_origin=True,
  points=points, integral=False)
  linrel_dict = {ii:i for i,ii in enumerate(indices)}
  linrel = np.array(linrel_rand[:,[linrel_dict[i] for i in range(linrel_rand.shape[1])]])
- basis_ind = np.array(sorted([i for i in range(linrel.shape[1]) if linrel_dict[i] not in basis_exc]), dtype=int)
- basis_exc = [i for i in range(linrel.shape[1]) if i not in basis_ind]
+ basis_ind = np.array([i for i in range(linrel.shape[1]) if linrel_dict[i] not in basis_exc], dtype=int)
+ basis_exc = np.array([indices[i] for i in basis_exc])
  glsm = np.zeros((linrel.shape[1]-linrel.shape[0],linrel.shape[1]), dtype=int)
- glsm[:,basis_ind] = np.eye(len(basis_ind),dtype=int)
- for nb in basis_exc:
+ glsm[:,basis_ind] = np.eye(len(basis_ind), dtype=int)
+ for nb in basis_exc[::-1]:
  tup = [(k,kk) for k,kk in enumerate(linrel[:,nb]) if kk]
- if len(tup) != 1 or abs(tup[0][1]) != 1:
+ if sublat_ind % tup[-1][1] != 0:
  raise Exception("Problem with linear relations")
- i,ii = tup[0]
- glsm[:,nb] = -ii*glsm.dot(linrel[i])
+ i,ii = tup[-1]
+ glsm[:,nb] = -glsm.dot(linrel[i])//ii
  # Check that everything was computed correctly
  if (np.linalg.matrix_rank(glsm[:,basis_ind]) != len(basis_ind)
  or any(glsm.dot(linrel.T).flat)

cytools/utils.py

@@ -598,12 +598,13 @@ def set_divisor_basis(tv_or_cy, basis, include_origin=True):
  linrels_new[:,b] = linrels_tmp[:,len(nobasis):]
  self._curve_basis_mat = np.zeros(glsm_cm.shape, dtype=int)
  self._curve_basis_mat[:,b] = np.eye(len(b),dtype=int)
- for nb in nobasis:
+ sublat_ind = int(round(np.linalg.det(np.array(fmpz_mat(linrels.tolist()).snf().tolist(), dtype=int)[:,:linrels.shape[0]])))
+ for nb in nobasis[::-1]:
  tup = [(k,kk) for k,kk in enumerate(linrels_new[:,nb]) if kk]
- if len(tup) != 1 or abs(tup[0][1]) != 1:
+ if sublat_ind % tup[-1][1] != 0:
  raise Exception("Problem with linear relations")
- i,ii = tup[0]
- self._curve_basis_mat[:,nb] = -ii*self._curve_basis_mat.dot(linrels_new[i])
+ i,ii = tup[-1]
+ self._curve_basis_mat[:,nb] = -self._curve_basis_mat.dot(linrels_new[i])//ii
  # Else if input is a matrix
  elif len(b.shape) == 2:
  if not config._exp_features_enabled:
@@ -638,11 +639,12 @@ def set_divisor_basis(tv_or_cy, basis, include_origin=True):
  points=self.polytope().points_to_indices(self.triangulation().points()))
  self._divisor_basis_mat = np.array(new_b)
  nobasis = np.array([i for i in range(glsm_cm.shape[1]) if i not in standard_basis])
- for nb in nobasis:
+ sublat_ind = int(round(np.linalg.det(np.array(fmpz_mat(linrels.tolist()).snf().tolist(), dtype=int)[:,:linrels.shape[0]])))
+ for nb in nobasis[::-1]:
  tup = [(k,kk) for k,kk in enumerate(linrels[:,nb]) if kk]
- if len(tup) != 1 or abs(tup[0][1]) != 1:
+ if sublat_ind % tup[-1][1] != 0:
  raise Exception("Problem with linear relations")
- i,ii = tup[0]
+ i,ii = tup[-1]
  for j in range(self._divisor_basis_mat.shape[0]):
  self._divisor_basis_mat[j] -= self._divisor_basis_mat[j,nb]*linrels[i]
  # Finally, we invert the matrix and construct the dual curve basis
@@ -655,12 +657,12 @@ def set_divisor_basis(tv_or_cy, basis, include_origin=True):
  inv_mat *= -1;
  self._curve_basis_mat = np.zeros(glsm_cm.shape, dtype=int)
  self._curve_basis_mat[:,standard_basis] = np.array(inv_mat).T
- for nb in nobasis:
+ for nb in nobasis[::-1]:
  tup = [(k,kk) for k,kk in enumerate(linrels[:,nb]) if kk]
- if len(tup) != 1 or abs(tup[0][1]) != 1:
+ if sublat_ind % tup[-1][1] != 0:
  raise Exception("Problem with linear relations")
- i,ii = tup[0]
- self._curve_basis_mat[:,nb] = -ii*self._curve_basis_mat.dot(linrels[i])
+ i,ii = tup[-1]
+ self._curve_basis_mat[:,nb] = -self._curve_basis_mat.dot(linrels[i])//ii
  self._curve_basis = np.array(self._curve_basis_mat)
  else:
  raise Exception("Input must be either a vector or a matrix.")

scripts/linux/cytools

@@ -106,7 +106,7 @@ cat << EOF
  ░░█████████ █████ █████ ░░██████ ░░██████ █████ ██████
  ░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░ ░░░░░░ ░░░░░ ░░░░░░
 
- Developed by Liam McAllister's Group | Version 0.3.3
+ Developed by Liam McAllister's Group | Version 0.3.4
  https://cytools.liammcallistergroup.com
 
 EOF

scripts/macos/cytools

@@ -100,7 +100,7 @@ cat << EOF
  ░░█████████ █████ █████ ░░██████ ░░██████ █████ ██████
  ░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░ ░░░░░░ ░░░░░ ░░░░░░
 
- Developed by Liam McAllister's Group | Version 0.3.3
+ Developed by Liam McAllister's Group | Version 0.3.4
  https://cytools.liammcallistergroup.com
 
 EOF

scripts/windows/launcher.ps1

@@ -38,7 +38,7 @@ $banner=@"
  ░░█████████ █████ █████ ░░██████ ░░██████ █████ ██████
  ░░░░░░░░░ ░░░░░ ░░░░░ ░░░░░░ ░░░░░░ ░░░░░ ░░░░░░
 
- Developed by Liam McAllister's Group | Version 0.3.3
+ Developed by Liam McAllister's Group | Version 0.3.4
  https://cytools.liammcallistergroup.com
 
 "@

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
  name="cytools",
- version="0.3.3",
+ version="0.3.4",
  author="Liam McAllister Group",
  author_email="",
  description="A software package for analyzing Calabi-Yau hypersurfaces in toric varieties.",