← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/sandbox/timeseries/src/c_tseries.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/sandbox/timeseries/src/c_tseries.c
 
1
#include "c_tdates.h"
 
2
#include "c_tseries.h"
 
3
 
 
4
/* Helper function for TimeSeries_convert:
 
5
    determine the size of the second dimension for the resulting
 
6
    converted array */
 
7
static long get_height(int fromFreq, int toFreq) {
 
8
 
 
9
    int maxBusDaysPerYear, maxBusDaysPerQuarter, maxBusDaysPerMonth;
 
10
    int maxDaysPerYear, maxDaysPerQuarter, maxDaysPerMonth;
 
11
 
 
12
    int fromGroup = get_freq_group(fromFreq);
 
13
    int toGroup = get_freq_group(toFreq);
 
14
 
 
15
    if (fromGroup == FR_UND) { fromGroup = FR_DAY; }
 
16
 
 
17
    maxBusDaysPerYear = 262;
 
18
    maxBusDaysPerQuarter = 66;
 
19
    maxBusDaysPerMonth = 23;
 
20
 
 
21
    maxDaysPerYear = 366;
 
22
    maxDaysPerQuarter = 92;
 
23
    maxDaysPerMonth = 31;
 
24
 
 
25
    switch(fromGroup)
 
26
    {
 
27
        case FR_ANN: return 1;
 
28
        case FR_QTR:
 
29
            switch(toGroup)
 
30
            {
 
31
                case FR_ANN: return 4;
 
32
                default: return 1;
 
33
            }
 
34
        case FR_MTH: //monthly
 
35
            switch(toGroup)
 
36
            {
 
37
                case FR_ANN: return 12;
 
38
                case FR_QTR: return 3;
 
39
                default: return 1;
 
40
            }
 
41
        case FR_WK: //weekly
 
42
            switch(toGroup)
 
43
            {
 
44
                case FR_ANN: return 53;
 
45
                case FR_QTR: return 13;
 
46
                case FR_MTH: return 4;
 
47
                default: return 1;
 
48
            }
 
49
        case FR_BUS: //business
 
50
            switch(toGroup)
 
51
            {
 
52
                case FR_ANN: return maxBusDaysPerYear;;
 
53
                case FR_QTR: return maxBusDaysPerQuarter;
 
54
                case FR_MTH: return maxBusDaysPerMonth;
 
55
                case FR_WK: return 5;
 
56
                default: return 1;
 
57
            }
 
58
        case FR_DAY: //daily
 
59
            switch(toGroup)
 
60
            {
 
61
                case FR_ANN: return maxDaysPerYear;;
 
62
                case FR_QTR: return maxDaysPerQuarter;
 
63
                case FR_MTH: return maxDaysPerMonth;
 
64
                case FR_WK: return 7;
 
65
                default: return 1;
 
66
            }
 
67
        case FR_HR: //hourly
 
68
            switch(toGroup)
 
69
            {
 
70
                case FR_ANN: return 24 * maxDaysPerYear;;
 
71
                case FR_QTR: return 24 * maxDaysPerQuarter;
 
72
                case FR_MTH: return 24 * maxDaysPerMonth;
 
73
                case FR_WK: return 24 * 7;
 
74
                case FR_DAY: return 24;
 
75
                case FR_BUS: return 24;
 
76
                default: return 1;
 
77
            }
 
78
        case FR_MIN: //minutely
 
79
            switch(toGroup)
 
80
            {
 
81
                case FR_ANN: return 24 * 60 * maxDaysPerYear;;
 
82
                case FR_QTR: return 24 * 60 * maxDaysPerQuarter;
 
83
                case FR_MTH: return 24 * 60 * maxDaysPerMonth;
 
84
                case FR_WK: return 24 * 60 * 7;
 
85
                case FR_DAY: return 24 * 60;
 
86
                case FR_BUS: return 24 * 60;
 
87
                case FR_HR: return 60;
 
88
                default: return 1;
 
89
            }
 
90
        case FR_SEC: //minutely
 
91
            switch(toGroup)
 
92
            {
 
93
                case FR_ANN: return 24 * 60 * 60 * maxDaysPerYear;;
 
94
                case FR_QTR: return 24 * 60 * 60 * maxDaysPerQuarter;
 
95
                case FR_MTH: return 24 * 60 * 60 * maxDaysPerMonth;
 
96
                case FR_WK: return 24 * 60 * 60 * 7;
 
97
                case FR_DAY: return 24 * 60 * 60;
 
98
                case FR_BUS: return 24 * 60 * 60;
 
99
                case FR_HR: return 60 * 60;
 
100
                case FR_MIN: return 60;
 
101
                default: return 1;
 
102
            }
 
103
        default: return 1;
 
104
    }
 
105
}
 
106
 
 
107
PyObject *
 
108
TimeSeries_convert(PyObject *self, PyObject *args)
 
109
{
 
110
    PyObject *arrayTest;
 
111
    PyArrayObject *array, *newArray;
 
112
    PyArrayObject *mask, *newMask;
 
113
 
 
114
    PyObject *returnVal = NULL;
 
115
    PyObject *start_index_retval;
 
116
 
 
117
    long startIndex;
 
118
    long newStart, newStartTemp;
 
119
    long newEnd, newEndTemp;
 
120
    long newLen, newHeight;
 
121
    long currIndex, prevIndex;
 
122
    long nd;
 
123
    npy_intp *dim, *newIdx;
 
124
    long currPerLen;
 
125
    char *position;
 
126
    PyObject *fromFreq_arg, *toFreq_arg;
 
127
    int fromFreq, toFreq;
 
128
    char relation;
 
129
    asfreq_info af_info;
 
130
    int i;
 
131
 
 
132
    PyObject *val, *valMask;
 
133
 
 
134
    long (*asfreq_main)(long, char, asfreq_info*) = NULL;
 
135
    long (*asfreq_endpoints)(long, char, asfreq_info*) = NULL;
 
136
    long (*asfreq_reverse)(long, char, asfreq_info*) = NULL;
 
137
 
 
138
    returnVal = PyDict_New();
 
139
 
 
140
    if (!PyArg_ParseTuple(args,
 
141
        "OOOslO:convert(array, fromfreq, tofreq, position, startIndex, mask)",
 
142
        &array, &fromFreq_arg, &toFreq_arg,
 
143
        &position, &startIndex, &mask)) return NULL;
 
144
 
 
145
    if((fromFreq = check_freq(fromFreq_arg)) == INT_ERR_CODE) return NULL;
 
146
    if((toFreq = check_freq(toFreq_arg)) == INT_ERR_CODE) return NULL;
 
147
 
 
148
    if (toFreq == fromFreq)
 
149
    {
 
150
        PyObject *sidx;
 
151
        newArray = (PyArrayObject *)PyArray_Copy(array);
 
152
        newMask = (PyArrayObject *)PyArray_Copy(mask);
 
153
        sidx = PyInt_FromLong(startIndex);
 
154
 
 
155
        PyDict_SetItemString(returnVal, "values", (PyObject*)newArray);
 
156
        PyDict_SetItemString(returnVal, "mask", (PyObject*)newMask);
 
157
        PyDict_SetItemString(returnVal, "startindex", sidx);
 
158
 
 
159
        Py_DECREF(newArray);
 
160
        Py_DECREF(newMask);
 
161
        Py_DECREF(sidx);
 
162
 
 
163
        return returnVal;
 
164
    }
 
165
 
 
166
    switch(position[0])
 
167
    {
 
168
        case 'S':
 
169
            // start -> before
 
170
            relation = 'B';
 
171
            break;
 
172
        case 'E':
 
173
            // end -> after
 
174
            relation = 'A';
 
175
            break;
 
176
        default:
 
177
            return NULL;
 
178
            break;
 
179
    }
 
180
 
 
181
    get_asfreq_info(fromFreq, toFreq, &af_info);
 
182
 
 
183
    asfreq_main = get_asfreq_func(fromFreq, toFreq, 1);
 
184
    asfreq_endpoints = get_asfreq_func(fromFreq, toFreq, 0);
 
185
 
 
186
    //convert start index to new frequency
 
187
    CHECK_ASFREQ(newStartTemp = asfreq_main(startIndex, 'B', &af_info));
 
188
    if (newStartTemp < 1) {
 
189
        CHECK_ASFREQ(newStart = asfreq_endpoints(startIndex, 'A', &af_info));
 
190
    }
 
191
    else { newStart = newStartTemp; }
 
192
 
 
193
    //convert end index to new frequency
 
194
    CHECK_ASFREQ(newEndTemp = asfreq_main(startIndex+array->dimensions[0]-1, 'A', &af_info));
 
195
    if (newEndTemp < 1) {
 
196
        CHECK_ASFREQ(newEnd = asfreq_endpoints(startIndex+array->dimensions[0]-1, 'B', &af_info));
 
197
    }
 
198
    else { newEnd = newEndTemp; }
 
199
 
 
200
    if (newStart < 1) {
 
201
        PyErr_SetString(PyExc_ValueError, "start_date outside allowable range for destination frequency");
 
202
        return NULL;
 
203
    }
 
204
 
 
205
    newLen = newEnd - newStart + 1;
 
206
    newHeight = get_height(fromFreq, toFreq);
 
207
 
 
208
    if (newHeight > 1) {
 
209
        long tempval;
 
210
        asfreq_info af_info_rev;
 
211
 
 
212
        get_asfreq_info(toFreq, fromFreq, &af_info_rev);
 
213
        asfreq_reverse = get_asfreq_func(toFreq, fromFreq, 0);
 
214
 
 
215
        CHECK_ASFREQ(tempval = asfreq_reverse(newStart, 'B', &af_info_rev));
 
216
        currPerLen = startIndex - tempval;
 
217
 
 
218
        nd = 2;
 
219
        dim = PyDimMem_NEW(nd);
 
220
        dim[0] = (npy_intp)newLen;
 
221
        dim[1] = (npy_intp)newHeight;
 
222
    } else {
 
223
        nd = 1;
 
224
        dim = PyDimMem_NEW(nd);
 
225
        dim[0] = (npy_intp)newLen;
 
226
    }
 
227
 
 
228
    newIdx = PyDimMem_NEW(nd);
 
229
    arrayTest = PyArray_SimpleNew(nd, dim, array->descr->type_num);
 
230
    if (arrayTest == NULL) { return NULL; }
 
231
    newArray = (PyArrayObject*)arrayTest;
 
232
    newMask  = (PyArrayObject*)PyArray_SimpleNew(nd, dim, mask->descr->type_num);
 
233
 
 
234
    PyDimMem_FREE(dim);
 
235
 
 
236
    PyArray_FILLWBYTE(newArray,0);
 
237
    PyArray_FILLWBYTE(newMask,1);
 
238
 
 
239
    prevIndex = newStart;
 
240
 
 
241
    //set values in the new array
 
242
    for (i = 0; i < array->dimensions[0]; i++) {
 
243
 
 
244
        npy_intp idx = (npy_intp)i;
 
245
 
 
246
        val = PyArray_GETITEM(array, PyArray_GetPtr(array, &idx));
 
247
        valMask = PyArray_GETITEM(mask, PyArray_GetPtr(mask, &idx));
 
248
 
 
249
        CHECK_ASFREQ(currIndex = asfreq_main(startIndex + i, relation, &af_info));
 
250
 
 
251
        newIdx[0] = (npy_intp)(currIndex-newStart);
 
252
 
 
253
        if (newHeight > 1) {
 
254
 
 
255
                if (currIndex != prevIndex)
 
256
                {
 
257
                    //reset period length
 
258
                    currPerLen = 0;
 
259
                    prevIndex = currIndex;
 
260
                }
 
261
 
 
262
                newIdx[1] = (npy_intp)currPerLen;
 
263
                currPerLen++;
 
264
        }
 
265
 
 
266
        if (newIdx[0] > -1) {
 
267
            PyArray_SETITEM(newArray, PyArray_GetPtr(newArray, newIdx), val);
 
268
            PyArray_SETITEM(newMask, PyArray_GetPtr(newMask, newIdx), valMask);
 
269
        }
 
270
 
 
271
        Py_DECREF(val);
 
272
        Py_DECREF(valMask);
 
273
 
 
274
    }
 
275
 
 
276
    PyDimMem_FREE(newIdx);
 
277
 
 
278
    start_index_retval = (PyObject*)PyInt_FromLong(newStart);
 
279
 
 
280
    PyDict_SetItemString(returnVal, "values", (PyObject*)newArray);
 
281
    PyDict_SetItemString(returnVal, "mask", (PyObject*)newMask);
 
282
    PyDict_SetItemString(returnVal, "startindex", start_index_retval);
 
283
 
 
284
    Py_DECREF(newArray);
 
285
    Py_DECREF(newMask);
 
286
    Py_DECREF(start_index_retval);
 
287
 
 
288
    return returnVal;
 
289
}
 
290
 
 
291
 
 
292
/* This function is directly copied from the numpy source  */
 
293
/* Return typenumber from dtype2 unless it is NULL, then return
 
294
   NPY_DOUBLE if dtype1->type_num is integer or bool
 
295
   and dtype1->type_num otherwise.
 
296
*/
 
297
static int
 
298
_get_type_num_double(PyArray_Descr *dtype1, PyArray_Descr *dtype2)
 
299
{
 
300
    if (dtype2 != NULL) {
 
301
        return dtype2->type_num;
 
302
    }
 
303
 
 
304
    /* For integer or bool data-types */
 
305
    if (dtype1->type_num < NPY_FLOAT) {
 
306
        return NPY_DOUBLE;
 
307
    }
 
308
    else {
 
309
        return dtype1->type_num;
 
310
    }
 
311
}
 
312
 
 
313
static int
 
314
_get_type_num(PyArray_Descr *dtype1, PyArray_Descr *dtype2)
 
315
{
 
316
    if (dtype2 != NULL) {
 
317
        return dtype2->type_num;
 
318
    } else {
 
319
        return dtype1->type_num;
 
320
    }
 
321
}
 
322
 
 
323
 
 
324
/* validates the standard arguments to moving functions and set the original
 
325
   mask, original ndarray, and mask for the result */
 
326
static PyObject *
 
327
check_mov_args(PyObject *orig_arrayobj, int span, int min_win_size,
 
328
               PyObject **orig_ndarray, PyObject **result_mask) {
 
329
 
 
330
    PyObject *orig_mask=NULL;
 
331
    PyArrayObject **orig_ndarray_tmp, **result_mask_tmp;
 
332
    int *raw_result_mask;
 
333
 
 
334
    if (!PyArray_Check(orig_arrayobj)) {
 
335
        PyErr_SetString(PyExc_ValueError, "array must be a valid subtype of ndarray");
 
336
        return NULL;
 
337
    }
 
338
 
 
339
    // check if array has a mask, and if that mask is an array
 
340
    if (PyObject_HasAttrString(orig_arrayobj, "_mask")) {
 
341
        PyObject *tempMask = PyObject_GetAttrString(orig_arrayobj, "_mask");
 
342
        if (PyArray_Check(tempMask)) {
 
343
            orig_mask = PyArray_EnsureArray(tempMask);
 
344
        } else {
 
345
            Py_DECREF(tempMask);
 
346
        }
 
347
    }
 
348
 
 
349
    *orig_ndarray = PyArray_EnsureArray(orig_arrayobj);
 
350
    orig_ndarray_tmp = (PyArrayObject**)orig_ndarray;
 
351
 
 
352
    if ((*orig_ndarray_tmp)->nd != 1) {
 
353
        PyErr_SetString(PyExc_ValueError, "array must be 1 dimensional");
 
354
        return NULL;
 
355
    }
 
356
 
 
357
    if (span < min_win_size) {
 
358
        char *error_str;
 
359
        error_str = malloc(60 * sizeof(char));
 
360
        MEM_CHECK(error_str)
 
361
        sprintf(error_str,
 
362
                "span must be greater than or equal to %i",
 
363
                min_win_size);
 
364
        PyErr_SetString(PyExc_ValueError, error_str);
 
365
        free(error_str);
 
366
        return NULL;
 
367
    }
 
368
 
 
369
    raw_result_mask = malloc((*orig_ndarray_tmp)->dimensions[0] * sizeof(int));
 
370
    MEM_CHECK(raw_result_mask)
 
371
 
 
372
    {
 
373
        PyArrayObject *orig_mask_tmp;
 
374
        int i, valid_points=0, is_masked;
 
375
 
 
376
        orig_mask_tmp = (PyArrayObject*)orig_mask;
 
377
 
 
378
        for (i=0; i<((*orig_ndarray_tmp)->dimensions[0]); i++) {
 
379
 
 
380
            npy_intp idx = (npy_intp)i;
 
381
            is_masked=0;
 
382
 
 
383
            if (orig_mask != NULL) {
 
384
                PyObject *valMask;
 
385
                valMask = PyArray_GETITEM(orig_mask_tmp,
 
386
                                          PyArray_GetPtr(orig_mask_tmp, &idx));
 
387
                is_masked = (int)PyInt_AsLong(valMask);
 
388
                Py_DECREF(valMask);
 
389
            }
 
390
 
 
391
            if (is_masked) {
 
392
                valid_points=0;
 
393
            } else {
 
394
                if (valid_points < span) { valid_points += 1; }
 
395
                if (valid_points < span) { is_masked = 1; }
 
396
            }
 
397
 
 
398
            raw_result_mask[i] = is_masked;
 
399
        }
 
400
    }
 
401
 
 
402
    *result_mask = PyArray_SimpleNewFromData(
 
403
                             1, (*orig_ndarray_tmp)->dimensions,
 
404
                             PyArray_INT32, raw_result_mask);
 
405
    MEM_CHECK(*result_mask)
 
406
    result_mask_tmp = (PyArrayObject**)result_mask;
 
407
    (*result_mask_tmp)->flags = ((*result_mask_tmp)->flags) | NPY_OWNDATA;
 
408
    return 0;
 
409
}
 
410
 
 
411
/* computation portion of moving sum. Appropriate mask is overlayed on top
 
412
   afterwards */
 
413
static PyObject*
 
414
calc_mov_sum(PyArrayObject *orig_ndarray, int span, int rtype)
 
415
{
 
416
    PyArrayObject *result_ndarray=NULL;
 
417
    int i;
 
418
 
 
419
    result_ndarray = (PyArrayObject*)PyArray_ZEROS(
 
420
                                       orig_ndarray->nd,
 
421
                                       orig_ndarray->dimensions,
 
422
                                       rtype, 0);
 
423
    ERR_CHECK(result_ndarray)
 
424
 
 
425
    for (i=0; i<orig_ndarray->dimensions[0]; i++) {
 
426
 
 
427
        PyObject *val=NULL, *mov_sum_val=NULL;
 
428
        npy_intp idx = (npy_intp)i;
 
429
 
 
430
        val = PyArray_GETITEM(orig_ndarray, PyArray_GetPtr(orig_ndarray, &idx));
 
431
 
 
432
        if (i == 0) {
 
433
            mov_sum_val = val;
 
434
        } else {
 
435
            idx = (npy_intp)(i-1);
 
436
            PyObject *mov_sum_prevval;
 
437
            mov_sum_prevval= PyArray_GETITEM(result_ndarray,
 
438
                                   PyArray_GetPtr(result_ndarray, &idx));
 
439
            mov_sum_val = np_add(val, mov_sum_prevval);
 
440
            Py_DECREF(mov_sum_prevval);
 
441
            ERR_CHECK(mov_sum_val)
 
442
 
 
443
            if (i >= span) {
 
444
                PyObject *temp_val, *rem_val;
 
445
                idx = (npy_intp)(i-span);
 
446
                temp_val = mov_sum_val;
 
447
                rem_val = PyArray_GETITEM(orig_ndarray,
 
448
                                   PyArray_GetPtr(orig_ndarray, &idx));
 
449
 
 
450
                mov_sum_val = np_subtract(temp_val, rem_val);
 
451
                ERR_CHECK(mov_sum_val)
 
452
 
 
453
                Py_DECREF(temp_val);
 
454
                Py_DECREF(rem_val);
 
455
            }
 
456
        }
 
457
 
 
458
        idx = (npy_intp)i;
 
459
 
 
460
        PyArray_SETITEM(result_ndarray,
 
461
                        PyArray_GetPtr(result_ndarray, &idx),
 
462
                        mov_sum_val);
 
463
 
 
464
        if (mov_sum_val != val) { Py_DECREF(val); }
 
465
 
 
466
        Py_DECREF(mov_sum_val);
 
467
    }
 
468
 
 
469
    return (PyObject*)result_ndarray;
 
470
 
 
471
}
 
472
 
 
473
PyObject *
 
474
MaskedArray_mov_sum(PyObject *self, PyObject *args, PyObject *kwds)
 
475
{
 
476
    PyObject *orig_arrayobj=NULL, *orig_ndarray=NULL,
 
477
             *result_ndarray=NULL, *result_mask=NULL,
 
478
             *result_dict=NULL;
 
479
    PyArray_Descr *dtype=NULL;
 
480
 
 
481
    int rtype, span;
 
482
 
 
483
    static char *kwlist[] = {"array", "span", "dtype", NULL};
 
484
 
 
485
    if (!PyArg_ParseTupleAndKeywords(args, kwds,
 
486
                "Oi|O&:mov_sum(array, span, dtype)", kwlist,
 
487
                &orig_arrayobj, &span,
 
488
                PyArray_DescrConverter2, &dtype)) return NULL;
 
489
 
 
490
    check_mov_args(orig_arrayobj, span, 1,
 
491
                   &orig_ndarray, &result_mask);
 
492
 
 
493
    rtype = _get_type_num(((PyArrayObject*)orig_ndarray)->descr, dtype);
 
494
 
 
495
    result_ndarray = calc_mov_sum((PyArrayObject*)orig_ndarray,
 
496
                                  span, rtype);
 
497
    ERR_CHECK(result_ndarray)
 
498
 
 
499
    result_dict = PyDict_New();
 
500
    MEM_CHECK(result_dict)
 
501
    PyDict_SetItemString(result_dict, "array", result_ndarray);
 
502
    PyDict_SetItemString(result_dict, "mask", result_mask);
 
503
 
 
504
    Py_DECREF(result_ndarray);
 
505
    Py_DECREF(result_mask);
 
506
    return result_dict;
 
507
}
 
508
 
 
509
PyObject *
 
510
MaskedArray_mov_average(PyObject *self, PyObject *args, PyObject *kwds)
 
511
{
 
512
    PyObject *orig_arrayobj=NULL, *orig_ndarray=NULL,
 
513
             *result_ndarray=NULL, *result_mask=NULL,
 
514
             *result_dict=NULL,
 
515
             *mov_sum=NULL, *denom=NULL;
 
516
    PyArray_Descr *dtype=NULL;
 
517
 
 
518
    int rtype, span;
 
519
 
 
520
    static char *kwlist[] = {"array", "span", "dtype", NULL};
 
521
 
 
522
    if (!PyArg_ParseTupleAndKeywords(args, kwds,
 
523
                "Oi|O&:mov_average(array, span, dtype)", kwlist,
 
524
                &orig_arrayobj, &span,
 
525
                PyArray_DescrConverter2, &dtype)) return NULL;
 
526
 
 
527
 
 
528
    check_mov_args(orig_arrayobj, span, 2,
 
529
                   &orig_ndarray, &result_mask);
 
530
 
 
531
    rtype = _get_type_num_double(((PyArrayObject*)orig_ndarray)->descr, dtype);
 
532
 
 
533
    mov_sum = calc_mov_sum((PyArrayObject*)orig_ndarray, span, rtype);
 
534
    ERR_CHECK(mov_sum)
 
535
 
 
536
    denom = PyFloat_FromDouble(1.0/(double)(span));
 
537
 
 
538
    result_ndarray = np_multiply(mov_sum, denom);
 
539
    ERR_CHECK(result_ndarray)
 
540
 
 
541
    Py_DECREF(mov_sum);
 
542
    Py_DECREF(denom);
 
543
 
 
544
    result_dict = PyDict_New();
 
545
    MEM_CHECK(result_dict)
 
546
    PyDict_SetItemString(result_dict, "array", result_ndarray);
 
547
    PyDict_SetItemString(result_dict, "mask", result_mask);
 
548
 
 
549
    Py_DECREF(result_ndarray);
 
550
    Py_DECREF(result_mask);
 
551
    return result_dict;
 
552
}
 
553
 
 
554
 
 
555
/* computation portion of moving median. Appropriate mask is overlayed on top
 
556
   afterwards.
 
557
 
 
558
   The algorithm used here is based on the code found at:
 
559
    http://cran.r-project.org/src/contrib/Devel/runStat_1.1.tar.gz
 
560
 
 
561
   This code was originally released under the GPL, but the author
 
562
   (David Brahm) has granted me (and scipy) permission to use it under the BSD
 
563
   license. */
 
564
PyObject*
 
565
calc_mov_median(PyArrayObject *orig_ndarray, int span, int rtype)
 
566
{
 
567
    PyArrayObject *result_ndarray=NULL;
 
568
    PyObject **result_array, **ref_array, **even_array=NULL;
 
569
    PyObject *new_val, *old_val;
 
570
    PyObject *temp_add, *one_half;
 
571
    int a, i, k, R, arr_size, z;
 
572
    int *r;
 
573
    npy_intp idx;
 
574
 
 
575
    arr_size = (int)(orig_ndarray->dimensions[0]);
 
576
 
 
577
    result_ndarray = (PyArrayObject*)PyArray_ZEROS(
 
578
                                       orig_ndarray->nd,
 
579
                                       orig_ndarray->dimensions,
 
580
                                       rtype, 0);
 
581
    ERR_CHECK(result_ndarray)
 
582
 
 
583
    if (arr_size >= span) {
 
584
        result_array = calloc(arr_size, sizeof(PyObject*));
 
585
        MEM_CHECK(result_array)
 
586
 
 
587
        /* this array will be used for quick access to the data in the original
 
588
           array (so PyArray_GETITEM doesn't have to be used over and over in the
 
589
           main loop) */
 
590
        ref_array = malloc(arr_size * sizeof(PyObject*));
 
591
        MEM_CHECK(ref_array)
 
592
 
 
593
        for (i=0; i<arr_size; i++) {
 
594
            idx = (npy_intp)i;
 
595
            ref_array[i] = PyArray_GETITEM(orig_ndarray, PyArray_GetPtr(orig_ndarray, &idx));
 
596
        }
 
597
 
 
598
        /* this array wll be used for keeping track of the "ranks" of the values
 
599
           in the current window */
 
600
        r = malloc(span * sizeof(int));
 
601
        MEM_CHECK(r)
 
602
 
 
603
        for (i=0; i < span; i++) {
 
604
            r[i] = 1;
 
605
        }
 
606
 
 
607
        if ((span % 2) == 0) {
 
608
            // array to store two median values when span is an even #
 
609
            even_array = calloc(2, sizeof(PyObject*));
 
610
            MEM_CHECK(even_array)
 
611
        }
 
612
 
 
613
        R = (span + 1)/2;
 
614
        one_half = PyFloat_FromDouble(0.5);
 
615
 
 
616
        z = arr_size - span;
 
617
 
 
618
        /* Calculate initial ranks "r" */
 
619
        for (i=0; i < span; i++) {
 
620
 
 
621
            for (k=0;   k < i;  k++) {
 
622
                if (np_greater_equal(ref_array[z+i], ref_array[z+k])) {
 
623
                    r[i]++;
 
624
                }
 
625
            }
 
626
            for (k=i+1; k < span; k++) {
 
627
                if (np_greater(ref_array[z+i], ref_array[z+k])) {
 
628
                    r[i]++;
 
629
                }
 
630
            }
 
631
 
 
632
            /* If rank=R, this is the median */
 
633
            if (even_array != NULL) {
 
634
                if (r[i]==R) {
 
635
                    even_array[0] = ref_array[z+i];
 
636
                } else if (r[i] == (R+1)) {
 
637
                    even_array[1] = ref_array[z+i];
 
638
                }
 
639
            } else {
 
640
                if (r[i]==R) {
 
641
                    result_array[arr_size-1] = ref_array[z+i];
 
642
                }
 
643
            }
 
644
        }
 
645
 
 
646
        if (even_array != NULL) {
 
647
            temp_add = np_add(even_array[0], even_array[1]);
 
648
            result_array[arr_size-1] = np_multiply(temp_add, one_half);
 
649
            Py_DECREF(temp_add);
 
650
        }
 
651
 
 
652
        for (i=arr_size-2; i >= span-1; i--) {
 
653
            a = span;
 
654
            z = i - span + 1;
 
655
            old_val = ref_array[i+1];
 
656
            new_val = ref_array[i-span+1];
 
657
 
 
658
            for (k=span-1; k > 0; k--) {
 
659
                r[k] = r[k-1]; /* Shift previous iteration's ranks */
 
660
                if (np_greater_equal(ref_array[z+k], new_val)) {r[k]++; a--;}
 
661
                if (np_greater(ref_array[z+k], old_val)) {r[k]--;}
 
662
 
 
663
                if (r[k]==R) {
 
664
                    result_array[i] = ref_array[z+k];
 
665
                }
 
666
 
 
667
                if (even_array != NULL) {
 
668
                    if (r[k]==R) {
 
669
                        even_array[0] = ref_array[z+k];
 
670
                    } else if (r[k] == (R+1)) {
 
671
                        even_array[1] = ref_array[z+k];
 
672
                    }
 
673
                } else {
 
674
                    if (r[k]==R) {
 
675
                        result_array[i] = ref_array[z+k];
 
676
                    }
 
677
                }
 
678
 
 
679
            }
 
680
 
 
681
            r[0] = a;
 
682
 
 
683
            if (even_array != NULL) {
 
684
                if (a==R) {
 
685
                    even_array[0] = new_val;
 
686
                } else if (a == (R+1)) {
 
687
                    even_array[1] = new_val;
 
688
                }
 
689
 
 
690
                temp_add = np_add(even_array[0], even_array[1]);
 
691
                result_array[i] = np_multiply(temp_add, one_half);;
 
692
                Py_DECREF(temp_add);
 
693
 
 
694
            } else {
 
695
                if (a==R) {
 
696
                    result_array[i] = new_val;
 
697
                }
 
698
            }
 
699
 
 
700
        }
 
701
 
 
702
        Py_DECREF(one_half);
 
703
 
 
704
        for (i=span-1; i<arr_size; i++) {
 
705
            idx = (npy_intp)i;
 
706
            PyArray_SETITEM(result_ndarray,
 
707
                            PyArray_GetPtr(result_ndarray, &idx),
 
708
                            result_array[i]);
 
709
        }
 
710
 
 
711
        for (i=0; i<arr_size; i++) {
 
712
            Py_DECREF(ref_array[i]);
 
713
        }
 
714
 
 
715
        if (even_array != NULL) {
 
716
            for (i=span-1; i<arr_size; i++) {
 
717
                Py_DECREF(result_array[i]);
 
718
            }
 
719
        }
 
720
 
 
721
        free(ref_array);
 
722
        free(result_array);
 
723
    }
 
724
 
 
725
    return (PyObject*)result_ndarray;
 
726
 
 
727
}
 
728
 
 
729
PyObject *
 
730
MaskedArray_mov_median(PyObject *self, PyObject *args, PyObject *kwds)
 
731
{
 
732
    PyObject *orig_arrayobj=NULL, *orig_ndarray=NULL,
 
733
             *result_ndarray=NULL, *result_mask=NULL, *result_dict=NULL;
 
734
    PyArray_Descr *dtype=NULL;
 
735
 
 
736
    int rtype, span;
 
737
 
 
738
    static char *kwlist[] = {"array", "span", "dtype", NULL};
 
739
 
 
740
    if (!PyArg_ParseTupleAndKeywords(args, kwds,
 
741
                "Oi|O&:mov_median(array, span, dtype)", kwlist,
 
742
                &orig_arrayobj, &span,
 
743
                PyArray_DescrConverter2, &dtype)) return NULL;
 
744
 
 
745
    check_mov_args(orig_arrayobj, span, 1,
 
746
                   &orig_ndarray, &result_mask);
 
747
 
 
748
    if ((span % 2) == 0) {
 
749
        rtype = _get_type_num_double(((PyArrayObject*)orig_ndarray)->descr, dtype);
 
750
    } else {
 
751
        rtype = _get_type_num(((PyArrayObject*)orig_ndarray)->descr, dtype);
 
752
    }
 
753
 
 
754
    result_ndarray = calc_mov_median((PyArrayObject*)orig_ndarray,
 
755
                                     span, rtype);
 
756
    ERR_CHECK(result_ndarray)
 
757
 
 
758
    result_dict = PyDict_New();
 
759
    MEM_CHECK(result_dict)
 
760
    PyDict_SetItemString(result_dict, "array", result_ndarray);
 
761
    PyDict_SetItemString(result_dict, "mask", result_mask);
 
762
 
 
763
    Py_DECREF(result_ndarray);
 
764
    Py_DECREF(result_mask);
 
765
    return result_dict;
 
766
}
 
767
 
 
768
 
 
769
PyObject *
 
770
MaskedArray_mov_stddev(PyObject *self, PyObject *args, PyObject *kwds)
 
771
{
 
772
 
 
773
    PyObject *orig_ndarray=NULL, *orig_arrayobj=NULL,
 
774
             *result_ndarray=NULL, *result_mask=NULL,
 
775
             *result_dict=NULL,
 
776
             *result_temp1=NULL, *result_temp2=NULL, *result_temp3=NULL,
 
777
             *mov_sum=NULL, *mov_sum_sq=NULL,
 
778
             *denom1=NULL, *denom2=NULL;
 
779
 
 
780
    PyArray_Descr *dtype=NULL;
 
781
 
 
782
    int rtype, span, is_variance, bias;
 
783
 
 
784
    static char *kwlist[] = {"array", "span", "is_variance", "bias", "dtype", NULL};
 
785
 
 
786
    if (!PyArg_ParseTupleAndKeywords(args, kwds,
 
787
                "Oiii|O&:mov_stddev(array, span, is_variance, bias, dtype)",
 
788
                kwlist, &orig_arrayobj, &span, &is_variance, &bias,
 
789
                PyArray_DescrConverter2, &dtype)) return NULL;
 
790
 
 
791
 
 
792
    check_mov_args(orig_arrayobj, span, 2,
 
793
                   &orig_ndarray, &result_mask);
 
794
 
 
795
    rtype = _get_type_num_double(((PyArrayObject*)orig_ndarray)->descr, dtype);
 
796
 
 
797
    mov_sum = calc_mov_sum((PyArrayObject*)orig_ndarray, span, rtype);
 
798
    ERR_CHECK(mov_sum)
 
799
 
 
800
    result_temp1 = np_multiply(orig_ndarray, orig_ndarray);
 
801
    ERR_CHECK(result_temp1)
 
802
 
 
803
    mov_sum_sq = calc_mov_sum((PyArrayObject*)result_temp1, span, rtype);
 
804
    Py_DECREF(result_temp1);
 
805
    ERR_CHECK(mov_sum_sq)
 
806
 
 
807
 
 
808
    /*
 
809
      formulas from:
 
810
      http://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods
 
811
     */
 
812
    if (bias == 0) {
 
813
        denom1 = PyFloat_FromDouble(1.0/(double)(span-1));
 
814
        denom2 = PyFloat_FromDouble(1.0/(double)(span*(span-1)));
 
815
    } else {
 
816
        denom1 = PyFloat_FromDouble(1.0/(double)span);
 
817
        denom2 = PyFloat_FromDouble(1.0/(double)(span*span));
 
818
    }
 
819
 
 
820
    result_temp1 = np_multiply(mov_sum_sq, denom1);
 
821
    ERR_CHECK(result_temp1)
 
822
    Py_DECREF(mov_sum_sq);
 
823
    Py_DECREF(denom1);
 
824
 
 
825
    result_temp3 = np_multiply(mov_sum, mov_sum);
 
826
    ERR_CHECK(result_temp3)
 
827
    Py_DECREF(mov_sum);
 
828
 
 
829
    result_temp2 = np_multiply(result_temp3, denom2);
 
830
    ERR_CHECK(result_temp2)
 
831
    Py_DECREF(result_temp3);
 
832
    Py_DECREF(denom2);
 
833
 
 
834
    result_temp3 = np_subtract(result_temp1, result_temp2);
 
835
    ERR_CHECK(result_temp3)
 
836
    Py_DECREF(result_temp1);
 
837
    Py_DECREF(result_temp2);
 
838
 
 
839
    if (is_variance) {
 
840
        result_ndarray = result_temp3;
 
841
    } else {
 
842
        result_temp1 = np_sqrt(result_temp3);
 
843
        ERR_CHECK(result_temp1)
 
844
        Py_DECREF(result_temp3);
 
845
        result_ndarray = result_temp1;
 
846
    }
 
847
 
 
848
    result_dict = PyDict_New();
 
849
    MEM_CHECK(result_dict)
 
850
    PyDict_SetItemString(result_dict, "array", result_ndarray);
 
851
    PyDict_SetItemString(result_dict, "mask", result_mask);
 
852
 
 
853
    Py_DECREF(result_ndarray);
 
854
    Py_DECREF(result_mask);
 
855
    return result_dict;
 
856
 
 
857
}
 
858
 
 
859
void import_c_tseries(PyObject *m) { import_array(); }