[Yt-svn] commit/yt: MatthewTurk: Remove this inv_bin_indices nonsense from BinnedProfile1D and use a hardcoded C
Bitbucket
commits-noreply at bitbucket.org
Fri Sep 2 21:13:37 PDT 2011
1 new changeset in yt:
http://bitbucket.org/yt_analysis/yt/changeset/22c102c7299f/
changeset: 22c102c7299f
branch: yt
user: MatthewTurk
date: 2011-09-03 06:12:39
summary: Remove this inv_bin_indices nonsense from BinnedProfile1D and use a hardcoded C
function. Why did it take so long to make this decision? Factor of ten
speedup for my profile run.
affected #: 2 files (4.2 KB)
--- a/yt/data_objects/profiles.py Fri Sep 02 10:20:42 2011 -0400
+++ b/yt/data_objects/profiles.py Sat Sep 03 00:12:39 2011 -0400
@@ -30,8 +30,8 @@
from yt.funcs import *
-from yt.utilities.data_point_utilities import Bin2DProfile, \
- Bin3DProfile
+from yt.utilities.data_point_utilities import \
+ Bin1DProfile, Bin2DProfile, Bin3DProfile
from yt.utilities.parallel_tools.parallel_analysis_interface import \
ParallelAnalysisInterface
@@ -239,22 +239,19 @@
mi, inv_bin_indices = args # Args has the indices to use as input
# check_cut is set if source != self._data_source
# (i.e., lazy_reader)
- source_data = self._get_field(source, field, check_cut)[mi]
- if weight: weight_data = self._get_field(source, weight, check_cut)[mi]
+ source_data = self._get_field(source, field, check_cut)
+ if weight: weight_data = self._get_field(source, weight, check_cut)
+ else: weight_data = na.ones(source_data.shape, dtype='float64')
+ self.total_stuff = source_data.sum()
binned_field = self._get_empty_field()
weight_field = self._get_empty_field()
used_field = self._get_empty_field()
- # Now we perform the actual binning
- for bin in inv_bin_indices.keys():
- # temp_field is *all* the points from source that go into this bin
- temp_field = source_data[inv_bin_indices[bin]]
- if weight:
- # now w_i * v_i and store sum(w_i)
- weight_field[bin] = weight_data[inv_bin_indices[bin]].sum()
- temp_field *= weight_data[inv_bin_indices[bin]]
- binned_field[bin] = temp_field.sum()
- # inv_bin_indices is a tuple of indices
- if inv_bin_indices[bin][0].size > 0: used_field[bin] = 1
+ mi = args[0]
+ bin_indices_x = args[1].ravel().astype('int64')
+ source_data = source_data[mi]
+ weight_data = weight_data[mi]
+ Bin1DProfile(bin_indices_x, weight_data, source_data,
+ weight_field, binned_field, used_field)
# Fix for laziness, because at the *end* we will be
# summing up all of the histograms and dividing by the
# weights. Accumulation likely doesn't work with weighted
@@ -270,26 +267,21 @@
raise EmptyProfileData()
# Truncate at boundaries.
if self.end_collect:
- mi = na.arange(source_data.size)
+ sd = source_data[:]
else:
- mi = na.where( (source_data > self._bins.min())
- & (source_data < self._bins.max()))
- sd = source_data[mi]
+ mi = ((source_data > self._bins.min())
+ & (source_data < self._bins.max()))
+ sd = source_data[mi]
if sd.size == 0:
raise EmptyProfileData()
# Stick the bins into our fixed bins, set at initialization
bin_indices = na.digitize(sd, self._bins)
if self.end_collect: #limit the range of values to 0 and n_bins-1
- bin_indices = na.minimum(na.maximum(1, bin_indices), self.n_bins) - 1
+ bin_indices = na.clip(bin_indices, 0, self.n_bins - 1)
else: #throw away outside values
bin_indices -= 1
- # Now we set up our inverse bin indices
- inv_bin_indices = {}
- for bin in range(self[self.bin_field].size):
- # Which fall into our bin?
- inv_bin_indices[bin] = na.where(bin_indices == bin)
- return (mi, inv_bin_indices)
+ return (mi, bin_indices)
def choose_bins(self, bin_style):
# Depending on the bin_style, choose from bin edges 0...N either:
--- a/yt/utilities/data_point_utilities.c Fri Sep 02 10:20:42 2011 -0400
+++ b/yt/utilities/data_point_utilities.c Sat Sep 03 00:12:39 2011 -0400
@@ -273,6 +273,111 @@
}
+static PyObject *_profile1DError;
+
+static PyObject *Py_Bin1DProfile(PyObject *obj, PyObject *args)
+{
+ int i, j;
+ PyObject *obins_x, *owsource, *obsource, *owresult, *obresult, *oused;
+ PyArrayObject *bins_x, *wsource, *bsource, *wresult, *bresult, *used;
+ bins_x = wsource = bsource = wresult = bresult = used = NULL;
+
+ if (!PyArg_ParseTuple(args, "OOOOOO",
+ &obins_x, &owsource, &obsource,
+ &owresult, &obresult, &oused))
+ return PyErr_Format(_profile1DError,
+ "Bin1DProfile: Invalid parameters.");
+ i = 0;
+
+ bins_x = (PyArrayObject *) PyArray_FromAny(obins_x,
+ PyArray_DescrFromType(NPY_INT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if(bins_x==NULL) {
+ PyErr_Format(_profile1DError,
+ "Bin1DProfile: One dimension required for bins_x.");
+ goto _fail;
+ }
+
+ wsource = (PyArrayObject *) PyArray_FromAny(owsource,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((wsource==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(wsource))) {
+ PyErr_Format(_profile1DError,
+ "Bin1DProfile: One dimension required for wsource, same size as bins_x.");
+ goto _fail;
+ }
+
+ bsource = (PyArrayObject *) PyArray_FromAny(obsource,
+ PyArray_DescrFromType(NPY_FLOAT64), 1, 1,
+ NPY_IN_ARRAY, NULL);
+ if((bsource==NULL) || (PyArray_SIZE(bins_x) != PyArray_SIZE(bsource))) {
+ PyErr_Format(_profile1DError,
+ "Bin1DProfile: One dimension required for bsource, same size as bins_x.");
+ goto _fail;
+ }
+
+ wresult = (PyArrayObject *) PyArray_FromAny(owresult,
+ PyArray_DescrFromType(NPY_FLOAT64), 1,1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if(wresult==NULL){
+ PyErr_Format(_profile1DError,
+ "Bin1DProfile: Two dimensions required for wresult.");
+ goto _fail;
+ }
+
+ bresult = (PyArrayObject *) PyArray_FromAny(obresult,
+ PyArray_DescrFromType(NPY_FLOAT64), 1,1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((bresult==NULL) ||(PyArray_SIZE(wresult) != PyArray_SIZE(bresult))
+ || (PyArray_DIM(bresult,0) != PyArray_DIM(wresult,0))){
+ PyErr_Format(_profile1DError,
+ "Bin1DProfile: Two dimensions required for bresult, same shape as wresult.");
+ goto _fail;
+ }
+
+ used = (PyArrayObject *) PyArray_FromAny(oused,
+ PyArray_DescrFromType(NPY_FLOAT64), 1,1,
+ NPY_INOUT_ARRAY | NPY_UPDATEIFCOPY, NULL);
+ if((used==NULL) ||(PyArray_SIZE(used) != PyArray_SIZE(wresult))
+ || (PyArray_DIM(used,0) != PyArray_DIM(wresult,0))){
+ PyErr_Format(_profile1DError,
+ "Bin1DProfile: Two dimensions required for used, same shape as wresult.");
+ goto _fail;
+ }
+
+ npy_float64 wval, bval;
+ int n;
+
+ for(n=0; n<bins_x->dimensions[0]; n++) {
+ i = *(npy_int64*)PyArray_GETPTR1(bins_x, n);
+ bval = *(npy_float64*)PyArray_GETPTR1(bsource, n);
+ wval = *(npy_float64*)PyArray_GETPTR1(wsource, n);
+ *(npy_float64*)PyArray_GETPTR1(wresult, i) += wval;
+ *(npy_float64*)PyArray_GETPTR1(bresult, i) += wval*bval;
+ *(npy_float64*)PyArray_GETPTR1(used, i) = 1.0;
+ }
+
+ Py_DECREF(bins_x);
+ Py_DECREF(wsource);
+ Py_DECREF(bsource);
+ Py_DECREF(wresult);
+ Py_DECREF(bresult);
+ Py_DECREF(used);
+
+ PyObject *onum_found = PyInt_FromLong((long)1);
+ return onum_found;
+
+ _fail:
+ Py_XDECREF(bins_x);
+ Py_XDECREF(wsource);
+ Py_XDECREF(bsource);
+ Py_XDECREF(wresult);
+ Py_XDECREF(bresult);
+ Py_XDECREF(used);
+ return NULL;
+
+}
+
static PyObject *_profile2DError;
static PyObject *Py_Bin2DProfile(PyObject *obj, PyObject *args)
@@ -1674,6 +1779,7 @@
{"Interpolate", Py_Interpolate, METH_VARARGS},
{"DataCubeRefine", Py_DataCubeRefine, METH_VARARGS},
{"DataCubeReplace", Py_DataCubeReplace, METH_VARARGS},
+ {"Bin1DProfile", Py_Bin1DProfile, METH_VARARGS},
{"Bin2DProfile", Py_Bin2DProfile, METH_VARARGS},
{"Bin3DProfile", Py_Bin3DProfile, METH_VARARGS},
{"FindContours", Py_FindContours, METH_VARARGS},
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list