Ignore:
Timestamp:
Jul 17, 2008, 1:46:09 PM (16 years ago)
Author:
nariman
Message:

Elimination of memory leaks
Code reformatting
Expansion of first for loop to include first muxfile.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/urs_ext.c

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