source: branches/numpy/anuga/shallow_water/urs_ext.c @ 6780

Last change on this file since 6780 was 6780, checked in by rwilson, 15 years ago

Changes to make legacy Numeric API work in numpy.

File size: 25.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#include "Python.h"
7#include "numpy/arrayobject.h"
8#include "structure.h"
9#include "math.h"
10#include <stdio.h>
11#include <string.h>
12#include <errno.h>
13#include <float.h>
14#include <time.h>
15
16#include "numpy_shim.h"
17
18#define MAX_FILE_NAME_LENGTH 128
19#define NODATA 99.0
20#define EPSILON  0.00001
21
22#define DEBUG 0
23
24#define POFFSET 5 //Number of site_params
25
26static int *fros=NULL;  // First recorded output step
27static int *lros=NULL;  // Last recorded output step
28static struct tgsrwg* mytgs0=NULL;
29
30static long numDataMax=0;
31
32
33/*The MUX file format
34
35
36 */
37
38
39
40/////////////////////////////////////////////////////////////////////////
41//Auxiliary functions
42void fillDataArray(int ista, int total_number_of_stations, int nt, int ig, int *nst, 
43                   int *nft, float *data, int *istart_p, 
44                   int *istop_p, float *muxData)
45{
46    int it, last_it, jsta;
47    long int offset=0;
48
49
50    last_it = -1;
51    /* Make arrays of starting and finishing time steps for the tide gauges */
52    /* and fill them from the file */
53
54    /* Update start and stop timesteps for this gauge */
55    if (nst[ista]!= -1)
56    {
57        if(*istart_p == -1)
58        {
59            *istart_p = nst[ista];
60        }
61        else
62        {
63            *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
64        }
65    }
66   
67    if (nft[ista] != -1)
68    {
69        if (*istop_p == -1)
70        {
71            *istop_p = nft[ista];
72        }
73        else
74        {
75            *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
76        }
77    }     
78   
79    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
80    {
81        /* gauge never started recording, or was outside of all grids,
82        fill array with 0 */
83        for(it = 0; it < nt; it++)
84        {
85            data[it] = 0.0;
86        }
87    }   
88    else
89    {
90        for(it = 0; it < nt; it++)
91        {
92            last_it = it;
93            /* skip t record of data block */
94            offset++;
95            /* skip records from earlier tide gauges */
96            for(jsta = 0; jsta < ista; jsta++)
97                if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
98                    offset++;
99
100            /* deal with the tide gauge at hand */
101            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
102            {
103                /* gauge is recording at this time */
104                memcpy(data + it, muxData + offset, sizeof(float));
105               
106                //printf("%d: muxdata=%f\n", it, muxData[offset]);                     
107                //printf("data[%d]=%f, offset=%d\n", it, data[it], offset);         
108                offset++;
109            }
110            else if (it + 1 < nst[ista])
111            {
112                /* gauge has not yet started recording */
113                data[it] = 0.0;
114            }   
115            else
116                /* gauge has finished recording */
117            {
118                data[it] = NODATA;
119                break;
120            }
121
122            /* skip records from later tide gauges */
123            for(jsta = ista + 1; jsta < total_number_of_stations; jsta++)
124                if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
125                    offset++;
126        }
127
128        if(last_it < nt - 1)
129            /* the loop was exited early because the gauge had
130            finished recording */
131            for(it = last_it+1; it < nt; it++)
132                data[it] = NODATA;
133    }
134} 
135
136
137char isdata(float x)
138{
139    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
140    {
141      return 0;
142    }
143    else
144    {
145      return 1; 
146    }
147}
148
149
150long getNumData(const int *fros, const int *lros, const int total_number_of_stations)
151/* calculates the number of data in the data block of a mux file */
152/* based on the first and last recorded output steps for each gauge */ 
153{
154    int ista, last_output_step;
155    long numData = 0;
156
157    last_output_step = 0;   
158    for(ista = 0; ista < total_number_of_stations; ista++)
159        if(*(fros + ista) != -1)
160        {
161            numData += *(lros + ista) - *(fros + ista) + 1;
162            last_output_step = (last_output_step < *(lros+ista) ? 
163                            *(lros+ista):last_output_step);
164        }   
165        numData += last_output_step*total_number_of_stations; /* these are the t records */
166        return numData;
167}
168
169/////////////////////////////////////////////////////////////////////////
170//Internal Functions
171int _read_mux2_headers(int numSrc, 
172                       char **muxFileNameArray, 
173                       int* total_number_of_stations,
174                       int* number_of_time_steps,
175                       double* delta_t,
176                       //long* numDataMax,
177                       int verbose)
178{
179    FILE *fp;
180    int numsta, i, j;
181    struct tgsrwg *mytgs=0;
182    char *muxFileName;                                                                 
183    char susMuxFileName;
184    long numData;
185    size_t elements_read; // fread return value
186    int block_size;
187
188    /* Allocate space for the names and the weights and pointers to the data*/
189
190    /* Check that the input files have mux2 extension*/
191    susMuxFileName = 0;
192    for(i = 0; i < numSrc; i++)
193    { 
194        muxFileName = muxFileNameArray[i];
195        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, 
196                     "mux2") != 0)
197        {
198            susMuxFileName = 1;
199            break;
200        }
201    }
202
203    if(susMuxFileName)
204    {
205        printf("\n**************************************************************************\n");
206        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n"); 
207        printf("   At least one input file name does not end with mux2\n");
208        printf("   Check your results carefully!\n");
209        printf("**************************************************************************\n\n");
210    }   
211
212    if (verbose)
213    {
214        printf("Reading mux header information\n");
215    }
216
217    // Loop over all sources, read headers and check compatibility
218    for (i = 0; i < numSrc; i++)
219    {
220        muxFileName = muxFileNameArray[i];
221
222        // Open the mux file
223        if((fp = fopen(muxFileName, "r")) == NULL)
224        {
225            char *err_msg = strerror(errno);
226
227            fprintf(stderr, "cannot open file '%s': %s\n", muxFileName, err_msg);
228            return -1; 
229        }
230       
231        if (!i)
232        {
233            elements_read = fread(total_number_of_stations, sizeof(int), 1, fp);
234            if ((int) elements_read == 0 && ferror(fp)){
235              fprintf(stderr, "Error reading total number of stations\n");
236              return -2;
237            }
238       
239            fros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int)); 
240            lros = (int*) malloc(*total_number_of_stations*numSrc*sizeof(int));
241     
242            mytgs0 = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg));
243            mytgs = (struct tgsrwg*) malloc(*total_number_of_stations*sizeof(struct tgsrwg));
244
245            block_size = *total_number_of_stations*sizeof(struct tgsrwg);
246            elements_read = fread(mytgs0, block_size , 1, fp);
247            if ((int) elements_read == 0 && ferror(fp)){
248              fprintf(stderr, "Error reading mytgs0\n");
249              return -2;
250            }
251        }
252        else
253        {
254            // Check that the mux files are compatible
255            elements_read = fread(&numsta, sizeof(int), 1, fp);
256            if ((int) elements_read == 0 && ferror(fp)){
257              fprintf(stderr, "Error reading numsta\n");
258              return -2;
259            }
260           
261            if(numsta != *total_number_of_stations)
262            {
263                fprintf(stderr,"%s has different number of stations to %s\n", 
264                muxFileName, 
265                muxFileNameArray[0]);
266                fclose(fp);
267                return -1;   
268            }
269
270            block_size = numsta*sizeof(struct tgsrwg);
271            elements_read = fread(mytgs, block_size, 1, fp); 
272            if ((int) elements_read == 0 && ferror(fp)){
273              fprintf(stderr, "Error reading mgtgs\n");
274              return -2;
275            }       
276           
277           
278            for (j = 0; j < numsta; j++)
279            {
280                if (mytgs[j].dt != mytgs0[j].dt)
281                {
282                    fprintf(stderr, "%s has different sampling rate to %s\n", 
283                    muxFileName, 
284                    muxFileNameArray[0]);
285                    fclose(fp);
286                    return -1;           
287                }   
288                if (mytgs[j].nt != mytgs0[j].nt)
289                {
290                    fprintf(stderr, "%s has different series length to %s\n", 
291                    muxFileName, 
292                    muxFileNameArray[0]);
293                    fclose(fp);
294                    return -1;           
295                }
296
297                if (mytgs[j].nt != mytgs0[0].nt)
298                {
299                    printf("Station 0 has different series length to Station %d\n", j); 
300                }
301            }
302        }
303
304        /* Read the start and stop times for this source */
305        elements_read = fread(fros + i*(*total_number_of_stations), 
306                              *total_number_of_stations*sizeof(int), 1, fp);
307        if ((int) elements_read == 0 && ferror(fp)){
308          fprintf(stderr, "Error reading start times\n");
309          return -3;
310        }           
311                           
312                           
313        elements_read = fread(lros + i*(*total_number_of_stations), 
314                              *total_number_of_stations*sizeof(int), 1, fp);
315        if ((int) elements_read == 0 && ferror(fp)){
316          fprintf(stderr, "Error reading stop times\n");       
317          return -3;
318        }                     
319
320        /* Compute the size of the data block for this source */
321        numData = getNumData(fros + i*(*total_number_of_stations), 
322                             lros + i*(*total_number_of_stations), 
323                             (*total_number_of_stations));
324
325        /* Sanity check */
326        if (numData < 0)
327        {
328            fprintf(stderr,"Size of data block appears to be negative!\n");
329            return -1;       
330        }
331
332        if (numDataMax < numData)
333        {
334            numDataMax = numData;
335        }
336
337        fclose(fp);         
338    }
339
340   
341    // Store time resolution and number of timesteps   
342    // These are the same for all stations as tested above, so
343    // we take the first one.
344    *delta_t = (double)mytgs0[0].dt;
345    *number_of_time_steps = mytgs0[0].nt;
346
347    free(mytgs);
348
349    return 0; // Succesful execution
350}
351
352
353float** _read_mux2(int numSrc, 
354                   char **muxFileNameArray, 
355                   float *weights, 
356                   double *params, 
357                   int *number_of_stations,
358                   long *permutation,
359                   int verbose)
360{
361    FILE *fp;
362    int total_number_of_stations, i, isrc, ista, k;
363    char *muxFileName;
364    int istart=-1, istop=-1;
365    int number_of_selected_stations;
366    float *muxData=NULL; // Suppress warning
367    long numData;
368
369    int len_sts_data, error_code;
370    float **sts_data;
371    float *temp_sts_data;
372
373    long int offset;
374
375    int number_of_time_steps, N;
376    double delta_t;
377   
378    size_t elements_read;
379   
380    // Shorthands pointing to memory blocks for each source
381    int *fros_per_source=NULL;     
382    int *lros_per_source=NULL;         
383
384   
385    error_code = _read_mux2_headers(numSrc, 
386                                    muxFileNameArray, 
387                                    &total_number_of_stations,
388                                    &number_of_time_steps,
389                                    &delta_t,
390                                    verbose);
391    if (error_code != 0) {
392      printf("urs_ext.c: Internal function _read_mux2_headers failed: Error code = %d\n", 
393             error_code);
394     
395      return NULL;
396    }
397   
398
399    // Apply rule that an empty permutation file means 'take all stations'
400    // We could change this later by passing in None instead of the empty
401    // permutation.
402    number_of_selected_stations = *number_of_stations; 
403    if (number_of_selected_stations == 0)
404    {
405        number_of_selected_stations = total_number_of_stations; 
406   
407        // Return possibly updated number of stations
408        *number_of_stations = total_number_of_stations;     
409     
410        // Create the Identity permutation vector
411        permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
412        for (i = 0; i < number_of_selected_stations; i++)
413        {
414            permutation[i] = (long) i; 
415        }
416    }
417   
418    // The params array is used only for passing data back to Python.
419    params[0] = (double) number_of_selected_stations;
420    params[1] = (double) delta_t;
421    params[2] = (double) number_of_time_steps;
422   
423   
424   
425    // Make array(s) to hold demuxed data for stations given in the
426    // permutation file
427    sts_data = (float**) malloc(number_of_selected_stations*sizeof(float*));
428    if (sts_data == NULL)
429    {
430        printf("ERROR: Memory for sts_data could not be allocated.\n");
431        return NULL;
432    }
433
434    // For each selected station, allocate space for its data
435    len_sts_data = number_of_time_steps + POFFSET; // Max length of each timeseries?
436    for (i = 0; i < number_of_selected_stations; i++)
437    {
438        // Initialise sts_data to zero
439        sts_data[i] = (float*) calloc(len_sts_data, sizeof(float));
440        if (sts_data[i] == NULL)
441        {
442            printf("ERROR: Memory for sts_data could not be allocated.\n");
443            return NULL;
444        }
445    }
446
447    temp_sts_data = (float*) calloc(len_sts_data, sizeof(float));
448
449    muxData = (float*) malloc(numDataMax*sizeof(float));
450   
451    // Loop over all sources
452    for (isrc = 0; isrc < numSrc; isrc++)
453    {
454   
455        // Shorthands to local memory
456        fros_per_source = (int*) fros + isrc*total_number_of_stations; 
457        lros_per_source = (int*) lros + isrc*total_number_of_stations;     
458           
459   
460        // Read in data block from mux2 file
461        muxFileName = muxFileNameArray[isrc];
462        if((fp = fopen(muxFileName, "r")) == NULL)
463        {
464            fprintf(stderr, "cannot open file %s\n", muxFileName);
465            return NULL;                   
466        }
467
468        if (verbose){
469            printf("Reading mux file %s\n", muxFileName);
470        }
471
472        offset = sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
473        fseek(fp, offset, 0);
474
475        numData = getNumData(fros_per_source, 
476                             lros_per_source, 
477                             total_number_of_stations);
478                             
479                             
480                                     
481        elements_read = fread(muxData, ((int) numData)*sizeof(float), 1, fp); 
482        if ((int) elements_read == 0 && ferror(fp)) {
483          fprintf(stderr, "Error reading mux data\n");
484          return NULL;
485        }
486               
487        fclose(fp);
488       
489
490        // FIXME (Ole): This is where Nariman and Ole traced the platform dependent
491        // difference on 11 November 2008. We don't think the problem lies in the
492        // C code. Maybe it is a problem with the MUX files written by the unit test
493        // that fails on Windows but works OK on Linux. JJ's test on 17th November shows
494        // that as far as Python is concerned this file should be OK on both platforms.
495       
496       
497       
498        // These printouts are enough to show the problem when compared
499        // on the two platforms
500        //printf("\nRead %d elements, ", (int) numData);
501        //printf("muxdata[%d]=%f\n", 39, muxData[39]);         
502       
503       
504        /*
505        In Windows we get
506       
507        Read 85 elements, muxdata[39]=0.999574
508        Read 85 elements, muxdata[39]=-0.087599
509        Read 85 elements, muxdata[39]=-0.087599
510       
511       
512        In Linux we get (the correct)
513       
514        Read 85 elements, muxdata[39]=0.999574
515        Read 85 elements, muxdata[39]=-0.087599
516        Read 85 elements, muxdata[39]=1.490349
517        */
518       
519       
520       
521       
522       
523        // loop over stations present in the permutation array
524        //     use ista with mux data
525        //     use i with the processed data to be returned         
526        for(i = 0; i < number_of_selected_stations; i++)
527        {               
528   
529            ista = (int) permutation[i]; // Get global index into mux data 
530       
531            // fill the data0 array from the mux file, and weight it
532            fillDataArray(ista, 
533                          total_number_of_stations, 
534                          number_of_time_steps,
535                          mytgs0[ista].ig, // Grid number (if -1 fill with zeros)
536                          fros_per_source, 
537                          lros_per_source, 
538                          temp_sts_data, 
539                          &istart, 
540                          &istop, 
541                          muxData);
542
543            // Weight appropriately and add
544            for(k = 0; k < mytgs0[ista].nt; k++)
545            {
546                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
547                {
548                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
549                }
550                else
551                {
552                    sts_data[i][k] = NODATA;
553                }
554                //printf("%d: temp_sts_data[%d]=%f\n", i, k, temp_sts_data[k]);     
555           
556            }
557           
558           
559            // Update metadata (e.g. start time and end time)
560            N = number_of_time_steps;
561           
562            if (isrc == 0) {
563                // Assign values for first source
564                sts_data[i][N] = (float)mytgs0[ista].geolat;
565                sts_data[i][N+1] = (float)mytgs0[ista].geolon;
566                sts_data[i][N+2] = (float)mytgs0[ista].z;
567                sts_data[i][N+3] = (float)fros_per_source[ista];
568                sts_data[i][N+4] = (float)lros_per_source[ista];
569            } else {
570                // Update first and last timesteps for subsequent sources
571                if (sts_data[i][N+3] > (float)fros_per_source[ista]) {         
572                    if (verbose) {
573                        printf("Adjusting start time for station %d and source %d",
574                               ista, isrc);
575                        printf(" from %f to %f\n", 
576                               sts_data[i][N+3], 
577                               (float) fros_per_source[ista]); 
578                    }
579                    sts_data[i][N+3] = (float) fros_per_source[ista];
580                }
581               
582                if (sts_data[i][N+4] < (float) lros_per_source[ista]) {         
583                    if (verbose) {
584                        printf("Adjusting end time for station %d and source %d",
585                               ista, isrc);
586                        printf(" from %f to %f\n", 
587                               sts_data[i][N+4], 
588                               (float) lros_per_source[ista]); 
589                    }
590                    sts_data[i][N+4] = (float) lros_per_source[ista];
591                }               
592            }
593        }
594    }
595
596    //printf("sts_data[1,8]=%f\n", sts_data[1][8]);
597   
598    free(muxData);
599    free(temp_sts_data);
600    free(fros);
601    free(lros);
602    free(mytgs0);
603
604    return sts_data;
605}
606
607/////////////////////////////////////////////////////////////////////////
608//Python gateways
609PyObject *read_mux2(PyObject *self, PyObject *args)
610{
611    /*Read in mux 2 file
612
613    Python call:
614    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
615
616    NOTE:
617    A Python int is equivalent to a C long
618    (this becomes really important on 64 bit architectures)
619   
620    A Python double corresponds to a C double
621    */
622   
623    PyObject *filenames;
624    PyArrayObject *pyweights;
625    PyArrayObject *file_params;
626    PyArrayObject *permutation;  // Ordering of selected stations   
627    PyArrayObject *pydata;
628    PyObject *fname;
629
630    char **muxFileNameArray;
631    float **cdata;
632    float *weights;
633    int dimensions[2];
634    int numSrc;
635    int verbose;
636    int total_number_of_stations;
637    int number_of_selected_stations;   
638    int nt;
639    double dt;
640    int i;
641    int j;
642    int start_tstep;
643    int finish_tstep;
644    int it;
645    int time;
646    int num_ts;
647   
648    // Convert Python arguments to C
649    if (!PyArg_ParseTuple(args, "iOOOOi",
650              &numSrc, &filenames, &pyweights, &file_params, 
651              &permutation, &verbose)) 
652    {
653            PyErr_SetString(PyExc_RuntimeError, 
654                "Input arguments to read_mux2 failed");
655            return NULL;
656    }
657
658    if(!PyList_Check(filenames)) 
659    {
660        PyErr_SetString(PyExc_TypeError, "get_first_elem expects a list");
661        return NULL;
662    }
663
664    if(PyList_Size(filenames) == 0)
665    {
666        PyErr_SetString(PyExc_ValueError, "empty lists not allowed");
667        return NULL;
668    }
669
670    if (pyweights->nd != 1 || pyweights->descr->type_num != PyArray_DOUBLE) 
671    {
672        PyErr_SetString(PyExc_ValueError,
673            "pyweights must be one-dimensional and of type double");
674        return NULL; 
675    }
676
677    if(PyList_Size(filenames) != pyweights->dimensions[0])
678    {
679        PyErr_SetString(PyExc_ValueError, 
680            "Must specify one weight for each filename");
681        return NULL;
682    }
683
684    muxFileNameArray = (char**)malloc(numSrc*sizeof(char*));
685    if (muxFileNameArray == NULL) 
686    {
687        PyErr_SetString(PyExc_ValueError, 
688                        "ERROR: Memory for muxFileNameArray could not be allocated.");
689        return NULL;
690    }
691
692    for (i = 0; i < numSrc; i++)
693    {
694
695        fname = PyList_GetItem(filenames, i);
696        if (!fname) 
697        {
698            PyErr_SetString(PyExc_ValueError, "filename not a string");
699            return NULL;
700        }       
701
702        muxFileNameArray[i] = PyString_AsString(fname);
703        if (muxFileNameArray[i] == NULL) 
704        {
705            PyErr_SetString(PyExc_ValueError, 
706          "ERROR: Memory for muxFileNameArray could not be allocated.\n");
707            return NULL;
708        }
709    }
710
711    if (file_params->nd != 1 || file_params->descr->type_num != PyArray_DOUBLE) 
712    {
713        PyErr_SetString(PyExc_ValueError,
714          "file_params must be one-dimensional and of type double");
715        return NULL; 
716    }
717
718
719    // Create array for weights which are passed to read_mux2
720    weights = (float*) malloc(numSrc*sizeof(float));
721    for (i = 0; i < numSrc; i++)
722    {
723        weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0]));
724    }
725   
726    // Desired number of stations
727    number_of_selected_stations = (int) permutation->dimensions[0];
728
729    // Read in mux2 data from file
730    cdata = _read_mux2(numSrc, 
731                       muxFileNameArray, 
732                       weights, 
733                       (double*)file_params->data,
734                       &number_of_selected_stations, 
735                       (long*) permutation->data,
736                       verbose);
737
738    if (!cdata) 
739    {
740        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
741        return NULL;
742    }       
743               
744               
745    // Allocate space for return vector
746    total_number_of_stations = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
747    dt = *(double*)(file_params->data + 1*file_params->strides[0]);
748    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
749
750   
751    // Find min and max start times of all gauges
752    start_tstep = nt + 1;
753    finish_tstep = -1;
754    for (i = 0; i < number_of_selected_stations; i++)
755    {
756        //printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
757        // printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]);   
758   
759        if ((int)cdata[i][nt + 3] < start_tstep)
760        {
761            start_tstep = (int)cdata[i][nt + 3];
762        }
763        if ((int)cdata[i][nt + 4] > finish_tstep)
764        {
765            finish_tstep = (int)cdata[i][nt + 4]; 
766        }
767    }
768
769    if ((start_tstep > nt) | (finish_tstep < 0))
770    {
771        printf("ERROR: Gauge data has incorrect start and finish times:\n");
772        printf("   start_tstep = %d, max_number_of_steps = %d\n", 
773               start_tstep, nt);
774        printf("   finish_tstep = %d, min_number_of_steps = %d\n", 
775               finish_tstep, 0);   
776           
777        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
778        return NULL;
779    }
780
781    if (start_tstep >= finish_tstep)
782    {
783        PyErr_SetString(PyExc_ValueError,
784                    "ERROR: Gauge data has non-postive_length");
785        return NULL;
786    }
787
788    num_ts = finish_tstep - start_tstep + 1;
789    dimensions[0] = number_of_selected_stations;
790    dimensions[1] = num_ts + POFFSET;
791   
792    pydata = (PyArrayObject*) anuga_FromDims(2, dimensions, PyArray_DOUBLE);
793    if(pydata == NULL)
794    {
795        PyErr_SetString(PyExc_ValueError, 
796          "ERROR: Memory for pydata array could not be allocated.");
797        return NULL;
798    }
799
800   
801    // Each gauge begins and ends recording at different times. When a gauge is
802    // not recording but at least one other gauge is.
803    // Pad the non-recording gauge array with zeros.
804    //printf("\nData put into the cdata array from C code\n");
805    for (i = 0; i < number_of_selected_stations; i++)
806    {
807        time = 0;
808        for (it = 0; it < finish_tstep; it++)
809        {
810            if ((it + 1 >= start_tstep) && (it + 1 <= finish_tstep))
811            {
812                if (it + 1 > (int)cdata[i][nt + 4])
813                {
814                    // This gauge has stopped recording but others are
815                    // still recording
816                    *(double*)(pydata->data + i*pydata->strides[0] 
817                               + time*pydata->strides[1]) = 
818                      0.0;
819                }
820                else
821                {
822                  //printf("cdata[%d][%d] = %f\n", i, it, cdata[i][it]);
823               
824                    *(double*)(pydata->data + i*pydata->strides[0] 
825                                    + time*pydata->strides[1]) = 
826                        cdata[i][it];
827                }
828                time++;
829            }
830        }
831        // Pass back lat,lon,elevation
832        for (j = 0; j < POFFSET; j++)
833        {
834          *(double*)(pydata->data + i*pydata->strides[0] 
835                     + (num_ts + j)*pydata->strides[1]) = 
836            cdata[i][nt + j];
837        }
838    }
839
840    free(weights);
841   
842    // Free filename array, but not independent Python strings
843    // FIXME(Ole): Do we need to update a reference counter in this case?
844    free(muxFileNameArray);
845   
846    for (i = 0; i < number_of_selected_stations; ++i)
847    {
848        free(cdata[i]);
849    }
850    free(cdata);
851
852    return PyArray_Return(pydata);
853}
854
855//-------------------------------
856// Method table for python module
857//-------------------------------
858static struct PyMethodDef MethodTable[] = {
859    {"read_mux2", read_mux2, METH_VARARGS, "Print out"},
860    {NULL, NULL}
861};
862
863// Module initialisation
864void initurs_ext(void){
865    Py_InitModule("urs_ext", MethodTable);
866
867    import_array(); // Necessary for handling of NumPY structures
868}
Note: See TracBrowser for help on using the repository browser.