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

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

Creation of static work arrays. Cleanup of code.

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