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

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

Attempt to fix pointer bug in urs_ext.c

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