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

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

A bit more cleanup in urs2sts - all tests pass, but I think we need a few more.

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