ECC-1467: Complex packing single-precision

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 */

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<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, 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

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

tests/grib_ecc-1467.c

@@ -9,84 +9,85 @@
  */
 #include <stdio.h>
 #include <stdlib.h>
+#include <math.h>
 #undef NDEBUG
 #include <assert.h>
 #include "eccodes.h"
+#include "grib_api_internal.h"
 
 int main(int argc, char** argv)
 {
-    int err           = 0;
+    int err            = 0;
-    float* fvalues    = NULL; /* data values as floats */
+    float* fvalues     = NULL; /* data values as floats */
-    double* dvalues   = NULL; /* data values as doubles */
+    double* dvalues    = NULL; /* data values as doubles */
-    size_t values_len = 0;  // number of data points
+    size_t values_len  = 0;    // number of data points
-    size_t cvalues_len = 0; // coded values excluding missing
+    size_t cvalues_len = 0;    // coded values excluding missing
-    size_t i = 0;
+    size_t i           = 0;
-    int mode = 2; // 1=single-precision, 2=double-precision
+    int mode           = 2;  // 1=single-precision, 2=double-precision
 
-    double daverage    = 0;
+    double daverage      = 0;
-    float  faverage    = 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));
-        fvalues = (float*)malloc(values_len * sizeof(float));
+    dvalues = (double*)malloc(values_len * sizeof(double));
-    if (mode==2)
+    CODES_CHECK(codes_get_float_array(h, "values", fvalues, &values_len), 0);
-        dvalues = (double*)malloc(values_len * sizeof(double));
+    CODES_CHECK(codes_get_double_array(h, "values", dvalues, &values_len), 0);
 
-    /* get data values*/
+    for (i = 0; i < values_len; i++) {
-    if(mode==1)
+        abs_error = fabs(dvalues[i] - (double)fvalues[i]);
-        CODES_CHECK(codes_get_float_array(h, "values", fvalues, &values_len), 0);
+        if (abs_error > max_abs_error) {
-    if(mode==2)
+            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);
-        CODES_CHECK(codes_get_double_array(h, "values", dvalues, &values_len), 0);
+            Assert(!"Absolute error test failed\n");
+        }
 
-    faverage = 0;
+        dmin = dvalues[i] >= 0 ? dvalues[i] / (1 + tolerance) : dvalues[i] * (1 + tolerance);
-    daverage = 0;
+        dmax = dvalues[i] >= 0 ? dvalues[i] * (1 + tolerance) : dvalues[i] / (1 + tolerance);
-    if(mode==1) {
+        fval = fvalues[i];
-        for (i = 0; i < values_len; i++) {
+
-            if (fvalues[i] != 9999) {
+        if (!((dmin <= fval) && (fval <= dmax))) {
-                //if(i<10)printf("%10.15f\n",fvalues[i]);
+            fprintf(stderr, "Error:\n");
-                faverage += fvalues[i];
+            fprintf(stderr, "dvalue: %f, fvalue: %f\n", dvalues[i], fvalues[i]);
-                cvalues_len++;
+            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;
 }

tests/grib_ecc-1467.sh

@@ -10,23 +10,40 @@
 
 . ./include.ctest.sh
 
-gfiles="
+# Constant Fields
-    sst_globus0083.grib
+gfiles+="constant_field.grib1 " # grid_simple, edition=1
-    constant_field.grib1
+gfiles+="constant_field.grib2 " # grid_simple, edition=2
-    constant_field.grib2
+
-    reduced_latlon_surface.grib1
+# Grid Simple
-    reduced_latlon_surface.grib2
+gfiles+="sst_globus0083.grib " # grid_simple, edition=1
-    regular_latlon_surface.grib1
+gfiles+="reduced_latlon_surface.grib1 " # grid_simple, edition=1
-    regular_latlon_surface.grib2
+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
+    ${tools_dir}/grib_ls -w count=1 -p numberOfDataPoints,numberOfCodedValues,numberOfMissing,avg $input
-    $EXEC ${test_dir}/grib_ecc-1467 double $input
+    $EXEC ${test_dir}/grib_ecc-1467 $input
-    $EXEC ${test_dir}/grib_ecc-1467 float  $input
 done
+
+
+