[Scipy-svn] r2936 - in trunk/Lib/sandbox/timeseries: lib src
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Apr 24 15:39:27 EDT 2007
Author: mattknox_ca
Date: 2007-04-24 14:39:22 -0500 (Tue, 24 Apr 2007)
New Revision: 2936
Modified:
trunk/Lib/sandbox/timeseries/lib/moving_funcs.py
trunk/Lib/sandbox/timeseries/src/cseries.c
Log:
added mov_corr and mov_covar functions. Also modified mov functions to work on 2D arrays
Modified: trunk/Lib/sandbox/timeseries/lib/moving_funcs.py
===================================================================
--- trunk/Lib/sandbox/timeseries/lib/moving_funcs.py 2007-04-19 19:57:23 UTC (rev 2935)
+++ trunk/Lib/sandbox/timeseries/lib/moving_funcs.py 2007-04-24 19:39:22 UTC (rev 2936)
@@ -26,7 +26,7 @@
__all__ = ['mov_sum', 'mov_median',
'mov_average', 'mov_mean', 'mov_average_expw',
- 'mov_stddev', 'mov_var', 'mov_sample_stddev', 'mov_sample_var',
+ 'mov_stddev', 'mov_var', 'mov_covar', 'mov_corr',
'cmov_average', 'cmov_mean', 'cmov_window'
]
@@ -40,160 +40,162 @@
data = orig_data.astype(rtype)
data[:] = result_dict['array']
- return marray(data, mask=rmask, copy=True, subok=True)
+ return marray(data, mask=rmask, copy=False, subok=True)
+def _moving_func(data, cfunc, kwargs):
-def mov_sum(data, window_size, dtype=None):
- kwargs = {'array':data,
- 'window_size':window_size}
+ if data.ndim == 1:
+ kwargs['array'] = data
- if dtype is not None:
- kwargs['dtype'] = dtype
-
- result_dict = MA_mov_sum(**kwargs)
- return _process_result_dict(data, result_dict)
+ result_dict = cfunc(**kwargs)
+ return _process_result_dict(data, result_dict)
-def mov_median(data, window_size, dtype=None):
- kwargs = {'array':data,
- 'window_size':window_size}
+ elif data.ndim == 2:
+ for i in range(data.shape[-1]):
+ kwargs['array'] = data[:,i]
+ result_dict = cfunc(**kwargs)
+
+ if i == 0:
+ rtype = result_dict['array'].dtype
+ result = data.astype(rtype)
+ print data.dtype, result.dtype
+
+ rmask = result_dict['mask']
- if dtype is not None:
- kwargs['dtype'] = dtype
-
- result_dict = MA_mov_median(**kwargs)
- return _process_result_dict(data, result_dict)
+ curr_col = marray(result_dict['array'], mask=rmask, copy=False)
+ result[:,i] = curr_col
-def mov_average(data, window_size, dtype=None):
- kwargs = {'array':data,
- 'window_size':window_size}
+ return result
+ else:
+ raise ValueError, "Data should be at most 2D"
+
+#...............................................................................
+def mov_sum(data, span, dtype=None):
+ """Calculates the moving sum of a series.
+
+:Parameters:
+ $$data$$
+ $$span$$
+ $$dtype$$"""
+
+ kwargs = {'span':span}
+ if dtype is None: dtype = data.dtype
+ kwargs['dtype'] = dtype
+
+
+ return _moving_func(data, MA_mov_sum, kwargs)
+#...............................................................................
+def mov_median(data, span, dtype=None):
+ """Calculates the moving median of a series.
+
+:Parameters:
+ $$data$$
+ $$span$$
+ $$dtype$$"""
+
+ kwargs = {'span':span}
+ if dtype is None: dtype = data.dtype
+ kwargs['dtype'] = dtype
+
+ return _moving_func(data, MA_mov_median, kwargs)
+#...............................................................................
+def mov_average(data, span, dtype=None):
+ """Calculates the moving average of a series.
+
+:Parameters:
+ $$data$$
+ $$span$$
+ $$dtype$$"""
+
+ kwargs = {'span':span}
if dtype is not None:
kwargs['dtype'] = dtype
-
- result_dict = MA_mov_average(**kwargs)
- return _process_result_dict(data, result_dict)
+
+ return _moving_func(data, MA_mov_average, kwargs)
mov_mean = mov_average
-
-def _mov_var_stddev(data, window_size, is_variance, is_sample, dtype):
+#...............................................................................
+def _mov_var_stddev(data, span, is_variance, bias, dtype):
"helper function for mov_var and mov_stddev functions"
- kwargs = {'array':data,
- 'window_size':window_size,
+ kwargs = {'span':span,
'is_variance':is_variance,
- 'is_sample':is_sample}
-
+ 'bias':bias}
if dtype is not None:
kwargs['dtype'] = dtype
-
- result_dict = MA_mov_stddev(**kwargs)
- return _process_result_dict(data, result_dict)
+ return _moving_func(data, MA_mov_stddev, kwargs)
+#...............................................................................
+def mov_var(data, span, bias=False, dtype=None):
+ """Calculates the moving variance of a 1-D array.
-def mov_var(data, window_size, dtype=None):
- """Calculates the moving variance of a 1-D array. This is the population
-variance. See "mov_sample_var" for moving sample variance.
-
:Parameters:
- data : ndarray
- Data as a valid (subclass of) ndarray or MaskedArray. In particular,
- TimeSeries objects are valid here.
- window_size : int
- Time periods to use for each calculation.
- dtype : numpy data type specification (*None*)
- Behaves the same as the dtype parameter for the numpy.var function.
-
-:Return value:
- The result is always a masked array (preserves subclass attributes). The
- result at index i uses values from [i-window_size:i+1], and will be masked
- for the first `window_size` values. The result will also be masked at i
- if any of the input values in the slice [i-window_size:i+1] are masked."""
-
+ $$data$$
+ $$span$$
+ $$bias$$
+ $$dtype$$"""
- return _mov_var_stddev(data=data, window_size=window_size,
- is_variance=1, is_sample=0, dtype=dtype)
+ return _mov_var_stddev(data=data, span=span,
+ is_variance=1, bias=int(bias), dtype=dtype)
+#...............................................................................
+def mov_stddev(data, span, bias=False, dtype=None):
+ """Calculates the moving standard deviation of a 1-D array.
-def mov_stddev(data, window_size, dtype=None):
- """Calculates the moving standard deviation of a 1-D array. This is the
-population standard deviation. See "mov_sample_stddev" for moving sample standard
-deviation.
-
:Parameters:
- data : ndarray
- Data as a valid (subclass of) ndarray or MaskedArray. In particular,
- TimeSeries objects are valid here.
- window_size : int
- Time periods to use for each calculation.
- dtype : numpy data type specification (*None*)
- Behaves the same as the dtype parameter for the numpy.std function.
-
-:Return value:
- The result is always a masked array (preserves subclass attributes). The
- result at index i uses values from [i-window_size:i+1], and will be masked
- for the first `window_size` values. The result will also be masked at i
- if any of the input values in the slice [i-window_size:i+1] are masked."""
+ $$data$$
+ $$span$$
+ $$bias$$
+ $$dtype$$"""
- return _mov_var_stddev(data=data, window_size=window_size,
- is_variance=0, is_sample=0, dtype=dtype)
+ return _mov_var_stddev(data=data, span=span,
+ is_variance=0, bias=int(bias), dtype=dtype)
+#...............................................................................
+def mov_covar(x, y, span, bias=False, dtype=None):
+ """Calculates the moving covariance of two 1-D arrays.
-
-def mov_sample_var(data, window_size, dtype=None):
- """Calculates the moving sample variance of a 1-D array.
-
:Parameters:
- data : ndarray
- Data as a valid (subclass of) ndarray or MaskedArray. In particular,
- TimeSeries objects are valid here.
- window_size : int
- Time periods to use for each calculation.
- dtype : numpy data type specification (*None*)
- Behaves the same as the dtype parameter for the numpy.var function.
-
-:Return value:
- The result is always a masked array (preserves subclass attributes). The
- result at index i uses values from [i-window_size:i+1], and will be masked
- for the first `window_size` values. The result will also be masked at i
- if any of the input values in the slice [i-window_size:i+1] are masked."""
-
+ $$x$$
+ $$y$$
+ $$span$$
+ $$bias$$
+ $$dtype$$"""
- return _mov_var_stddev(data=data, window_size=window_size,
- is_variance=1, is_sample=1, dtype=dtype)
+ result = x - mov_average(x, span, dtype=dtype)
+ result = result * (y - mov_average(y, span, dtype=dtype))
+
+ if bias: denom = span
+ else: denom = span - 1
+
+ return result/denom
+#...............................................................................
+def mov_corr(x, y, span, dtype=None):
+ """Calculates the moving correlation of two 1-D arrays.
-def mov_sample_stddev(data, window_size, dtype=None):
- """Calculates the moving sample standard deviation of a 1-D array.
-
:Parameters:
- data : ndarray
- Data as a valid (subclass of) ndarray or MaskedArray. In particular,
- TimeSeries objects are valid here.
- window_size : int
- Time periods to use for each calculation.
- dtype : numpy data type specification (*None*)
- Behaves the same as the dtype parameter for the numpy.std function.
-
-:Return value:
- The result is always a masked array (preserves subclass attributes). The
- result at index i uses values from [i-window_size:i+1], and will be masked
- for the first `window_size` values. The result will also be masked at i
- if any of the input values in the slice [i-window_size:i+1] are masked."""
-
- return _mov_var_stddev(data=data, window_size=window_size,
- is_variance=0, is_sample=1, dtype=dtype)
+ $$x$$
+ $$y$$
+ $$span$$
+ $$dtype$$"""
+ result = mov_covar(x, y, span, bias=True, dtype=dtype)
+ result = result / mov_stddev(x, span, bias=True, dtype=dtype)
+ result = result / mov_stddev(y, span, bias=True, dtype=dtype)
+
+ return result
+#...............................................................................
def mov_average_expw(data, span, tol=1e-6):
"""Calculates the exponentially weighted moving average of a series.
:Parameters:
- data : ndarray
- Data as a valid (subclass of) ndarray or MaskedArray. In particular,
- TimeSeries objects are valid here.
+ $$data$$
span : int
Time periods. The smoothing factor is 2/(span + 1)
tol : float, *[1e-6]*
Tolerance for the definition of the mask. When data contains masked
values, this parameter determinea what points in the result should be masked.
Values in the result that would not be "significantly" impacted (as
- determined by this parameter) by the masked values are left unmasked.
-"""
+ determined by this parameter) by the masked values are left unmasked."""
+
data = marray(data, copy=True, subok=True)
ismasked = (data._mask is not nomask)
data._mask = N.zeros(data.shape, bool_)
@@ -211,40 +213,16 @@
data._mask[0] = True
#
return data
-
-"""
-def weightmave(data, span):
- data = marray(data, subok=True, copy=True)
- data._mask = N.zeros(data.shape, bool_)
- # Set the data
- _data = data._data
- tmp = N.empty_like(_data)
- tmp[:span] = _data[:span]
- s = 0
- for i in range(span, len(data)):
- s += _data[i] - _data[i-span]
- tmp[i] = span*_data[i] + tmp[i-1] - s
- tmp *= 2./(span*(n+1))
- data._data.flat = tmp
- # Set the mask
- if data._mask is not nomask:
- msk = data._mask.nonzero()[0].repeat(span).reshape(-1,span)
- msk += range(span)
- data._mask[msk.ravel()] = True
- data._mask[:span] = True
- return data
-"""
-
#...............................................................................
def cmov_window(data, span, window_type):
- """Applies a centered moving window of type window_type and size span on the
- data.
+ """Applies a centered moving window of type window_type and size span on the
+data.
+
+Returns a (subclass of) MaskedArray. The k first and k last data are always
+masked (with k=span//2). When data has a missing value at position i, the
+result has missing values in the interval [i-k:i+k+1].
- Returns a (subclass of) MaskedArray. The k first and k last data are always
- masked (with k=span//2). When data has a missing value at position i,
- the result has missing values in the interval [i-k:i+k+1].
-
:Parameters:
data : ndarray
Data to process. The array should be at most 2D. On 2D arrays, the window
@@ -265,9 +243,8 @@
the needed parameters. If window_type is a floating point number, it is interpreted
as the beta parameter of the kaiser window.
-Note also that only boxcar has been thoroughly tested.
- """
- #
+Note also that only boxcar has been thoroughly tested."""
+
data = marray(data, copy=True, subok=True)
if data._mask is nomask:
data._mask = N.zeros(data.shape, bool_)
@@ -299,8 +276,55 @@
Data to process. The array should be at most 2D. On 2D arrays, the window
is applied recursively on each column.
span : integer
- The width of the window.
- """
+ The width of the window."""
return cmov_window(data, span, 'boxcar')
cmov_mean = cmov_average
+
+param_doc = {}
+param_doc['data'] = \
+"""data : ndarray
+ Data must be an ndarray (or subclass). In particular, note that
+ TimeSeries objects are valid here."""
+
+param_doc['x'] = \
+"""x : ndarray
+ First array to be included in the calculation. x must be an ndarray (or
+ subclass). In particular, note that TimeSeries objects are valid here."""
+
+param_doc['y'] = \
+"""y : ndarray
+ Second array to be included in the calculation. y must be an ndarray (or
+ subclass). In particular, note that TimeSeries objects are valid here."""
+
+param_doc['span'] = \
+"""span : int
+ Time periods to use for each calculation."""
+
+param_doc['bias'] = \
+"""bias : boolean (*False*)
+ If False, Normalization is by (N-1) where N == span (unbiased
+ estimate). If True then normalization is by N."""
+
+param_doc['dtype'] = \
+"""dtype : numpy data type specification (*None*)
+ dtype for the result"""
+
+mov_result_doc = \
+"""
+
+:Return value:
+ The result is always a masked array (preserves subclass attributes). The
+ result at index i uses values from [i-span:i+1], and will be masked for the
+ first `span` values. The result will also be masked at i if any of the
+ input values in the slice [i-span:i+1] are masked."""
+
+_g = globals()
+
+# generate function doc strings
+for fn in (x for x in __all__ if x[:4] == 'mov_' and x[4:] != 'mean'):
+ fdoc = _g[fn].func_doc
+ for prm, dc in param_doc.iteritems():
+ fdoc = fdoc.replace('$$'+prm+'$$', dc)
+ fdoc += mov_result_doc
+ _g[fn].func_doc = fdoc
Modified: trunk/Lib/sandbox/timeseries/src/cseries.c
===================================================================
--- trunk/Lib/sandbox/timeseries/src/cseries.c 2007-04-19 19:57:23 UTC (rev 2935)
+++ trunk/Lib/sandbox/timeseries/src/cseries.c 2007-04-24 19:39:22 UTC (rev 2936)
@@ -2887,7 +2887,81 @@
return returnVal;
}
+static PyObject *NP_ADD, *NP_MULTIPLY, *NP_SUBTRACT, *NP_SQRT,
+ *NP_GREATER, *NP_GREATER_EQUAL;
+static PyObject*
+np_add(PyObject *left_val, PyObject *right_val) {
+
+ PyObject *result;
+
+ result = PyObject_CallFunction(
+ NP_ADD, "OO",
+ (PyArrayObject*)left_val,
+ right_val);
+ return result;
+}
+
+static PyObject*
+np_subtract(PyObject *left_val, PyObject *right_val) {
+
+ PyObject *result;
+
+ result = PyObject_CallFunction(
+ NP_SUBTRACT, "OO",
+ (PyArrayObject*)left_val,
+ right_val);
+ return result;
+}
+
+static PyObject*
+np_multiply(PyObject *left_val, PyObject *right_val) {
+
+ PyObject *result;
+
+ result = PyObject_CallFunction(
+ NP_MULTIPLY, "OO",
+ (PyArrayObject*)left_val,
+ right_val);
+ return result;
+}
+
+static PyObject*
+np_sqrt(PyObject *val) {
+ return PyObject_CallFunction(NP_SQRT, "(O)", val);
+}
+
+static int np_greater(PyObject *left_val, PyObject *right_val) {
+
+ PyObject *temp;
+ int result;
+
+ temp = PyObject_CallFunction(
+ NP_GREATER, "OO",
+ (PyArrayObject*)left_val,
+ right_val);
+
+ result = (int)PyInt_AsLong(temp);
+ Py_DECREF(temp);
+ return result;
+}
+
+static int np_greater_equal(PyObject *left_val, PyObject *right_val) {
+
+ PyObject *temp;
+ int result;
+
+ temp = PyObject_CallFunction(
+ NP_GREATER_EQUAL, "OO",
+ (PyArrayObject*)left_val,
+ right_val);
+
+ result = (int)PyInt_AsLong(temp);
+ Py_DECREF(temp);
+ return result;
+}
+
+
/* This function is directly copied from direct copy of function in */
/* Return typenumber from dtype2 unless it is NULL, then return
NPY_DOUBLE if dtype1->type_num is integer or bool
@@ -2913,10 +2987,11 @@
/* validates the standard arguments to moving functions and set the original
mask, original ndarray, and mask for the result */
static PyObject *
-check_mov_args(PyObject *orig_arrayobj, int window_size, int min_win_size,
- PyArrayObject **orig_ndarray, PyArrayObject **result_mask) {
+check_mov_args(PyObject *orig_arrayobj, int span, int min_win_size,
+ PyObject **orig_ndarray, PyObject **result_mask) {
- PyArrayObject *orig_mask=NULL;
+ PyObject *orig_mask=NULL;
+ PyArrayObject **orig_ndarray_tmp, **result_mask_tmp;
int *raw_result_mask;
if (!PyArray_Check(orig_arrayobj)) {
@@ -2928,44 +3003,49 @@
if (PyObject_HasAttrString(orig_arrayobj, "_mask")) {
PyObject *tempMask = PyObject_GetAttrString(orig_arrayobj, "_mask");
if (PyArray_Check(tempMask)) {
- orig_mask = (PyArrayObject*)PyArray_EnsureArray(tempMask);
+ orig_mask = PyArray_EnsureArray(tempMask);
} else {
Py_DECREF(tempMask);
}
}
- *orig_ndarray = (PyArrayObject*)PyArray_EnsureArray(orig_arrayobj);
+ *orig_ndarray = PyArray_EnsureArray(orig_arrayobj);
+ orig_ndarray_tmp = (PyArrayObject**)orig_ndarray;
- if ((*orig_ndarray)->nd != 1) {
+ if ((*orig_ndarray_tmp)->nd != 1) {
PyErr_SetString(PyExc_ValueError, "array must be 1 dimensional");
return NULL;
}
- if (window_size < min_win_size) {
+ if (span < min_win_size) {
char *error_str;
error_str = malloc(60 * sizeof(char));
MEM_CHECK(error_str)
sprintf(error_str,
- "window_size must be greater than or equal to %i",
+ "span must be greater than or equal to %i",
min_win_size);
PyErr_SetString(PyExc_ValueError, error_str);
free(error_str);
return NULL;
}
- raw_result_mask = malloc((*orig_ndarray)->dimensions[0] * sizeof(int));
+ raw_result_mask = malloc((*orig_ndarray_tmp)->dimensions[0] * sizeof(int));
MEM_CHECK(raw_result_mask)
{
+ PyArrayObject *orig_mask_tmp;
int i, valid_points=0, is_masked;
- for (i=0; i<((*orig_ndarray)->dimensions[0]); i++) {
+ orig_mask_tmp = (PyArrayObject*)orig_mask;
+ for (i=0; i<((*orig_ndarray_tmp)->dimensions[0]); i++) {
+
is_masked=0;
if (orig_mask != NULL) {
PyObject *valMask;
- valMask = PyArray_GETITEM(orig_mask, PyArray_GetPtr(orig_mask, &i));
+ valMask = PyArray_GETITEM(orig_mask_tmp,
+ PyArray_GetPtr(orig_mask_tmp, &i));
is_masked = (int)PyInt_AsLong(valMask);
Py_DECREF(valMask);
}
@@ -2973,29 +3053,26 @@
if (is_masked) {
valid_points=0;
} else {
- if (valid_points < window_size) { valid_points += 1; }
- if (valid_points < window_size) { is_masked = 1; }
+ if (valid_points < span) { valid_points += 1; }
+ if (valid_points < span) { is_masked = 1; }
}
raw_result_mask[i] = is_masked;
}
}
- *result_mask = (PyArrayObject*)PyArray_SimpleNewFromData(
- 1, (*orig_ndarray)->dimensions,
- PyArray_INT32, raw_result_mask);
+ *result_mask = PyArray_SimpleNewFromData(
+ 1, (*orig_ndarray_tmp)->dimensions,
+ PyArray_INT32, raw_result_mask);
MEM_CHECK(*result_mask)
- (*result_mask)->flags = ((*result_mask)->flags) | NPY_OWNDATA;
+ result_mask_tmp = (PyArrayObject**)result_mask;
+ (*result_mask_tmp)->flags = ((*result_mask_tmp)->flags) | NPY_OWNDATA;
}
-
-static PyObject *NP_ADD, *NP_MULTIPLY, *NP_SUBTRACT, *NP_SQRT,
- *NP_GREATER, *NP_GREATER_EQUAL;
-
/* computation portion of moving sum. Appropriate mask is overlayed on top
afterwards */
static PyObject*
-calc_mov_sum(PyArrayObject *orig_ndarray, int window_size, int rtype)
+calc_mov_sum(PyArrayObject *orig_ndarray, int span, int rtype)
{
PyArrayObject *result_ndarray=NULL;
int i;
@@ -3019,22 +3096,18 @@
PyObject *mov_sum_prevval;
mov_sum_prevval= PyArray_GETITEM(result_ndarray,
PyArray_GetPtr(result_ndarray, &prev_idx));
- mov_sum_val = PyObject_CallFunction(NP_ADD, "OO", (PyArrayObject*)val,
- mov_sum_prevval);
+ mov_sum_val = np_add(val, mov_sum_prevval);
Py_DECREF(mov_sum_prevval);
ERR_CHECK(mov_sum_val)
- if (i >= (window_size - 1)) {
+ if (i >= span) {
PyObject *temp_val, *rem_val;
- int rem_idx = i-window_size;
+ int rem_idx = i-span;
temp_val = mov_sum_val;
rem_val = PyArray_GETITEM(orig_ndarray,
PyArray_GetPtr(orig_ndarray, &rem_idx));
- mov_sum_val = PyObject_CallFunction(
- NP_SUBTRACT,
- "OO", (PyArrayObject*)temp_val,
- rem_val);
+ mov_sum_val = np_subtract(temp_val, rem_val);
ERR_CHECK(mov_sum_val)
Py_DECREF(temp_val);
@@ -3042,7 +3115,9 @@
}
}
- PyArray_SETITEM(result_ndarray, PyArray_GetPtr(result_ndarray, &i), mov_sum_val);
+ PyArray_SETITEM(result_ndarray,
+ PyArray_GetPtr(result_ndarray, &i),
+ mov_sum_val);
if (mov_sum_val != val) { Py_DECREF(val); }
@@ -3057,31 +3132,33 @@
static PyObject *
MaskedArray_mov_sum(PyObject *self, PyObject *args, PyObject *kwds)
{
- PyObject *orig_arrayobj=NULL, *result_dict=NULL;
- PyArrayObject *orig_ndarray=NULL, *result_ndarray=NULL, *result_mask=NULL;
+ PyObject *orig_arrayobj=NULL, *orig_ndarray=NULL,
+ *result_ndarray=NULL, *result_mask=NULL,
+ *result_dict=NULL;
PyArray_Descr *dtype=NULL;
- int rtype, window_size;
+ int rtype, span;
- static char *kwlist[] = {"array", "window_size", "dtype", NULL};
+ static char *kwlist[] = {"array", "span", "dtype", NULL};
if (!PyArg_ParseTupleAndKeywords(args, kwds,
- "Oi|O&:mov_sum(array, window_size, dtype)", kwlist,
- &orig_arrayobj, &window_size,
+ "Oi|O&:mov_sum(array, span, dtype)", kwlist,
+ &orig_arrayobj, &span,
PyArray_DescrConverter2, &dtype)) return NULL;
- check_mov_args(orig_arrayobj, window_size, 1,
+ check_mov_args(orig_arrayobj, span, 1,
&orig_ndarray, &result_mask);
rtype = _CHKTYPENUM(dtype);
- result_ndarray = (PyArrayObject*)calc_mov_sum(orig_ndarray, window_size, rtype);
+ result_ndarray = calc_mov_sum((PyArrayObject*)orig_ndarray,
+ span, rtype);
ERR_CHECK(result_ndarray)
result_dict = PyDict_New();
MEM_CHECK(result_dict)
- PyDict_SetItemString(result_dict, "array", (PyObject*)result_ndarray);
- PyDict_SetItemString(result_dict, "mask", (PyObject*)result_mask);
+ PyDict_SetItemString(result_dict, "array", result_ndarray);
+ PyDict_SetItemString(result_dict, "mask", result_mask);
Py_DECREF(result_ndarray);
Py_DECREF(result_mask);
@@ -3092,117 +3169,49 @@
static PyObject *
MaskedArray_mov_average(PyObject *self, PyObject *args, PyObject *kwds)
{
- PyObject *orig_arrayobj=NULL, *result_dict=NULL;
- PyArrayObject *orig_ndarray=NULL, *result_ndarray=NULL, *result_mask=NULL,
- *mov_sum=NULL;
- PyObject *denom=NULL;
-
+ PyObject *orig_arrayobj=NULL, *orig_ndarray=NULL,
+ *result_ndarray=NULL, *result_mask=NULL,
+ *result_dict=NULL,
+ *mov_sum=NULL, *denom=NULL;
PyArray_Descr *dtype=NULL;
- int rtype, window_size;
+ int rtype, span;
- static char *kwlist[] = {"array", "window_size", "dtype", NULL};
+ static char *kwlist[] = {"array", "span", "dtype", NULL};
if (!PyArg_ParseTupleAndKeywords(args, kwds,
- "Oi|O&:mov_average(array, window_size, dtype)", kwlist,
- &orig_arrayobj, &window_size,
+ "Oi|O&:mov_average(array, span, dtype)", kwlist,
+ &orig_arrayobj, &span,
PyArray_DescrConverter2, &dtype)) return NULL;
- check_mov_args(orig_arrayobj, window_size, 2,
+ check_mov_args(orig_arrayobj, span, 2,
&orig_ndarray, &result_mask);
- rtype = _get_type_num_double(orig_ndarray->descr, dtype);
+ rtype = _get_type_num_double(((PyArrayObject*)orig_ndarray)->descr, dtype);
- mov_sum = (PyArrayObject*)calc_mov_sum(orig_ndarray, window_size, rtype);
+ mov_sum = calc_mov_sum((PyArrayObject*)orig_ndarray, span, rtype);
ERR_CHECK(mov_sum)
- denom = PyFloat_FromDouble(1.0/(double)(window_size));
+ denom = PyFloat_FromDouble(1.0/(double)(span));
- result_ndarray = (PyArrayObject*)PyObject_CallFunction(
- NP_MULTIPLY,
- "OO", mov_sum,
- denom);
+ result_ndarray = np_multiply(mov_sum, denom);
ERR_CHECK(result_ndarray)
+
Py_DECREF(mov_sum);
Py_DECREF(denom);
result_dict = PyDict_New();
MEM_CHECK(result_dict)
- PyDict_SetItemString(result_dict, "array", (PyObject*)result_ndarray);
- PyDict_SetItemString(result_dict, "mask", (PyObject*)result_mask);
+ PyDict_SetItemString(result_dict, "array", result_ndarray);
+ PyDict_SetItemString(result_dict, "mask", result_mask);
Py_DECREF(result_ndarray);
Py_DECREF(result_mask);
return result_dict;
}
-static PyObject*
-np_add(PyObject *left_val, PyObject *right_val) {
- PyObject *result;
-
- result = PyObject_CallFunction(
- NP_ADD, "OO",
- (PyArrayObject*)left_val,
- right_val);
- return result;
-}
-
-static PyObject*
-np_subtract(PyObject *left_val, PyObject *right_val) {
-
- PyObject *result;
-
- result = PyObject_CallFunction(
- NP_SUBTRACT, "OO",
- (PyArrayObject*)left_val,
- right_val);
- return result;
-}
-
-static PyObject*
-np_multiply(PyObject *left_val, PyObject *right_val) {
-
- PyObject *result;
-
- result = PyObject_CallFunction(
- NP_MULTIPLY, "OO",
- (PyArrayObject*)left_val,
- right_val);
- return result;
-}
-
-static int np_greater(PyObject *left_val, PyObject *right_val) {
-
- PyObject *temp;
- int result;
-
- temp = PyObject_CallFunction(
- NP_GREATER, "OO",
- (PyArrayObject*)left_val,
- right_val);
-
- result = (int)PyInt_AsLong(temp);
- Py_DECREF(temp);
- return result;
-}
-
-static int np_greater_equal(PyObject *left_val, PyObject *right_val) {
-
- PyObject *temp;
- int result;
-
- temp = PyObject_CallFunction(
- NP_GREATER_EQUAL, "OO",
- (PyArrayObject*)left_val,
- right_val);
-
- result = (int)PyInt_AsLong(temp);
- Py_DECREF(temp);
- return result;
-}
-
/* computation portion of moving median. Appropriate mask is overlayed on top
afterwards.
@@ -3213,13 +3222,13 @@
(David Brahm) has granted me (and scipy) permission to use it under the BSD
license. */
static PyObject*
-calc_mov_median(PyArrayObject *orig_ndarray, int window_size, int rtype)
+calc_mov_median(PyArrayObject *orig_ndarray, int span, int rtype)
{
PyArrayObject *result_ndarray=NULL;
PyObject **result_array, **ref_array, **even_array=NULL;
- PyObject *new_val, *old_val, *temp;
+ PyObject *new_val, *old_val;
PyObject *temp_add, *one_half;
- int a, i, k, R, arr_size, z, is_odd;
+ int a, i, k, R, arr_size, z;
int *r;
arr_size = orig_ndarray->dimensions[0];
@@ -3230,7 +3239,7 @@
rtype, 0);
ERR_CHECK(result_ndarray)
- if (arr_size >= window_size) {
+ if (arr_size >= span) {
result_array = calloc(arr_size, sizeof(PyObject*));
MEM_CHECK(result_array)
@@ -3246,35 +3255,33 @@
/* this array wll be used for keeping track of the "ranks" of the values
in the current window */
- r = malloc(window_size * sizeof(int));
+ r = malloc(span * sizeof(int));
MEM_CHECK(r)
- for (i=0; i < window_size; i++) {
+ for (i=0; i < span; i++) {
r[i] = 1;
}
- if ((window_size % 2) == 0) {
- // array to store two median values when window_size is an even #
+ if ((span % 2) == 0) {
+ // array to store two median values when span is an even #
even_array = calloc(2, sizeof(PyObject*));
MEM_CHECK(even_array)
}
- R = (window_size + 1)/2;
+ R = (span + 1)/2;
one_half = PyFloat_FromDouble(0.5);
- z = arr_size - window_size;
+ z = arr_size - span;
- //printf("yep 1: %f, %f\n", PyFloat_AsDouble(data_array[z+i]), PyFloat_AsDouble(data_array[z+k]));
-
/* Calculate initial ranks "r" */
- for (i=0; i < window_size; i++) {
+ for (i=0; i < span; i++) {
for (k=0; k < i; k++) {
if (np_greater_equal(ref_array[z+i], ref_array[z+k])) {
r[i]++;
}
}
- for (k=i+1; k < window_size; k++) {
+ for (k=i+1; k < span; k++) {
if (np_greater(ref_array[z+i], ref_array[z+k])) {
r[i]++;
}
@@ -3300,13 +3307,13 @@
Py_DECREF(temp_add);
}
- for (i=arr_size-2; i >= window_size-1; i--) {
- a = window_size;
- z = i - window_size + 1;
+ for (i=arr_size-2; i >= span-1; i--) {
+ a = span;
+ z = i - span + 1;
old_val = ref_array[i+1];
- new_val = ref_array[i-window_size+1];
+ new_val = ref_array[i-span+1];
- for (k=window_size-1; k > 0; k--) {
+ for (k=span-1; k > 0; k--) {
r[k] = r[k-1]; /* Shift previous iteration's ranks */
if (np_greater_equal(ref_array[z+k], new_val)) {r[k]++; a--;}
if (np_greater(ref_array[z+k], old_val)) {r[k]--;}
@@ -3352,7 +3359,7 @@
Py_DECREF(one_half);
- for (i=window_size-1; i<arr_size; i++) {
+ for (i=span-1; i<arr_size; i++) {
PyArray_SETITEM(result_ndarray,
PyArray_GetPtr(result_ndarray, &i),
result_array[i]);
@@ -3363,7 +3370,7 @@
}
if (even_array != NULL) {
- for (i=window_size-1; i<arr_size; i++) {
+ for (i=span-1; i<arr_size; i++) {
Py_DECREF(result_array[i]);
}
}
@@ -3380,31 +3387,32 @@
static PyObject *
MaskedArray_mov_median(PyObject *self, PyObject *args, PyObject *kwds)
{
- PyObject *orig_arrayobj=NULL, *result_dict=NULL;
- PyArrayObject *orig_ndarray=NULL, *result_ndarray=NULL, *result_mask=NULL;
+ PyObject *orig_arrayobj=NULL, *orig_ndarray=NULL,
+ *result_ndarray=NULL, *result_mask=NULL, *result_dict=NULL;
PyArray_Descr *dtype=NULL;
- int rtype, window_size;
+ int rtype, span;
- static char *kwlist[] = {"array", "window_size", "dtype", NULL};
+ static char *kwlist[] = {"array", "span", "dtype", NULL};
if (!PyArg_ParseTupleAndKeywords(args, kwds,
- "Oi|O&:mov_median(array, window_size, dtype)", kwlist,
- &orig_arrayobj, &window_size,
+ "Oi|O&:mov_median(array, span, dtype)", kwlist,
+ &orig_arrayobj, &span,
PyArray_DescrConverter2, &dtype)) return NULL;
- check_mov_args(orig_arrayobj, window_size, 1,
+ check_mov_args(orig_arrayobj, span, 1,
&orig_ndarray, &result_mask);
rtype = _CHKTYPENUM(dtype);
- result_ndarray = (PyArrayObject*)calc_mov_median(orig_ndarray, window_size, rtype);
+ result_ndarray = calc_mov_median((PyArrayObject*)orig_ndarray,
+ span, rtype);
ERR_CHECK(result_ndarray)
result_dict = PyDict_New();
MEM_CHECK(result_dict)
- PyDict_SetItemString(result_dict, "array", (PyObject*)result_ndarray);
- PyDict_SetItemString(result_dict, "mask", (PyObject*)result_mask);
+ PyDict_SetItemString(result_dict, "array", result_ndarray);
+ PyDict_SetItemString(result_dict, "mask", result_mask);
Py_DECREF(result_ndarray);
Py_DECREF(result_mask);
@@ -3416,39 +3424,38 @@
static PyObject *
MaskedArray_mov_stddev(PyObject *self, PyObject *args, PyObject *kwds)
{
- PyObject *orig_arrayobj=NULL, *result_dict=NULL;
- PyArrayObject *orig_ndarray=NULL, *result_ndarray=NULL, *result_mask=NULL,
- *result_temp1=NULL, *result_temp2=NULL, *result_temp3=NULL;
- PyArrayObject *mov_sum=NULL, *mov_sum_sq=NULL;
- PyObject *denom1=NULL, *denom2=NULL;
+ PyObject *orig_ndarray=NULL, *orig_arrayobj=NULL,
+ *result_ndarray=NULL, *result_mask=NULL,
+ *result_dict=NULL,
+ *result_temp1=NULL, *result_temp2=NULL, *result_temp3=NULL,
+ *mov_sum=NULL, *mov_sum_sq=NULL,
+ *denom1=NULL, *denom2=NULL;
+
PyArray_Descr *dtype=NULL;
- int rtype, window_size, is_variance, is_sample;
+ int rtype, span, is_variance, bias;
- static char *kwlist[] = {"array", "window_size", "is_variance", "is_sample", "dtype", NULL};
+ static char *kwlist[] = {"array", "span", "is_variance", "bias", "dtype", NULL};
if (!PyArg_ParseTupleAndKeywords(args, kwds,
- "Oiii|O&:mov_stddev(array, window_size, is_variance, is_sample, dtype)", kwlist,
- &orig_arrayobj, &window_size, &is_variance, &is_sample,
+ "Oiii|O&:mov_stddev(array, span, is_variance, bias, dtype)",
+ kwlist, &orig_arrayobj, &span, &is_variance, &bias,
PyArray_DescrConverter2, &dtype)) return NULL;
- check_mov_args(orig_arrayobj, window_size, 2,
+ check_mov_args(orig_arrayobj, span, 2,
&orig_ndarray, &result_mask);
- rtype = _get_type_num_double(orig_ndarray->descr, dtype);
+ rtype = _get_type_num_double(((PyArrayObject*)orig_ndarray)->descr, dtype);
- mov_sum = (PyArrayObject*)calc_mov_sum(orig_ndarray, window_size, rtype);
+ mov_sum = calc_mov_sum((PyArrayObject*)orig_ndarray, span, rtype);
ERR_CHECK(mov_sum)
- result_temp1 = (PyArrayObject*)PyObject_CallFunction(
- NP_MULTIPLY, "OO",
- orig_ndarray, (PyObject*)orig_ndarray);
+ result_temp1 = np_multiply(orig_ndarray, orig_ndarray);
ERR_CHECK(result_temp1)
- mov_sum_sq = (PyArrayObject*)calc_mov_sum(result_temp1, window_size, rtype);
-
+ mov_sum_sq = calc_mov_sum((PyArrayObject*)result_temp1, span, rtype);
Py_DECREF(result_temp1);
ERR_CHECK(mov_sum_sq)
@@ -3457,40 +3464,29 @@
formulas from:
http://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods
*/
- if (is_sample) {
- denom1 = PyFloat_FromDouble(1.0/(double)(window_size-1));
- denom2 = PyFloat_FromDouble(1.0/(double)(window_size*(window_size-1)));
+ if (bias == 0) {
+ denom1 = PyFloat_FromDouble(1.0/(double)(span-1));
+ denom2 = PyFloat_FromDouble(1.0/(double)(span*(span-1)));
} else {
- denom1 = PyFloat_FromDouble(1.0/(double)window_size);
- denom2 = PyFloat_FromDouble(1.0/(double)(window_size*window_size));
+ denom1 = PyFloat_FromDouble(1.0/(double)span);
+ denom2 = PyFloat_FromDouble(1.0/(double)(span*span));
}
- result_temp1 = (PyArrayObject*)PyObject_CallFunction(
- NP_MULTIPLY,
- "OO", mov_sum_sq,
- denom1);
+ result_temp1 = np_multiply(mov_sum_sq, denom1);
ERR_CHECK(result_temp1)
Py_DECREF(mov_sum_sq);
Py_DECREF(denom1);
- result_temp3 = (PyArrayObject*)PyObject_CallFunction(
- NP_MULTIPLY,
- "OO", mov_sum,
- (PyObject*)mov_sum);
+ result_temp3 = np_multiply(mov_sum, mov_sum);
ERR_CHECK(result_temp3)
Py_DECREF(mov_sum);
- result_temp2 = (PyArrayObject*)PyObject_CallFunction(
- NP_MULTIPLY,
- "OO", result_temp3,
- denom2);
+
+ result_temp2 = np_multiply(result_temp3, denom2);
ERR_CHECK(result_temp2)
Py_DECREF(result_temp3);
Py_DECREF(denom2);
- result_temp3 = (PyArrayObject*)PyObject_CallFunction(
- NP_SUBTRACT,
- "OO", result_temp1,
- (PyObject*)result_temp2);
+ result_temp3 = np_subtract(result_temp1, result_temp2);
ERR_CHECK(result_temp3)
Py_DECREF(result_temp1);
Py_DECREF(result_temp2);
@@ -3498,8 +3494,7 @@
if (is_variance) {
result_ndarray = result_temp3;
} else {
- result_temp1 = (PyArrayObject*)PyObject_CallFunction(
- NP_SQRT, "(O)", result_temp3);
+ result_temp1 = np_sqrt(result_temp3);
ERR_CHECK(result_temp1)
Py_DECREF(result_temp3);
result_ndarray = result_temp1;
@@ -3507,14 +3502,15 @@
result_dict = PyDict_New();
MEM_CHECK(result_dict)
- PyDict_SetItemString(result_dict, "array", (PyObject*)result_ndarray);
- PyDict_SetItemString(result_dict, "mask", (PyObject*)result_mask);
+ PyDict_SetItemString(result_dict, "array", result_ndarray);
+ PyDict_SetItemString(result_dict, "mask", result_mask);
Py_DECREF(result_ndarray);
Py_DECREF(result_mask);
return result_dict;
}
+
static char DateArray_asfreq_doc[] = "";
static PyObject *
DateArray_asfreq(PyObject *self, PyObject *args)
More information about the Scipy-svn
mailing list