Ignore:
Timestamp:
Jul 18, 2008, 4:59:46 PM (14 years ago)
Author:
ole
Message:

Work on permutations in urs_ext.c
It does not work yet, but the array is being passed in and the algorithm
has been described.

File:
1 edited

Legend:

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

    r5531 r5537  
    2828long getNumData(const int *fros, const int *lros, const int nsta);
    2929char isdata(float);
    30 float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose);
     30float** _read_mux2(int numSrc,
     31                   char **muxFileNameArray,
     32                   float *weights,
     33                   double *params,
     34                   int number_of_selected_stations,
     35                   int *permutation,
     36                   int verbose);
    3137
    3238PyObject *read_mux2(PyObject *self, PyObject *args){
     
    3440
    3541    Python call:
    36     read_mux2(numSrc,filenames,weights,file_params,verbose)
     42    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
    3743
    3844    NOTE:
     
    4349    PyArrayObject *pyweights;
    4450    PyArrayObject *file_params;
     51    PyArrayObject *permutation;   
    4552    PyArrayObject *pydata;
    4653    PyObject *fname;
     
    5360    int verbose;
    5461    int nsta0;
     62    int number_of_selected_stations;   
    5563    int nt;
    5664    double dt;
     
    6674
    6775    // Convert Python arguments to C
    68     if (!PyArg_ParseTuple(args, "iOOOi",
    69         &numSrc, &filenames, &pyweights, &file_params, &verbose))
     76    if (!PyArg_ParseTuple(args, "iOOOOi",
     77                          &numSrc, &filenames, &pyweights, &file_params, &permutation, &verbose))
    7078    {
    7179            PyErr_SetString(PyExc_RuntimeError,
     
    106114    }
    107115
    108     //printf("Checkpoint 2 for urs2sts_ext.c\n");   
     116    printf("Checkpoint 2 for urs2sts_ext.c\n");   
    109117    //printf("numSrc = %d\n", numSrc);
    110118    for (i = 0; i < numSrc; i++)
     
    128136        if (muxFileNameArray[i] == NULL)
    129137        {
    130             printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
     138            PyErr_SetString(PyExc_ValueError,
     139                            "ERROR: Memory for muxFileNameArray could not be allocated.\n");
    131140            return NULL;
    132141        }
     
    137146    {
    138147        PyErr_SetString(PyExc_ValueError,
    139             "file_params must be one-dimensional and of type double");
     148                        "file_params must be one-dimensional and of type double");
    140149        return NULL;
    141150    }
    142151
    143     //printf("Checkpoint 4 for urs2sts_ext.c\n");       
     152    printf("Checkpoint 4 for urs2sts_ext.c\n");       
    144153    // Create array for weights which are passed to read_mux2
    145154    weights = (float*) malloc(numSrc*sizeof(float));
     
    149158    }
    150159   
    151     //printf("Checkpoint 5 for urs2sts_ext.c\n");           
     160   
     161    number_of_selected_stations = (int) permutation->dimensions[0];
     162    printf("Checkpoint 5 for urs2sts_ext.c\n");       
    152163    // Read in mux2 data from file
    153164    cdata = _read_mux2(numSrc,
    154165                       muxFileNameArray,
    155166                       weights,
    156                        (double*)file_params->data,
     167                       (double*)file_params->data,
     168                       number_of_selected_stations, // Desired number of stations
     169                       (int*) permutation->data, // Ordering of selected stations
    157170                       verbose);
    158171
     172    if (!cdata)
     173    {
     174        PyErr_SetString(PyExc_ValueError, "No STS_DATA returned");
     175        return NULL;
     176    }       
     177                       
     178                       
    159179    //printf("Checkpoint 6 for urs2sts_ext.c\n");               
    160180    // Allocate space for return vector
     
    163183    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
    164184
    165     //printf("Checkpoint 7 for urs2sts_ext.c\n");                   
     185    printf("Checkpoint 7 for urs2sts_ext.c\n");                   
     186   
    166187    // Find min and max start times of all gauges
    167188    start_tstep = nt + 1;
    168189    finish_tstep = -1;
    169     for (i = 0; i < nsta0; i++)
    170     {
     190    for (i = 0; i < number_of_selected_stations; i++)
     191    {
     192        printf("cdata[%d] start = %f\n", i, (double) cdata[i][nt+3]);
     193        printf("cdata[%d] finish = %f\n", i, (double) cdata[i][nt+4]); 
     194       
    171195        if ((int)cdata[i][nt + 3] < start_tstep)
    172196        {
     
    179203    }
    180204
    181     //printf("Checkpoint 8 for urs2sts_ext.c\n");                       
     205    printf("Checkpoint 8 for urs2sts_ext.c\n");                       
    182206    if ((start_tstep > nt) | (finish_tstep < 0))
    183207    {
    184         printf("ERROR: Gauge data has incorrect start and finsh times\n");
     208        printf("ERROR: Gauge data has incorrect start and finish times:\n");
     209        printf("   start_tstep = %d, max_number_of_steps = %d\n", start_tstep, nt);
     210        printf("   finish_tstep = %d, min_number_of_steps = %d\n", finish_tstep, 0);   
     211        PyErr_SetString(PyExc_ValueError, "Incorrect start and finish times"); 
    185212        return NULL;
    186213    }
     
    188215    if (start_tstep >= finish_tstep)
    189216    {
    190         printf("ERROR: Gauge data has non-postive_length\n");
     217        PyErr_SetString(PyExc_ValueError, "ERROR: Gauge data has non-postive_length");
    191218        return NULL;
    192219    }
    193220
    194221    num_ts = finish_tstep - start_tstep + 1;
    195     dimensions[0] = nsta0;
     222    dimensions[0] = number_of_selected_stations;
    196223    dimensions[1] = num_ts + POFFSET;
    197224    pydata = (PyArrayObject*)PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
    198225    if(pydata == NULL)
    199226    {
    200         printf("ERROR: Memory for pydata array could not be allocated.\n");
    201         return NULL;
    202     }
    203 
    204     //printf("Checkpoint 9 for urs2sts_ext.c\n");                           
     227        PyErr_SetString(PyExc_ValueError, "ERROR: Memory for pydata array could not be allocated.");
     228        return NULL;
     229    }
     230
     231    printf("Checkpoint 9 for urs2sts_ext.c\n");                           
    205232    // Each gauge begins and ends recording at different times. When a gauge is
    206233    // not recording but at least one other gauge is. Pad the non-recording gauge
    207234    // array with zeros.
    208     for (i = 0; i < nsta0; i++)
     235    for (i = 0; i < number_of_selected_stations; i++)
    209236    {
    210237        time = 0;
     
    232259    }
    233260
    234     //printf("Checkpoint 10 for urs2sts_ext.c\n");                               
     261    printf("Checkpoint 10 for urs2sts_ext.c\n");                               
    235262    free(weights);
    236263   
     
    239266    free(muxFileNameArray);
    240267   
    241     for (i = 0; i < nsta0; ++i)
     268    for (i = 0; i < number_of_selected_stations; ++i)
    242269    {
    243270        free(cdata[i]);
     
    245272    free(cdata);
    246273
    247     //printf("Checkpoint 11 for urs2sts_ext.c\n");                                   
     274    printf("Checkpoint 11 for urs2sts_ext.c\n");                                   
    248275    return  PyArray_Return(pydata);
    249276}
     
    253280                   float *weights,
    254281                   double *params,
     282                   int number_of_selected_stations,
     283                   int *permutation,
    255284                   int verbose)
    256285{
     
    299328    }
    300329
    301     /* loop over all sources, read headers and check compatibility */
     330    printf("P1\n");
     331    /* Loop over all sources, read headers and check compatibility */
    302332    for (isrc = 0; isrc < numSrc; isrc++)
    303333    {
    304334        muxFileName = muxFileNameArray[isrc];
    305335
    306                 /* open the mux file */
     336        /* open the mux file */
    307337        if((fp = fopen(muxFileName, "r")) == NULL)
    308338        {
    309339            fprintf(stderr, "cannot open file %s\n", muxFileName);
    310             exit(-1)
     340            return NULL
    311341        }
    312342       
     
    329359            if(nsta != nsta0)
    330360            {
    331                 fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
     361                fprintf(stderr,"%s has different number of stations to %s\n",
     362                        muxFileName,
     363                        muxFileNameArray[0]);
    332364                fclose(fp);
    333                 exit(-1);   
     365                return NULL;   
    334366            }
    335367
     
    340372                if (mytgs[ista].dt != mytgs0[ista].dt)
    341373                {
    342                     fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
     374                    fprintf(stderr,"%s has different sampling rate to %s\n",
     375                            muxFileName,
     376                            muxFileNameArray[0]);
    343377                    fclose(fp);
    344                     exit(-1);                
     378                    return NULL;                   
    345379                }   
    346380                if (mytgs[ista].nt != mytgs0[ista].nt)
    347381                {
    348                     fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
     382                    fprintf(stderr,"%s has different series length to %s\n",
     383                            muxFileName,
     384                            muxFileNameArray[0]);
    349385                    fclose(fp);
    350                     exit(-1);                
     386                    return NULL;                   
    351387                }
    352388            }
    353389        }
    354390
    355         /* read the start and stop times for this source */
     391        /* Read the start and stop times for this source */
    356392        fread(fros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
    357393        fread(lros + isrc*nsta0, nsta0*sizeof(int), 1, fp);
    358394
    359         /* compute the size of the data block for this source */
     395        /* Compute the size of the data block for this source */
    360396        numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
    361397
    362         /* Burbidge: Added a sanity check here */
     398        /* Sanity check */
    363399        if (numData < 0)
    364400        {
    365401            fprintf(stderr,"Size of data block appears to be negative!\n");
    366             exit(-1);
     402            return NULL;           
    367403        }
    368404
     
    370406    }
    371407
     408    printf("P2\n");   
    372409    params[0] = (double)nsta0;
    373410    params[1] = (double)mytgs0[0].dt;
    374411    params[2] = (double)mytgs0[0].nt;
    375412
    376     /* make array(s) to hold the demuxed data */
    377     //FIXME: This can be reduced to only contain stations given in permutation file
    378     sts_data = (float**)malloc(nsta0*sizeof(float *));
     413    // Apply rule that an empty permutation file means 'take all stations'
     414    // We can change this later by passing in None instead of the empty permutation.
     415   
     416    printf("number_of_selected_stations B4 = %d\n", number_of_selected_stations);   
     417    if (number_of_selected_stations == 0)
     418    {
     419        number_of_selected_stations = nsta0; 
     420     
     421        printf("Creating identity permutation vector of length = %d\n", nsta0);
     422        // Create the Identity permutation vector
     423        permutation = (int *) malloc(number_of_selected_stations*sizeof(int));
     424        for (i = 0; i < number_of_selected_stations; i++)
     425        {
     426          permutation[i] = i; 
     427        }
     428     
     429    }
     430   
     431    printf("P3\n");       
     432    printf("Number of selected stations = %d\n", number_of_selected_stations);
     433    /* Make array(s) to hold demuxed data for stations given in the permutation file */
     434    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
    379435    if (sts_data == NULL)
    380436    {
    381437        printf("ERROR: Memory for sts_data could not be allocated.\n");
    382         exit(-1);
    383     }
    384 
    385     len_sts_data = mytgs0[0].nt + POFFSET;
    386     for (ista = 0; ista < nsta0; ista++)
    387     {
     438        return NULL;
     439    }
     440
     441    printf("P4\n");           
     442    // For each selected station, allocate space for its data
     443    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries (I think)
     444    for (i = 0; i < number_of_selected_stations; i++)
     445    {
     446        printf("P4.1, i=%d\n", i);               
     447   
    388448        // Initialise sts_data to zero
    389         sts_data[ista] = (float*)calloc(len_sts_data, sizeof(float));
    390         if (sts_data[ista] == NULL)
     449        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
     450        if (sts_data[i] == NULL)
    391451        {
    392452            printf("ERROR: Memory for sts_data could not be allocated.\n");
    393             exit(-1);
    394         }
    395 
    396         sts_data[ista][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
    397         sts_data[ista][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
    398         sts_data[ista][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
    399         sts_data[ista][mytgs0[0].nt + 3] = (float)fros[ista];
    400         sts_data[ista][mytgs0[0].nt + 4] = (float)lros[ista];
    401     }
    402 
     453            return NULL;
     454        }
     455
     456       
     457        printf("P4.2, i=%d\n", i);                     
     458        ista = permutation[i]; // Get global index into mux data
     459       
     460
     461        printf("P4.3, i=%d\n", i);                     
     462        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
     463        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
     464        sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
     465        sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
     466        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
     467    }
     468
     469    printf("P5\n");               
    403470    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
    404471
     
    415482            //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
    416483            fprintf(stderr, "cannot open file %s\n", muxFileName);
    417             exit(-1); 
     484            return NULL;                           
    418485        }
    419486
    420487        if (verbose){
    421             printf("Reading mux file %s\n",muxFileName);
     488            printf("Reading mux file %s\n", muxFileName);
    422489        }
    423490
     
    430497        fclose(fp);
    431498
    432         /* loop over stations */
    433         /* This is where we should only take stations with indices in the permutation array
    434            (and in that order) */
    435        
    436         //for i in range(len(permution)):
    437         //     ista = permutation[i]
     499        // loop over stations present in the permutation array
    438500        //     use ista with mux data
    439501        //     use i with the processed data to be returned         
    440         i=0;
    441502       
    442         for(ista = 0; ista < nsta0; ista++)
     503        for(i = 0; i < number_of_selected_stations; i++)
    443504        {               
     505       
     506            ista = permutation[i]; // Get global index into mux data   
     507           
    444508            /* fill the data0 array from the mux file, and weight it */
    445             fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros + isrc*nsta0, lros + isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
     509            fillDataArray(ista,
     510                          nsta0,
     511                          mytgs0[ista].nt,
     512                          mytgs0[ista].ig,
     513                          fros + isrc*nsta0,
     514                          lros + isrc*nsta0,
     515                          temp_sts_data,
     516                          &istart,
     517                          &istop,
     518                          muxData);
    446519
    447520            /* weight appropriately and add */
    448521            for(k = 0; k < mytgs0[ista].nt; k++)
    449522            {
    450                 if((isdata(sts_data[ista][k])) && isdata(temp_sts_data[k]))
     523                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
    451524                {
    452                     sts_data[ista][k] += temp_sts_data[k] * weights[isrc];
     525                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
    453526                    //printf("%d,%d,%f\n",ista,k,sts_data[ista][k]);
    454527                }
    455528                else
    456529                {
    457                     sts_data[ista][k] = NODATA;
     530                    sts_data[i][k] = NODATA;
    458531                }
    459532            }
     
    461534    }
    462535
     536    printf("P6\n");                   
    463537    free(muxData);
    464538    free(temp_sts_data);
Note: See TracChangeset for help on using the changeset viewer.