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

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

'added check to ensure that permutation indices are not out of bounds'

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