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

Last change on this file since 5542 was 5542, checked in by ole, 11 years ago

Suppressed muxdata warning and did some formatting

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