Merge pull request #34 from chazeon/master

Cleanup deprecated and removed NumPy APIs in C extension

c/correlation.c

@@ -109,26 +109,25 @@ static PyObject* correlation_par (PyObject* self, PyObject *arg, PyObject *keywo
     static char *kwlist[] = {"frequency", "velocity", "time_step", "step", "integration_method", NULL};
     if (!PyArg_ParseTupleAndKeywords(arg, keywords, "OOd|ii", kwlist, &frequency_obj, &velocity_obj, &TimeStep, &Increment, &IntMethod))  return NULL;
 
-    PyObject *velocity_array = PyArray_FROM_OTF(velocity_obj, NPY_CDOUBLE, NPY_IN_ARRAY);
-    PyObject *frequency_array = PyArray_FROM_OTF(frequency_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *velocity_array = PyArray_FROM_OTF(velocity_obj, NPY_CDOUBLE, NPY_ARRAY_IN_ARRAY);
+    PyObject *frequency_array = PyArray_FROM_OTF(frequency_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
 
     if (velocity_array == NULL || frequency_array == NULL ) {
         Py_XDECREF(velocity_array);
         Py_XDECREF(frequency_array);
         return NULL;
     }
 
-    _Dcomplex  *Velocity = (_Dcomplex *)PyArray_DATA(velocity_array);
-    double *Frequency    = (double*)PyArray_DATA(frequency_array);
-    int     NumberOfData = (int)PyArray_DIM(velocity_array, 0);
-    int     NumberOfFrequencies = (int)PyArray_DIM(frequency_array, 0);
+    _Dcomplex   *Velocity           = (_Dcomplex *)PyArray_DATA((PyArrayObject *)velocity_array);
+    double      *Frequency          = (double*)PyArray_DATA((PyArrayObject *)frequency_array);
+    npy_intp    NumberOfData        = PyArray_DIM((PyArrayObject *)velocity_array, 0);
+    npy_intp    NumberOfFrequencies = PyArray_DIM((PyArrayObject *)frequency_array, 0);
 
-
-    //Create new numpy array for storing result
+    // Create new numpy array for storing result
     PyArrayObject *PowerSpectrum_object;
-    int dims[1]={NumberOfFrequencies};
-    PowerSpectrum_object = (PyArrayObject *) PyArray_FromDims(1,dims,NPY_DOUBLE);
-    double *PowerSpectrum  = (double*)PyArray_DATA(PowerSpectrum_object);
+    npy_intp dims[] = {NumberOfFrequencies};
+    PowerSpectrum_object = (PyArrayObject *)PyArray_SimpleNew(1, dims, NPY_DOUBLE);
+    double *PowerSpectrum = (double *)PyArray_DATA(PowerSpectrum_object);
 
     // Maximum Entropy Method Algorithm
     if (IntMethod < 0 || IntMethod > 1) {

c/displacements.c

@@ -99,9 +99,9 @@ static PyObject *atomic_displacements(PyObject *self, PyObject *arg, PyObject *k
     static char *kwlist[] = {"trajectory", "positions", "cell", NULL};
     if (!PyArg_ParseTupleAndKeywords(arg, keywords, "OOO", kwlist,  &Trajectory_obj, &Positions_obj, &Cell_obj))  return NULL;
 
-    PyObject *Trajectory_array = PyArray_FROM_OTF(Trajectory_obj, NPY_CDOUBLE, NPY_IN_ARRAY);
-    PyObject *Positions_array = PyArray_FROM_OTF(Positions_obj, NPY_DOUBLE, NPY_IN_ARRAY);
-    PyObject *Cell_array = PyArray_FROM_OTF(Cell_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *Trajectory_array = PyArray_FROM_OTF(Trajectory_obj, NPY_CDOUBLE, NPY_ARRAY_IN_ARRAY);
+    PyObject *Positions_array = PyArray_FROM_OTF(Positions_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    PyObject *Cell_array = PyArray_FROM_OTF(Cell_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
 
     if (Cell_array == NULL || Trajectory_array == NULL) {
          Py_XDECREF(Cell_array);
@@ -112,22 +112,21 @@ static PyObject *atomic_displacements(PyObject *self, PyObject *arg, PyObject *k
 
 
 //  double _Complex *Cell           = (double _Complex*)PyArray_DATA(Cell_array);
-    _Dcomplex *Trajectory     = (_Dcomplex*)PyArray_DATA(Trajectory_array);
-    double *Positions               = (double*)PyArray_DATA(Positions_array);
+    _Dcomplex *Trajectory           = (_Dcomplex*)PyArray_DATA((PyArrayObject *)Trajectory_array);
+    double *Positions               = (double*)PyArray_DATA((PyArrayObject *)Positions_array);
 
-    int NumberOfData       = (int)PyArray_DIM(Trajectory_array, 0);
-    int NumberOfDimensions = (int)PyArray_DIM(Cell_array, 0);
+    npy_intp NumberOfData           = PyArray_DIM((PyArrayObject *)Trajectory_array, 0);
+    npy_intp NumberOfDimensions     = PyArray_DIM((PyArrayObject *)Cell_array, 0);
 
 //  Create new Numpy array to store the result
     _Dcomplex **Displacement;
     PyArrayObject *Displacement_object;
 
-    npy_intp dims[]={NumberOfData, NumberOfDimensions};
-    Displacement_object=(PyArrayObject *) PyArray_SimpleNew(2, dims, NPY_CDOUBLE);
-    Displacement=pymatrix_to_c_array_complex( Displacement_object);
+    npy_intp dims[] = {NumberOfData, NumberOfDimensions};
+    Displacement_object = (PyArrayObject *)PyArray_SimpleNew(2, dims, NPY_CDOUBLE);
+    Displacement = pymatrix_to_c_array_complex(Displacement_object);
 
-
-//  Create a pointer array from cell matrix (to be improved)
+    //  Create a pointer array from cell matrix (to be improved)
     double  **Cell_c = pymatrix_to_c_array_real((PyArrayObject *) Cell_array);
 
 /*
@@ -196,11 +195,12 @@ static PyObject *atomic_displacements(PyObject *self, PyObject *arg, PyObject *k
 
 static _Dcomplex **pymatrix_to_c_array_complex(PyArrayObject *array)  {
 
-      long n=(*array).dimensions[0];
-      long m=(*array).dimensions[1];
+      npy_intp n = PyArray_DIM(array, 0);
+      npy_intp m = PyArray_DIM(array, 1);
+
       _Dcomplex ** c = malloc(n*sizeof(_Dcomplex));
 
-      _Dcomplex *a = (_Dcomplex *) (*array).data;  /* pointer to array data as double _Complex */
+      _Dcomplex *a = (_Dcomplex *)PyArray_DATA(array);  /* pointer to array data as double _Complex */
       for ( int i=0; i<n; i++)  {
           c[i]=a+i*m;
       }
@@ -211,12 +211,12 @@ static _Dcomplex **pymatrix_to_c_array_complex(PyArrayObject *array)  {
 
 static double  **pymatrix_to_c_array_real(PyArrayObject *array)  {
 
-      long n=(*array).dimensions[0];
-      long m=(*array).dimensions[1];
+      npy_intp n = PyArray_DIM(array, 0);
+      npy_intp m = PyArray_DIM(array, 1);
       //PyObject *transpose_array = PyArray_Transpose(array, dims);
 
       double  ** c = malloc(n*sizeof(double));
-      double  *a = (double  *) (*array).data;
+      double *a = (double *)PyArray_DATA(array);
 
       for ( int i=0; i<n; i++)  {
           double  *b = malloc(m*sizeof(double));

c/mem.c

@@ -107,24 +107,24 @@ static PyObject* MaximumEntropyMethod (PyObject* self, PyObject *arg, PyObject *
     static char *kwlist[] = {"frequency", "velocity", "time_step", "coefficients", NULL};
     if (!PyArg_ParseTupleAndKeywords(arg, keywords, "OOd|i", kwlist, &frequency_obj, &velocity_obj, &TimeStep, &NumberOfCoefficients))  return NULL;
 
-    PyObject *velocity_array = PyArray_FROM_OTF(velocity_obj, NPY_CDOUBLE, NPY_IN_ARRAY);
-    PyObject *frequency_array = PyArray_FROM_OTF(frequency_obj, NPY_DOUBLE, NPY_IN_ARRAY);
+    PyObject *velocity_array = PyArray_FROM_OTF(velocity_obj, NPY_CDOUBLE, NPY_ARRAY_IN_ARRAY);
+    PyObject *frequency_array = PyArray_FROM_OTF(frequency_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
 
     if (velocity_array == NULL || frequency_array == NULL ) {
         Py_XDECREF(velocity_array);
         Py_XDECREF(frequency_array);
         return NULL;
     }
 
-    _Dcomplex  *Velocity = (_Dcomplex *)PyArray_DATA(velocity_array);
-    double *Frequency    = (double*)PyArray_DATA(frequency_array);
-    int    NumberOfData = (int)PyArray_DIM(velocity_array, 0);
-    int     NumberOfFrequencies = (int)PyArray_DIM(frequency_array, 0);
+    _Dcomplex   *Velocity               = (_Dcomplex *)PyArray_DATA((PyArrayObject *)velocity_array);
+    double      *Frequency              = (double*)PyArray_DATA((PyArrayObject *)frequency_array);
+    npy_intp     NumberOfData           = PyArray_DIM((PyArrayObject *)velocity_array, 0);
+    npy_intp     NumberOfFrequencies    = PyArray_DIM((PyArrayObject *)frequency_array, 0);
 
 
     //Create new numpy array for storing result
     PyArrayObject *PowerSpectrum_object;
-    npy_intp dims[]={NumberOfFrequencies};
+    npy_intp dims[] = {NumberOfFrequencies};
     PowerSpectrum_object = (PyArrayObject *) PyArray_SimpleNew(1,dims,NPY_DOUBLE);
 
 

dynaphopy/__init__.py

@@ -910,7 +910,7 @@ def plot_power_spectrum_full(self):
         # plt.legend()
         plt.show()
 
-        total_integral = integrate.simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
+        total_integral = integrate.simpson(self.get_power_spectrum_full(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def plot_power_spectrum_wave_vector(self):
@@ -920,7 +920,7 @@ def plot_power_spectrum_wave_vector(self):
         plt.ylabel('eV * ps')
         plt.axhline(y=0, color='k', ls='dashed')
         plt.show()
-        total_integral = integrate.simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
+        total_integral = integrate.simpson(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def plot_power_spectrum_phonon(self):
@@ -1031,14 +1031,14 @@ def write_power_spectrum_full(self, file_name):
         reading.write_curve_to_file(self.get_frequency_range(),
                                     self.get_power_spectrum_full()[None].T,
                                     file_name)
-        total_integral = integrate.simps(self.get_power_spectrum_full(), x=self.get_frequency_range())
+        total_integral = integrate.simpson(self.get_power_spectrum_full(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def write_power_spectrum_wave_vector(self, file_name):
         reading.write_curve_to_file(self.get_frequency_range(),
                                     self.get_power_spectrum_wave_vector()[None].T,
                                     file_name)
-        total_integral = integrate.simps(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
+        total_integral = integrate.simpson(self.get_power_spectrum_wave_vector(), x=self.get_frequency_range())
         print("Total Area (Kinetic energy <K>): {0} eV".format(total_integral))
 
     def write_power_spectrum_phonon(self, file_name):
@@ -1260,7 +1260,7 @@ def get_thermal_properties(self, force_constants=None):
         free_energy = thm.get_free_energy(temperature, phonopy_dos[0], phonopy_dos[1])
         entropy = thm.get_entropy(temperature, phonopy_dos[0], phonopy_dos[1])
         c_v = thm.get_cv(temperature, phonopy_dos[0], phonopy_dos[1])
-        integration = integrate.simps(phonopy_dos[1], x=phonopy_dos[0]) / (
+        integration = integrate.simpson(phonopy_dos[1], x=phonopy_dos[0]) / (
             self.dynamic.structure.get_number_of_atoms() *
             self.dynamic.structure.get_number_of_dimensions())
         total_energy = thm.get_total_energy(temperature, phonopy_dos[0], phonopy_dos[1])
@@ -1319,7 +1319,7 @@ def display_thermal_properties(self, from_power_spectrum=False, normalize_dos=Fa
 
             power_spectrum_dos = thm.get_dos(temperature, frequency_range, self.get_power_spectrum_full(),
                                              normalization)
-            integration = integrate.simps(power_spectrum_dos, x=frequency_range) / (
+            integration = integrate.simpson(power_spectrum_dos, x=frequency_range) / (
                 self.dynamic.structure.get_number_of_atoms() *
                 self.dynamic.structure.get_number_of_dimensions())
 

dynaphopy/analysis/fitting/__init__.py

@@ -1,6 +1,6 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.integrate import simps
+from scipy.integrate import simpson
 from dynaphopy.analysis.fitting import fitting_functions
 
 h_planck = 4.135667662e-3  # eV/ps
@@ -79,7 +79,7 @@ def phonon_fitting_analysis(original, ps_frequencies,
         maximum = fitting_parameters['maximum']
         error = fitting_parameters['global_error']
 
-        total_integral = simps(power_spectrum, x=ps_frequencies)
+        total_integral = simpson(power_spectrum, x=ps_frequencies)
 
         # Calculated properties
         dt_Q2_lor = 2 * area

dynaphopy/analysis/thermal_properties.py

@@ -37,7 +37,7 @@ def n(temp, freq):
     total_energy = np.nan_to_num([dos[i] * h_bar * freq * (0.5 + n(temperature, freq))
                                  for i, freq in enumerate(frequency)])
 
-    total_energy = integrate.simps(total_energy, frequency) * N_a / 1000  # KJ/K/mol
+    total_energy = integrate.simpson(total_energy, x=frequency) * N_a / 1000  # KJ/K/mol
     return total_energy
 
 
@@ -47,7 +47,7 @@ def get_free_energy(temperature, frequency, dos):
                                  for i, freq in enumerate(frequency)])
 
     free_energy[0] = 0
-    free_energy = integrate.simps(free_energy, frequency) * N_a / 1000  # KJ/K/mol
+    free_energy = integrate.simpson(free_energy, x=frequency) * N_a / 1000  # KJ/K/mol
     return free_energy
 
 
@@ -59,7 +59,7 @@ def n(temp, freq):
     free_energy_c = np.nan_to_num([dos[i] * -h_bar/2 *shift*(n(temperature, freq) + 1 / 2.)
                                    for i, freq in enumerate(frequency)])
 
-    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
+    free_energy_c = integrate.simpson(free_energy_c, x=frequency) * N_a / 1000 # KJ/K/mol
     return free_energy_c
 
 
@@ -76,7 +76,7 @@ def n(temp, freq):
 
     free_energy_c = free_energy_1 - free_energy_2
 
-    free_energy_c = integrate.simps(free_energy_c, frequency) * N_a / 1000 # KJ/K/mol
+    free_energy_c = integrate.simpson(free_energy_c, x=frequency) * N_a / 1000 # KJ/K/mol
     return free_energy_c
 
 
@@ -88,7 +88,7 @@ def coth(x):
     entropy = np.nan_to_num([dos[i]*(1.0 / (2. * temperature) * h_bar * freq * coth(h_bar * freq / (2 * k_b * temperature))
                                      - k_b * np.log(2 * np.sinh(h_bar * freq / (2 * k_b * temperature))))
                              for i, freq in enumerate(frequency)])
-    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
+    entropy = integrate.simpson(entropy, x=frequency) * N_a # J/K/mol
     return entropy
 
 # Alternative way to calculate entropy (not used)
@@ -100,7 +100,7 @@ def n(temp, freq):
     entropy = np.nan_to_num([dos[i] * k_b * ((n(temperature, freq) + 1) * np.log(n(temperature, freq) + 1)
                                              - n(temperature, freq) * np.log(n(temperature, freq)))
                          for i, freq in enumerate(frequency)])
-    entropy = integrate.simps(entropy, frequency) * N_a # J/K/mol
+    entropy = integrate.simpson(entropy, x=frequency) * N_a # J/K/mol
     return entropy
 
 
@@ -111,7 +111,7 @@ def z(temp, freq):
 
     c_v = np.nan_to_num([dos[i] * k_b * pow(z(temperature, freq), 2) * np.exp(z(temperature, freq)) / pow(np.exp(z(temperature, freq)) - 1, 2)
                          for i, freq in enumerate(frequency)])
-    c_v = integrate.simps(c_v, frequency) * N_a # J/K/mol
+    c_v = integrate.simpson(c_v, x=frequency) * N_a # J/K/mol
 
     return c_v
 
@@ -171,7 +171,7 @@ def z(temp, freq):
     print ('Entropy: {0} J/K/mol'.format(entropy))
     print ('Cv: {0} J/K/mol'.format(c_v))
     print (np.trapz(dos_r, x=frequency_r))/(8*3)
-    print (integrate.simps(dos_r,x=frequency_r)/(8*3))
+    print (integrate.simpson(dos_r,x=frequency_r)/(8*3))
 
     print ('\nFrom MD')
     print ('-------------------------')
@@ -183,7 +183,7 @@ def z(temp, freq):
     print ('Entropy: {0} J/K/mol'.format(entropy))
     print ('Cv: {0} J/K/mol'.format(c_v))
     print (np.trapz(power_spectrum, x=frequency_p))/(8*3)
-    print (integrate.simps(power_spectrum, x=frequency_p))/(8*3)
+    print (integrate.simpson(power_spectrum, x=frequency_p))/(8*3)
 
     print ('\nHARMONIC')
     print ('-------------------------')
@@ -195,4 +195,4 @@ def z(temp, freq):
     print ('Entropy: {0} J/K/mol'.format(entropy))
     print ('Cv: {0} J/K/mol'.format(c_v))
     print (np.trapz(dos, x=frequency)/(8*3))
-    print (integrate.simps(dos, x=frequency)/(8*3))
+    print (integrate.simpson(dos, x=frequency)/(8*3))

dynaphopy/atoms.py

@@ -187,7 +187,7 @@ def get_force_sets(self):
 
         if not isinstance(self._force_sets,type(None)):
             force_atoms_file = self._force_sets.get_dict()['natom']
-            force_atoms_input = np.product(np.diagonal(self._force_sets.get_supercell())) * self.get_number_of_atoms()
+            force_atoms_input = np.prod(np.diagonal(self._force_sets.get_supercell())) * self.get_number_of_atoms()
 
             if force_atoms_file != force_atoms_input:
                 print("Error: FORCE_SETS file does not match with SUPERCELL MATRIX")

dynaphopy/dynamics.py

@@ -108,7 +108,7 @@ def crop_trajectory(self, last_steps):
 
     def get_number_of_atoms(self):
         if self._number_of_atoms is None:
-            self._number_of_atoms = self.structure.get_number_of_atoms()*np.product(self.get_supercell_matrix())
+            self._number_of_atoms = self.structure.get_number_of_atoms()*np.prod(self.get_supercell_matrix())
         return self._number_of_atoms
 
     def set_time(self, time):

setup.py

@@ -77,7 +77,7 @@ def get_version_number():
                'scripts/fitdata',
                'scripts/qha_extract',
                'scripts/rfc_calc'],
-      install_requires=['phonopy', 'numpy', 'scipy', 'matplotlib'] + ["windows-curses"] if sys.platform in ["win32", "cygwin"] else [],
+      install_requires=['phonopy', 'numpy', 'scipy', 'matplotlib'] + (["windows-curses"] if sys.platform in ["win32", "cygwin"] else []),
       license='MIT License',
       long_description=open('README.md').read(),
       long_description_content_type='text/markdown',