Merge pull request #600 from bobmyhill/pycddlib3

enable pycddlib3

Author: bobmyhill

.github/workflows/main.yml

@@ -6,7 +6,7 @@ on: [push, pull_request, workflow_dispatch]
 jobs:  
   indent:
     name: indent
-    runs-on: [ubuntu-20.04]
+    runs-on: [ubuntu-24.04]
 
     strategy:
       matrix:
@@ -34,7 +34,7 @@ jobs:
 
   linux:
     name: test
-    runs-on: [ubuntu-20.04]
+    runs-on: [ubuntu-24.04]
 
     strategy:
       matrix:
@@ -48,14 +48,11 @@ jobs:
         python-version: ${{ matrix.python-versions }}
     - name: setup
       run: |
-        # pycddlib requires libgmp3-dev
         sudo apt update && sudo apt install --yes \
          numdiff \
-         libgmp3-dev \
          texlive \
          texlive-latex-extra
         python -m pip install --upgrade pip
-        pip install pycddlib==2.1.7 # this is optional
         pip install autograd # this is optional
         python --version
     - name: test

burnman/classes/polytope.py

@@ -12,18 +12,19 @@
 from scipy.spatial import Delaunay
 from scipy.special import comb
 from copy import copy
+import cdd as cdd_float
 
 from .material import cached_property
 
 from ..utils.math import independent_row_indices
 
+
 try:
-    cdd = importlib.import_module("cdd")
-except ImportError as err:
-    print(
-        f"Warning: {err}. "
-        "For full functionality of BurnMan, please install pycddlib."
-    )
+    cdd_fraction = importlib.import_module("cdd.gmp")
+    cdd_gmp_loaded = True
+except ImportError:
+    cdd_fraction = importlib.import_module("cdd")
+    cdd_gmp_loaded = False
 
 
 class SimplexGrid(object):
@@ -111,9 +112,7 @@ def n_points(self):
 
 class MaterialPolytope(object):
     """
-    A class that can be instantiated to create pycddlib polytope objects.
-    These objects can be interrogated to provide the vertices satisfying the
-    input constraints.
+    A class for creating and manipulating polytope objects using pycddlib.
 
     This class is available as :class:`burnman.polytope.MaterialPolytope`.
     """
@@ -122,7 +121,6 @@ def __init__(
         self,
         equalities,
         inequalities,
-        number_type="fraction",
         return_fractions=False,
         independent_endmember_occupancies=None,
     ):
@@ -132,31 +130,45 @@ def __init__(
 
         :param equalities: A numpy array containing all the
             equalities defining the polytope. Each row should evaluate to 0.
-        :type equalities: numpy.array (2D)
+        :type equalities: numpy.array (2D) of floats or Fractions
         :param inequalities: A numpy array containing all the inequalities
             defining the polytope. Each row should evaluate to <= 0.
-        :type inequalities: numpy.array (2D)
-        :param number_type: Whether pycddlib should read the input arrays as
-            fractions or floats. Valid options are 'fraction' or 'float'.
-        :type number_type: str
-        :param return_fractions: Choose whether the generated polytope object
-            should return fractions or floats.
+        :type inequalities: numpy.array (2D) of the same type as equalities
+        :param return_fractions: Whether or not to return fractions.
         :type return_fractions: bool
-        :param independent_endmember_occupancies: If specified, this array provides
-            the independent endmember set against which the dependent endmembers
-            are defined.
+        :param independent_endmember_occupancies: If specified, this array
+            provides the independent endmember set against which the
+            dependent endmembers are defined.
         :type independent_endmember_occupancies: numpy.array (2D) or None
         """
+        if equalities.dtype != inequalities.dtype:
+            raise Exception(
+                f"The equalities and inequalities arrays should have the same type ({equalities.dtype} != {inequalities.dtype})."
+            )
+
         self.set_return_type(return_fractions)
         self.equality_matrix = equalities[:, 1:]
         self.equality_vector = -equalities[:, 0]
 
-        self.polytope_matrix = cdd.Matrix(
-            equalities, linear=True, number_type=number_type
+        if equalities.dtype == Fraction and cdd_gmp_loaded is True:
+            self.number_type = Fraction
+            self.cdd = cdd_fraction
+        elif equalities.dtype == float or cdd_gmp_loaded is False:
+            self.number_type = float
+            self.cdd = cdd_float
+        else:
+            raise Exception(
+                "equalities should be arrays of either floats or Fractions."
+            )
+
+        self.polytope_matrix = self.cdd.matrix_from_array(equalities)
+        self.polytope_matrix.lin_set = set(range(len(equalities)))
+        self.polytope_matrix.rep_type = self.cdd.RepType.INEQUALITY
+
+        self.cdd.matrix_append_to(
+            self.polytope_matrix, self.cdd.matrix_from_array(inequalities)
         )
-        self.polytope_matrix.rep_type = cdd.RepType.INEQUALITY
-        self.polytope_matrix.extend(inequalities, linear=False)
-        self.polytope = cdd.Polyhedron(self.polytope_matrix)
+        self.polytope = self.cdd.polyhedron_from_matrix(self.polytope_matrix)
 
         if independent_endmember_occupancies is not None:
             self.independent_endmember_occupancies = independent_endmember_occupancies
@@ -182,14 +194,14 @@ def raw_vertices(self):
         Returns a list of the vertices of the polytope without any
         postprocessing. See also endmember_occupancies.
         """
-        return self.polytope.get_generators()[:]
+        return self.cdd.copy_generators(self.polytope).array
 
     @cached_property
     def limits(self):
         """
         Return the limits of the polytope (the set of bounding inequalities).
         """
-        return np.array(self.polytope.get_inequalities(), dtype=float)
+        return np.array(self.cdd.copy_inequalities(self.polytope).array, dtype=float)
 
     @cached_property
     def n_endmembers(self):
@@ -205,27 +217,23 @@ def endmember_occupancies(self):
         Return the endmember occupancies
         (a processed list of all of the vertex locations).
         """
-        if self.return_fractions:
-            if self.polytope.number_type == "fraction":
-                v = np.array(
-                    [[Fraction(value) for value in v] for v in self.raw_vertices]
-                )
-            else:
-                v = np.array(
-                    [
-                        [Rational(value).limit_denominator(1000000) for value in v]
-                        for v in self.raw_vertices
-                    ]
-                )
+        if self.number_type == float and self.return_fractions:
+            v = np.array(
+                [
+                    [Fraction(value).limit_denominator(1000000) for value in v]
+                    for v in self.raw_vertices
+                ]
+            )
+        elif self.number_type == Fraction and not self.return_fractions:
+            v = np.array(self.raw_vertices).astype(float)
         else:
-            v = np.array([[float(value) for value in v] for v in self.raw_vertices])
+            v = np.array(self.raw_vertices)
 
         if len(v.shape) == 1:
             raise ValueError(
                 "The combined equality and positivity "
                 "constraints result in a null polytope."
             )
-
         return v[:, 1:] / v[:, 0, np.newaxis]
 
     @cached_property
@@ -281,37 +289,41 @@ def independent_endmember_polytope(self):
         """
         arr = self.endmembers_as_independent_endmember_amounts
         arr = np.hstack((np.ones((len(arr), 1)), arr[:, :-1]))
-        M = cdd.Matrix(arr, number_type="fraction")
-        M.rep_type = cdd.RepType.GENERATOR
-        return cdd.Polyhedron(M)
+        arr = [[Fraction(value).limit_denominator(1000000) for value in v] for v in arr]
+        M = cdd_fraction.matrix_from_array(arr)
+        M.rep_type = cdd_fraction.RepType.GENERATOR
+        return cdd_fraction.polyhedron_from_matrix(M)
 
     @cached_property
     def independent_endmember_limits(self):
         """
         Gets the limits of the polytope as a function of the independent
         endmembers.
         """
-        return np.array(
-            self.independent_endmember_polytope.get_inequalities(), dtype=float
-        )
+        ind_poly = self.independent_endmember_polytope
+        inequalities = cdd_fraction.copy_inequalities(ind_poly).array
+        return np.array(inequalities, dtype=float)
 
     def subpolytope_from_independent_endmember_limits(self, limits):
         """
         Returns a smaller polytope by applying additional limits to the amounts
         of the independent endmembers.
         """
-        modified_limits = self.independent_endmember_polytope.get_inequalities().copy()
-        modified_limits.extend(limits, linear=False)
-        return cdd.Polyhedron(modified_limits)
+        ind_poly = self.independent_endmember_polytope
+        modified_limits = cdd_fraction.copy_inequalities(ind_poly)
+        cdd_fraction.matrix_append_to(
+            modified_limits, cdd_fraction.matrix_from_array(limits)
+        )
+        return cdd_fraction.polyhedron_from_matrix(modified_limits)
 
     def subpolytope_from_site_occupancy_limits(self, limits):
         """
         Returns a smaller polytope by applying additional limits to the
         individual site occupancies.
         """
-        modified_limits = self.polytope_matrix.copy()
-        modified_limits.extend(limits, linear=False)
-        return cdd.Polyhedron(modified_limits)
+        modified_limits = copy(self.polytope_matrix)
+        self.cdd.matrix_append_to(modified_limits, self.cdd.matrix_from_array(limits))
+        return self.cdd.polyhedron_from_matrix(modified_limits)
 
     def grid(
         self,
@@ -325,15 +337,16 @@ def grid(
 
         :param points_per_edge: Number of points per edge of the polytope.
         :type points_per_edge: int
-        :param unique_sorted: The gridding is done by splitting the polytope into
-            a set of simplices. This means that points will be duplicated along
-            vertices, faces etc. If unique_sorted is True, this function
+        :param unique_sorted: The gridding is done by splitting the polytope
+            into a set of simplices. This means that points will be duplicated
+            along vertices, faces etc. If unique_sorted is True, this function
             will sort and make the points unique. This is an expensive
             operation for large polytopes, and may not always be necessary.
         :type unique_sorted: bool
         :param grid_type: Whether to grid the polytope in terms of
             independent endmember proportions or site occupancies.
-            Choices are 'independent endmember proportions' or 'site occupancies'
+            Choices are 'independent endmember proportions' or
+            'site occupancies'
         :type grid_type: str
         :param limits: Additional inequalities restricting the
             gridded area of the polytope.
@@ -361,11 +374,8 @@ def grid(
             )
         else:
             if grid_type == "independent endmember proportions":
-                ppns = np.array(
-                    self.subpolytope_from_independent_endmember_limits(
-                        limits
-                    ).get_generators()[:]
-                )[:, 1:]
+                plims = self.subpolytope_from_site_occupancy_limits(limits)
+                ppns = np.array(self.cdd.copy_generators(plims).array)[:, 1:]
                 last_ppn = np.array([1.0 - sum(p) for p in ppns]).reshape(
                     (len(ppns), 1)
                 )
@@ -377,11 +387,8 @@ def grid(
                 )
 
             elif grid_type == "site occupancies":
-                occ = np.array(
-                    self.subpolytope_from_site_occupancy_limits(
-                        limits
-                    ).get_generators()[:]
-                )[:, 1:]
+                plims = self.subpolytope_from_site_occupancy_limits(limits)
+                occ = np.array(self.cdd.copy_generators(plims).array)[:, 1:]
                 f_occ = occ / (points_per_edge - 1)
 
                 ind = self.independent_endmember_occupancies

burnman/tools/chemistry.py

@@ -276,21 +276,28 @@ def reactions_from_stoichiometric_matrix(stoichiometric_matrix):
     """
     n_components = stoichiometric_matrix.shape[1]
 
-    equalities = np.concatenate(([np.zeros(n_components)], stoichiometric_matrix)).T
+    equalities = np.concatenate(
+        ([np.zeros(n_components)], np.array(stoichiometric_matrix).astype(float))
+    ).T
 
     polys = [
-        MaterialPolytope(equalities, np.diag(v))
+        MaterialPolytope(equalities, np.diag(v).astype(float))
         for v in itertools.product(*[[-1, 1]] * len(equalities[0]))
     ]
     reactions = []
     for p in polys:
-        v = np.array([[value for value in v] for v in p.raw_vertices])
+        v = np.array(
+            [
+                [Rational(value).limit_denominator(1000000) for value in v]
+                for v in p.raw_vertices
+            ],
+            dtype=float,
+        )
 
         if v is not []:
             reactions.extend(v)
 
-    reactions = np.unique(np.array(reactions, dtype=float), axis=0)
-
+    reactions = np.unique(np.array(reactions), axis=0)
     reactions = np.array(
         [[Rational(value).limit_denominator(1000000) for value in v] for v in reactions]
     )

burnman/tools/polytope.py

@@ -164,9 +164,7 @@ def composite_polytope_at_constrained_composition(
     eoccs = block_diag(*eoccs)
     inequalities = np.concatenate((np.zeros((len(eoccs), 1)), eoccs), axis=1)
 
-    return MaterialPolytope(
-        equalities, inequalities, number_type="float", return_fractions=return_fractions
-    )
+    return MaterialPolytope(equalities, inequalities, return_fractions=return_fractions)
 
 
 def simplify_composite_with_composition(composite, composition):

examples/example_polytopetools.py

@@ -35,6 +35,7 @@
 """
 from __future__ import absolute_import
 import numpy as np
+from sympy import Rational
 
 import burnman
 from burnman.minerals import SLB_2011
@@ -197,14 +198,21 @@
     )
     print([formula_to_string(f) for f in assemblage.endmember_formulae])
 
+    # The conversions from Fraction to Rational below
+    # are solely for pretty-printing purposes.
     print(
         "Which led to the following extreme vertices satifying "
         "the bulk compositional constraints:"
     )
-    print(old_polytope.endmember_occupancies)
+
+    print(
+        np.array([[Rational(v) for v in i] for i in old_polytope.endmember_occupancies])
+    )
 
     print("\nThe new assemblage has fewer endmembers, with compositions:")
     print([formula_to_string(f) for f in new_assemblage.endmember_formulae])
 
     print("The same vertices now look like this:")
-    print(new_polytope.endmember_occupancies)
+    print(
+        np.array([[Rational(v) for v in i] for i in new_polytope.endmember_occupancies])
+    )

pyproject.toml

@@ -25,6 +25,7 @@ sympy = "^1.12"
 cvxpy = "^1.3"
 matplotlib = "^3.7"
 numba = "^0.59"
+pycddlib-standalone = "^3.0"
 
 [tool.poetry.dev-dependencies]
 ipython = "^8.5"

test.sh

@@ -29,8 +29,8 @@ $PYTHON -m pipdeptree -p burnman -d 1 2> /dev/null
 echo ""
 
 # Quietly install optional modules after burnman
-echo "Installing optional cvxpy, pycddlib, autograd and jupyter modules ..."
-$PYTHON -m pip install -q cvxpy pycddlib==2.1.7 autograd jupyter
+echo "Installing optional autograd and jupyter modules ..."
+$PYTHON -m pip install -q autograd jupyter
 echo ""
 
 function testit {