Merge 25dc9f7 into 94a33b5

pycalphad/core/calculate.py

@@ -120,8 +120,9 @@ def _sample_phase_constitution(phase_name, phase_constituents, sublattice_dof, c
 
 
 def _compute_phase_values(components, statevar_dict,
-                          points, phase_record, output, maximum_internal_dof, broadcast=True, fake_points=False,
-                          largest_energy=None):
+                          points, phase_record, output, maximum_internal_dof,
+                          broadcast=True, fake_points=False,
+                          largest_energy=None, parallel=True,):
     """
     Calculate output values for a particular phase.
 
@@ -147,6 +148,8 @@ def _compute_phase_values(components, statevar_dict,
         If True, the first few points of the output surface will be fictitious
         points used to define an equilibrium hyperplane guaranteed to be above
         all the other points. This is used for convex hull computations.
+    parallel : bool, optional
+        If True, will calculate the objective function values in parallel. Defaults to True.
 
     Returns
     -------
@@ -183,7 +186,10 @@ def _compute_phase_values(components, statevar_dict,
     dof = np.ascontiguousarray(np.concatenate((bc_statevars.T, pts), axis=1))
     phase_output = np.zeros(dof.shape[0], order='C')
     phase_compositions = np.zeros((dof.shape[0], len(pure_elements)), order='F')
-    phase_record.obj(phase_output, dof)
+    if parallel:
+        phase_record.obj_parallel(phase_output, dof)
+    else:
+        phase_record.obj(phase_output, dof)
     for el_idx in range(len(pure_elements)):
         phase_record.mass_obj(phase_compositions[:,el_idx], dof, el_idx)
 
@@ -244,7 +250,8 @@ def _compute_phase_values(components, statevar_dict,
     return Dataset(data_arrays, coords=coordinate_dict)
 
 
-def calculate(dbf, comps, phases, mode=None, output='GM', fake_points=False, broadcast=True, parameters=None, **kwargs):
+def calculate(dbf, comps, phases, mode=None, output='GM', fake_points=False,
+              broadcast=True, parameters=None, parallel=True, **kwargs):
     """
     Sample the property surface of 'output' containing the specified
     components and phases. Model parameters are taken from 'dbf' and any
@@ -285,6 +292,8 @@ def calculate(dbf, comps, phases, mode=None, output='GM', fake_points=False, bro
         The density of points is determined by 'pdens'
     parameters : dict, optional
         Maps SymPy Symbol to numbers, for overriding the values of parameters in the Database.
+    parallel : bool, optional
+        If True, will calculate phase values in parallel. Defaults to True.
 
     Returns
     -------
@@ -410,7 +419,8 @@ def calculate(dbf, comps, phases, mode=None, output='GM', fake_points=False, bro
         phase_ds = _compute_phase_values(nonvacant_components, str_statevar_dict,
                                          points, phase_record, output,
                                          maximum_internal_dof, broadcast=broadcast,
-                                         largest_energy=float(largest_energy), fake_points=fp)
+                                         largest_energy=float(largest_energy),
+                                         fake_points=fp, parallel=parallel)
         all_phase_data.append(phase_ds)
 
     # speedup for single-phase case (found by profiling)

pycalphad/core/phase_rec.pxd

@@ -25,6 +25,7 @@ cdef public class PhaseRecord(object)[type PhaseRecordType, object PhaseRecordOb
     cdef public double[:,:,::1] composition_matrices
     cdef public int vacancy_index
     cpdef void obj(self, double[::1] out, double[:,::1] dof) nogil
+    cpdef void obj_parallel(self, double[::1] out, double[:,::1] dof) nogil
     cpdef void grad(self, double[::1] out, double[::1] dof) nogil
     cpdef void hess(self, double[:,::1] out, double[::1] dof) nogil
     cpdef void mass_obj(self, double[::1] out, double[:, ::1] dof, int comp_idx) nogil
@@ -33,4 +34,4 @@ cdef public class PhaseRecord(object)[type PhaseRecordType, object PhaseRecordOb
 
 cpdef PhaseRecord PhaseRecord_from_cython(object comps, object variables, double[::1] num_sites, double[::1] parameters,
                                           object ofunc, object gfunc, object hfunc,
-                                          object massfuncs, object massgradfuncs)
+                                          object massfuncs, object massgradfuncs)

pycalphad/core/phase_rec.pyx

@@ -1,4 +1,9 @@
+# distutils: extra_compile_args = -fopenmp
+# distutils: extra_link_args = -fopenmp
+
+cimport openmp
 cimport cython
+from cython.parallel cimport prange
 from cpython.mem cimport PyMem_Malloc, PyMem_Free
 import numpy as np
 cimport numpy as np
@@ -33,6 +38,34 @@ cdef public class PhaseRecord(object)[type PhaseRecordType, object PhaseRecordOb
         if self._obj != NULL:
             self._obj(&out[0], &dof[0,0], &self.parameters[0], <int>out.shape[0])
 
+    @cython.boundscheck(False)
+    @cython.wraparound(False)
+    cpdef void obj_parallel(self, double[::1] out, double[:,::1] dof) nogil:
+        """Wrapper around PhaseRecord.obj to compute the phase values in parallel"""
+        # we define the total number of chunks to be one greater than what we will parallelize
+        # this ensures that we can do a final pass and get any remaining rows left by integer division
+        # the whole reason we are doing this, rather than letting prange sort it out,
+        # is because prx.obj specifically expects that the inputs are vectorized
+        # e.g. out must be 1d and dof must be 2d.
+        cdef int nprocs = openmp.omp_get_num_procs()
+        cdef int chunks = nprocs + 1
+        cdef int chunksize = out.shape[0] // chunks
+        cdef int final_chunk_idx = (chunks - 1)*chunksize
+        cdef int chunk_idx
+        cdef int i
+        cdef int chunk_start_idx
+        cdef int chunk_end_idx
+
+        # TODO: num_threads should be num_procs, at least on a Mac, for optimal use.
+        # check for other systems
+        with nogil:
+            for i in prange(chunks-1, num_threads=nprocs):
+                chunk_start_idx = i*chunksize
+                chunk_end_idx = (i+1)*chunksize
+                self.obj(out[chunk_start_idx:chunk_end_idx], dof[chunk_start_idx:chunk_end_idx, :])
+            # go from the last parallel chunk through to the end
+            self.obj(out[final_chunk_idx:], dof[final_chunk_idx:, :])
+
     @cython.boundscheck(False)
     @cython.wraparound(False)
     cpdef void grad(self, double[::1] out, double[::1] dof) nogil:

pycalphad/tests/test_calculate.py

@@ -23,6 +23,12 @@ def test_surface():
     calculate(DBF, ['AL', 'CR', 'NI'], 'L12_FCC',
                 T=1273., mode='numpy')
 
+def test_parallel_serial_calculate():
+    "Serial and parallel objective functions should produce the same result."
+    serial = calculate(DBF, ['AL', 'CR', 'NI'], 'L12_FCC', T=1273., parallel=False)
+    parallel = calculate(DBF, ['AL', 'CR', 'NI'], 'L12_FCC', T=1273., parallel=True)
+    np.testing.assert_array_equal(serial.GM.values, parallel.GM.values)
+
 @nose.tools.raises(AttributeError)
 def test_unknown_model_attribute():
     "Sampling an unknown model attribute raises exception."

setup.py

@@ -27,7 +27,8 @@ def read(fname):
     ext_modules=cythonize(['pycalphad/core/hyperplane.pyx', 'pycalphad/core/eqsolver.pyx',
                            'pycalphad/core/phase_rec.pyx',
                            'pycalphad/core/composition_set.pyx',
-                           'pycalphad/core/problem.pyx']),
+                           'pycalphad/core/problem.pyx',
+                           ]),
     package_data={
         'pycalphad/core': ['*.pxd'],
     },