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

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

Removed duplication of sts data originally used to pass back to python

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