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

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

Checked in debug information for urs_ext.c

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