From 663a2f1e9c2a0152b346aabec003ce59c5af2eba Mon Sep 17 00:00:00 2001 From: Shahram Najm Date: Wed, 19 Apr 2017 17:36:16 +0100 Subject: [PATCH] ECC-458: MeteoFrance contribution: GRIB spectral complex packing (Part 2) --- ...grib_accessor_class_data_complex_packing.c | 609 +----------------- 1 file changed, 2 insertions(+), 607 deletions(-) diff --git a/src/grib_accessor_class_data_complex_packing.c b/src/grib_accessor_class_data_complex_packing.c index 4dd2d16be..7380176a8 100644 --- a/src/grib_accessor_class_data_complex_packing.c +++ b/src/grib_accessor_class_data_complex_packing.c @@ -221,231 +221,7 @@ static int value_count(grib_accessor* a,long* count) return ret; } -static int unpack_double_standard(grib_accessor* a, double* val, size_t *len) -{ - grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a; - grib_handle* gh = grib_handle_of_accessor(a); - - size_t i = 0; - int ret = GRIB_SUCCESS; - long hcount = 0; - long lcount = 0; - long hpos = 0; - long lup = 0; - long mmax = 0; - long n_vals = 0; - double *scals = NULL; - double *pscals=NULL,*pval=NULL; - - double s = 0; - double d = 0; - double laplacianOperator = 0; - unsigned char* buf = NULL; - unsigned char* hres = NULL; - unsigned char* lres = NULL; - unsigned long packed_offset; - long lpos = 0; - - long maxv = 0; - long GRIBEX_sh_bug_present =0; - long ieee_floats = 0; - - long offsetdata = 0; - long bits_per_value = 0; - double reference_value = 0; - long binary_scale_factor = 0; - long decimal_scale_factor = 0; - - long sub_j= 0; - long sub_k= 0; - long sub_m= 0; - long pen_j= 0; - long pen_k= 0; - long pen_m= 0; - - double operat= 0; - int bytes; - int err=0; - - decode_float_proc decode_float = NULL; - - err=grib_value_count(a,&n_vals); - if (err) return err; - - if(*len < n_vals){ - *len = n_vals; - return GRIB_ARRAY_TOO_SMALL; - } - - if((ret = grib_get_long_internal(gh,self->offsetdata,&offsetdata)) - != GRIB_SUCCESS) return ret; - if((ret = grib_get_long_internal(gh,self->bits_per_value,&bits_per_value)) - != GRIB_SUCCESS) return ret; - if((ret = grib_get_double_internal(gh,self->reference_value,&reference_value)) - != GRIB_SUCCESS) return ret; - if((ret = grib_get_long_internal(gh,self->binary_scale_factor,&binary_scale_factor)) - != GRIB_SUCCESS) return ret; - - if((ret = grib_get_long_internal(gh,self->decimal_scale_factor,&decimal_scale_factor)) - != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->GRIBEX_sh_bug_present,&GRIBEX_sh_bug_present)) - != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->ieee_floats,&ieee_floats)) != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_double_internal(gh,self->laplacianOperator,&laplacianOperator)) - != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->sub_j,&sub_j)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->sub_k,&sub_k)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->sub_m,&sub_m)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->pen_j,&pen_j)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->pen_k,&pen_k)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->pen_m,&pen_m)) != GRIB_SUCCESS) - return ret; - - self->dirty=0; - - switch (ieee_floats) { - case 0: - decode_float=grib_long_to_ibm; - bytes=4; - break; - case 1: - decode_float=grib_long_to_ieee; - bytes=4; - break; - case 2: - decode_float=grib_long_to_ieee64; - bytes=8; - break; - default: - return GRIB_NOT_IMPLEMENTED; - } - - Assert (sub_j == sub_k); - Assert (sub_j == sub_m); - Assert (pen_j == pen_k); - Assert (pen_j == pen_m); - - buf = (unsigned char*)gh->buffer->data; - - maxv = pen_j+1; - - buf += grib_byte_offset(a); - hres = buf; - lres = buf; - - if (pen_j == sub_j) { - n_vals = (pen_j+1)*(pen_j+2); - d = grib_power(-decimal_scale_factor,10) ; - grib_ieee_decode_array(a->context,buf,n_vals,bytes,val); - if (d) { - for (i=0;icontext,maxv*sizeof(double)); - Assert(scals); - - scals[0] = 0; - for(i=1;icontext,GRIB_LOG_WARNING, - "COMPLEX_PACKING : problem with operator div by zero at index %d of %d \n", - i , maxv); - scals[i] = 0; - } - } - - /* - printf("UNPACKING LAPLACE=%.20f\n",laplacianOperator); - - printf("packed offset=%ld\n",packed_offset); - for(i=0;i0) - { - lup=mmax; - if(sub_k>=0) - { - for(hcount=0;hcount= i); - *len = i; - - if(d != 1) { - for(i=0;i<*len;i++) - val[i++] *= d; - } - - grib_context_free(a->context,scals); - - return ret; -} - -static int unpack_double_optimised(grib_accessor* a, double* val, size_t *len) +static int unpack_double(grib_accessor* a, double* val, size_t *len) { grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a; grib_handle* gh = grib_handle_of_accessor(a); @@ -769,353 +545,7 @@ static double calculate_pfactor(grib_context *ctx,const double* spectralField, l return pFactor; } -static int pack_double_standard(grib_accessor* a, const double* val, size_t *len) -{ - grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a; - grib_handle* gh = grib_handle_of_accessor(a); - - size_t i = 0; - int ret = GRIB_SUCCESS; - long hcount = 0; - long lcount = 0; - long hpos = 0; - long lup = 0; - long mmax = 0; - long n_vals = 0; - double *scals = NULL; - - double s = 0; - double d = 0; - - unsigned char* buf = NULL; - - size_t buflen = 0; - size_t hsize = 0; - size_t lsize = 0; - - unsigned char* hres = NULL; - unsigned char* lres = NULL; - - long lpos = 0; - long maxv = 0; - - long offsetdata = 0; - long bits_per_value = 0; - double reference_value = 0; - long binary_scale_factor = 0; - long decimal_scale_factor = 0; - long optimize_scaling_factor = 0; - long laplacianOperatorIsSet = 0; - - double laplacianOperator = 0; - long sub_j= 0; - long sub_k= 0; - long sub_m= 0; - long pen_j= 0; - long pen_k= 0; - long pen_m= 0; - long GRIBEX_sh_bug_present =0; - long ieee_floats =0; - double min = 0; - double max = 0; - double current_val = 0; - short mixmax_unset = 0; - int bytes; - - encode_float_proc encode_float = NULL; - - if (*len ==0) return GRIB_NO_VALUES; - - if((ret = grib_get_long_internal(gh,self->offsetdata,&offsetdata)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->bits_per_value,&bits_per_value)) != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->decimal_scale_factor,&decimal_scale_factor)) != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->optimize_scaling_factor,&optimize_scaling_factor)) != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->GRIBEX_sh_bug_present,&GRIBEX_sh_bug_present)) != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->ieee_floats,&ieee_floats)) != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->laplacianOperatorIsSet,&laplacianOperatorIsSet)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_double_internal(gh,self->laplacianOperator,&laplacianOperator)) != GRIB_SUCCESS) - return ret; - - if((ret = grib_get_long_internal(gh,self->sub_j,&sub_j)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->sub_k,&sub_k)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->sub_m,&sub_m)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->pen_j,&pen_j)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->pen_k,&pen_k)) != GRIB_SUCCESS) - return ret; - if((ret = grib_get_long_internal(gh,self->pen_m,&pen_m)) != GRIB_SUCCESS) - return ret; - - self->dirty=1; - - switch (ieee_floats) { - case 0: - encode_float =grib_ibm_to_long; - bytes=4; - break; - case 1: - encode_float =grib_ieee_to_long; - bytes=4; - break; - case 2: - encode_float =grib_ieee64_to_long; - bytes=8; - break; - default: - return GRIB_NOT_IMPLEMENTED; - } - - Assert (sub_j == sub_k); Assert( sub_j == sub_m); - Assert (pen_j == pen_k); Assert( pen_j == pen_m); - - n_vals = (pen_j+1)*(pen_j+2); - d = grib_power(decimal_scale_factor,10) ; - - if(*len != n_vals){ - grib_context_log(a->context,GRIB_LOG_ERROR,"COMPLEX_PACKING : wrong number of values, expected %d - got %d",n_vals,*len); - return GRIB_INTERNAL_ERROR; - } - - if (pen_j == sub_j) { - double* values; - if (d) { - values=(double*)grib_context_malloc_clear(a->context,sizeof(double)*n_vals); - for (i=0;icontext,buflen); - grib_ieee_encode_array(a->context,values,n_vals,bytes,buf); - if (d) grib_context_free(a->context,values); - grib_buffer_replace(a, buf, buflen,1,1); - grib_context_free(a->context,buf); - return 0; - } - - if(!laplacianOperatorIsSet) { - laplacianOperator = calculate_pfactor(a->context,val,pen_j,sub_j); - if((ret = grib_set_double_internal(gh,self->laplacianOperator,laplacianOperator)) - != GRIB_SUCCESS) return ret; - grib_get_double_internal(gh,self->laplacianOperator,&laplacianOperator); - } - - /* - printf("PACKING LAPLACE set=%ld value=%.20f\n",laplacianOperatorIsSet,laplacianOperator); - */ - - hsize = 4*(sub_k+1)*(sub_k+2); - lsize = ((n_vals - ((sub_k+1)*(sub_k+2)))*bits_per_value)/8; - - buflen = hsize+lsize; - - buf = (unsigned char*)grib_context_malloc(a->context,buflen); - hres = buf; - lres = buf+hsize; - - maxv = pen_j+1; - - lpos = 0; - hpos = 0; - - scals = (double*) grib_context_malloc(a->context,maxv*sizeof(double)); - Assert(scals); - - scals[0] =0; - for(i=1;i0) - { - lup=mmax; - - if(sub_k>=0) - { - i += 2*(sub_k+1); - lup += sub_k+1 ; - hcount += sub_k+1 ; - sub_k--; - } - - for(lcount=hcount; lcount < maxv ; lcount++) - { - current_val = ((val[i++]*d) * scals[lup]); - if(mixmax_unset == 0){ - max = current_val; - min = current_val; - mixmax_unset = 1; - } - - if(current_val > max) max = current_val; - if(current_val < min) min = current_val; - - current_val = ((val[i++]*d) * scals[lup]); - if(current_val > max) max = current_val; - if(current_val < min) min = current_val; - - lup++; - } - maxv--; - hcount=0; - mmax++; - } - - if (grib_get_nearest_smaller_value(gh,self->reference_value,min,&reference_value) - !=GRIB_SUCCESS) { - grib_context_log(a->context,GRIB_LOG_ERROR, - "unable to find nearest_smaller_value of %g for %s",min,self->reference_value); - return GRIB_INTERNAL_ERROR; - } - binary_scale_factor = grib_get_binary_scale_fact(max,reference_value,bits_per_value,&ret); - - if (ret==GRIB_UNDERFLOW) { - d=0; - binary_scale_factor = 0; - reference_value=0; - } else { - if (ret!=GRIB_SUCCESS) { - grib_context_log(a->context,GRIB_LOG_ERROR,"COMPLEX_PACKING : Cannot compute binary_scale_factor"); - return ret; - } - } - s = grib_power(-binary_scale_factor,2); - - /* printf("D : %.30f\n",d); */ - - i=0; - - mmax = 0; - maxv = pen_j+1; - i=0; - lcount=0; - hcount=0; - sub_k = sub_j; - - while(maxv>0) - { - lup=mmax; - - if(sub_k>=0) - { - for(hcount=0;hcountcontext,GRIB_LOG_ERROR, - "COMPLEX_PACKING : negative coput before packing (%g)", current_val); - grib_encode_unsigned_longb(lres, current_val, &lpos, bits_per_value); - - current_val = (((((val[i++]*d) * scals[lup])-reference_value)*s)+0.5); - if(current_val < 0) - grib_context_log(a->context,GRIB_LOG_ERROR, - "COMPLEX_PACKING : negative coput before packing (%g)", current_val); - grib_encode_unsigned_longb(lres, current_val, &lpos, bits_per_value); - lup++; - } - } else { - for(lcount=hcount; lcount < maxv ; lcount++) - { - current_val = (((((val[i++]*d) * scals[lup])-reference_value)*s)+0.5); - if(current_val < 0) - grib_context_log(a->context,GRIB_LOG_ERROR, - "COMPLEX_PACKING : negative coput before packing (%g)", current_val); - grib_encode_unsigned_long(lres, current_val, &lpos, bits_per_value); - - current_val = (((((val[i++]*d) * scals[lup])-reference_value)*s)+0.5); - if(current_val < 0) - grib_context_log(a->context,GRIB_LOG_ERROR, - "COMPLEX_PACKING : negative coput before packing (%g)", current_val); - grib_encode_unsigned_long(lres, current_val, &lpos, bits_per_value); - lup++; - } - } -#endif - - maxv--; - hcount=0; - mmax++; - } - - if(((hpos/8) != hsize) &&((lpos/8) != lsize)) - { - grib_context_log(a->context,GRIB_LOG_ERROR, - "COMPLEX_PACKING : Mismatch in packing between high resolution and low resolution part"); - grib_context_free(a->context,buf); - grib_context_free(a->context,scals); - return GRIB_INTERNAL_ERROR; - } - - buflen = ((hpos + lpos)/8); - - if((ret = grib_set_double_internal(gh,self->reference_value, reference_value)) != GRIB_SUCCESS) - return ret; - { - /* Make sure we can decode it again */ - double ref = 1e-100; - grib_get_double_internal(gh,self->reference_value,&ref); - Assert(ref == reference_value); - } - - if((ret = grib_set_long_internal(gh,self->binary_scale_factor, binary_scale_factor)) != GRIB_SUCCESS) - return ret; - - grib_buffer_replace(a, buf, buflen,1,1); - grib_context_free(a->context,buf); - grib_context_free(a->context,scals); - - return ret; -} - -static int pack_double_optimised(grib_accessor* a, const double* val, size_t *len) +static int pack_double(grib_accessor* a, const double* val, size_t *len) { grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a; grib_handle* gh = grib_handle_of_accessor(a); @@ -1474,38 +904,3 @@ static int pack_double_optimised(grib_accessor* a, const double* val, size_t *le return ret; } - -/* The driver unpack and pack routines. See ECC-261 */ -static int unpack_double(grib_accessor* a, double* val, size_t *len) -{ - grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a; - grib_handle* gh = grib_handle_of_accessor(a); - int ret = GRIB_SUCCESS; - long optimize_scaling_factor = 0; - - if ((ret = grib_get_long_internal(gh, self->optimize_scaling_factor, &optimize_scaling_factor)) != GRIB_SUCCESS) - return ret; - - if (optimize_scaling_factor) { - return unpack_double_optimised(a, val, len); - } else { - return unpack_double_standard(a, val, len); - } -} - -static int pack_double(grib_accessor* a, const double* val, size_t *len) -{ - grib_accessor_data_complex_packing* self = (grib_accessor_data_complex_packing*)a; - grib_handle* gh = grib_handle_of_accessor(a); - int ret = GRIB_SUCCESS; - long optimize_scaling_factor = 0; - - if ((ret = grib_get_long_internal(gh, self->optimize_scaling_factor, &optimize_scaling_factor)) != GRIB_SUCCESS) - return ret; - - if (optimize_scaling_factor) { - return pack_double_optimised(a, val, len); - } else { - return pack_double_standard(a, val, len); - } -}