From 2cad3be88bdca45650e55bb952c395123c4ed916 Mon Sep 17 00:00:00 2001 From: Eugen Betke Date: Mon, 30 Jan 2023 13:38:13 +0000 Subject: [PATCH] ECC-1467: Complex packing single-precision --- src/eccodes_prototypes.h | 1 + ...grib_accessor_class_data_complex_packing.c | 227 +++++++++++++++++- src/grib_ieeefloat.c | 40 +++ tests/grib_ecc-1467.c | 105 ++++---- tests/grib_ecc-1467.sh | 43 +++- 5 files changed, 350 insertions(+), 66 deletions(-) diff --git a/src/eccodes_prototypes.h b/src/eccodes_prototypes.h index ac0793765..a4215546c 100644 --- a/src/eccodes_prototypes.h +++ b/src/eccodes_prototypes.h @@ -840,6 +840,7 @@ int grib_nearest_smaller_ieee_float(double a, double* ret); unsigned long grib_ieee64_to_long(double x); double grib_long_to_ieee64(unsigned long x); int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, double* val); +int grib_ieee_decode_array_float(grib_context* c, unsigned char* buf, size_t nvals, int bytes, float* val); int grib_ieee_encode_array(grib_context* c, double* val, size_t nvals, int bytes, unsigned char* buf); /* grib_accessor_class_reference_value_error.c */ diff --git a/src/grib_accessor_class_data_complex_packing.c b/src/grib_accessor_class_data_complex_packing.c index bf3eec285..652b7f3bf 100644 --- a/src/grib_accessor_class_data_complex_packing.c +++ b/src/grib_accessor_class_data_complex_packing.c @@ -442,6 +442,12 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len) mmax++; } + + for (int i = n_vals-10; i < n_vals; ++i) { + printf("d: %f\n", val[i]); + } + + Assert(*len >= i); *len = i; @@ -449,10 +455,229 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len) return ret; } + + + // TODO(maee): ECC-1467 static int unpack_float(grib_accessor* a, float* val, size_t* len) { - return GRIB_NOT_IMPLEMENTED; + 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; + float* scals = NULL; + float *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; + + float 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; + + /* ECC-774: don't use grib_get_long_internal */ + if ((ret = grib_get_long(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_float(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) + bytes * (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 = (float*)grib_context_malloc(a->context, maxv * sizeof(float)); + 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 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, 8 * bytes)); + val[i++] = decode_float(grib_decode_unsigned_long(hres, &hpos, 8 * bytes)); + + 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 * (float)((grib_decode_unsigned_long(lres, &lpos, bits_per_value) * s) + reference_value) * scals[lup]; + val[i++] = d * (float)((grib_decode_unsigned_long(lres, &lpos, bits_per_value) * s) + reference_value) * scals[lup]; + /* These values should always be zero, but as they are packed, + it is necessary to force them back to zero */ + if (mmax == 0) + val[i - 1] = 0; + lup++; + } +#endif + + maxv--; + hcount = 0; + mmax++; + } + + for (int i = n_vals-10; i < n_vals; ++i) { + printf("f: %f\n", val[i]); + } + + Assert(*len >= i); + *len = i; + + grib_context_free(a->context, scals); + + return ret; } #define MAXVAL(a, b) a > b ? a : b diff --git a/src/grib_ieeefloat.c b/src/grib_ieeefloat.c index 2c88d7630..b8bb60fd4 100644 --- a/src/grib_ieeefloat.c +++ b/src/grib_ieeefloat.c @@ -246,6 +246,7 @@ double grib_long_to_ieee(unsigned long x) return val; } + unsigned long grib_ieee_nearest_smaller_to_long(double x) { unsigned long l; @@ -386,6 +387,34 @@ int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, in return err; } +int grib_ieee_decode_array_float(grib_context* c, unsigned char* buf, size_t nvals, int bytes, float* val) +{ + int err = 0, i = 0, j = 0; + unsigned char s[4] = {0,}; + float fval; + + switch (bytes) { + case 4: + for (i = 0; i < nvals; i++) { +#if IEEE_LE + for (j = 3; j >= 0; j--) + s[j] = *(buf++); + memcpy(&val[i], s, 4); +#elif IEEE_BE + memcpy(&val[i], buf, 4); + buf += 4; +#endif + } + break; + default: + grib_context_log(c, GRIB_LOG_ERROR, + "grib_ieee_decode_array_float: %d bits not implemented", bytes * 8); + return GRIB_NOT_IMPLEMENTED; + } + + return err; +} + #else int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, double* val) @@ -399,6 +428,17 @@ int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, in return err; } +int grib_ieee_decode_array_float(grib_context* c, unsigned char* buf, size_t nvals, int bytes, float* val) +{ + int err = 0, i = 0; + long bitr = 0; + + for (i = 0; i < nvals; i++) + val[i] = (float) grib_long_to_ieee(grib_decode_unsigned_long(buf, &bitr, bytes * 8)); + + return err; +} + #endif #ifdef IEEE diff --git a/tests/grib_ecc-1467.c b/tests/grib_ecc-1467.c index 809c6c346..887964c8e 100644 --- a/tests/grib_ecc-1467.c +++ b/tests/grib_ecc-1467.c @@ -9,84 +9,85 @@ */ #include #include +#include #undef NDEBUG #include #include "eccodes.h" +#include "grib_api_internal.h" int main(int argc, char** argv) { - int err = 0; - float* fvalues = NULL; /* data values as floats */ - double* dvalues = NULL; /* data values as doubles */ - size_t values_len = 0; // number of data points - size_t cvalues_len = 0; // coded values excluding missing - size_t i = 0; - int mode = 2; // 1=single-precision, 2=double-precision + int err = 0; + float* fvalues = NULL; /* data values as floats */ + double* dvalues = NULL; /* data values as doubles */ + size_t values_len = 0; // number of data points + size_t cvalues_len = 0; // coded values excluding missing + size_t i = 0; + int mode = 2; // 1=single-precision, 2=double-precision - double daverage = 0; - float faverage = 0; + double daverage = 0; + float faverage = 0; + double abs_error = 0; + double max_abs_error = 1e-04; + double tolerance = 1e-04; + double dmin; + double dmax; + float fval; FILE* in = NULL; const char* filename = 0; codes_handle* h = NULL; - if (argc!=3) {fprintf(stderr,"usage: %s mode file\n",argv[0]); return 1;} - if (strcmp(argv[1], "double")==0) mode=2; - else if (strcmp(argv[1], "float")==0) mode=1; - else { fprintf(stderr,"Invalid mode: Use float or double\n");return 1; } - filename = argv[2]; - printf( "Opening %s, mode=%s\n",filename, (mode==1?"float":"double") ); + if (argc != 2) { + fprintf(stderr, "usage: %s file\n", argv[0]); + return 1; + } + filename = argv[1]; + + printf("Opening %s\n", filename); in = fopen(filename, "rb"); assert(in); - /* create new handle from the first message in the file*/ h = codes_handle_new_from_file(0, in, PRODUCT_GRIB, &err); assert(h); - fclose(in); - /* get the size of the values array*/ CODES_CHECK(codes_get_size(h, "values", &values_len), 0); - if (mode==1) - fvalues = (float*)malloc(values_len * sizeof(float)); - if (mode==2) - dvalues = (double*)malloc(values_len * sizeof(double)); + fvalues = (float*)malloc(values_len * sizeof(float)); + dvalues = (double*)malloc(values_len * sizeof(double)); + CODES_CHECK(codes_get_float_array(h, "values", fvalues, &values_len), 0); + CODES_CHECK(codes_get_double_array(h, "values", dvalues, &values_len), 0); - /* get data values*/ - if(mode==1) - CODES_CHECK(codes_get_float_array(h, "values", fvalues, &values_len), 0); - if(mode==2) - CODES_CHECK(codes_get_double_array(h, "values", dvalues, &values_len), 0); + for (i = 0; i < values_len; i++) { + abs_error = fabs(dvalues[i] - (double)fvalues[i]); + if (abs_error > max_abs_error) { + fprintf(stderr, "ERROR:\n\tfvalue %e\n\tdvalue %e\n\terror %e\n\tmax_abs_error %e\n", fvalues[i], dvalues[i], abs_error, max_abs_error); + Assert(!"Absolute error test failed\n"); + } - faverage = 0; - daverage = 0; - if(mode==1) { - for (i = 0; i < values_len; i++) { - if (fvalues[i] != 9999) { - //if(i<10)printf("%10.15f\n",fvalues[i]); - faverage += fvalues[i]; - cvalues_len++; - } + dmin = dvalues[i] >= 0 ? dvalues[i] / (1 + tolerance) : dvalues[i] * (1 + tolerance); + dmax = dvalues[i] >= 0 ? dvalues[i] * (1 + tolerance) : dvalues[i] / (1 + tolerance); + fval = fvalues[i]; + + if (!((dmin <= fval) && (fval <= dmax))) { + fprintf(stderr, "Error:\n"); + fprintf(stderr, "dvalue: %f, fvalue: %f\n", dvalues[i], fvalues[i]); + fprintf(stderr, "\tmin < fvalue < max = %.20e < %.20e < %.20e FAILED\n", + dmin, fvalues[i], dmax); + fprintf(stderr, "\tfvalue - min = %.20e (%s)\n", + fvalues[i] - dmin, fvalues[i] - dmin >= 0 ? "OK" : "FAILED (should be positive)"); + fprintf(stderr, "\tmax - fvalue = %.20e (%s)\n", + dmax - fvalues[i], dmax - fvalues[i] >= 0 ? "OK" : "FAILED (should be positive)"); + + Assert(!"Relative tolerance test failed\n"); } - faverage /= (float)cvalues_len; - free(fvalues); - printf("\tThere are %zu total values, %zu coded, float average = %10.15f\n", values_len, cvalues_len, faverage); - } - if (mode==2){ - for (i = 0; i < values_len; i++) { - if (dvalues[i] != 9999) { - //if(i<10)printf("%10.15f\n",dvalues[i]); - daverage += dvalues[i]; - cvalues_len++; - } - } - daverage /= (double)cvalues_len; - free(dvalues); - printf("\tThere are %zu total values, %zu coded, double average = %10.15g\n", values_len, cvalues_len, daverage); } + free(fvalues); + free(dvalues); + codes_handle_delete(h); - + fclose(in); return 0; } diff --git a/tests/grib_ecc-1467.sh b/tests/grib_ecc-1467.sh index a38df98ca..15f6a3157 100755 --- a/tests/grib_ecc-1467.sh +++ b/tests/grib_ecc-1467.sh @@ -10,23 +10,40 @@ . ./include.ctest.sh -gfiles=" - sst_globus0083.grib - constant_field.grib1 - constant_field.grib2 - reduced_latlon_surface.grib1 - reduced_latlon_surface.grib2 - regular_latlon_surface.grib1 - regular_latlon_surface.grib2 -" +# Constant Fields +gfiles+="constant_field.grib1 " # grid_simple, edition=1 +gfiles+="constant_field.grib2 " # grid_simple, edition=2 + +# Grid Simple +gfiles+="sst_globus0083.grib " # grid_simple, edition=1 +gfiles+="reduced_latlon_surface.grib1 " # grid_simple, edition=1 +gfiles+="reduced_latlon_surface.grib2 " # grid_simple, edition=2 +gfiles+="regular_latlon_surface.grib1 " # grid_simple, edition=1 +gfiles+="regular_latlon_surface.grib2 " # grid_simple, edition=2 + +# Spherical Complex +gfiles+="spherical_pressure_level.grib1 " # spectral_complex, edition=1 +gfiles+="spherical_pressure_level.grib2 " # spectral_complex, edition=2 + +# Complex +#gfiles+=" "# grid_complex, edition=1 +#gfiles+="gfs.complex.mvmu.grib2 " # grid_complex, edition=2, g22order_packing + +# Second Order +#gfiles+="lfpw.grib1 " # grid_second_order, edition=1 +#gfiles+=" "# grid_second_order, edition=2 + +# CCSDS if [ $HAVE_AEC -eq 1 ]; then echo "Adding extra files (HAVE_AEC=1)" - gfiles="ccsds.grib2 "$gfiles + gfiles+="ccsds.grib2 " fi for f in $gfiles; do input=${data_dir}/$f - ${tools_dir}/grib_ls -p numberOfDataPoints,numberOfCodedValues,numberOfMissing,avg $input - $EXEC ${test_dir}/grib_ecc-1467 double $input - $EXEC ${test_dir}/grib_ecc-1467 float $input + ${tools_dir}/grib_ls -w count=1 -p numberOfDataPoints,numberOfCodedValues,numberOfMissing,avg $input + $EXEC ${test_dir}/grib_ecc-1467 $input done + + +