source: anuga_core/source/anuga/shallow_water/urs_ext.c @ 5552

Last change on this file since 5552 was 5552, checked in by nariman, 11 years ago

Refactoring of code and creation of _read_mux2_headers.
Memory leak caused by muxData.

File size: 21.8 KB
Line 
1/*
2gcc -fPIC -c urs_ext.c -I/usr/include/python2.5 -o urs_ext.o -Wall -O
3gcc -shared urs_ext.o  -o urs_ext.so
4*/
5
6/*
7This file was reverted from changeset:5484 to changeset:5470 on 10th July
8by Ole.
9*/
10
11#include "Python.h"
12#include "Numeric/arrayobject.h"
13#include "structure.h"
14#include "math.h"
15#include <stdio.h>
16#include <float.h>
17#include <time.h>
18
19#define MAX_FILE_NAME_LENGTH 128
20#define NODATA 99.0
21#define EPSILON  0.00001
22
23#define DEBUG 0
24
25#define POFFSET 5 //Number of site_params
26
27void fillDataArray(int, int, int, int, int *, int *, float *, 
28           int *, int *, float *);
29long getNumData(const int *fros, const int *lros, const int nsta);
30char isdata(float);
31void _read_mux2_headers(int numSrc, 
32                           char **muxFileNameArray, 
33                   double *params,
34                                   int** fros, int** lros, struct tgsrwg ** mytgs0,
35                   int* nsta0,
36                   int verbose);
37
38float** _read_mux2(int numSrc, 
39           char **muxFileNameArray, 
40           float *weights, 
41           double *params, 
42           int *number_of_stations,
43           long *permutation,
44           int verbose);
45
46    PyObject *read_mux2(PyObject *self, PyObject *args){
47    /*Read in mux 2 file
48
49    Python call:
50    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
51
52    NOTE:
53    A Python int is equivalent to a C long
54    (this becomes really important on 64 bit architectures)
55   
56    A Python double corresponds to a C double
57    */
58   
59    PyObject *filenames;
60    PyArrayObject *pyweights;
61    PyArrayObject *file_params;
62    PyArrayObject *permutation;  // Ordering of selected stations   
63    PyArrayObject *pydata;
64    PyObject *fname;
65
66    char **muxFileNameArray;
67    float **cdata;
68    float *weights;
69    int dimensions[2];
70    int numSrc;
71    int verbose;
72    int nsta0;
73    int number_of_selected_stations;   
74    int nt;
75    double dt;
76    int i;
77    int j;
78    int start_tstep;
79    int finish_tstep;
80    int it;
81    int time;
82    int num_ts;
83   
84    // Convert Python arguments to C
85    if (!PyArg_ParseTuple(args, "iOOOOi",
86              &numSrc, &filenames, &pyweights, &file_params, 
87              &permutation, &verbose)) 
88    {
89            PyErr_SetString(PyExc_RuntimeError, 
90                "Input arguments to read_mux2 failed");
91            return NULL;
92    }
93
94    if(!PyList_Check(filenames)) 
95    {
96        PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list");
97        return NULL;
98    }
99
100    if(PyList_Size(filenames) == 0)
101    {
102        PyErr_SetString(PyExc_ValueError, "empty lists not allowed");
103        return NULL;
104    }
105
106    if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) 
107    {
108        PyErr_SetString(PyExc_ValueError,
109            "pyweights must be one-dimensional and of type double");
110        return NULL; 
111    }
112
113    if(PyList_Size(filenames) != pyweights->dimensions[0])
114    {
115        PyErr_SetString(PyExc_ValueError, 
116            "Must specify one weight for each filename");
117        return NULL;
118    }
119
120    muxFileNameArray = (char**)malloc(numSrc*sizeof(char*));
121    if (muxFileNameArray == NULL) 
122    {
123        PyErr_SetString(PyExc_ValueError, 
124      "ERROR: Memory for muxFileNameArray could not be allocated.");
125        return NULL;
126    }
127
128    for (i = 0; i < numSrc; i++)
129    {
130
131        fname = PyList_GetItem(filenames, i);
132        if (!fname) 
133        {
134            PyErr_SetString(PyExc_ValueError, "filename not a string");
135            return NULL;
136        }       
137
138        muxFileNameArray[i] = PyString_AsString(fname);
139        if (muxFileNameArray[i] == NULL) 
140        {
141            PyErr_SetString(PyExc_ValueError, 
142          "ERROR: Memory for muxFileNameArray could not be allocated.\n");
143            return NULL;
144        }
145    }
146
147    if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) 
148    {
149        PyErr_SetString(PyExc_ValueError,
150          "file_params must be one-dimensional and of type double");
151        return NULL; 
152    }
153
154
155    // Create array for weights which are passed to read_mux2
156    weights = (float*) malloc(numSrc*sizeof(float));
157    for (i = 0; i < numSrc; i++)
158    {
159        weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0]));
160    }
161   
162    // Desired number of stations
163    number_of_selected_stations = (int) permutation->dimensions[0];
164
165    // Read in mux2 data from file
166    cdata = _read_mux2(numSrc, 
167                       muxFileNameArray, 
168                       weights, 
169                       (double*)file_params->data,
170                       &number_of_selected_stations, 
171                       (long*) permutation->data,
172                       verbose);
173
174    if (!cdata) 
175    {
176        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
177        return NULL;
178    }       
179               
180               
181    // Allocate space for return vector
182    nsta0 = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
183    dt = *(double*)(file_params->data + 1*file_params->strides[0]);
184    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
185
186   
187    // Find min and max start times of all gauges
188    start_tstep = nt + 1;
189    finish_tstep = -1;
190    for (i = 0; i < number_of_selected_stations; i++)
191    {
192        //printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
193        // printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]);   
194   
195        if ((int)cdata[i][nt + 3] < start_tstep)
196        {
197            start_tstep = (int)cdata[i][nt + 3];
198        }
199        if ((int)cdata[i][nt + 4] > finish_tstep)
200        {
201            finish_tstep = (int)cdata[i][nt + 4]; 
202        }
203    }
204
205    if ((start_tstep > nt) | (finish_tstep < 0))
206    {
207        printf("ERROR: Gauge data has incorrect start and finish times:\n");
208        printf("   start_tstep = %d, max_number_of_steps = %d\n", 
209               start_tstep, nt);
210        printf("   finish_tstep = %d, min_number_of_steps = %d\n", 
211               finish_tstep, 0);   
212           
213        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
214        return NULL;
215    }
216
217    if (start_tstep >= finish_tstep)
218    {
219        PyErr_SetString(PyExc_ValueError,
220                    "ERROR: Gauge data has non-postive_length");
221        return NULL;
222    }
223
224    num_ts = finish_tstep - start_tstep + 1;
225    dimensions[0] = number_of_selected_stations;
226    dimensions[1] = num_ts + POFFSET;
227   
228    pydata = (PyArrayObject*)PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
229    if(pydata == NULL)
230    {
231        PyErr_SetString(PyExc_ValueError, 
232          "ERROR: Memory for pydata array could not be allocated.");
233        return NULL;
234    }
235
236   
237    // Each gauge begins and ends recording at different times. When a gauge is
238    // not recording but at least one other gauge is.
239    // Pad the non-recording gauge array with zeros.
240    for (i = 0; i < number_of_selected_stations; i++)
241    {
242        time = 0;
243        for (it = 0; it < finish_tstep; it++)
244        {
245            if ((it + 1 >= start_tstep) && (it + 1 <= finish_tstep))
246            {
247                if (it + 1 > (int)cdata[i][nt + 4])
248                {
249                    // This gauge has stopped recording but others are
250            // still recording
251                    *(double*)(pydata->data + i*pydata->strides[0] 
252                                    + time*pydata->strides[1]) = 
253                        0.0;
254                }
255                else
256                {
257                    *(double*)(pydata->data + i*pydata->strides[0] 
258                                    + time*pydata->strides[1]) = 
259                        cdata[i][it];
260                }
261                time++;
262            }
263        }
264        // Pass back lat,lon,elevation
265        for (j = 0; j < POFFSET; j++)
266        {
267            *(double*)(pydata->data + i*pydata->strides[0] 
268                                + (num_ts + j)*pydata->strides[1]) = 
269                    cdata[i][nt + j];
270        }
271    }
272
273    free(weights);
274   
275    // Free filename array, but not independent Python strings
276    // FIXME(Ole): Do we need to update a reference counter in this case?
277    free(muxFileNameArray);
278   
279    for (i = 0; i < number_of_selected_stations; ++i)
280    {
281        free(cdata[i]);
282    }
283    free(cdata);
284
285    return  PyArray_Return(pydata);
286}
287
288
289float** _read_mux2(int numSrc, 
290                   char **muxFileNameArray, 
291                   float *weights, 
292                   double *params, 
293                   int *number_of_stations,
294                   long *permutation,
295                   int verbose)
296{
297    FILE *fp;
298    int nsta0, i, isrc, ista, k;
299    struct tgsrwg* mytgs0=NULL;
300    char *muxFileName;
301    int istart, istop;
302    int *fros=0, *lros=0;
303    int number_of_selected_stations;
304    float *muxData=NULL; // Suppress warning
305    long numData;
306
307    int len_sts_data;
308    float **sts_data;
309    float *temp_sts_data;
310
311    long int offset;
312
313    _read_mux2_headers(numSrc, 
314                       muxFileNameArray, 
315                       params,
316                       &fros, 
317                       &lros, 
318                       &mytgs0,
319                       &nsta0,
320                       verbose);
321    // Apply rule that an empty permutation file means 'take all stations'
322    // We could change this later by passing in None instead of the empty
323    // permutation.
324    number_of_selected_stations = *number_of_stations; 
325    if (number_of_selected_stations == 0)
326    {
327        number_of_selected_stations = nsta0; 
328   
329        // Return possibly updated number of stations
330        *number_of_stations = nsta0;     
331     
332        // Create the Identity permutation vector
333        permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
334        for (i = 0; i < number_of_selected_stations; i++)
335        {
336            permutation[i] = (long) i; 
337        }
338    }
339   
340    /*
341    printf("number_of_selected_stations = %d\n", number_of_selected_stations);
342    for (i = 0; i < number_of_selected_stations; i++) {
343      printf("permutation[%d] = %d\n", i, (int) permutation[i]);
344    }   
345    */
346   
347       
348   
349    // Make array(s) to hold demuxed data for stations given in the
350    // permutation file
351    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
352    if (sts_data == NULL)
353    {
354        printf("ERROR: Memory for sts_data could not be allocated.\n");
355        return NULL;
356    }
357
358    // For each selected station, allocate space for its data
359
360    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
361    for (i = 0; i < number_of_selected_stations; i++)
362    {
363        // Initialise sts_data to zero
364        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
365        if (sts_data[i] == NULL)
366        {
367            printf("ERROR: Memory for sts_data could not be allocated.\n");
368            return NULL;
369        }
370
371        ista = (int) permutation[i]; // Get global index into mux data
372   
373        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
374        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
375        sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
376        sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
377        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
378    }
379
380    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
381
382    /* Loop over all sources */
383    //FIXME: remove istart and istop they are not used.
384    istart = -1;
385    istop = -1;
386    for (isrc = 0; isrc < numSrc; isrc++)
387    {
388        /* Read in data block from mux2 file */
389        muxFileName = muxFileNameArray[isrc];
390        if((fp = fopen(muxFileName, "r")) == NULL)
391        {
392            fprintf(stderr, "cannot open file %s\n", muxFileName);
393            return NULL;                   
394        }
395
396        if (verbose){
397            printf("Reading mux file %s\n", muxFileName);
398        }
399
400        offset = sizeof(int) + nsta0*(sizeof(struct tgsrwg) + 2*sizeof(int));
401        fseek(fp, offset, 0);
402
403        numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
404        muxData = (float*)malloc(numData*sizeof(float));
405        fread(muxData, numData*sizeof(float), 1, fp); 
406        fclose(fp);
407
408        // loop over stations present in the permutation array
409        //     use ista with mux data
410        //     use i with the processed data to be returned         
411       
412        for(i = 0; i < number_of_selected_stations; i++)
413        {               
414   
415            ista = (int) permutation[i]; // Get global index into mux data 
416       
417            /* fill the data0 array from the mux file, and weight it */
418            fillDataArray(ista, 
419                          nsta0, 
420                          mytgs0[ista].nt, 
421                          mytgs0[ista].ig, 
422                          fros + isrc*nsta0, 
423                          lros + isrc*nsta0, 
424                          temp_sts_data, 
425                          &istart, 
426                          &istop, 
427                          muxData);
428
429            /* weight appropriately and add */
430            for(k = 0; k < mytgs0[ista].nt; k++)
431            {
432                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
433                {
434                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
435                }
436                else
437                {
438                    sts_data[i][k] = NODATA;
439                }
440            }
441        }
442        free(muxData);
443    }
444
445    free(temp_sts_data);
446    free(fros);
447    free(lros);
448    free(mytgs0);
449
450    return sts_data;
451}
452
453void _read_mux2_headers(int numSrc, 
454                        char **muxFileNameArray, 
455                        double *params,
456                        int** fros, 
457                        int** lros, 
458                        struct tgsrwg ** mytgs0,
459                        int* nsta0,
460                        int verbose)
461{
462    FILE *fp;
463    int nsta, isrc, ista;
464    struct tgsrwg *mytgs=0;
465    char *muxFileName;                                                                 
466    char susMuxFileName;
467    long numData;
468
469    /* Allocate space for the names and the weights and pointers to the data*/
470
471    /* Check that the input files have mux2 extension*/
472    susMuxFileName = 0;
473    for(isrc = 0; isrc < numSrc; isrc++)
474    { 
475        muxFileName = muxFileNameArray[isrc];
476        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, 
477                     "mux2") != 0)
478        {
479            susMuxFileName = 1;
480            break;
481        }
482    }
483
484    if(susMuxFileName)
485    {
486        printf("\n**************************************************************************\n");
487        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n"); 
488        printf("   At least one input file name does not end with mux2\n");
489        printf("   Check your results carefully!\n");
490        printf("**************************************************************************\n\n");
491    }   
492
493    if (verbose)
494    {
495        printf("Reading mux header information\n");
496    }
497
498    /* Loop over all sources, read headers and check compatibility */
499    for (isrc = 0; isrc < numSrc; isrc++)
500    {
501        muxFileName = muxFileNameArray[isrc];
502
503        /* open the mux file */
504        if((fp = fopen(muxFileName, "r")) == NULL)
505        {
506            fprintf(stderr, "cannot open file %s\n", muxFileName);
507            exit(-1); 
508        }
509       
510        if (!isrc)
511        {
512            fread(nsta0, sizeof(int), 1, fp);
513       
514            *fros = (int*)malloc(*nsta0*numSrc*sizeof(int)); 
515            *lros = (int*)malloc(*nsta0*numSrc*sizeof(int));
516     
517            *mytgs0 = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
518            mytgs = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
519
520            fread(*mytgs0, *nsta0*sizeof(struct tgsrwg), 1, fp);
521        }
522        else
523        {
524            /* check that the mux files are compatible */     
525            fread(&nsta, sizeof(int), 1, fp);
526            if(nsta != *nsta0)
527            {
528                fprintf(stderr,"%s has different number of stations to %s\n", 
529                muxFileName, 
530                muxFileNameArray[0]);
531                fclose(fp);
532                exit(-1);   
533            }
534
535            fread(mytgs, nsta*sizeof(struct tgsrwg), 1, fp); 
536           
537            for (ista = 0; ista < nsta; ista++)
538            {
539                if (mytgs[ista].dt != (*mytgs0)[ista].dt)
540                {
541                    fprintf(stderr,"%s has different sampling rate to %s\n", 
542                    muxFileName, 
543                    muxFileNameArray[0]);
544                    fclose(fp);
545                    exit(-1);           
546                }   
547                if (mytgs[ista].nt != (*mytgs0)[ista].nt)
548                {
549                    fprintf(stderr,"%s has different series length to %s\n", 
550                    muxFileName, 
551                    muxFileNameArray[0]);
552                    fclose(fp);
553                    exit(-1);           
554                }
555            }
556        }
557
558        /* Read the start and stop times for this source */
559        fread(*fros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
560        fread(*lros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
561
562        /* Compute the size of the data block for this source */
563        numData = getNumData(*fros + isrc*(*nsta0), *lros + isrc*(*nsta0), (*nsta0));
564
565        /* Sanity check */
566        if (numData < 0)
567        {
568            fprintf(stderr,"Size of data block appears to be negative!\n");
569            exit(-1);       
570        }
571
572        fclose(fp);         
573    }
574
575    params[0] = (double)*nsta0;
576    params[1] = (double)(*mytgs0)[0].dt;
577    params[2] = (double)(*mytgs0)[0].nt;
578
579        free(mytgs);
580}
581
582
583/* thomas */
584void fillDataArray(int ista, int nsta, int nt, int ig, int *nst, 
585                   int *nft, float *data, int *istart_p, 
586           int *istop_p, float *muxData)
587{
588    int it, last_it, jsta;
589    long int offset=0;
590
591
592    last_it = -1;
593    /* make arrays of starting and finishing time steps for the tide gauges */
594    /* and fill them from the file */
595
596    /* update start and stop timesteps for this gauge */
597    if (nst[ista]!= -1)
598    {
599        if(*istart_p == -1)
600        {
601            *istart_p = nst[ista];
602        }
603        else
604        {
605            *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
606        }
607    }
608    if (nft[ista] != -1)
609    {
610        if (*istop_p == -1)
611        {
612            *istop_p = nft[ista];
613        }
614        else
615        {
616            *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
617        }
618    }     
619    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
620    {
621        /* gauge never started recording, or was outside of all grids,
622    fill array with 0 */
623        for(it = 0; it < nt; it++)
624        {
625            data[it] = 0.0;
626        }
627    }   
628    else
629    {
630        for(it = 0; it < nt; it++)
631        {
632            last_it = it;
633            /* skip t record of data block */
634            offset++;
635            /* skip records from earlier tide gauges */
636            for(jsta = 0; jsta < ista; jsta++)
637                if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
638                    offset++;
639
640            /* deal with the tide gauge at hand */
641            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
642                /* gauge is recording at this time */
643            {
644                memcpy(data + it, muxData + offset, sizeof(float));
645                offset++;
646            }
647            else if (it + 1 < nst[ista])
648            {
649                /* gauge has not yet started recording */
650                data[it] = 0.0;
651            }   
652            else
653                /* gauge has finished recording */
654            {
655                data[it] = NODATA;
656                break;
657            }
658
659            /* skip records from later tide gauges */
660            for(jsta = ista + 1; jsta < nsta; jsta++)
661                if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
662                    offset++;
663        }
664
665        if(last_it < nt - 1)
666            /* the loop was exited early because the gauge had
667        finished recording */
668            for(it = last_it+1; it < nt; it++)
669                data[it] = NODATA;
670    }
671} 
672
673char isdata(float x)
674{
675    //char value;
676    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
677    {
678       return 0;
679    }
680    else
681    {
682        return 1; 
683    }
684}
685
686
687long getNumData(const int *fros, const int *lros, const int nsta)
688/* calculates the number of data in the data block of a mux file */
689/* based on the first and last recorded output steps for each gauge */ 
690{
691    int ista, last_output_step;
692    long numData = 0;
693
694    last_output_step = 0;   
695    for(ista = 0; ista < nsta; ista++)
696        if(*(fros + ista) != -1)
697        {
698            numData += *(lros + ista) - *(fros + ista) + 1;
699            last_output_step = (last_output_step < *(lros+ista) ? 
700                            *(lros+ista):last_output_step);
701        }   
702        numData += last_output_step*nsta; /* these are the t records */
703        return numData;
704}   
705
706
707
708//-------------------------------
709// Method table for python module
710//-------------------------------
711static struct PyMethodDef MethodTable[] = {
712    {"read_mux2", read_mux2, METH_VARARGS, "Print out"},
713    {NULL, NULL}
714};
715
716// Module initialisation
717void initurs_ext(void){
718    Py_InitModule("urs_ext", MethodTable);
719
720    import_array(); // Necessary for handling of NumPY structures
721}
Note: See TracBrowser for help on using the repository browser.