Source code for burnman.classes.polytope

# This file is part of BurnMan - a thermoelastic and thermodynamic toolkit
# for the Earth and Planetary Sciences
# Copyright (C) 2012 - 2021 by the BurnMan team, released under the GNU
# GPL v2 or later.

from __future__ import absolute_import

import importlib
import numpy as np
from sympy import Rational
from fractions import Fraction
from scipy.spatial import Delaunay
from scipy.special import comb
from copy import copy

from .material import cached_property

from ..utils.math import independent_row_indices

    cdd = importlib.import_module('cdd')
except ImportError as err:
    print(f'Warning: {err}. '
          'For full functionality of BurnMan, please install pycddlib.')

class SimplexGrid(object):
    A class that creates objects that can efficiently generate a set of points
    that grid a simplex with a user-defined number of vertices. The class
    contains both a generator method and a grid method. It also contains
    an n_points attribute that returns the number of points in the gridded

    This class is available as :class:`burnman.polytope.SimplexGrid`.

    def __init__(self, vertices, points_per_edge):
        Initialize SimplexGrid object with the desired number of vertices
        and points per edge.
        assert vertices >= 2, 'need at least two vertices'
        assert points_per_edge >= 2, 'need at least 2 points per edge'

        self.vertices = vertices
        self.points_per_edge = points_per_edge

    def generate(self, generate_type='list'):
        Generates the grid points of the simplex in lexicographic order.

        generate_type : 'list' or 'array'
            Determines whether the generator returns lists or arrays
            corresponding to each point in the simplex grid.

        generator of lists or ndarrays (int, ndim=1)
            Grid points of the simplex.

        if generate_type == 'list':
            x = [0]*self.vertices
        elif generate_type == 'array':
            x = np.zeros(self.vertices, dtype=int)
            raise Exception('generate_type must be of type list or array.')

        x[self.vertices-1] = self.points_per_edge-1

        h = self.vertices
        while True:
            yield copy(x)

            h -= 1
            if h == 0:

            val = x[h]
            x[h] = 0
            x[self.vertices-1] = val - 1
            x[h-1] += 1
            if val != 1:
                h = self.vertices

    def grid(self, generate_type='list'):
        Returns either a list or a numpy array
        corresponding the the points in the simplex grid, depending on
        whether the user chooses 'list' (default) or 'array' as
        the generate_type parameter.
        if generate_type == 'list':
            return list(self.generate(generate_type))
        elif generate_type == 'array':
            return np.array(list(self.generate(generate_type)))
            raise Exception('generate_type must be of type list or array.')

    def n_points(self):
        The number of points corresponding to the number of vertices and
        points per edge chosen by the user.
        return comb(self.vertices+self.points_per_edge-2,
                    self.vertices-1, exact=True)

[docs]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. This class is available as :class:`burnman.polytope.MaterialPolytope`. """ def __init__(self, equalities, inequalities, number_type='fraction', return_fractions=False, independent_endmember_occupancies=None): """ Initialization function for the MaterialPolytope class. Declares basis attributes of the class. Parameters ---------- equalities: 2D numpy array A numpy array containing all the equalities of the polytope. Each row should evaluate to 0. inequalities: 2D numpy array A numpy array containing all the inequalities of the polytope. Each row should evaluate to <= 0. number_type: 'fraction' or 'float' (default is 'fraction') Whether pycddlib should read the input arrays as fractions or floats. return_fractions : boolean (default is False) Whether the generated polytope object should return fractions or floats. independent_endmember_occupancies : 2D numpy array (or None) If specified, this array provides the independent endmember set against which the dependent endmembers are defined. """ 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) self.polytope_matrix.rep_type = cdd.RepType.INEQUALITY self.polytope_matrix.extend(inequalities, linear=False) self.polytope = cdd.Polyhedron(self.polytope_matrix) if independent_endmember_occupancies is not None: self.independent_endmember_occupancies = independent_endmember_occupancies
[docs] def set_return_type(self, return_fractions=False): """ Sets the return_type for the polytope object. Also deletes the cached endmember_occupancies property. Parameters ---------- return_fractions : boolean (default is False) Whether the generated polytope object should return fractions or floats. """ try: del self.__dict__['endmember_occupancies'] except KeyError: pass self.return_fractions = return_fractions
[docs] @cached_property 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()[:]
[docs] @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)
[docs] @cached_property def n_endmembers(self): """ Return the number of endmembers (the number of vertices of the polytope). """ return len(self.raw_vertices)
[docs] @cached_property 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]) else: v = np.array([[float(value) for value in v] for v in 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]
[docs] @cached_property def independent_endmember_occupancies(self): """ Return an independent set of endmember occupancies (a linearly-independent set of vertex locations) """ arr = self.endmember_occupancies return arr[independent_row_indices(arr)]
[docs] @cached_property def endmembers_as_independent_endmember_amounts(self): """ Return a list of all the endmembers as a linear sum of the independent endmembers. """ ind = self.independent_endmember_occupancies sol = np.linalg.lstsq(np.array(ind.T).astype(float), np.array(self.endmember_occupancies.T).astype( float), rcond=0)[0].round(decimals=12).T return sol
def _decompose_vertices_into_simplices(self, vertices): """ Decomposes a set of vertices into simplices by Delaunay triangulation. """ # Delaunay triangulation only works in dimensions > 1 # and we remove the nullspace (sum(fractions) = 1) if len(vertices) > 2: nulls = np.repeat(vertices[:, -1], vertices.shape[1]).reshape(vertices.shape) tri = Delaunay((vertices - nulls)[:, :-1]) return tri.simplices else: return [[0, 1]]
[docs] @cached_property def independent_endmember_polytope(self): """ Returns the polytope expressed in terms of proportions of the independent endmembers. The polytope involves the first n-1 independent endmembers. The last endmember proportion makes the sum equal to one. """ 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)
[docs] @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)
[docs] 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)
[docs] 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)
[docs] def grid(self, points_per_edge=2, unique_sorted=True, grid_type='independent endmember proportions', limits=None): """ Create a grid of points which span the polytope. Parameters ---------- points_per_edge : integer (default is 2) Number of points per edge of the polytope. unique_sorted : boolean (default is True) 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. grid_type : 'independent endmember proportions' (default) or 'site occupancies' Whether to grid the polytope in terms of independent endmember proportions or site occupancies. limits : 2D numpy array Additional inequalities restricting the gridded area of the polytope. Returns ------- points : 2D numpy array A list of points gridding the polytope. """ if limits is None: if grid_type == 'independent endmember proportions': f_occ = (self.endmembers_as_independent_endmember_amounts / (points_per_edge - 1)) elif grid_type == 'site occupancies': f_occ = self.endmember_occupancies/(points_per_edge-1) else: raise Exception('grid type not recognised. Should be one of ' 'independent endmember proportions ' 'or site occupancies') simplices = self._decompose_vertices_into_simplices( self.endmembers_as_independent_endmember_amounts) else: if grid_type == 'independent endmember proportions': ppns = np.array(self.subpolytope_from_independent_endmember_limits( limits).get_generators()[:])[:, 1:] last_ppn = np.array([1. - sum(p) for p in ppns]).reshape((len(ppns), 1)) vertices_as_independent_endmember_proportions = np.hstack( (ppns, last_ppn)) f_occ = vertices_as_independent_endmember_proportions / \ (points_per_edge-1) elif grid_type == 'site occupancies': occ = np.array(self.subpolytope_from_site_occupancy_limits( limits).get_generators()[:])[:, 1:] f_occ = occ/(points_per_edge-1) ind = self.independent_endmember_occupancies vertices_as_independent_endmember_proportions = np.linalg.lstsq(np.array(ind.T).astype(float), np.array(occ.T).astype( float), rcond=None)[0].round(decimals=12).T else: raise Exception('grid_type not recognised. ' 'Should be one of ' 'independent endmember proportions ' 'or site occupancies') simplices = self._decompose_vertices_into_simplices( vertices_as_independent_endmember_proportions) n_ind = f_occ.shape[1] n_simplices = len(simplices) dim = len(simplices[0]) simplex_grid = SimplexGrid(dim, points_per_edge) grid = simplex_grid.grid('array') points_per_simplex = simplex_grid.n_points() n_points = n_simplices*points_per_simplex points = np.empty((n_points, n_ind)) idx = 0 for i in range(0, n_simplices): points[idx:idx+points_per_simplex] =[simplices[i]]) idx += points_per_simplex if unique_sorted: points = np.unique(points, axis=0) return points