Changeset 5473 for anuga_core/source/anuga/shallow_water/urs_ext.c
- Timestamp:
- Jul 7, 2008, 3:34:26 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/urs_ext.c
r5471 r5473 23 23 char isdata(float); 24 24 void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int); 25 PyArrayObject *_read_mux2(int, char **, float *, double *, int);25 PyArrayObject *_read_mux2(int, char **, float *, double *, PyObject *, int); 26 26 27 27 PyObject *read_mux2(PyObject *self, PyObject *args){ … … 29 29 30 30 Python call: 31 read_mux2(numSrc,filenames,weights,file_params, verbose)31 read_mux2(numSrc,filenames,weights,file_params,permutation,verbose) 32 32 33 33 NOTE: … … 38 38 PyArrayObject *pyweights,*file_params; 39 39 PyArrayObject *sts_data; 40 PyObject *permutation; 40 41 PyObject *fname; 41 42 … … 51 52 52 53 // Convert Python arguments to C 53 if (!PyArg_ParseTuple(args, "iOOO i",54 &numSrc,&filenames,&pyweights,&file_params,& verbose)) {54 if (!PyArg_ParseTuple(args, "iOOOOi", 55 &numSrc,&filenames,&pyweights,&file_params,&permutation,&verbose)) { 55 56 56 57 PyErr_SetString(PyExc_RuntimeError, … … 82 83 muxFileNameArray = (char **) malloc((int)numSrc*sizeof(char *)); 83 84 if (muxFileNameArray == NULL) { 84 printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");85 exit(-1);85 PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated"); 86 return NULL; 86 87 } 87 88 for (i=0;i<PyList_Size(filenames);i++){ 88 89 muxFileNameArray[i] = (char *) malloc((MAX_FILE_NAME_LENGTH+1)*sizeof(char)); 89 90 if (muxFileNameArray[i] == NULL) { 90 printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");91 exit(-1);91 PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated"); 92 return NULL; 92 93 } 93 94 fname=PyList_GetItem(filenames, i); 94 95 if (!PyString_Check(fname)) { 95 96 PyErr_SetString(PyExc_ValueError, "filename not a string"); 97 return NULL; 96 98 } 97 99 muxFileNameArray[i]=PyString_AsString(fname); … … 111 113 112 114 // Read in mux2 data from file 113 sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose); 114 115 sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,permutation,(int)verbose); 115 116 116 117 // Allocate space for return vector … … 125 126 } 126 127 127 PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)128 PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, PyObject *permutation, int verbose) 128 129 { 129 130 FILE *fp; 130 int nsta, nsta0, i, t, isrc, ista;131 int nsta, nsta0, i, j, t, isrc, ista, N; 131 132 struct tgsrwg *mytgs, *mytgs0; 132 133 char *muxFileName; … … 144 145 float *temp_sts_data; 145 146 PyArrayObject *sts_data; 147 148 int permuation_dimensions[1]; 149 PyArrayObject *permutation_array; 146 150 147 151 long int offset; … … 174 178 /* Open the first muxfile */ 175 179 if((fp=fopen(muxFileNameArray[0],"r"))==NULL){ 176 fprintf(stderr, "cannot open file %s\n", muxFileNameArray[0]); 177 exit(-1); 180 fprintf(stderr,"Cannot open file %s\n", muxFileNameArray[0]); 181 PyErr_SetString(PyExc_RuntimeError,""); 182 return NULL; 178 183 } 179 184 … … 181 186 /* First read the number of stations*/ 182 187 fread(&nsta0,sizeof(int),1,fp); 188 189 if (permutation == Py_None){ 190 // if no permutation is specified return the times series of all the gauges in the mux file 191 permuation_dimensions[0]=nsta0; 192 permutation_array=(PyArrayObject *) PyArray_FromDims(1, permuation_dimensions, PyArray_INT); 193 for (ista=0; ista<nsta0; ista++){ 194 *(long *) (permutation_array -> data+ista*permutation_array->strides[0])=ista; 195 } 196 }else{ 197 // Specifies the gauge numbers that for which data is to be extracted 198 permutation_array=(PyArrayObject *) PyArray_ContiguousFromObject(permutation,PyArray_INT,1,1); 199 } 200 if(permutation_array == NULL){ 201 PyErr_SetString(PyExc_RuntimeError,"ERROR: Memory for permutation_array array could not be allocated"); 202 return NULL; 203 } 183 204 184 205 /*now allocate space for, and read in, the structures for each station*/ … … 203 224 if (numData < 0) { 204 225 fprintf(stderr,"Size of data block appears to be negative!\n"); 205 exit(-1); 226 PyErr_SetString(PyExc_RuntimeError,""); 227 return NULL; 206 228 } 207 229 fclose(fp); … … 229 251 { 230 252 fprintf(stderr, "cannot open file %s\n", muxFileName); 231 exit(-1); 253 PyErr_SetString(PyExc_RuntimeError,""); 254 return NULL; 232 255 } 233 256 … … 236 259 if(nsta != nsta0){ 237 260 fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]); 261 PyErr_SetString(PyExc_RuntimeError,""); 238 262 fclose(fp); 239 exit(-1);263 return NULL; 240 264 } 241 265 fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp); … … 244 268 fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]); 245 269 fclose(fp); 246 exit(-1);270 return NULL; 247 271 } 248 272 if(mytgs[ista].nt != num_ts){ 249 273 fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]); 274 PyErr_SetString(PyExc_RuntimeError,""); 250 275 fclose(fp); 251 exit(-1);276 return NULL; 252 277 } 253 278 } … … 263 288 if (numData < 0){ 264 289 fprintf(stderr,"Size of data block appears to be negative!\n"); 265 exit(-1); 290 PyErr_SetString(PyExc_RuntimeError,""); 291 return NULL; 266 292 } 267 293 fclose(fp); 268 294 } 269 params[0]=(double)nsta0; 295 N = permutation_array->dimensions[0]; 296 297 params[0]=(double)N; 270 298 params[1]=(double)mytgs0[0].dt; 271 params[2]=(double) mytgs0[0].nt;299 params[2]=(double)num_ts; 272 300 273 301 274 302 // Find min and max start times of all gauges 303 // FIXME: ONLY CHECK gauges in permutation for start and finsh times 275 304 start_tstep=num_ts+1; 276 305 finish_tstep=-1; … … 285 314 286 315 if ((start_tstep>num_ts) | (finish_tstep < 0)){ 287 printf("ERROR: Gauge data has incorrect start and finsh times\n");316 PyErr_SetString(PyExc_RuntimeError,"Gauge data has incorrect start and finsh times"); 288 317 return NULL; 289 318 } 290 319 291 320 if (start_tstep>=finish_tstep){ 292 printf("ERROR: Gauge data has non-postive_length\n");321 PyErr_SetString(PyExc_RuntimeError,"Gauge data has non-postive_length"); 293 322 return NULL; 294 323 } … … 296 325 /* Make array(s) to hold the demuxed data */ 297 326 len_sts_data=num_ts+POFFSET; 298 dimensions[0]= nsta0;327 dimensions[0]=N; 299 328 dimensions[1]=len_sts_data; 300 329 sts_data = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE); 301 330 if(sts_data == NULL){ 302 printf("ERROR: Memory for pydata array could not be allocated.\n");331 PyErr_SetString(PyExc_RuntimeError,"Memory for pydata array could not be allocated"); 303 332 return NULL; 304 333 } 305 334 306 335 /* Initialise sts data to zero */ 307 for (ista=0;ista<nsta0;ista++){ 336 for(i=0; i<N; i++){ 337 ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]); 308 338 for (t=0;t<num_ts;t++){ 309 *(double *) (sts_data -> data+i sta*sts_data->strides[0]+t*sts_data->strides[1])=0.0;310 } 311 *(double *) (sts_data -> data+i sta*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;312 *(double *) (sts_data -> data+i sta*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;313 *(double *) (sts_data -> data+i sta*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;314 *(double *) (sts_data -> data+i sta*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];315 *(double *) (sts_data -> data+i sta*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];339 *(double *) (sts_data -> data+i*sts_data->strides[0]+t*sts_data->strides[1])=0.0; 340 } 341 *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat; 342 *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon; 343 *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z; 344 *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista]; 345 *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista]; 316 346 } 317 347 … … 327 357 if((fp=fopen(muxFileName,"r"))==NULL) 328 358 { 329 //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);330 fprintf(stderr, "cannot open file %s\n", muxFileName);331 exit(-1);359 fprintf(stderr, "Cannot open file %s\n", muxFileName); 360 PyErr_SetString(PyExc_RuntimeError,""); 361 return NULL; 332 362 } 333 363 … … 345 375 346 376 /* loop over stations */ 347 for(ista=0; ista<nsta0; ista++){ 377 for(j=0; j<N; j++){ 378 ista = *(long *) (permutation_array -> data+j*permutation_array->strides[0]); 379 //printf("%d\n",(int)*(long *) (permutation_array -> data+j*permutation_array->strides[0])); 348 380 349 381 /* fill the data0 array from the mux file, and weight it */ … … 351 383 352 384 /* 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];385 for(t=0; t<num_ts; t++){ 386 if((isdata(*(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]))) && isdata(temp_sts_data[t])){ 387 *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) += temp_sts_data[t] * weights[isrc]; 356 388 }else{ 357 *(double *) (sts_data -> data+ ista*sts_data->strides[0]+i*sts_data->strides[1]) = 0.0;389 *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) = 0.0; 358 390 } 359 391 } 360 392 } 361 } 393 } 362 394 free(muxData); 363 395 free(temp_sts_data); 396 free(fros); 397 free(lros); 364 398 return sts_data; 365 399 }
Note: See TracChangeset
for help on using the changeset viewer.