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

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

Work on permutations in urs_ext.c
It does not work yet, but the array is being passed in and the algorithm
has been described.

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