diff --git a/src/grib_accessor_class_data_complex_packing.c b/src/grib_accessor_class_data_complex_packing.c index 5fcd7eb06..0b3fb999e 100644 --- a/src/grib_accessor_class_data_complex_packing.c +++ b/src/grib_accessor_class_data_complex_packing.c @@ -221,8 +221,7 @@ static int value_count(grib_accessor* a,long* count) return ret; } - -static int unpack_double(grib_accessor* a, double* val, size_t *len) +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); @@ -446,6 +445,224 @@ static int unpack_double(grib_accessor* a, double* val, size_t *len) return ret; } +static int unpack_double_optimised(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; + + grib_context_free(a->context,scals); + + return ret; +} #define MAXVAL(a,b) a>b?a:b @@ -552,9 +769,8 @@ static double calculate_pfactor(grib_context *ctx,const double* spectralField, l return pFactor; } - -static int pack_double(grib_accessor* a, const double* val, size_t *len) { - +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); @@ -647,7 +863,6 @@ static int pack_double(grib_accessor* a, const double* val, size_t *len) { self->dirty=1; - switch (ieee_floats) { case 0: encode_float =grib_ibm_to_long; @@ -725,7 +940,6 @@ static int pack_double(grib_accessor* a, const double* val, size_t *len) { for(i=1;icontext,scals); return ret; - +} + +static int pack_double_optimised(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); + + 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; + d = grib_power(decimal_scale_factor,10) ; + 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); + } + + 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++]) * 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++]) * scals[lup]); + if(current_val > max) max = current_val; + if(current_val < min) min = current_val; + + lup++; + } + maxv--; + hcount=0; + mmax++; + } + + Assert(optimize_scaling_factor); + ret = grib_optimize_decimal_factor (a, self->reference_value, + max, min, bits_per_value, 0, 1, + &decimal_scale_factor, + &binary_scale_factor, + &reference_value); + if (ret !=GRIB_SUCCESS) { + grib_context_log(gh->context,GRIB_LOG_ERROR, + "unable to find nearest_smaller_value of %g for %s",min,self->reference_value); + return GRIB_INTERNAL_ERROR; + } + + d = grib_power(+decimal_scale_factor,10); + s = grib_power(- binary_scale_factor, 2); + + 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; + if((ret = grib_set_long_internal(gh,self->decimal_scale_factor,decimal_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; +} + +/* 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); + } } diff --git a/src/grib_accessor_class_data_g1second_order_general_extended_packing.c b/src/grib_accessor_class_data_g1second_order_general_extended_packing.c index 6b05dba22..b42f296b0 100644 --- a/src/grib_accessor_class_data_g1second_order_general_extended_packing.c +++ b/src/grib_accessor_class_data_g1second_order_general_extended_packing.c @@ -204,13 +204,23 @@ static void init_class(grib_accessor_class* c) #define MAX_NUMBER_OF_GROUPS 65534 #define EFDEBUG 0 -static unsigned long nbits[32]={ - 0x1, 0x2, 0x4, 0x8, 0x10, 0x20, - 0x40, 0x80, 0x100, 0x200, 0x400, 0x800, - 0x1000, 0x2000, 0x4000, 0x8000, 0x10000, 0x20000, - 0x40000, 0x80000, 0x100000, 0x200000, 0x400000, 0x800000, - 0x1000000, 0x2000000, 0x4000000, 0x8000000, 0x10000000, 0x20000000, - 0x40000000, 0x80000000 +static unsigned long nbits[64]={ + 0x1, 0x2, 0x4, 0x8, + 0x10, 0x20, 0x40, 0x80, + 0x100, 0x200, 0x400, 0x800, + 0x1000, 0x2000, 0x4000, 0x8000, + 0x10000, 0x20000, 0x40000, 0x80000, + 0x100000, 0x200000, 0x400000, 0x800000, + 0x1000000, 0x2000000, 0x4000000, 0x8000000, + 0x10000000, 0x20000000, 0x40000000, 0x80000000, + 0x100000000, 0x200000000, 0x400000000, 0x800000000, + 0x1000000000, 0x2000000000, 0x4000000000, 0x8000000000, + 0x10000000000, 0x20000000000, 0x40000000000, 0x80000000000, + 0x100000000000, 0x200000000000, 0x400000000000, 0x800000000000, + 0x1000000000000, 0x2000000000000, 0x4000000000000, 0x8000000000000, + 0x10000000000000, 0x20000000000000, 0x40000000000000, 0x80000000000000, + 0x100000000000000, 0x200000000000000, 0x400000000000000, 0x800000000000000, + 0x1000000000000000, 0x2000000000000000, 0x4000000000000000, 0x8000000000000000 }; static long number_of_bits(grib_handle*h, unsigned long x) @@ -600,7 +610,7 @@ static void grib_split_long_groups(grib_handle* hand, grib_context* c,long* numb grib_context_free(c,localFirstOrderValues); } -static int pack_double(grib_accessor* a, const double* val, size_t *len) +static int pack_double_standard(grib_accessor* a, const double* val, size_t *len) { grib_accessor_data_g1second_order_general_extended_packing* self = (grib_accessor_data_g1second_order_general_extended_packing*)a; int ret=0; @@ -643,9 +653,6 @@ static int pack_double(grib_accessor* a, const double* val, size_t *len) long dataHeadersLength,widthsLength,lengthsLength,firstOrderValuesLength; long decimal_scale_factor; grib_handle* handle = grib_handle_of_accessor(a); - long optimize_scaling_factor = 0; - grib_context* c=handle->context; - int compat_gribex = c->gribex_mode_on && self->edition==1; self->dirty=1; @@ -1207,6 +1214,628 @@ static int pack_double(grib_accessor* a, const double* val, size_t *len) return ret; } +static int pack_double_optimised(grib_accessor* a, const double* val, size_t *len) +{ + grib_accessor_data_g1second_order_general_extended_packing* self = (grib_accessor_data_g1second_order_general_extended_packing*)a; + int ret=0; + int grib2=0; + long bits_per_value,orderOfSPD,binary_scale_factor; + long numberOfValues; + double max,min; + double decimal,divisor; + double reference_value; + size_t size,sizebits; + long half_byte; + long* X; + long* Xp; + long i; + long incrementGroupLengthA,groupWidthA,prevGroupLength,offsetD,remainingValuesB,groupLengthB; + long maxB,minB,maxAB,minAB; + long offsetBeforeData,offsetSection4; + unsigned char* buffer = NULL; + long maxWidth,maxLength,widthOfWidths,NL,widthOfLengths,N1,N2,extraValues,codedNumberOfGroups,numberOfSecondOrderPackedValues; + long pos; + + long numberOfGroups; + long groupLengthC,groupLengthA,remainingValues,count; + long maxA=0,minA=0; + long maxC,minC,offsetC; + long maxAC,minAC; + long range,bias=0,maxSPD; + long firstOrderValuesMax,offset,groupLength,j,groupWidth,firstOrderValue,lengthOfSecondOrderValues; + long *groupLengths,*groupWidths,*firstOrderValues; + /* long groupLengths[MAX_NUMBER_OF_GROUPS],groupWidths[MAX_NUMBER_OF_GROUPS],firstOrderValues[MAX_NUMBER_OF_GROUPS]; */ + + /* TODO put these parameters in def file */ + long startGroupLength=15; + long incrementGroupLength=3; + long minGroupLength=3; + long widthOfSPD=0,widthOfBias=0; + long *offsets; + long widthOfFirstOrderValues; + int computeGroupA=1; + long dataHeadersLength,widthsLength,lengthsLength,firstOrderValuesLength; + long decimal_scale_factor; + grib_handle* handle = grib_handle_of_accessor(a); + long optimize_scaling_factor = 0; + grib_context* c=handle->context; + int compat_gribex = c->gribex_mode_on && self->edition==1; + + self->dirty=1; + + numberOfValues=*len; + + min = max = val[0]; + for(i=1;i< numberOfValues;i++) { + if (val[i] > max ) max = val[i]; + if (val[i] < min ) min = val[i]; + } + + if((ret=grib_get_long_internal(handle,self->bits_per_value,&bits_per_value)) != GRIB_SUCCESS) + return ret; + + if((ret = grib_get_long_internal(handle,self->optimize_scaling_factor, &optimize_scaling_factor)) + != GRIB_SUCCESS) + return ret; + + Assert(optimize_scaling_factor); + if ((ret=grib_optimize_decimal_factor (a, self->reference_value, + max, min, bits_per_value, + compat_gribex, 1, + &decimal_scale_factor, &binary_scale_factor, &reference_value)) != GRIB_SUCCESS) + return ret; + + decimal = grib_power(decimal_scale_factor,10); + divisor = grib_power(-binary_scale_factor,2); + + min = min * decimal; + max = max * decimal; + + if((ret = grib_set_long_internal(handle,self->decimal_scale_factor, decimal_scale_factor)) != + GRIB_SUCCESS) + return ret; + + if((ret = grib_set_long_internal(handle,self->binary_scale_factor, binary_scale_factor)) != + GRIB_SUCCESS) + return ret; + + if((ret = grib_set_double_internal(handle,self->reference_value, reference_value)) != + GRIB_SUCCESS) + return ret; + + if((ret=grib_get_long_internal(handle,self->offsetdata,&offsetBeforeData)) != GRIB_SUCCESS) + return ret; + + if((ret=grib_get_long_internal(handle,self->offsetsection,&offsetSection4)) != GRIB_SUCCESS) + return ret; + + if((ret=grib_get_long_internal(handle,self->orderOfSPD,&orderOfSPD)) != GRIB_SUCCESS) + return ret; + + X=(long*)grib_context_malloc_clear(a->context,sizeof(long)*numberOfValues); + for(i=0;i< numberOfValues;i++){ + X[i] = (((val[i]*decimal)-reference_value)*divisor)+0.5; + } + + groupLengths=(long*)grib_context_malloc_clear(a->context,sizeof(long)*numberOfValues); + groupWidths=(long*)grib_context_malloc_clear(a->context,sizeof(long)*numberOfValues); + firstOrderValues=(long*)grib_context_malloc_clear(a->context,sizeof(long)*numberOfValues); + + /* spatial differencing */ + switch (orderOfSPD) { + case 1: + for (i=numberOfValues-1;i>0;i--) { + X[i]-=X[i-1]; + } + break; + case 2: + for (i=numberOfValues-1;i>1;i--) { + X[i]-=2*X[i-1]-X[i-2]; + } + break; + case 3: + for (i=numberOfValues-1;i>2;i--) { + X[i]-=3*(X[i-1]-X[i-2])+X[i-3]; + } + break; + } + if (orderOfSPD) { + Assert(orderOfSPD >=0 && orderOfSPD < numberOfValues); + bias=X[orderOfSPD]; + for (i=orderOfSPD+1;i X[i] ) bias=X[i]; + } + for (i=orderOfSPD;iX[count+i]) minA=X[count+i]; + } + } + groupWidthA=number_of_bits(handle, maxA-minA); + range=(long)grib_power(groupWidthA,2)-1; + + offsetC=count+groupLengthA; + if (offsetC==numberOfValues) { + /* no more values close group A and end loop */ + groupLengths[numberOfGroups]=groupLengthA; + groupWidths[numberOfGroups]=groupWidthA; + /* firstOrderValues[numberOfGroups]=minA; */ + /* to optimise the width of first order variable */ + firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0; + numberOfGroups++; + break; + } + + /* group C created with length=incrementGroupLength (fixed) + or remaining values if close to end + */ + groupLengthC=incrementGroupLength; + if ( groupLengthC + offsetC > numberOfValues - startGroupLength/2) { + groupLengthC=numberOfValues-offsetC; + } + maxC=X[offsetC]; + minC=X[offsetC]; + for (i=1;iX[offsetC+i]) minC=X[offsetC+i]; + } + + maxAC= maxA > maxC ? maxA : maxC; + minAC= minA < minC ? minA : minC; + + /* check if A+C can be represented with the same width as A*/ + if (maxAC-minAC > range) { + /* A could not be expanded adding C. Check if A could be expanded taking + some elements from preceding group. The condition is always that width of + A doesn't increase. + */ + if (numberOfGroups>0 && groupWidths[numberOfGroups-1] > groupWidthA ) { + prevGroupLength=groupLengths[numberOfGroups-1]-incrementGroupLength; + offsetC=count-incrementGroupLength; + /* preceding group length cannot be less than a minimum value */ + while (prevGroupLength >= minGroupLength) { + maxAC=maxA; + minAC=minA; + for (i=0;iX[offsetC+i]) minAC=X[offsetC+i]; + } + + /* no more elements can be transfered, exit loop*/ + if (maxAC-minAC > range) break; + + maxA=maxAC; + minA=minAC; + groupLengths[numberOfGroups-1]-=incrementGroupLength; + groupLengthA+=incrementGroupLength; + count-=incrementGroupLength; + remainingValues+=incrementGroupLength; + + offsetC-=incrementGroupLength; + prevGroupLength-=incrementGroupLength; + } + } + /* close group A*/ + groupLengths[numberOfGroups]=groupLengthA; + groupWidths[numberOfGroups]=groupWidthA; + /* firstOrderValues[numberOfGroups]=minA; */ + /* to optimise the width of first order variable */ + firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0; + count+=groupLengthA; + remainingValues-=groupLengthA; + numberOfGroups++; + /* incrementGroupLengthA is reset to the fixed startGroupLength as it + could have been changed after the A+C or A+B ok condition. + */ + incrementGroupLengthA=startGroupLength; + computeGroupA=1; +#if 0 + if (numberOfGroups==MAX_NUMBER_OF_GROUPS) { + groupLengthA= remainingValues ; + maxA=X[count]; + minA=X[count]; + for (i=1;iX[count+i]) minA=X[count+i]; + } + groupWidthA=number_of_bits(maxA-minA); + range=(long)grib_power(groupWidthA,2)-1; + groupLengths[numberOfGroups]=groupLengthA; + groupWidths[numberOfGroups]=groupWidthA; + firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0; + break; + } +#endif + continue; + } + + /* A+C could be coded with the same width as A*/ + offsetD=offsetC+groupLengthC; + if (offsetD==numberOfValues) { + groupLengths[numberOfGroups]=groupLengthA+groupLengthC; + groupWidths[numberOfGroups]=groupWidthA; + + /* range of AC is the same as A*/ + /* firstOrderValues[numberOfGroups]=minAC; */ + /* to optimise the width of first order variable */ + firstOrderValues[numberOfGroups] = maxAC-range > 0 ? maxAC-range : 0; + numberOfGroups++; + break; + } + + /* group B is created with length startGroupLength, starting at the + same offset as C. + */ + remainingValuesB=numberOfValues-offsetC; + groupLengthB= startGroupLength < remainingValuesB ? startGroupLength : remainingValuesB ; + maxB=maxC; + minB=minC; + for (i=groupLengthC;iX[offsetC+i]) minB=X[offsetC+i]; + } + + /* check if group B can be coded with a smaller width than A */ + if (maxB-minB <= range/2 && range>0 ) { + + /* TODO Add code to try if A can be expanded taking some elements + from the left (preceding) group. + A possible variation is to do this left check (and the previous one) + in the final loop when checking that the width of each group. + */ + + /* close group A and continue loop*/ + groupLengths[numberOfGroups]=groupLengthA; + groupWidths[numberOfGroups]=groupWidthA; + /* firstOrderValues[numberOfGroups]=minA; */ + /* to optimise the width of first order variable */ + firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0; + count+=groupLengthA; + remainingValues-=groupLengthA; + numberOfGroups++; +#if 0 + if (numberOfGroups==MAX_NUMBER_OF_GROUPS) { + groupLengthA= remainingValues ; + maxA=X[count]; + minA=X[count]; + for (i=1;iX[count+i]) minA=X[count+i]; + } + groupWidthA=number_of_bits(maxA-minA); + range=(long)grib_power(groupWidthA,2)-1; + groupLengths[numberOfGroups]=groupLengthA; + groupWidths[numberOfGroups]=groupWidthA; + firstOrderValues[numberOfGroups] = maxA-range > 0 ? maxA-range : 0; + break; + } +#endif + incrementGroupLengthA=startGroupLength; + computeGroupA=1; + continue; + } + + /* check if A+B can be coded with same with as A */ + maxAB= maxA > maxB ? maxA : maxB; + minAB= minA < minB ? minA : minB; + if (maxAB-minAB <= range) { + /* A+B can be merged. The increment used at the beginning of the loop to + build group C is increased to the size of group B + */ + incrementGroupLengthA+=groupLengthB; + maxA=maxAB; + minA=minAB; + computeGroupA=0; + continue; + } + + /* A+B cannot be merged, A+C can be merged*/ + incrementGroupLengthA+=groupLengthC; + computeGroupA=1; + + } /* end of the while*/ + + /* computing bitsPerValue as the number of bits needed to represent + the firstOrderValues. + */ + max=firstOrderValues[0]; + min=firstOrderValues[0]; + for (i=1;ifirstOrderValues[i]) min=firstOrderValues[i]; + } + widthOfFirstOrderValues=number_of_bits(handle, max-min); + firstOrderValuesMax=(long)grib_power(widthOfFirstOrderValues,2)-1; + + if (numberOfGroups>2) { + /* loop through all the groups except the last in reverse order to + check if each group width is still appropriate for the group. + Focus on groups which have been shrank as left groups of an A group taking + some of their elements. + */ + offsets=(long*)grib_context_malloc_clear(a->context,sizeof(long)*numberOfGroups); + offsets[0]=orderOfSPD; + for (i=1;i=0;i--) { + offset=offsets[i]; + groupLength=groupLengths[i]; + + if (groupLength >= startGroupLength) continue; + + max=X[offset]; + min=X[offset]; + for (j=1;jX[offset+j]) min=X[offset+j]; + } + groupWidth=number_of_bits(handle, max-min); + range=(long)grib_power(groupWidth,2)-1; + + /* width of first order values has to be unchanged.*/ + for (j=groupWidth;jrange ? max-range : 0; + if (firstOrderValue <= firstOrderValuesMax ) { + groupWidths[i]=j; + firstOrderValues[i]=firstOrderValue; + break; + } + } + + offsetC=offset; + /* group width of the current group (i) can have been reduced + and it is worth to try to expand the group to get some elements + from the left group if it has bigger width. + */ + if (i>0 && (groupWidths[i-1] > groupWidths[i]) ) { + prevGroupLength=groupLengths[i-1]-incrementGroupLength; + offsetC-=incrementGroupLength; + while (prevGroupLength >= minGroupLength) { + for (j=0;jX[offsetC+j]) min=X[offsetC+j]; + } + + /* width of first order values has to be unchanged*/ + firstOrderValue=max>range ? max-range : 0; + if (max-min > range || firstOrderValue > firstOrderValuesMax ) break; + + groupLengths[i-1]-=incrementGroupLength; + groupLengths[i]+=incrementGroupLength; + firstOrderValues[i]=firstOrderValue; + + offsetC-=incrementGroupLength; + prevGroupLength-=incrementGroupLength; + } + } + + } + grib_context_free(a->context,offsets); + } + + maxWidth=groupWidths[0]; + maxLength=groupLengths[0]; + for (i=1;iparent->h->context, GRIB_LOG_ERROR, "Cannot compute parameters for second order packing."); + return GRIB_ENCODING_ERROR; + } + widthOfWidths=number_of_bits(handle, maxWidth); + widthOfLengths=number_of_bits(handle, maxLength); + + lengthOfSecondOrderValues=0; + for ( i=0; icontext->no_big_group_split) { + grib_split_long_groups(handle, a->context,&numberOfGroups,&lengthOfSecondOrderValues, + groupLengths,&widthOfLengths,groupWidths,widthOfWidths, + firstOrderValues,widthOfFirstOrderValues); + } + + Xp=X+orderOfSPD; + for ( i=0; iwidthOfSPD, widthOfSPD)) + != GRIB_SUCCESS) + return ret; + } + + /* end writing SPD */ + if((ret = grib_set_long_internal(handle,self->widthOfFirstOrderValues, widthOfFirstOrderValues)) + != GRIB_SUCCESS) + return ret; + + dataHeadersLength=25; + if (orderOfSPD) dataHeadersLength+=1+((orderOfSPD+1)*widthOfSPD+7)/8; + widthsLength=(widthOfWidths*numberOfGroups+7)/8; + lengthsLength=(widthOfLengths*numberOfGroups+7)/8; + firstOrderValuesLength=(widthOfFirstOrderValues*numberOfGroups+7)/8; + + NL=widthsLength+dataHeadersLength+1; + N1=NL+lengthsLength; + N2=N1+firstOrderValuesLength; + + NL= NL > 65535 ? 65535 : NL; + N2= N2 > 65535 ? 65535 : N2; + N1= N1 > 65535 ? 65535 : N1; + + grib_set_long(handle,self->NL, NL); + grib_set_long(handle,self->N1, N1); + grib_set_long(handle,self->N2, N2); + + if (numberOfGroups > 65535 ) { + extraValues=numberOfGroups/65536; + codedNumberOfGroups=numberOfGroups%65536; + } else { + extraValues=0; + codedNumberOfGroups=numberOfGroups; + } + + /* if no extraValues key present it is a GRIB2*/ + grib2=0; + if((ret = grib_set_long(handle,self->extraValues, extraValues)) != GRIB_SUCCESS) { + codedNumberOfGroups=numberOfGroups; + grib2=1; + } + + if((ret = grib_set_long_internal(handle,self->codedNumberOfGroups, codedNumberOfGroups)) != GRIB_SUCCESS) + return ret; + + numberOfSecondOrderPackedValues=numberOfValues-orderOfSPD; + if (!grib2 && numberOfSecondOrderPackedValues > 65535 ) + numberOfSecondOrderPackedValues= 65535; + + if((ret = grib_set_long_internal(handle,self->numberOfSecondOrderPackedValues, numberOfSecondOrderPackedValues)) + != GRIB_SUCCESS) + return ret; + + if (grib2) { + if((ret = grib_set_long_internal(handle,self->bits_per_value, bits_per_value)) != GRIB_SUCCESS) + return ret; + } else { + if((ret = grib_set_long_internal(handle,self->bits_per_value, 0)) != GRIB_SUCCESS) + return ret; + } + + if((ret = grib_set_long_internal(handle,self->widthOfWidths, widthOfWidths)) != GRIB_SUCCESS) + return ret; + + if((ret = grib_set_long_internal(handle,self->widthOfLengths, widthOfLengths)) != GRIB_SUCCESS) + return ret; + + lengthOfSecondOrderValues=0; + for ( i=0; ihalf_byte, half_byte)) != GRIB_SUCCESS) + return ret; + + buffer=(unsigned char*)grib_context_malloc_clear(a->context,size); + + pos=0; + if (orderOfSPD) { + long SPD[4]={0,}; + size_t nSPD=orderOfSPD+1; + Assert(orderOfSPD<=3); + for (i=0;iSPD,SPD,nSPD); + if(ret) return ret; + } + + ret=grib_set_long_array_internal(handle,self->groupWidths,groupWidths,(size_t)numberOfGroups); + if(ret) return ret; + + ret=grib_set_long_array_internal(handle,self->groupLengths,groupLengths,(size_t)numberOfGroups); + if(ret) return ret; + + ret=grib_set_long_array_internal(handle,self->firstOrderValues,firstOrderValues,(size_t)numberOfGroups); + if(ret) return ret; + + Xp=X+orderOfSPD; + pos=0; + count=0; + for (i=0;i0) { + for (j=0;jparent->h,self->number_of_values, *len); + if(ret) return ret; + + grib_buffer_replace(a, buffer, size,1,1); + + grib_context_free(a->context,buffer); + grib_context_free(a->context,X); + grib_context_free(a->context,groupLengths); + grib_context_free(a->context,groupWidths); + grib_context_free(a->context,firstOrderValues); + + return ret; +} + +/* The driver pack routine. See ECC-261 */ +static int pack_double(grib_accessor* a, const double* val, size_t *len) +{ + grib_accessor_data_g1second_order_general_extended_packing* self = (grib_accessor_data_g1second_order_general_extended_packing*)a; + int ret = GRIB_SUCCESS; + grib_handle* handle = grib_handle_of_accessor(a); + long optimize_scaling_factor = 0; + + if((ret = grib_get_long_internal(handle,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); + } +} + static void destroy(grib_context* context,grib_accessor* a) { grib_accessor_data_g1second_order_general_extended_packing *self =(grib_accessor_data_g1second_order_general_extended_packing*)a;