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

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

Resolved 64 bit issue for permutation in urs_ext.c

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