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

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

Reverted MUX reader to revision 5470 - the last working version.
This one still has servere memory leaks and needs to be redeveloped.

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