Changeset 5552


Ignore:
Timestamp:
Jul 22, 2008, 2:27:54 PM (15 years ago)
Author:
nariman
Message:

Refactoring of code and creation of _read_mux2_headers.
Memory leak caused by muxData.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5542 r5552  
    2626
    2727void fillDataArray(int, int, int, int, int *, int *, float *,
    28                    int *, int *, float *);
     28           int *, int *, float *);
    2929long getNumData(const int *fros, const int *lros, const int nsta);
    3030char isdata(float);
     31void _read_mux2_headers(int numSrc,
     32                           char **muxFileNameArray,
     33                   double *params,
     34                                   int** fros, int** lros, struct tgsrwg ** mytgs0,
     35                   int* nsta0,
     36                   int verbose);
     37
    3138float** _read_mux2(int numSrc,
    32                    char **muxFileNameArray,
    33                    float *weights,
    34                    double *params,
    35                    int *number_of_stations,
    36                    long *permutation,
    37                    int verbose);
    38 
    39 PyObject *read_mux2(PyObject *self, PyObject *args){
     39           char **muxFileNameArray,
     40           float *weights,
     41           double *params,
     42           int *number_of_stations,
     43           long *permutation,
     44           int verbose);
     45
     46    PyObject *read_mux2(PyObject *self, PyObject *args){
    4047    /*Read in mux 2 file
    4148
    … …  
    7784    // Convert Python arguments to C
    7885    if (!PyArg_ParseTuple(args, "iOOOOi",
    79                           &numSrc, &filenames, &pyweights, &file_params,
    80                           &permutation, &verbose))
     86              &numSrc, &filenames, &pyweights, &file_params,
     87              &permutation, &verbose))
    8188    {
    8289            PyErr_SetString(PyExc_RuntimeError,
    … …  
    107114    {
    108115        PyErr_SetString(PyExc_ValueError,
    109                         "Must specify one weight for each filename");
     116            "Must specify one weight for each filename");
    110117        return NULL;
    111118    }
    … …  
    115122    {
    116123        PyErr_SetString(PyExc_ValueError,
    117           "ERROR: Memory for muxFileNameArray could not be allocated.");
     124      "ERROR: Memory for muxFileNameArray could not be allocated.");
    118125        return NULL;
    119126    }
    … …  
    133140        {
    134141            PyErr_SetString(PyExc_ValueError,
    135               "ERROR: Memory for muxFileNameArray could not be allocated.\n");
     142          "ERROR: Memory for muxFileNameArray could not be allocated.\n");
    136143            return NULL;
    137144        }
    … …  
    141148    {
    142149        PyErr_SetString(PyExc_ValueError,
    143               "file_params must be one-dimensional and of type double");
     150          "file_params must be one-dimensional and of type double");
    144151        return NULL;
    145152    }
    … …  
    150157    for (i = 0; i < numSrc; i++)
    151158    {
    152         weights[i] = (float)(*(double*)
    153                      (pyweights->data + i*pyweights->strides[0]));
     159        weights[i] = (float)(*(double*)(pyweights->data + i*pyweights->strides[0]));
    154160    }
    155161   
    … …  
    162168                       weights,
    163169                       (double*)file_params->data,
    164                        &number_of_selected_stations,
    165                        (long*) permutation->data,
     170                       &number_of_selected_stations,
     171                       (long*) permutation->data,
    166172                       verbose);
    167173
    … …  
    169175    {
    170176        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
    171         return NULL;
     177        return NULL;
    172178    }       
    173                        
    174                        
     179               
     180               
    175181    // Allocate space for return vector
    176182    nsta0 = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
    … …  
    185191    {
    186192        //printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
    187         // printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]);       
    188        
     193        // printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]);   
     194   
    189195        if ((int)cdata[i][nt + 3] < start_tstep)
    190196        {
    … …  
    200206    {
    201207        printf("ERROR: Gauge data has incorrect start and finish times:\n");
    202         printf("   start_tstep = %d, max_number_of_steps = %d\n",
    203                    start_tstep, nt);
    204         printf("   finish_tstep = %d, min_number_of_steps = %d\n",
    205                    finish_tstep, 0);   
    206                    
    207         PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
     208        printf("   start_tstep = %d, max_number_of_steps = %d\n",
     209               start_tstep, nt);
     210        printf("   finish_tstep = %d, min_number_of_steps = %d\n",
     211               finish_tstep, 0);   
     212           
     213        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
    208214        return NULL;
    209215    }
    … …  
    212218    {
    213219        PyErr_SetString(PyExc_ValueError,
    214                         "ERROR: Gauge data has non-postive_length");
     220                    "ERROR: Gauge data has non-postive_length");
    215221        return NULL;
    216222    }
    … …  
    224230    {
    225231        PyErr_SetString(PyExc_ValueError,
    226               "ERROR: Memory for pydata array could not be allocated.");
     232          "ERROR: Memory for pydata array could not be allocated.");
    227233        return NULL;
    228234    }
    … …  
    242248                {
    243249                    // This gauge has stopped recording but others are
    244                     // still recording
     250            // still recording
    245251                    *(double*)(pydata->data + i*pydata->strides[0]
    246                                             + time*pydata->strides[1]) =
    247                                             0.0;
     252                                    + time*pydata->strides[1]) =
     253                        0.0;
    248254                }
    249255                else
    250256                {
    251257                    *(double*)(pydata->data + i*pydata->strides[0]
    252                                             + time*pydata->strides[1]) =
    253                                             cdata[i][it];
     258                                    + time*pydata->strides[1]) =
     259                        cdata[i][it];
    254260                }
    255261                time++;
    … …  
    260266        {
    261267            *(double*)(pydata->data + i*pydata->strides[0]
    262                                     + (num_ts + j)*pydata->strides[1]) =
    263                                     cdata[i][nt + j];
     268                                + (num_ts + j)*pydata->strides[1]) =
     269                    cdata[i][nt + j];
    264270        }
    265271    }
    … …  
    279285    return  PyArray_Return(pydata);
    280286}
     287
    281288
    282289float** _read_mux2(int numSrc,
    … …  
    284291                   float *weights,
    285292                   double *params,
    286                    int *number_of_stations,
    287                    long *permutation,
     293                   int *number_of_stations,
     294                   long *permutation,
    288295                   int verbose)
    289296{
    290297    FILE *fp;
    291     int nsta, nsta0, i, isrc, ista, k;
    292     struct tgsrwg *mytgs=0, *mytgs0=0;
    293     char *muxFileName;                                                                 
     298    int nsta0, i, isrc, ista, k;
     299    struct tgsrwg* mytgs0=NULL;
     300    char *muxFileName;
    294301    int istart, istop;
    295302    int *fros=0, *lros=0;
    296303    int number_of_selected_stations;
    297     char susMuxFileName;
    298304    float *muxData=NULL; // Suppress warning
    299305    long numData;
    … …  
    305311    long int offset;
    306312
    307     /* Allocate space for the names and the weights and pointers to the data*/
    308 
    309     /* Check that the input files have mux2 extension*/
    310     susMuxFileName = 0;
    311     for(isrc = 0; isrc < numSrc; isrc++)
    312     {
    313         muxFileName = muxFileNameArray[isrc];
    314         if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4,
    315                                      "mux2") != 0)
    316         {
    317             susMuxFileName = 1;
    318             break;
    319         }
    320     }
    321 
    322     if(susMuxFileName)
    323     {
    324         printf("\n**************************************************************************\n");
    325         printf("   WARNING: This program operates only on multiplexed files in mux2 format\n");
    326         printf("   At least one input file name does not end with mux2\n");
    327         printf("   Check your results carefully!\n");
    328         printf("**************************************************************************\n\n");
    329     }   
    330 
    331     if (verbose)
    332     {
    333         printf("Reading mux header information\n");
    334     }
    335 
    336     /* Loop over all sources, read headers and check compatibility */
    337     for (isrc = 0; isrc < numSrc; isrc++)
    338     {
    339         muxFileName = muxFileNameArray[isrc];
    340 
    341         /* open the mux file */
    342         if((fp = fopen(muxFileName, "r")) == NULL)
    343         {
    344             fprintf(stderr, "cannot open file %s\n", muxFileName);
    345             return NULL; 
    346         }
    347        
    348         if (!isrc)
    349         {
    350             fread(&nsta0, sizeof(int), 1, fp);
    351        
    352             fros = (int*)malloc(nsta0*numSrc*sizeof(int));
    353             lros = (int*)malloc(nsta0*numSrc*sizeof(int));
    354      
    355             mytgs0 = (struct tgsrwg*)malloc(nsta0*sizeof(struct tgsrwg));
    356             mytgs = (struct tgsrwg*)malloc(nsta0*sizeof(struct tgsrwg));
    357 
    358             fread(mytgs0, nsta0*sizeof(struct tgsrwg), 1, fp);
    359         }
    360         else
    361         {
    362             /* check that the mux files are compatible */     
    363             fread(&nsta, sizeof(int), 1, fp);
    364             if(nsta != nsta0)
    365             {
    366                 fprintf(stderr,"%s has different number of stations to %s\n",
    367                         muxFileName,
    368                         muxFileNameArray[0]);
    369                 fclose(fp);
    370                 return NULL;   
    371             }
    372 
    373             fread(mytgs, nsta*sizeof(struct tgsrwg), 1, fp);
    374            
    375             for (ista = 0; ista < nsta; ista++)
    376             {
    377                 if (mytgs[ista].dt != mytgs0[ista].dt)
    378                 {
    379                     fprintf(stderr,"%s has different sampling rate to %s\n",
    380                             muxFileName,
    381                             muxFileNameArray[0]);
    382                     fclose(fp);
    383                     return NULL;                   
    384                 }   
    385                 if (mytgs[ista].nt != mytgs0[ista].nt)
    386                 {
    387                     fprintf(stderr,"%s has different series length to %s\n",
    388                             muxFileName,
    389                             muxFileNameArray[0]);
    390                     fclose(fp);
    391                     return NULL;                   
    392                 }
    393             }
    394         }
    395 
    396         /* Read the start and stop times for this source */
    397         fread(fros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
    398         fread(lros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
    399 
    400         /* Compute the size of the data block for this source */
    401         numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
    402 
    403         /* Sanity check */
    404         if (numData < 0)
    405         {
    406             fprintf(stderr,"Size of data block appears to be negative!\n");
    407             return NULL;           
    408         }
    409 
    410         fclose(fp);         
    411     }
    412 
    413     params[0] = (double)nsta0;
    414     params[1] = (double)mytgs0[0].dt;
    415     params[2] = (double)mytgs0[0].nt;
    416 
     313    _read_mux2_headers(numSrc,
     314                       muxFileNameArray,
     315                       params,
     316                       &fros,
     317                       &lros,
     318                       &mytgs0,
     319                       &nsta0,
     320                       verbose);
    417321    // Apply rule that an empty permutation file means 'take all stations'
    418322    // We could change this later by passing in None instead of the empty
    419323    // permutation.
    420     number_of_selected_stations = *number_of_stations; 
     324    number_of_selected_stations = *number_of_stations; 
    421325    if (number_of_selected_stations == 0)
    422326    {
    423327        number_of_selected_stations = nsta0; 
    424        
    425         // Return possibly updated number of stations
    426         *number_of_stations = nsta0;     
     328   
     329        // Return possibly updated number of stations
     330        *number_of_stations = nsta0;     
    427331     
    428332        // Create the Identity permutation vector
    … …  
    430334        for (i = 0; i < number_of_selected_stations; i++)
    431335        {
    432           permutation[i] = (long) i; 
     336            permutation[i] = (long) i; 
    433337        }
    434338    }
    … …  
    453357
    454358    // For each selected station, allocate space for its data
     359
    455360    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
    456361    for (i = 0; i < number_of_selected_stations; i++)
    … …  
    461366        {
    462367            printf("ERROR: Memory for sts_data could not be allocated.\n");
    463             return NULL;
    464         }
    465 
    466         ista = (int) permutation[i]; // Get global index into mux data
    467        
     368            return NULL;
     369        }
     370
     371        ista = (int) permutation[i]; // Get global index into mux data
     372   
    468373        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
    469374        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
    … …  
    472377        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
    473378    }
    474    
     379
    475380    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
    476381
    … …  
    486391        {
    487392            fprintf(stderr, "cannot open file %s\n", muxFileName);
    488             return NULL;                           
     393            return NULL;                   
    489394        }
    490395
    … …  
    507412        for(i = 0; i < number_of_selected_stations; i++)
    508413        {               
    509        
    510             ista = (int) permutation[i]; // Get global index into mux data     
    511            
     414   
     415            ista = (int) permutation[i]; // Get global index into mux data 
     416       
    512417            /* fill the data0 array from the mux file, and weight it */
    513418            fillDataArray(ista,
    514                           nsta0,
    515                           mytgs0[ista].nt,
    516                           mytgs0[ista].ig,
    517                           fros + isrc*nsta0,
    518                           lros + isrc*nsta0,
    519                           temp_sts_data,
    520                           &istart,
    521                           &istop,
    522                           muxData);
     419                          nsta0,
     420                          mytgs0[ista].nt,
     421                          mytgs0[ista].ig,
     422                          fros + isrc*nsta0,
     423                          lros + isrc*nsta0,
     424                          temp_sts_data,
     425                          &istart,
     426                          &istop,
     427                          muxData);
    523428
    524429            /* weight appropriately and add */
    … …  
    535440            }
    536441        }
    537     }
    538 
    539     free(muxData);
     442        free(muxData);
     443    }
     444
    540445    free(temp_sts_data);
    541446    free(fros);
    542447    free(lros);
    543448    free(mytgs0);
    544     free(mytgs);
    545449
    546450    return sts_data;
     451}
     452
     453void _read_mux2_headers(int numSrc,
     454                        char **muxFileNameArray,
     455                        double *params,
     456                        int** fros,
     457                        int** lros,
     458                        struct tgsrwg ** mytgs0,
     459                        int* nsta0,
     460                        int verbose)
     461{
     462    FILE *fp;
     463    int nsta, isrc, ista;
     464    struct tgsrwg *mytgs=0;
     465    char *muxFileName;                                                                 
     466    char susMuxFileName;
     467    long numData;
     468
     469    /* Allocate space for the names and the weights and pointers to the data*/
     470
     471    /* Check that the input files have mux2 extension*/
     472    susMuxFileName = 0;
     473    for(isrc = 0; isrc < numSrc; isrc++)
     474    {
     475        muxFileName = muxFileNameArray[isrc];
     476        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4,
     477                     "mux2") != 0)
     478        {
     479            susMuxFileName = 1;
     480            break;
     481        }
     482    }
     483
     484    if(susMuxFileName)
     485    {
     486        printf("\n**************************************************************************\n");
     487        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n");
     488        printf("   At least one input file name does not end with mux2\n");
     489        printf("   Check your results carefully!\n");
     490        printf("**************************************************************************\n\n");
     491    }   
     492
     493    if (verbose)
     494    {
     495        printf("Reading mux header information\n");
     496    }
     497
     498    /* Loop over all sources, read headers and check compatibility */
     499    for (isrc = 0; isrc < numSrc; isrc++)
     500    {
     501        muxFileName = muxFileNameArray[isrc];
     502
     503        /* open the mux file */
     504        if((fp = fopen(muxFileName, "r")) == NULL)
     505        {
     506            fprintf(stderr, "cannot open file %s\n", muxFileName);
     507            exit(-1); 
     508        }
     509       
     510        if (!isrc)
     511        {
     512            fread(nsta0, sizeof(int), 1, fp);
     513       
     514            *fros = (int*)malloc(*nsta0*numSrc*sizeof(int));
     515            *lros = (int*)malloc(*nsta0*numSrc*sizeof(int));
     516     
     517            *mytgs0 = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
     518            mytgs = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
     519
     520            fread(*mytgs0, *nsta0*sizeof(struct tgsrwg), 1, fp);
     521        }
     522        else
     523        {
     524            /* check that the mux files are compatible */     
     525            fread(&nsta, sizeof(int), 1, fp);
     526            if(nsta != *nsta0)
     527            {
     528                fprintf(stderr,"%s has different number of stations to %s\n",
     529                muxFileName,
     530                muxFileNameArray[0]);
     531                fclose(fp);
     532                exit(-1);   
     533            }
     534
     535            fread(mytgs, nsta*sizeof(struct tgsrwg), 1, fp);
     536           
     537            for (ista = 0; ista < nsta; ista++)
     538            {
     539                if (mytgs[ista].dt != (*mytgs0)[ista].dt)
     540                {
     541                    fprintf(stderr,"%s has different sampling rate to %s\n",
     542                    muxFileName,
     543                    muxFileNameArray[0]);
     544                    fclose(fp);
     545                    exit(-1);           
     546                }   
     547                if (mytgs[ista].nt != (*mytgs0)[ista].nt)
     548                {
     549                    fprintf(stderr,"%s has different series length to %s\n",
     550                    muxFileName,
     551                    muxFileNameArray[0]);
     552                    fclose(fp);
     553                    exit(-1);           
     554                }
     555            }
     556        }
     557
     558        /* Read the start and stop times for this source */
     559        fread(*fros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
     560        fread(*lros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
     561
     562        /* Compute the size of the data block for this source */
     563        numData = getNumData(*fros + isrc*(*nsta0), *lros + isrc*(*nsta0), (*nsta0));
     564
     565        /* Sanity check */
     566        if (numData < 0)
     567        {
     568            fprintf(stderr,"Size of data block appears to be negative!\n");
     569            exit(-1);       
     570        }
     571
     572        fclose(fp);         
     573    }
     574
     575    params[0] = (double)*nsta0;
     576    params[1] = (double)(*mytgs0)[0].dt;
     577    params[2] = (double)(*mytgs0)[0].nt;
     578
     579        free(mytgs);
    547580}
    548581
    … …  
    551584void fillDataArray(int ista, int nsta, int nt, int ig, int *nst,
    552585                   int *nft, float *data, int *istart_p,
    553                    int *istop_p, float *muxData)
     586           int *istop_p, float *muxData)
    554587{
    555588    int it, last_it, jsta;
    … …  
    587620    {
    588621        /* gauge never started recording, or was outside of all grids,
    589         fill array with 0 */
     622    fill array with 0 */
    590623        for(it = 0; it < nt; it++)
    591624        {
    … …  
    632665        if(last_it < nt - 1)
    633666            /* the loop was exited early because the gauge had
    634             finished recording */
     667        finished recording */
    635668            for(it = last_it+1; it < nt; it++)
    636669                data[it] = NODATA;
    … …  
    665698            numData += *(lros + ista) - *(fros + ista) + 1;
    666699            last_output_step = (last_output_step < *(lros+ista) ?
    667                                 *(lros+ista):last_output_step);
     700                            *(lros+ista):last_output_step);
    668701        }   
    669702        numData += last_output_step*nsta; /* these are the t records */
Note: See TracChangeset for help on using the changeset viewer.