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

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

Fixed urs2sts such that start and stop times are updated across mux sources

File size: 23.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
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    //FIXME: The params can be removed after the python interface is modified.
345    params[0] = (double)total_number_of_stations;
346    params[1] = (double)delta_t;
347    params[2] = (double)number_of_time_steps;
348   
349    // Apply rule that an empty permutation file means 'take all stations'
350    // We could change this later by passing in None instead of the empty
351    // permutation.
352    number_of_selected_stations = *number_of_stations; 
353    if (number_of_selected_stations == 0)
354    {
355        number_of_selected_stations = total_number_of_stations; 
356   
357        // Return possibly updated number of stations
358        *number_of_stations = total_number_of_stations;     
359     
360        // Create the Identity permutation vector
361        permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
362        for (i = 0; i < number_of_selected_stations; i++)
363        {
364            permutation[i] = (long) i; 
365        }
366    }
367   
368    // Make array(s) to hold demuxed data for stations given in the
369    // permutation file
370    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
371    if (sts_data == NULL)
372    {
373        printf("ERROR: Memory for sts_data could not be allocated.\n");
374        return NULL;
375    }
376
377    // For each selected station, allocate space for its data
378    len_sts_data = number_of_time_steps + POFFSET; // Max length of each timeseries?
379    for (i = 0; i < number_of_selected_stations; i++)
380    {
381        // Initialise sts_data to zero
382        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
383        if (sts_data[i] == NULL)
384        {
385            printf("ERROR: Memory for sts_data could not be allocated.\n");
386            return NULL;
387        }
388    }
389
390    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
391
392    muxData = (float*)malloc(numDataMax*sizeof(float));
393   
394    /* Loop over all sources */
395    for (isrc = 0; isrc < numSrc; isrc++)
396    {
397   
398        // Shorthands to local memory
399        fros_per_source = (int*) fros + isrc*total_number_of_stations; 
400        lros_per_source = (int*) lros + isrc*total_number_of_stations;     
401           
402   
403   
404        /* Read in data block from mux2 file */
405        muxFileName = muxFileNameArray[isrc];
406        if((fp = fopen(muxFileName, "r")) == NULL)
407        {
408            fprintf(stderr, "cannot open file %s\n", muxFileName);
409            return NULL;                   
410        }
411
412        if (verbose){
413            printf("Reading mux file %s\n", muxFileName);
414        }
415
416        offset = sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
417        fseek(fp, offset, 0);
418
419        numData = getNumData(fros_per_source, 
420                             lros_per_source, 
421                             total_number_of_stations);
422                             
423        fread(muxData, numData*sizeof(float), 1, fp); 
424        fclose(fp);
425
426        // loop over stations present in the permutation array
427        //     use ista with mux data
428        //     use i with the processed data to be returned         
429       
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           
464            /*
465            printf("station =  %d, source = %d, fros = %f, lros = %f\n",
466                   ista,
467                   isrc,
468                   (float) fros_per_source[ista],
469                   (float) lros_per_source[ista]);                 
470            */
471                   
472                   
473                   
474            if (isrc == 0) {
475                // Assign values for first source
476                sts_data[i][N] = (float)mytgs0[ista].geolat;
477                sts_data[i][N+1] = (float)mytgs0[ista].geolon;
478                sts_data[i][N+2] = (float)mytgs0[ista].z;
479                sts_data[i][N+3] = (float)fros_per_source[ista];
480                sts_data[i][N+4] = (float)lros_per_source[ista];
481            } else {
482                // Update first and last timesteps for subsequent sources
483                if (sts_data[i][N+3] > (float)fros_per_source[ista]) {         
484                    if (verbose) {
485                        printf("Adjusting start time for station %d and source %d",
486                               ista, isrc);
487                        printf(" from %f to %f\n", sts_data[i][N+3], (float)fros_per_source[ista]); 
488                    }
489                    sts_data[i][N+3] = (float)fros_per_source[ista];
490                }
491               
492               
493                if (sts_data[i][N+4] < (float)lros_per_source[ista]) {         
494                    if (verbose) {
495                        printf("Adjusting end time for station %d and source %d",
496                               ista, isrc);
497                        printf(" from %f to %f\n", sts_data[i][N+4], (float)lros_per_source[ista]); 
498                    }
499                    sts_data[i][N+4] = (float)lros_per_source[ista];
500                }               
501            }
502        }
503    }
504
505    free(muxData);
506    free(temp_sts_data);
507    free(fros);
508    free(lros);
509    free(mytgs0);
510
511    return sts_data;
512}
513
514/////////////////////////////////////////////////////////////////////////
515//Python gateways
516PyObject *read_mux2(PyObject *self, PyObject *args)
517{
518    /*Read in mux 2 file
519
520    Python call:
521    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
522
523    NOTE:
524    A Python int is equivalent to a C long
525    (this becomes really important on 64 bit architectures)
526   
527    A Python double corresponds to a C double
528    */
529   
530    PyObject *filenames;
531    PyArrayObject *pyweights;
532    PyArrayObject *file_params;
533    PyArrayObject *permutation;  // Ordering of selected stations   
534    PyArrayObject *pydata;
535    PyObject *fname;
536
537    char **muxFileNameArray;
538    float **cdata;
539    float *weights;
540    int dimensions[2];
541    int numSrc;
542    int verbose;
543    int total_number_of_stations;
544    int number_of_selected_stations;   
545    int nt;
546    double dt;
547    int i;
548    int j;
549    int start_tstep;
550    int finish_tstep;
551    int it;
552    int time;
553    int num_ts;
554   
555    // Convert Python arguments to C
556    if (!PyArg_ParseTuple(args, "iOOOOi",
557              &numSrc, &filenames, &pyweights, &file_params, 
558              &permutation, &verbose)) 
559    {
560            PyErr_SetString(PyExc_RuntimeError, 
561                "Input arguments to read_mux2 failed");
562            return NULL;
563    }
564
565    if(!PyList_Check(filenames)) 
566    {
567        PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list");
568        return NULL;
569    }
570
571    if(PyList_Size(filenames) == 0)
572    {
573        PyErr_SetString(PyExc_ValueError, "empty lists not allowed");
574        return NULL;
575    }
576
577    if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) 
578    {
579        PyErr_SetString(PyExc_ValueError,
580            "pyweights must be one-dimensional and of type double");
581        return NULL; 
582    }
583
584    if(PyList_Size(filenames) != pyweights->dimensions[0])
585    {
586        PyErr_SetString(PyExc_ValueError, 
587            "Must specify one weight for each filename");
588        return NULL;
589    }
590
591    muxFileNameArray = (char**)malloc(numSrc*sizeof(char*));
592    if (muxFileNameArray == NULL) 
593    {
594        PyErr_SetString(PyExc_ValueError, 
595      "ERROR: Memory for muxFileNameArray could not be allocated.");
596        return NULL;
597    }
598
599    for (i = 0; i < numSrc; i++)
600    {
601
602        fname = PyList_GetItem(filenames, i);
603        if (!fname) 
604        {
605            PyErr_SetString(PyExc_ValueError, "filename not a string");
606            return NULL;
607        }       
608
609        muxFileNameArray[i] = PyString_AsString(fname);
610        if (muxFileNameArray[i] == NULL) 
611        {
612            PyErr_SetString(PyExc_ValueError, 
613          "ERROR: Memory for muxFileNameArray could not be allocated.\n");
614            return NULL;
615        }
616    }
617
618    if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) 
619    {
620        PyErr_SetString(PyExc_ValueError,
621          "file_params must be one-dimensional and of type double");
622        return NULL; 
623    }
624
625
626    // Create array for weights which are passed to read_mux2
627    weights = (float*) malloc(numSrc*sizeof(float));
628    for (i = 0; i < numSrc; i++)
629    {
630        weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0]));
631    }
632   
633    // Desired number of stations
634    number_of_selected_stations = (int) permutation->dimensions[0];
635
636    // Read in mux2 data from file
637    cdata = _read_mux2(numSrc, 
638                       muxFileNameArray, 
639                       weights, 
640                       (double*)file_params->data,
641                       &number_of_selected_stations, 
642                       (long*) permutation->data,
643                       verbose);
644
645    if (!cdata) 
646    {
647        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
648        return NULL;
649    }       
650               
651               
652    // Allocate space for return vector
653    total_number_of_stations = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
654    dt = *(double*)(file_params->data + 1*file_params->strides[0]);
655    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
656
657   
658    // Find min and max start times of all gauges
659    start_tstep = nt + 1;
660    finish_tstep = -1;
661    for (i = 0; i < number_of_selected_stations; i++)
662    {
663        //printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
664        // printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]);   
665   
666        if ((int)cdata[i][nt + 3] < start_tstep)
667        {
668            start_tstep = (int)cdata[i][nt + 3];
669        }
670        if ((int)cdata[i][nt + 4] > finish_tstep)
671        {
672            finish_tstep = (int)cdata[i][nt + 4]; 
673        }
674    }
675
676    if ((start_tstep > nt) | (finish_tstep < 0))
677    {
678        printf("ERROR: Gauge data has incorrect start and finish times:\n");
679        printf("   start_tstep = %d, max_number_of_steps = %d\n", 
680               start_tstep, nt);
681        printf("   finish_tstep = %d, min_number_of_steps = %d\n", 
682               finish_tstep, 0);   
683           
684        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
685        return NULL;
686    }
687
688    if (start_tstep >= finish_tstep)
689    {
690        PyErr_SetString(PyExc_ValueError,
691                    "ERROR: Gauge data has non-postive_length");
692        return NULL;
693    }
694
695    num_ts = finish_tstep - start_tstep + 1;
696    dimensions[0] = number_of_selected_stations;
697    dimensions[1] = num_ts + POFFSET;
698   
699    pydata = (PyArrayObject*)PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
700    if(pydata == NULL)
701    {
702        PyErr_SetString(PyExc_ValueError, 
703          "ERROR: Memory for pydata array could not be allocated.");
704        return NULL;
705    }
706
707   
708    // Each gauge begins and ends recording at different times. When a gauge is
709    // not recording but at least one other gauge is.
710    // Pad the non-recording gauge array with zeros.
711    for (i = 0; i < number_of_selected_stations; i++)
712    {
713        time = 0;
714        for (it = 0; it < finish_tstep; it++)
715        {
716            if ((it + 1 >= start_tstep) && (it + 1 <= finish_tstep))
717            {
718                if (it + 1 > (int)cdata[i][nt + 4])
719                {
720                    // This gauge has stopped recording but others are
721            // still recording
722                    *(double*)(pydata->data + i*pydata->strides[0] 
723                                    + time*pydata->strides[1]) = 
724                        0.0;
725                }
726                else
727                {
728                    *(double*)(pydata->data + i*pydata->strides[0] 
729                                    + time*pydata->strides[1]) = 
730                        cdata[i][it];
731                }
732                time++;
733            }
734        }
735        // Pass back lat,lon,elevation
736        for (j = 0; j < POFFSET; j++)
737        {
738            *(double*)(pydata->data + i*pydata->strides[0] 
739                                + (num_ts + j)*pydata->strides[1]) = 
740                    cdata[i][nt + j];
741        }
742    }
743
744    free(weights);
745   
746    // Free filename array, but not independent Python strings
747    // FIXME(Ole): Do we need to update a reference counter in this case?
748    free(muxFileNameArray);
749   
750    for (i = 0; i < number_of_selected_stations; ++i)
751    {
752        free(cdata[i]);
753    }
754    free(cdata);
755
756    return  PyArray_Return(pydata);
757}
758
759//-------------------------------
760// Method table for python module
761//-------------------------------
762static struct PyMethodDef MethodTable[] = {
763    {"read_mux2", read_mux2, METH_VARARGS, "Print out"},
764    {NULL, NULL}
765};
766
767// Module initialisation
768void initurs_ext(void){
769    Py_InitModule("urs_ext", MethodTable);
770
771    import_array(); // Necessary for handling of NumPY structures
772}
Note: See TracBrowser for help on using the repository browser.