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

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

Added verbose flag to c code

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