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

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

Changed time vector in urs2sts to start at 0.0 insteade of delta_t

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