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

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

Modifications by John J.

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   /*now allocate space for, and read in, the structures for each station*/
213   mytgs0 = (struct tgsrwg *) malloc(nsta0*sizeof(struct tgsrwg));
214   fread(&mytgs0[0], nsta0*sizeof(struct tgsrwg), 1, fp);
215
216   /*make an array to hold the start and stop steps for each station for each
217     source*/   
218   //FIXME: Should only store start and stop times for one source
219   fros = (int *) malloc(nsta0*numSrc*sizeof(int));
220   lros = (int *) malloc(nsta0*numSrc*sizeof(int));
221   
222   /* read the start and stop steps for source 0 into the array */   
223   fread(fros,nsta0*sizeof(int),1,fp);
224   fread(lros,nsta0*sizeof(int),1,fp);
225
226   /* compute the size of the data block for source 0 */   
227   numData = getNumData(fros, lros, nsta0);
228   num_ts = mytgs0[0].nt;
229   
230   /* Burbidge: Added a sanity check here */
231   if (numData < 0) {
232     fprintf(stderr,"Size of data block appears to be negative!\n");
233     PyErr_SetString(PyExc_RuntimeError,"");
234     return NULL; 
235   }
236   fclose(fp); 
237   // Allocate header array for each remaining source file.
238   // FIXME: only need to store for one source and which is compared to all subsequent
239   // source headers
240   if(numSrc > 1){
241      /* allocate space for tgsrwg for the other sources */
242      mytgs = (struct tgsrwg *)malloc( nsta0*sizeof(struct tgsrwg) );
243   } else {
244     /* FIXME (JJ): What should happen in case the are no source files?*/
245     /* If we exit here, tests will break */
246     // fprintf(stderr, "No source file specified\n");
247     // exit(-1);       
248   }
249   
250   /* loop over remaining sources, check compatibility, and read them into
251    *muxData */
252   for(isrc=1; isrc<numSrc; isrc++){
253     muxFileName = muxFileNameArray[isrc];
254     
255     /* open the mux file */
256     if((fp=fopen(muxFileName,"r"))==NULL)
257       {
258         fprintf(stderr, "cannot open file %s\n", muxFileName);
259         PyErr_SetString(PyExc_RuntimeError,"");
260         return NULL; 
261       }
262     
263     /* check that the mux files are compatible */     
264     fread(&nsta,sizeof(int),1,fp);
265     if(nsta != nsta0){
266       fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
267       PyErr_SetString(PyExc_RuntimeError,"");
268       fclose(fp);
269       return NULL; 
270     }
271     fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp);
272     for(ista=0; ista < nsta; ista++){
273       if(mytgs[ista].dt != mytgs0[ista].dt){
274         fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
275         fclose(fp);
276         return NULL;               
277       }   
278       if(mytgs[ista].nt != num_ts){
279         fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
280         PyErr_SetString(PyExc_RuntimeError,"");
281         fclose(fp);
282         return NULL;             
283       }
284     }
285
286      /* read the start and stop times for this source */
287      fread(fros+isrc*nsta0,nsta0*sizeof(int),1,fp);
288      fread(lros+isrc*nsta0,nsta0*sizeof(int),1,fp);
289     
290      /* compute the size of the data block for this source */
291      numData = getNumData(fros+isrc*nsta0, lros+isrc*nsta0, nsta0);
292
293      /* Burbidge: Added a sanity check here */
294      if (numData < 0){
295          fprintf(stderr,"Size of data block appears to be negative!\n");
296          PyErr_SetString(PyExc_RuntimeError,"");
297          return NULL;
298      }
299      fclose(fp);             
300   }
301   N = permutation_array->dimensions[0];
302   
303   params[0]=(double)N;
304   params[1]=(double)mytgs0[0].dt;
305   params[2]=(double)num_ts;
306 
307   
308   // Find min and max start times of all gauges
309   // FIXME: ONLY CHECK gauges in permutation for start and finsh times
310   start_tstep=num_ts+1;
311   finish_tstep=-1;
312   for (ista=0;ista<nsta0;ista++){
313     if (fros[ista]< start_tstep){
314       start_tstep=fros[ista];
315     }
316     if (lros[0+nsta0*ista] > finish_tstep){
317       finish_tstep=lros[ista]; 
318     }
319   }
320   
321   if ((start_tstep>num_ts) | (finish_tstep < 0)){
322     PyErr_SetString(PyExc_RuntimeError,"Gauge data has incorrect start and finsh times");
323     return NULL;
324   }
325   
326   if (start_tstep>=finish_tstep){
327     PyErr_SetString(PyExc_RuntimeError,"Gauge data has non-postive_length");
328     return NULL;
329   }
330   
331   /* Make array(s) to hold the demuxed data */
332   len_sts_data=num_ts+POFFSET;
333   dimensions[0]=N;
334   dimensions[1]=len_sts_data;
335   sts_data = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
336   if(sts_data == NULL){
337    PyErr_SetString(PyExc_RuntimeError,"Memory for pydata array could not be allocated");
338    return NULL;
339   }
340   
341   /* Initialise sts data to zero */
342   for(i=0; i<N; i++){ 
343     ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]);
344     for (t=0;t<num_ts;t++){
345       *(double *) (sts_data -> data+i*sts_data->strides[0]+t*sts_data->strides[1])=0.0;
346     }
347     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;
348     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;
349     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;
350     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];
351     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];
352   } 
353
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 */
362     muxFileName = muxFileNameArray[isrc];
363     if((fp=fopen(muxFileName,"r"))==NULL)
364       {
365         fprintf(stderr, "Cannot open file %s\n", muxFileName);
366         PyErr_SetString(PyExc_RuntimeError,"");
367         return NULL;
368       }
369       
370     if (verbose){
371       printf("Reading mux file %s\n",muxFileName);
372     }
373
374     offset=sizeof(int)+nsta0*(sizeof(struct tgsrwg)+2*sizeof(int));
375     fseek(fp,offset,0);
376     
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(j=0; j<N; j++){ 
384       ista = *(long *) (permutation_array -> data+j*permutation_array->strides[0]);
385       //printf("%d\n",(int)*(long *) (permutation_array -> data+j*permutation_array->strides[0]));
386       
387       /* fill the data0 array from the mux file, and weight it */
388       fillDataArray(ista, nsta0, num_ts, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
389       
390       /* weight appropriately and add */
391       for(t=0; t<num_ts; t++){
392         if((isdata(*(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]))) && isdata(temp_sts_data[t])){
393           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) += temp_sts_data[t] * weights[isrc];
394         }else{
395           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) = 0.0;
396         }
397       }
398     }
399   }   
400   free(muxData);
401   free(temp_sts_data);
402   free(fros);
403   free(lros);
404   return sts_data;
405}   
406
407
408/* thomas */
409void fillDataArray(int ista, int nsta, int nt, int ig, int *nst, int *nft, float *data, int *istart_p, int *istop_p, float *muxData)
410{
411   int it, last_it, jsta;
412   long int offset=0;
413
414
415   last_it=-1;
416   /* make arrays of starting and finishing time steps for the tide gauges */
417   /* and fill them from the file */
418     
419   /* update start and stop timesteps for this gauge */
420   if(nst[ista]!=-1){
421     if(*istart_p==-1){
422         *istart_p=nst[ista];
423     }else{
424         *istart_p=((nst[ista]<*istart_p)?nst[ista]:*istart_p);
425     }
426   }
427   if(nft[ista]!=-1){
428     if(*istop_p==-1){
429         *istop_p=nft[ista];
430     }else{
431         *istop_p=((nft[ista]<*istop_p)?nft[ista]:*istop_p);
432     }
433   }     
434   if(ig==-1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
435   {
436      /* gauge never started recording, or was outside of all grids, fill array with 0 */
437      for(it=0; it<nt; it++)
438         data[it] = 0.0;
439   }   
440   else
441   {
442      for(it=0; it<nt; it++)
443      {
444         last_it = it;
445         /* skip t record of data block */
446         offset++;
447         /* skip records from earlier tide gauges */
448         for(jsta=0; jsta<ista; jsta++)
449            if(it+1>=nst[jsta]&&it+1<=nft[jsta])
450               offset++;
451               
452         /* deal with the tide gauge at hand */
453         if(it+1>=nst[ista]&&it+1<=nft[ista])
454         /* gauge is recording at this time */
455         {
456            memcpy(data+it,muxData+offset,sizeof(float));
457            offset++;
458         }
459         else if (it+1<nst[ista])
460         {
461            /* gauge has not yet started recording */
462            data[it] = 0.0;
463         }   
464         else
465         /* gauge has finished recording */                                           
466         {
467            data[it] = NODATA;
468            break;
469         }
470   
471         /* skip records from later tide gauges */
472         for(jsta=ista+1; jsta<nsta; jsta++)
473            if(it+1>=nst[jsta]&&it+1<=nft[jsta])
474               offset++;
475      }
476   
477      if(last_it < nt - 1)
478         /* the loop was exited early because the gauge had finished recording */
479         for(it=last_it+1; it < nt; it++)
480            data[it] = NODATA;
481   }
482} 
483
484/* Burbidge: No "offset" is sent. Replace with max. Added grid_id */
485void wrttxt( fname, dt, nt, x, beg, lat, lon, max, depth, grid_id, istart, istop )
486char *fname;
487float dt, *x, beg, max, lat, lon, depth;
488int grid_id;
489int nt;
490int istart, istop;
491{
492   int it;
493   float t;
494   FILE *fp;
495
496   fp = fopen(fname,"w");
497   fprintf(fp,"> lat %g lon %g depth %g max %g start_time %g stop_time %g grid_id %d\n", lat, lon, depth, max, 
498            istart*dt, istop*dt, grid_id );
499   for(it=0;it<nt;it++)
500   {
501      t=beg+it*dt;
502      fprintf(fp,"%9.3f %g\n", t, x[it]);
503   }
504   fclose(fp);
505}
506 
507char isdata(float x)
508{
509  //char value;
510   if(x < NODATA + EPSILON && NODATA < x + EPSILON)
511      return 0;
512   else
513      return 1; 
514}
515
516long getNumData(int *fros, int *lros, int nsta)
517/* calculates the number of data in the data block of a mux file */
518/* based on the first and last recorded output steps for each gauge */ 
519{
520   int ista, last_output_step;
521   long numData = 0;
522
523   last_output_step = 0;   
524   for(ista=0; ista < nsta; ista++)
525      if(*(fros + ista) != -1)
526      {
527         numData += *(lros + ista) - *(fros + ista) + 1;
528         last_output_step = (last_output_step < *(lros+ista) ? *(lros+ista):last_output_step);
529      }   
530   numData += last_output_step*nsta; /* these are the t records */
531   return numData;
532}   
533
534
535
536//-------------------------------
537// Method table for python module
538//-------------------------------
539static struct PyMethodDef MethodTable[] = {
540  {"read_mux2", read_mux2, METH_VARARGS, "Print out"},
541  {NULL, NULL}
542};
543
544// Module initialisation
545void initurs_ext(void){
546  Py_InitModule("urs_ext", MethodTable);
547
548  import_array(); // Necessary for handling of NumPY structures
549}
Note: See TracBrowser for help on using the repository browser.