mirror of https://github.com/ecmwf/eccodes.git
ECC-261: Reduce packing error by optimizing scaling factor (Part 2)
This commit is contained in:
Shahram Najm
parent
ea7ca7c180
commit
9acee6d01f
|
|
@ -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;i<n_vals;i++) val[i]*=d;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
packed_offset = grib_byte_offset(a) + 4*(sub_k+1)*(sub_k+2);
|
||||
|
||||
lpos = 8*(packed_offset-offsetdata);
|
||||
|
||||
s = grib_power(binary_scale_factor,2);
|
||||
d = grib_power(-decimal_scale_factor,10) ;
|
||||
|
||||
scals = (double*)grib_context_malloc(a->context,maxv*sizeof(double));
|
||||
Assert(scals);
|
||||
|
||||
scals[0] = 0;
|
||||
for(i=1;i<maxv;i++){
|
||||
operat = pow(i*(i+1),laplacianOperator);
|
||||
if(operat != 0)
|
||||
scals[i] = (1.0/operat);
|
||||
else{
|
||||
grib_context_log(a->context,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;i<maxv;i++)
|
||||
printf("scals[%d]=%g\n",i,scals[i]);*/
|
||||
|
||||
i=0;
|
||||
|
||||
while(maxv>0)
|
||||
{
|
||||
lup=mmax;
|
||||
if(sub_k>=0)
|
||||
{
|
||||
for(hcount=0;hcount<sub_k+1;hcount++)
|
||||
{
|
||||
val[i++] = decode_float(grib_decode_unsigned_long(hres,&hpos,32));
|
||||
val[i++] = decode_float(grib_decode_unsigned_long(hres,&hpos,32));
|
||||
|
||||
if (GRIBEX_sh_bug_present && hcount==sub_k){
|
||||
/* bug in ecmwf data, last row (K+1)is scaled but should not */
|
||||
val[i-2] *= scals[lup];
|
||||
val[i-1] *= scals[lup];
|
||||
}
|
||||
lup++;
|
||||
}
|
||||
sub_k--;
|
||||
}
|
||||
|
||||
pscals=scals+lup;
|
||||
pval=val+i;
|
||||
#if FAST_BIG_ENDIAN
|
||||
grib_decode_double_array_complex(lres,
|
||||
&lpos,bits_per_value,
|
||||
reference_value,s,pscals,(maxv-hcount)*2,pval);
|
||||
i+=(maxv-hcount)*2;
|
||||
#else
|
||||
(void)pscals; /* suppress gcc warning */
|
||||
(void)pval; /* suppress gcc warning */
|
||||
for(lcount=hcount; lcount < maxv ; lcount++)
|
||||
{
|
||||
val[i++] = d * (double) ((grib_decode_unsigned_long(lres, &lpos,
|
||||
bits_per_value)*s)+reference_value)*scals[lup];
|
||||
val[i++] = d * (double) ((grib_decode_unsigned_long(lres, &lpos,
|
||||
bits_per_value)*s)+reference_value)*scals[lup];
|
||||
lup++;
|
||||
}
|
||||
#endif
|
||||
|
||||
maxv--;
|
||||
hcount=0;
|
||||
mmax++;
|
||||
}
|
||||
|
||||
Assert(*len >= 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;i<maxv;i++)
|
||||
scals[i] = ((double)pow(i*(i+1),laplacianOperator));
|
||||
|
||||
|
||||
i=0;
|
||||
|
||||
mmax = 0;
|
||||
|
|
@ -899,5 +1113,376 @@ static int pack_double(grib_accessor* a, const double* val, size_t *len) {
|
|||
grib_context_free(a->context,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;i<n_vals;i++) values[i]=val[i]*d;
|
||||
} else {
|
||||
values=(double*)val;
|
||||
}
|
||||
buflen=n_vals*bytes;
|
||||
buf = (unsigned char*)grib_context_malloc_clear(a->context,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;i<maxv;i++)
|
||||
scals[i] = ((double)pow(i*(i+1),laplacianOperator));
|
||||
|
||||
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)
|
||||
{
|
||||
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;hcount<sub_k+1;hcount++)
|
||||
{
|
||||
if ( GRIBEX_sh_bug_present && hcount==sub_k ) {
|
||||
/* _test(val[i]*d*scals[lup],1); */
|
||||
grib_encode_unsigned_long(hres, encode_float((val[i++]*d)*scals[lup]) , &hpos, 32);
|
||||
/* _test(val[i]*d*scals[lup],1); */
|
||||
grib_encode_unsigned_long(hres, encode_float((val[i++]*d)*scals[lup]) , &hpos, 32);
|
||||
}else{
|
||||
|
||||
/* _test(val[i]*d,0); */
|
||||
|
||||
grib_encode_unsigned_long(hres, encode_float(val[i++]) , &hpos, 32);
|
||||
/* _test(val[i]*d,0); */
|
||||
grib_encode_unsigned_long(hres, encode_float(val[i++]) , &hpos, 32);
|
||||
}
|
||||
lup++;
|
||||
}
|
||||
sub_k--;
|
||||
}
|
||||
|
||||
#if FAST_BIG_ENDIAN
|
||||
grib_encode_double_array_complex((maxv-hcount)*2,&(val[i]),bits_per_value,reference_value,&(scals[lup]),d,s,lres,&lpos);
|
||||
i+=(maxv-hcount)*2;
|
||||
#else
|
||||
if (bits_per_value % 8) {
|
||||
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_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);
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<numberOfValues;i++) {
|
||||
if ( bias > X[i] ) bias=X[i];
|
||||
}
|
||||
for (i=orderOfSPD;i<numberOfValues;i++) {
|
||||
X[i]-=bias;
|
||||
}
|
||||
maxSPD=X[0];
|
||||
for (i=1;i<orderOfSPD;i++) {
|
||||
if ( maxSPD < X[i] ) maxSPD=X[i];
|
||||
}
|
||||
/* widthOfSPD=(long)ceil(log((double)(maxSPD+1))/log(2.0)); */
|
||||
widthOfSPD=number_of_bits(handle, maxSPD);
|
||||
widthOfBias=number_of_bits(handle, labs(bias))+1;
|
||||
|
||||
if ( widthOfSPD < widthOfBias ) widthOfSPD=widthOfBias;
|
||||
|
||||
}
|
||||
/* end of spatial differencing */
|
||||
|
||||
count=orderOfSPD;
|
||||
remainingValues=numberOfValues-count;
|
||||
numberOfGroups=0;
|
||||
incrementGroupLengthA=startGroupLength;
|
||||
|
||||
computeGroupA=1;
|
||||
while (remainingValues) {
|
||||
/* group A created with length=incrementGroupLengthA (if enough values remain)
|
||||
incrementGroupLengthA=startGroupLength always except when coming from an A+C or A+B ok branch
|
||||
*/
|
||||
groupLengthA= incrementGroupLengthA < remainingValues ? incrementGroupLengthA : remainingValues ;
|
||||
if (computeGroupA) {
|
||||
maxA=X[count];
|
||||
minA=X[count];
|
||||
for (i=1;i<groupLengthA;i++) {
|
||||
DebugAssertAccess(X, count+i, numberOfValues);
|
||||
if (maxA<X[count+i]) maxA=X[count+i];
|
||||
if (minA>X[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;i<groupLengthC;i++) {
|
||||
DebugAssertAccess(X, offsetC+i, numberOfValues);
|
||||
if (maxC<X[offsetC+i]) maxC=X[offsetC+i];
|
||||
if (minC>X[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;i<incrementGroupLength;i++) {
|
||||
if (maxAC<X[offsetC+i]) maxAC=X[offsetC+i];
|
||||
if (minAC>X[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;i<groupLengthA;i++) {
|
||||
if (maxA<X[count+i]) maxA=X[count+i];
|
||||
if (minA>X[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;i<groupLengthB;i++) {
|
||||
if (maxB<X[offsetC+i]) maxB=X[offsetC+i];
|
||||
if (minB>X[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;i<groupLengthA;i++) {
|
||||
if (maxA<X[count+i]) maxA=X[count+i];
|
||||
if (minA>X[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;i<numberOfGroups;i++) {
|
||||
if (max<firstOrderValues[i]) max=firstOrderValues[i];
|
||||
if (min>firstOrderValues[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<numberOfGroups;i++) offsets[i]=offsets[i-1]+groupLengths[i-1];
|
||||
for (i=numberOfGroups-2;i>=0;i--) {
|
||||
offset=offsets[i];
|
||||
groupLength=groupLengths[i];
|
||||
|
||||
if (groupLength >= startGroupLength) continue;
|
||||
|
||||
max=X[offset];
|
||||
min=X[offset];
|
||||
for (j=1;j<groupLength;j++) {
|
||||
if (max<X[offset+j]) max=X[offset+j];
|
||||
if (min>X[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;j<groupWidths[i];j++) {
|
||||
firstOrderValue= max>range ? 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;j<incrementGroupLength;j++) {
|
||||
if (max<X[offsetC+j]) max=X[offsetC+j];
|
||||
if (min>X[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;i<numberOfGroups;i++) {
|
||||
if (maxWidth<groupWidths[i]) maxWidth=groupWidths[i];
|
||||
if (maxLength<groupLengths[i]) maxLength=groupLengths[i];
|
||||
}
|
||||
|
||||
if (maxWidth < 0 || maxLength < 0) {
|
||||
grib_context_log(a->parent->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; i<numberOfGroups;i++) {
|
||||
lengthOfSecondOrderValues+=groupLengths[i]*groupWidths[i];
|
||||
}
|
||||
|
||||
if (!a->context->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; i<numberOfGroups;i++) {
|
||||
for (j=0; j<groupLengths[i]; j++) {
|
||||
*(Xp++)-=firstOrderValues[i];
|
||||
}
|
||||
}
|
||||
|
||||
/* start writing to message */
|
||||
|
||||
/* writing SPD */
|
||||
if (orderOfSPD) {
|
||||
if((ret = grib_set_long_internal(handle,self->widthOfSPD, 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; i<numberOfGroups;i++) {
|
||||
lengthOfSecondOrderValues+=groupLengths[i]*groupWidths[i];
|
||||
}
|
||||
|
||||
size=(lengthOfSecondOrderValues+7)/8;
|
||||
sizebits=lengthOfSecondOrderValues;
|
||||
|
||||
/* padding section 4 to an even number of octets*/
|
||||
size = (size+offsetBeforeData-offsetSection4) % 2 ? size+1 : size;
|
||||
half_byte=8*size-sizebits;
|
||||
if((ret = grib_set_long_internal(handle,self->half_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;i<orderOfSPD;i++) SPD[i]=X[i];
|
||||
SPD[orderOfSPD]=bias;
|
||||
ret=grib_set_long_array_internal(handle,self->SPD,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;i<numberOfGroups;i++) {
|
||||
if (groupWidths[i]>0) {
|
||||
for (j=0;j<groupLengths[i];j++) {
|
||||
#if EFDEBUG
|
||||
printf("CXXXXX %ld %ld %ld %ld\n",count,*Xp,groupWidths[i],groupLengths[i]);
|
||||
count++;
|
||||
#endif
|
||||
grib_encode_unsigned_longb(buffer,*(Xp++),&pos,groupWidths[i]);
|
||||
}
|
||||
} else {
|
||||
Xp+=groupLengths[i];
|
||||
#if EFDEBUG
|
||||
count+=groupLengths[i];
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
/* ECC-259: Set correct number of values */
|
||||
ret=grib_set_long_internal(a->parent->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;
|
||||
|
|
|
|||
Loading…
Reference in New Issue