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

Last change on this file since 5516 was 5516, checked in by nariman, 14 years ago

Elimination of memory leaks
Code reformatting
Expansion of first for loop to include first muxfile.

File size: 17.2 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 *, int *, int *, float *);
28long getNumData(const int *fros, const int *lros, const int nsta);
29char isdata(float);
30float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose);
31
32PyObject *read_mux2(PyObject *self, PyObject *args){
33    /*Read in mux 2 file
34
35    Python call:
36    read_mux2(numSrc,filenames,weights,file_params,verbose)
37
38    NOTE:
39    A Python int is equivalent to a C long
40    A Python double corresponds to a C double
41    */
42    PyObject *filenames;
43    PyArrayObject *pyweights;
44    PyArrayObject *file_params;
45    PyArrayObject *pydata;
46    PyObject *fname;
47
48    char **muxFileNameArray;
49    float **cdata;
50    float *weights;
51    int dimensions[2];
52    long numSrc;
53    long verbose;
54    int nsta0;
55    int nt;
56    double dt;
57    int i;
58    int j;
59    int start_tstep;
60    int finish_tstep;
61    int it;
62    int time;
63    int num_ts;
64
65    // Convert Python arguments to C
66    if (!PyArg_ParseTuple(args, "iOOOi",
67        &numSrc, &filenames, &pyweights, &file_params, &verbose)) 
68    {
69            PyErr_SetString(PyExc_RuntimeError, 
70                "Input arguments to read_mux2 failed");
71            return NULL;
72    }
73
74    if(!PyList_Check(filenames)) 
75    {
76        PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list");
77        return NULL;
78    }
79
80    if(PyList_Size(filenames) == 0)
81    {
82        PyErr_SetString(PyExc_ValueError, "empty lists not allowed");
83        return NULL;
84    }
85
86    if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) 
87    {
88        PyErr_SetString(PyExc_ValueError,
89            "pyweights must be one-dimensional and of type double");
90        return NULL; 
91    }
92
93    if(PyList_Size(filenames) != pyweights->dimensions[0])
94    {
95        PyErr_SetString(PyExc_ValueError, "Must specify one weight for each filename");
96        return NULL;
97    }
98
99    muxFileNameArray = (char**)malloc((int)numSrc*sizeof(char*));
100    if (muxFileNameArray == NULL) 
101    {
102        printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
103        exit(-1);
104    }
105
106    for (i = 0; i < numSrc; i++)
107    {
108        muxFileNameArray[i] = (char*)malloc((MAX_FILE_NAME_LENGTH + 1)*sizeof(char));
109       
110        if (muxFileNameArray[i] == NULL) 
111        {
112            printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
113            exit(-1);
114        }
115
116        fname = PyList_GetItem(filenames, i);
117        if (!PyString_Check(fname)) 
118        {
119            PyErr_SetString(PyExc_ValueError, "filename not a string");
120        }
121
122        muxFileNameArray[i] = PyString_AsString(fname);
123    }
124
125    if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) 
126    {
127        PyErr_SetString(PyExc_ValueError,
128            "file_params must be one-dimensional and of type double");
129        return NULL; 
130    }
131
132    // Create array for weights which are passed to read_mux2
133    weights = (float*) malloc((int)numSrc*sizeof(float));
134    for (i = 0; i < (int)numSrc; i++)
135    {
136        weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0]));
137    }
138   
139    // Read in mux2 data from file
140    cdata = _read_mux2((int)numSrc, muxFileNameArray, weights, (double*)file_params->data, (int)verbose);
141
142    // Allocate space for return vector
143    nsta0 = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
144    dt = *(double*)(file_params->data + 1*file_params->strides[0]);
145    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
146
147    // Find min and max start times of all gauges
148    start_tstep = nt + 1;
149    finish_tstep = -1;
150    for (i = 0; i < nsta0; i++)
151    {
152        if ((int)cdata[i][nt + 3] < start_tstep)
153        {
154            start_tstep = (int)cdata[i][nt + 3];
155        }
156        if ((int)cdata[i][nt + 4] > finish_tstep)
157        {
158            finish_tstep = (int)cdata[i][nt + 4]; 
159        }
160    }
161
162    if ((start_tstep > nt) | (finish_tstep < 0))
163    {
164        printf("ERROR: Gauge data has incorrect start and finsh times\n");
165        return NULL;
166    }
167
168    if (start_tstep >= finish_tstep)
169    {
170        printf("ERROR: Gauge data has non-postive_length\n");
171        return NULL;
172    }
173
174    num_ts = finish_tstep - start_tstep + 1;
175    dimensions[0] = nsta0;
176    dimensions[1] = num_ts + POFFSET;
177    pydata = (PyArrayObject*)PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
178    if(pydata == NULL)
179    {
180        printf("ERROR: Memory for pydata array could not be allocated.\n");
181        return NULL;
182    }
183
184    // Each gauge begins and ends recording at different times. When a gauge is
185    // not recording but at least one other gauge is. Pad the non-recording gauge
186    // array with zeros.
187    for (i = 0; i < nsta0; i++)
188    {
189        time = 0;
190        for (it = 0; it < finish_tstep; it++)
191        {
192            if ((it + 1 >= start_tstep) && (it + 1 <= finish_tstep))
193            {
194                if (it + 1 > (int)cdata[i][nt + 4])
195                {
196                    //This gauge has stopped recording but others are still recording
197                    *(double*)(pydata->data + i*pydata->strides[0] + time*pydata->strides[1]) = 0.0;
198                }
199                else
200                {
201                    *(double*)(pydata->data + i*pydata->strides[0] + time*pydata->strides[1]) = cdata[i][it];
202                }
203                time++;
204            }
205        }
206        // Pass back lat,lon,elevation
207        for (j = 0; j < POFFSET; j++)
208        {
209            *(double*)(pydata->data + i*pydata->strides[0] + (num_ts + j)*pydata->strides[1]) = cdata[i][nt + j];
210        }
211    }
212
213    free(weights);
214   
215    for (i = 0; i < numSrc; ++i)
216    {
217        free(muxFileNameArray[i]);
218    }
219    free(muxFileNameArray);
220   
221    for (i = 0; i < nsta0; ++i)
222    {
223        free(cdata[i]);
224    }
225    free(cdata);
226
227    return  PyArray_Return(pydata);
228}
229
230float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)
231{
232    FILE *fp;
233    int nsta, nsta0, i, isrc, ista;
234    struct tgsrwg *mytgs=0, *mytgs0=0;
235    char *muxFileName;                                                                 
236    int istart, istop;
237    int *fros=0, *lros=0;
238    char susMuxFileName;
239    float *muxData;
240    long numData;
241
242    int len_sts_data;
243    float **sts_data;
244    float *temp_sts_data;
245
246    long int offset;
247
248    /* Allocate space for the names and the weights and pointers to the data*/   
249
250    /* Check that the input files have mux2 extension*/
251    susMuxFileName = 0;
252    for(isrc = 0; isrc < numSrc; isrc++)
253    { 
254        muxFileName = muxFileNameArray[isrc];
255        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, "mux2") != 0)
256        {
257            susMuxFileName = 1;
258            break;
259        }
260    }
261
262    if(susMuxFileName)
263    {
264        printf("\n**************************************************************************\n");
265        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n"); 
266        printf("   At least one input file name does not end with mux2\n");
267        printf("   Check your results carefully!\n");
268        printf("**************************************************************************\n\n");
269    }   
270
271    if (verbose)
272    {
273        printf("Reading mux header information\n");
274    }
275
276    /* loop over remaining sources, check compatibility, and read them into
277    *muxData */ 
278        for (isrc = 0; isrc < numSrc; isrc++)
279    {
280        muxFileName = muxFileNameArray[isrc];
281
282                /* open the mux file */
283        if((fp = fopen(muxFileName, "r")) == NULL)
284        {
285            fprintf(stderr, "cannot open file %s\n", muxFileName);
286            exit(-1); 
287        }
288       
289        if (!isrc)
290        {
291                        fread(&nsta0, sizeof(int), 1, fp);
292                       
293                        fros = (int*)malloc(nsta0*numSrc*sizeof(int)); 
294                        lros = (int*)malloc(nsta0*numSrc*sizeof(int));
295     
296                        mytgs0 = (struct tgsrwg*)malloc(nsta0*sizeof(struct tgsrwg));
297                        mytgs = (struct tgsrwg*)malloc(nsta0*sizeof(struct tgsrwg));
298
299                        fread(mytgs0, nsta0*sizeof(struct tgsrwg), 1, fp);
300        }
301        else
302        {
303            /* check that the mux files are compatible */     
304            fread(&nsta, sizeof(int), 1, fp);
305            if(nsta != nsta0)
306            {
307                fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
308                fclose(fp);
309                exit(-1);   
310            }
311
312            fread(mytgs, nsta*sizeof(struct tgsrwg), 1, fp); 
313           
314            for (ista = 0; ista < nsta; ista++)
315            {
316                if (mytgs[ista].dt != mytgs0[ista].dt)
317                {
318                    fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
319                    fclose(fp);
320                    exit(-1);                 
321                }   
322                if (mytgs[ista].nt != mytgs0[ista].nt)
323                {
324                    fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
325                    fclose(fp);
326                    exit(-1);                 
327                }
328            }
329        }
330
331        /* read the start and stop times for this source */
332        fread(fros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
333        fread(lros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
334
335        /* compute the size of the data block for this source */
336        numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
337
338        /* Burbidge: Added a sanity check here */
339        if (numData < 0)
340        {
341            fprintf(stderr,"Size of data block appears to be negative!\n");
342            exit(-1);
343        }
344
345        fclose(fp);         
346    }
347
348    params[0] = (double)nsta0;
349    params[1] = (double)mytgs0[0].dt;
350    params[2] = (double)mytgs0[0].nt;
351
352    /* make array(s) to hold the demuxed data */
353    //FIXME: This can be reduced to only contain stations in order file
354    sts_data = (float**)malloc(nsta0*sizeof(float *));
355    if (sts_data == NULL)
356    {
357        printf("ERROR: Memory for sts_data could not be allocated.\n");
358        exit(-1);
359    }
360
361    len_sts_data = mytgs0[0].nt + POFFSET;
362    for (ista = 0; ista < nsta0; ista++)
363    {
364        // Initialise sts_data to zero
365        sts_data[ista] = (float*)calloc(len_sts_data, sizeof(float));
366        if (sts_data[ista] == NULL)
367        {
368            printf("ERROR: Memory for sts_data could not be allocated.\n");
369            exit(-1);
370        }
371
372        sts_data[ista][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
373        sts_data[ista][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
374        sts_data[ista][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
375        sts_data[ista][mytgs0[0].nt + 3] = (float)fros[ista];
376        sts_data[ista][mytgs0[0].nt + 4] = (float)lros[ista];
377    }
378
379    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
380
381    /* Loop over all sources */
382    //FIXME: remove istart and istop they are not used.
383    istart = -1;
384    istop = -1;
385    for (isrc = 0; isrc < numSrc; isrc++)
386    {
387        /* Read in data block from mux2 file */
388        muxFileName = muxFileNameArray[isrc];
389        if((fp = fopen(muxFileName, "r")) == NULL)
390        {
391            //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
392            fprintf(stderr, "cannot open file %s\n", muxFileName);
393            exit(-1); 
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 */
409        for(ista = 0; ista < nsta0; ista++)
410        {               
411            /* fill the data0 array from the mux file, and weight it */
412            fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros + isrc*nsta0, lros + isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
413
414            /* weight appropriately and add */
415            for(i = 0; i < mytgs0[ista].nt; i++)
416            {
417                if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i]))
418                {
419                    sts_data[ista][i] += temp_sts_data[i] * weights[isrc];
420                    //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]);
421                }
422                else
423                {
424                    sts_data[ista][i] = NODATA;
425                }
426            }
427        }
428    }
429
430    free(muxData);
431    free(temp_sts_data);
432    free(fros);
433    free(lros);
434    free(mytgs0);
435    free(mytgs);
436
437    return sts_data;
438}
439
440
441/* thomas */
442void fillDataArray(int ista, int nsta, int nt, int ig, int *nst, int *nft, float *data, int *istart_p, int *istop_p, float *muxData)
443{
444    int it, last_it, jsta;
445    long int offset=0;
446
447
448    last_it = -1;
449    /* make arrays of starting and finishing time steps for the tide gauges */
450    /* and fill them from the file */
451
452    /* update start and stop timesteps for this gauge */
453    if (nst[ista]!= -1)
454    {
455        if(*istart_p == -1)
456        {
457            *istart_p = nst[ista];
458        }
459        else
460        {
461            *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
462        }
463    }
464    if (nft[ista] != -1)
465    {
466        if (*istop_p == -1)
467        {
468            *istop_p = nft[ista];
469        }
470        else
471        {
472            *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
473        }
474    }     
475    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
476    {
477        /* gauge never started recording, or was outside of all grids, fill array with 0 */
478        for(it = 0; it < nt; it++)
479        {
480            data[it] = 0.0;
481        }
482    }   
483    else
484    {
485        for(it = 0; it < nt; it++)
486        {
487            last_it = it;
488            /* skip t record of data block */
489            offset++;
490            /* skip records from earlier tide gauges */
491            for(jsta = 0; jsta < ista; jsta++)
492                if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
493                    offset++;
494
495            /* deal with the tide gauge at hand */
496            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
497                /* gauge is recording at this time */
498            {
499                memcpy(data + it, muxData + offset, sizeof(float));
500                offset++;
501            }
502            else if (it + 1 < nst[ista])
503            {
504                /* gauge has not yet started recording */
505                data[it] = 0.0;
506            }   
507            else
508                /* gauge has finished recording */                                           
509            {
510                data[it] = NODATA;
511                break;
512            }
513
514            /* skip records from later tide gauges */
515            for(jsta = ista + 1; jsta < nsta; jsta++)
516                if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
517                    offset++;
518        }
519
520        if(last_it < nt - 1)
521            /* the loop was exited early because the gauge had finished recording */
522            for(it = last_it+1; it < nt; it++)
523                data[it] = NODATA;
524    }
525} 
526
527char isdata(float x)
528{
529    //char value;
530    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
531    {
532       return 0;
533    }
534    else
535    {
536        return 1; 
537    }
538}
539
540
541long getNumData(const int *fros, const int *lros, const int nsta)
542/* calculates the number of data in the data block of a mux file */
543/* based on the first and last recorded output steps for each gauge */ 
544{
545    int ista, last_output_step;
546    long numData = 0;
547
548    last_output_step = 0;   
549    for(ista = 0; ista < nsta; ista++)
550        if(*(fros + ista) != -1)
551        {
552            numData += *(lros + ista) - *(fros + ista) + 1;
553            last_output_step = (last_output_step < *(lros+ista) ? *(lros+ista):last_output_step);
554        }   
555        numData += last_output_step*nsta; /* these are the t records */
556        return numData;
557}   
558
559
560
561//-------------------------------
562// Method table for python module
563//-------------------------------
564static struct PyMethodDef MethodTable[] = {
565    {"read_mux2", read_mux2, METH_VARARGS, "Print out"},
566    {NULL, NULL}
567};
568
569// Module initialisation
570void initurs_ext(void){
571    Py_InitModule("urs_ext", MethodTable);
572
573    import_array(); // Necessary for handling of NumPY structures
574}
Note: See TracBrowser for help on using the repository browser.