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

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

Fixed urs_ext with permutation vector.

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