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

Last change on this file since 5529 was 5529, checked in by ole, 14 years ago

Verified that pointer fix works on Linux. All tests pass and boxing day example appears to work (inasmuch as it doesn't segfault anymore).
I have not looked for memory leaks.

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