FIX: Model/calculate: Piecewise branching performance and extrapolation

pycalphad/core/calculate.py

@@ -19,7 +19,7 @@
 from pycalphad.core.utils import endmember_matrix, extract_parameters, \
     get_pure_elements, filter_phases, instantiate_models, point_sample, \
     unpack_components, unpack_condition, unpack_kwarg
-from pycalphad.core.constants import MIN_SITE_FRACTION
+from pycalphad.core.constants import MIN_SITE_FRACTION, MAX_ENDMEMBER_PAIRS, MAX_EXTRA_POINTS
 
 
 def hr_point_sample(constraint_jac, constraint_rhs, initial_point, num_points):
@@ -150,12 +150,13 @@ def _sample_phase_constitution(model, sampler, fixed_grid, pdens):
     if (fixed_grid is True) and not linearly_constrained_space:
         # Sample along the edges of the endmembers
         # These constitution space edges are often the equilibrium points!
-        em_pairs = list(itertools.combinations(points, 2))
-        lingrid = np.linspace(0, 1, pdens)
-        extra_points = [first_em * lingrid[np.newaxis].T +
-                        second_em * lingrid[::-1][np.newaxis].T
-                        for first_em, second_em in em_pairs]
-        points = np.concatenate(list(itertools.chain([points], extra_points)))
+        em_pairs = list(itertools.combinations(points, 2))[:MAX_ENDMEMBER_PAIRS]
+        if len(em_pairs) > 0:
+            lingrid = np.linspace(0, 1, int(min(pdens, MAX_EXTRA_POINTS/len(em_pairs))))
+            extra_points = [first_em * lingrid[np.newaxis].T +
+                            second_em * lingrid[::-1][np.newaxis].T
+                            for first_em, second_em in em_pairs]
+            points = np.concatenate(list(itertools.chain([points], extra_points)))
     # Sample composition space for more points
     if sum(sublattice_dof) > len(sublattice_dof):
         if linearly_constrained_space:
@@ -477,6 +478,7 @@ def calculate(dbf, comps, phases, mode=None, output='GM', fake_points=False, bro
             raise ValueError(f"model must contain a Model instance for every active phase. Missing Model objects for {sorted(active_phases_without_models)}")
 
     maximum_internal_dof = max(len(models[phase_name].site_fractions) for phase_name in active_phases)
+
     for phase_name in sorted(active_phases):
         mod = models[phase_name]
         phase_record = phase_records[phase_name]

pycalphad/core/constants.py

@@ -12,3 +12,8 @@
 
 # Constraint scaling factors, for numerical stability
 INTERNAL_CONSTRAINT_SCALING = 1.0
+
+# Prevent excessive sampling for very complex phase models
+# This avoids running out of RAM
+MAX_ENDMEMBER_PAIRS = 5000 # ~100 endmembers
+MAX_EXTRA_POINTS = 90000

pycalphad/model.py

@@ -97,6 +97,12 @@ def __repr__(self):
             s = "ReferenceState('{}', '{}')".format(self.species.name, self.phase_name)
         return s
 
+class classproperty(object):
+    # https://stackoverflow.com/questions/5189699/how-to-make-a-class-property
+    def __init__(self, f):
+        self.f = f
+    def __get__(self, obj, owner):
+        return self.f(owner)
 
 class Model(object):
     """
@@ -154,6 +160,13 @@ class Model(object):
                      ('2st', 'twostate_energy'), ('ein', 'einstein_energy'),
                      ('vol', 'volume_energy'), ('ord', 'atomic_ordering_energy')]
 
+    # Behave as if piecewise temperature bounds extend to +-inf (i.e., ignore lower/upper T limits for parameters)
+    # This follows the behavior of most commercial codes, but may be undesirable in some circumstances
+    # Designed to be readonly on instances, because mutation after initialization will not work
+    @classproperty
+    def extrapolate_temperature_bounds(cls):
+        return True
+
     def __new__(cls, *args, **kwargs):
         target_cls = cls._dispatch_on(*args, **kwargs)
         instance = object.__new__(target_cls)
@@ -261,8 +274,47 @@ def __init__(self, dbe, comps, phase_name, parameters=None):
         self.site_fractions = sorted([x for x in self.variables if isinstance(x, v.SiteFraction)], key=str)
         self.state_variables = sorted([x for x in self.variables if not isinstance(x, v.SiteFraction)], key=str)
 
-    @staticmethod
-    def symbol_replace(obj, symbols):
+    @classmethod
+    def unwrap_piecewise(cls, graph):
+        from pycalphad.io.tdb import to_interval
+        replace_dict = {}
+        for atom in graph.atoms(Piecewise):
+            args = atom.args
+            # Unwrap temperature-dependent piecewise with zero-defaults
+            if len(args) == 4 and args[2] == 0 and args[3] == True and args[1].free_symbols == {v.T}:
+                replace_dict[atom] = args[0]
+            elif cls.extrapolate_temperature_bounds:
+                # Set lower and upper temperature limits to -+infinity
+                # First filter out default zero-branches
+                filtered_args = [(x, cond) for x, cond in zip(*[iter(args)]*2) if not ((cond == S.true) and (x == S.Zero))]
+                if len(filtered_args) == 0:
+                    continue
+                if not all([cond.free_symbols == {v.T} for _, cond in filtered_args]):
+                    # Only temperature-dependent piecewise conditions are supported for extrapolation
+                    continue
+                intervals = [to_interval(cond) for _, cond in filtered_args]
+                sortindices = [i[0] for i in sorted(enumerate(intervals), key=lambda x:x[1].args[0])]
+                if (intervals[sortindices[0]].args[0] == S.NegativeInfinity) and \
+                   (intervals[sortindices[-1]].args[1] == S.Infinity):
+                    # Nothing to do, temperature range already extrapolated
+                    continue
+                # First branch is special-cased to negative infinity
+                exprcondpairs = [(filtered_args[sortindices[0]][0], v.T < intervals[sortindices[0]].args[1])]
+                for idx in sortindices[1:-1]:
+                    exprcondpairs.append((filtered_args[sortindices[idx]][0],
+                                         And(v.T >= intervals[sortindices[idx]].args[0], v.T < intervals[sortindices[idx]].args[1])
+                    ))
+                # Last branch is special-cased to positive infinity
+                exprcondpairs.append((filtered_args[sortindices[-1]][0],
+                                      v.T >= intervals[sortindices[-1]].args[0]
+                ))
+                # Catch-all branch required for LLVM (should never hit in this formulation)
+                exprcondpairs.append((0, True))
+                replace_dict[atom] = Piecewise(*exprcondpairs)
+        return graph.xreplace(replace_dict)
+
+    @classmethod
+    def symbol_replace(cls, obj, symbols):
         """
         Substitute values of symbols into 'obj'.
 
@@ -280,6 +332,7 @@ def symbol_replace(obj, symbols):
             # of other symbols
             for iteration in range(_MAX_PARAM_NESTING):
                 obj = obj.xreplace(symbols)
+                obj = cls.unwrap_piecewise(obj)
                 undefs = [x for x in obj.free_symbols if not isinstance(x, v.StateVariable)]
                 if len(undefs) == 0:
                     break
@@ -1117,7 +1170,7 @@ def _partitioned_expr(disord_expr, ord_expr, disordered_mole_fraction_dict, orde
         return disord_expr + ordering_expr
 
     def atomic_ordering_energy(self, dbe):
-        """
+        r"""
         Return the atomic ordering contribution in symbolic form.
 
         If the current phase is anything other than the ordered phase in a