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

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

Nariman and Ole fixed 64 bit issue with urs_ext.c.
NumSrc? was declared as a long, but wasn't casted to int in a loop where the iterator i was an int.
Solution: NumSrc? now declared as an int.
32 bit compatibility yet to be verified.

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