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

Last change on this file since 5375 was 5375, checked in by ole, 16 years ago

Comments from Ole in regard to uninitialised pointer

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