core
/
eccodes
mirror of https://github.com/ecmwf/eccodes.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ECC-1467: Single-precision implementation with the use of templates

This commit is contained in:
Eugen Betke 2023-02-20 08:41:32 +00:00
parent 0a9d31197d
commit 687e3f4cf6
22 changed files with 940 additions and 1250 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
src/eccodes_prototypes.h
View File

@ -1,3 +1,5 @@
#pragma once
#ifdef ECCODES_ON_WINDOWS
#include <stdint.h>
#endif
@ -836,8 +838,9 @@ 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);
// ECC-1467
//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.cc*/
@ -1548,8 +1551,8 @@ int grib_encode_size_tb(unsigned char* p, size_t val, long* bitp, long nb);
/* grib_bits_any_endian_simple.cc*/
int grib_decode_long_array(const unsigned char* p, long* bitp, long bitsPerValue, size_t n_vals, long* val);
int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerValue, double reference_value, double s, double d, size_t n_vals, double* val);
int grib_decode_float_array(const unsigned char* p, long* bitp, long bitsPerValue, double reference_value, double s, double d, size_t n_vals, float* val);
//int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerValue, double reference_value, double s, double d, size_t n_vals, double* val);
//int grib_decode_float_array(const unsigned char* p, long* bitp, long bitsPerValue, double reference_value, double s, double d, size_t n_vals, float* val);
int grib_decode_double_array_complex(const unsigned char* p, long* bitp, long nbits, double reference_value, double s, double* d, size_t size, double* val);
int grib_encode_long_array(size_t n_vals, const long* val, long bits_per_value, unsigned char* p, long* off);
int grib_encode_double_array(size_t n_vals, const double* val, long bits_per_value, double reference_value, double d, double divisor, unsigned char* p, long* off);

61
src/grib_accessor_class_bitmap.cc
View File

@ -8,7 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_accessor_class_bitmap.h"
/*
This is used by make_class.pl
@ -40,8 +40,8 @@ or edit "accessor.class" and rerun ./make_class.pl
*/
static int unpack_double(grib_accessor*, double* val, size_t* len);
static int unpack_float(grib_accessor*, float* val, size_t* len);
//static int unpack_double(grib_accessor*, double* val, size_t* len);
//static int unpack_float(grib_accessor*, float* val, size_t* len);
static int unpack_long(grib_accessor*, long* val, size_t* len);
static int unpack_string(grib_accessor*, char*, size_t* len);
static long next_offset(grib_accessor*);
@ -89,8 +89,8 @@ static grib_accessor_class _grib_accessor_class_bitmap = {
&unpack_long, /* grib_unpack procedures long */
0, /* grib_pack procedures double */
0, /* grib_pack procedures float */
&unpack_double, /* grib_unpack procedures double */
&unpack_float, /* grib_unpack procedures float */
&GribAccessorClassBitmap<double>::unpack, /* grib_unpack procedures double */
&GribAccessorClassBitmap<float>::unpack, /* grib_unpack procedures float */
0, /* grib_pack procedures string */
&unpack_string, /* grib_unpack procedures string */
0, /* grib_pack array procedures string */
@ -247,57 +247,6 @@ static int unpack_long(grib_accessor* a, long* val, size_t* len)
return GRIB_SUCCESS;
}
static int unpack_double(grib_accessor* a, double* val, size_t* len)
{
long pos = a->offset * 8;
long tlen;
long i;
int err = 0;
grib_handle* hand = grib_handle_of_accessor(a);
err = grib_value_count(a, &tlen);
if (err)
return err;
if (*len < tlen) {
grib_context_log(a->context, GRIB_LOG_ERROR, "Wrong size for %s it contains %ld values", a->name, tlen);
*len = 0;
return GRIB_ARRAY_TOO_SMALL;
}
for (i = 0; i < tlen; i++) {
val[i] = (double)grib_decode_unsigned_long(hand->buffer->data, &pos, 1);
}
*len = tlen;
return GRIB_SUCCESS;
}
// TODO(maee): ECC-1467: Copied the 'double' version and reused by copy/paste!
static int unpack_float(grib_accessor* a, float* val, size_t* len)
{
long pos = a->offset * 8;
long tlen;
long i;
int err = 0;
grib_handle* hand = grib_handle_of_accessor(a);
err = grib_value_count(a, &tlen);
if (err)
return err;
if (*len < tlen) {
grib_context_log(a->context, GRIB_LOG_ERROR, "Wrong size for %s it contains %ld values", a->name, tlen);
*len = 0;
return GRIB_ARRAY_TOO_SMALL;
}
for (i = 0; i < tlen; i++) {
val[i] = (float)grib_decode_unsigned_long(hand->buffer->data, &pos, 1);
}
*len = tlen;
return GRIB_SUCCESS;
}
//
static int unpack_double_element(grib_accessor* a, size_t idx, double* val)
{
long pos = a->offset * 8;

39
src/grib_accessor_class_bitmap.h Normal file
View File

@ -0,0 +1,39 @@
// ECC-1467
#pragma once
#include "grib_api_internal_cpp.h"
#include <type_traits>
template <typename T>
class GribAccessorClassBitmap {
public:
static int unpack(grib_accessor* a, T* val, size_t* len);
};
template <typename T>
int GribAccessorClassBitmap<T>::unpack(grib_accessor* a, T* val, size_t* len)
{
static_assert(std::is_floating_point<T>::value, "Requires floating points numbers");
long pos = a->offset * 8;
long tlen;
long i;
int err = 0;
grib_handle* hand = grib_handle_of_accessor(a);
err = grib_value_count(a, &tlen);
if (err)
return err;
if (*len < tlen) {
grib_context_log(a->context, GRIB_LOG_ERROR, "Wrong size for %s it contains %ld values", a->name, tlen);
*len = 0;
return GRIB_ARRAY_TOO_SMALL;
}
for (i = 0; i < tlen; i++) {
val[i] = (T)grib_decode_unsigned_long(hand->buffer->data, &pos, 1);
}
*len = tlen;
return GRIB_SUCCESS;
}

181
src/grib_accessor_class_data_apply_bitmap.cc
View File

@ -9,6 +9,7 @@
*/
#include "grib_api_internal.h"
#include "grib_accessor_class_data_apply_bitmap.h"
/*
This is used by make_class.pl
@ -43,8 +44,6 @@ or edit "accessor.class" and rerun ./make_class.pl
static int get_native_type(grib_accessor*);
static int pack_double(grib_accessor*, const double* val, size_t* len);
static int unpack_double(grib_accessor*, double* val, size_t* len);
static int unpack_float(grib_accessor*, float* val, size_t* len);
static int value_count(grib_accessor*, long*);
static void dump(grib_accessor*, grib_dumper*);
static void init(grib_accessor*, const long, grib_arguments*);
@ -52,18 +51,6 @@ static void init_class(grib_accessor_class*);
static int unpack_double_element(grib_accessor*, size_t i, double* val);
static int unpack_double_element_set(grib_accessor*, const size_t* index_array, size_t len, double* val_array);
typedef struct grib_accessor_data_apply_bitmap
{
grib_accessor att;
/* Members defined in gen */
/* Members defined in data_apply_bitmap */
const char* coded_values;
const char* bitmap;
const char* missing_value;
const char* number_of_data_points;
const char* number_of_values;
const char* binary_scale_factor;
} grib_accessor_data_apply_bitmap;
extern grib_accessor_class* grib_accessor_class_gen;
@ -90,8 +77,8 @@ static grib_accessor_class _grib_accessor_class_data_apply_bitmap = {
0, /* grib_unpack procedures long */
&pack_double, /* grib_pack procedures double */
0, /* grib_pack procedures float */
&unpack_double, /* grib_unpack procedures double */
&unpack_float, /* grib_unpack procedures float */
&GribAccessorDataApplyBitmap<double>::unpack, /* grib_unpack procedures double */
&GribAccessorDataApplyBitmap<float>::unpack, /* grib_unpack procedures float */
0, /* grib_pack procedures string */
0, /* grib_unpack procedures string */
0, /* grib_pack array procedures string */
@ -189,168 +176,6 @@ static int value_count(grib_accessor* a, long* count)
return ret;
}
static int unpack_double(grib_accessor* a, double* val, size_t* len)
{
grib_accessor_data_apply_bitmap* self = (grib_accessor_data_apply_bitmap*)a;
size_t i = 0;
size_t j = 0;
size_t n_vals = 0;
long nn = 0;
int err = 0;
size_t coded_n_vals = 0;
double* coded_vals = NULL;
double missing_value = 0;
err = grib_value_count(a, &nn);
n_vals = nn;
if (err)
return err;
if (!grib_find_accessor(grib_handle_of_accessor(a), self->bitmap))
return grib_get_double_array(grib_handle_of_accessor(a), self->coded_values, val, len);
if ((err = grib_get_size(grib_handle_of_accessor(a), self->coded_values, &coded_n_vals)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(grib_handle_of_accessor(a), self->missing_value, &missing_value)) != GRIB_SUCCESS)
return err;
if (*len < n_vals) {
*len = n_vals;
return GRIB_ARRAY_TOO_SMALL;
}
if (coded_n_vals == 0) {
for (i = 0; i < n_vals; i++)
val[i] = missing_value;
*len = n_vals;
return GRIB_SUCCESS;
}
if ((err = grib_get_double_array_internal(grib_handle_of_accessor(a), self->bitmap, val, &n_vals)) != GRIB_SUCCESS)
return err;
coded_vals = (double*)grib_context_malloc(a->context, coded_n_vals * sizeof(double));
if (coded_vals == NULL)
return GRIB_OUT_OF_MEMORY;
if ((err = grib_get_double_array(grib_handle_of_accessor(a), self->coded_values, coded_vals, &coded_n_vals)) != GRIB_SUCCESS) {
grib_context_free(a->context, coded_vals);
return err;
}
grib_context_log(a->context, GRIB_LOG_DEBUG,
"grib_accessor_class_data_apply_bitmap: unpack_double : creating %s, %d values",
a->name, n_vals);
for (i = 0; i < n_vals; i++) {
if (val[i] == 0) {
val[i] = missing_value;
}
else {
val[i] = coded_vals[j++];
if (j > coded_n_vals) {
grib_context_free(a->context, coded_vals);
grib_context_log(a->context, GRIB_LOG_ERROR,
"grib_accessor_class_data_apply_bitmap [%s]:"
" unpack_double : number of coded values does not match bitmap %ld %ld",
a->name, coded_n_vals, n_vals);
return GRIB_ARRAY_TOO_SMALL;
}
}
}
*len = n_vals;
grib_context_free(a->context, coded_vals);
return err;
}
// TODO(maee): ECC-1467: Copied the 'double' version and reused by copy/paste!
// Should be merged with unpack_double and refactored! Most probably using C++ templates
static int unpack_float(grib_accessor* a, float* val, size_t* len)
{
grib_accessor_data_apply_bitmap* self = (grib_accessor_data_apply_bitmap*)a;
size_t i = 0;
size_t j = 0;
size_t n_vals = 0;
long nn = 0;
int err = 0;
size_t coded_n_vals = 0;
float* coded_vals = NULL;
double missing_value = 0;
err = grib_value_count(a, &nn);
n_vals = nn;
if (err)
return err;
if (!grib_find_accessor(grib_handle_of_accessor(a), self->bitmap))
return grib_get_float_array(grib_handle_of_accessor(a), self->coded_values, val, len);
if ((err = grib_get_size(grib_handle_of_accessor(a), self->coded_values, &coded_n_vals)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(grib_handle_of_accessor(a), self->missing_value, &missing_value)) != GRIB_SUCCESS)
return err;
if (*len < n_vals) {
*len = n_vals;
return GRIB_ARRAY_TOO_SMALL;
}
if (coded_n_vals == 0) {
for (i = 0; i < n_vals; i++)
val[i] = missing_value;
*len = n_vals;
return GRIB_SUCCESS;
}
if ((err = grib_get_float_array_internal(grib_handle_of_accessor(a), self->bitmap, val, &n_vals)) != GRIB_SUCCESS)
return err;
coded_vals = (float*)grib_context_malloc(a->context, coded_n_vals * sizeof(float));
if (coded_vals == NULL)
return GRIB_OUT_OF_MEMORY;
if ((err = grib_get_float_array(grib_handle_of_accessor(a), self->coded_values, coded_vals, &coded_n_vals)) != GRIB_SUCCESS) {
grib_context_free(a->context, coded_vals);
return err;
}
grib_context_log(a->context, GRIB_LOG_DEBUG,
"grib_accessor_class_data_apply_bitmap: unpack_float : creating %s, %d values",
a->name, n_vals);
for (i = 0; i < n_vals; i++) {
if (val[i] == 0) {
val[i] = missing_value;
}
else {
val[i] = coded_vals[j++];
if (j > coded_n_vals) {
grib_context_free(a->context, coded_vals);
grib_context_log(a->context, GRIB_LOG_ERROR,
"grib_accessor_class_data_apply_bitmap [%s]:"
" unpack_float : number of coded values does not match bitmap %ld %ld",
a->name, coded_n_vals, n_vals);
return GRIB_ARRAY_TOO_SMALL;
}
}
}
*len = n_vals;
grib_context_free(a->context, coded_vals);
return err;
}
static int unpack_double_element(grib_accessor* a, size_t idx, double* val)
{
grib_accessor_data_apply_bitmap* self = (grib_accessor_data_apply_bitmap*)a;

109
src/grib_accessor_class_data_apply_bitmap.h Normal file
View File

@ -0,0 +1,109 @@
// ECC-1467
#pragma once
#include "grib_api_internal_cpp.h"
#include <typeinfo>
#include <type_traits>
typedef struct grib_accessor_data_apply_bitmap
{
grib_accessor att;
/* Members defined in gen */
/* Members defined in data_apply_bitmap */
const char* coded_values;
const char* bitmap;
const char* missing_value;
const char* number_of_data_points;
const char* number_of_values;
const char* binary_scale_factor;
} grib_accessor_data_apply_bitmap;
template <typename T>
class GribAccessorDataApplyBitmap {
public:
static int unpack(grib_accessor* a, T* val, size_t* len);
};
template <typename T>
int GribAccessorDataApplyBitmap<T>::unpack(grib_accessor* a, T* val, size_t* len)
{
static_assert(std::is_floating_point<T>::value, "Requires floating point numbers");
grib_accessor_data_apply_bitmap* self = (grib_accessor_data_apply_bitmap*)a;
size_t i = 0;
size_t j = 0;
size_t n_vals = 0;
long nn = 0;
int err = 0;
size_t coded_n_vals = 0;
T* coded_vals = NULL;
double missing_value = 0;
err = grib_value_count(a, &nn);
n_vals = nn;
if (err)
return err;
if (!grib_find_accessor(grib_handle_of_accessor(a), self->bitmap))
return grib_get_array<T>(grib_handle_of_accessor(a), self->coded_values, val, len);
if ((err = grib_get_size(grib_handle_of_accessor(a), self->coded_values, &coded_n_vals)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(grib_handle_of_accessor(a), self->missing_value, &missing_value)) != GRIB_SUCCESS)
return err;
if (*len < n_vals) {
*len = n_vals;
return GRIB_ARRAY_TOO_SMALL;
}
if (coded_n_vals == 0) {
for (i = 0; i < n_vals; i++)
val[i] = missing_value;
*len = n_vals;
return GRIB_SUCCESS;
}
if ((err = grib_get_array_internal<T>(grib_handle_of_accessor(a), self->bitmap, val, &n_vals)) != GRIB_SUCCESS)
return err;
coded_vals = (T*)grib_context_malloc(a->context, coded_n_vals * sizeof(T));
if (coded_vals == NULL)
return GRIB_OUT_OF_MEMORY;
if ((err = grib_get_array<T>(grib_handle_of_accessor(a), self->coded_values, coded_vals, &coded_n_vals)) != GRIB_SUCCESS) {
grib_context_free(a->context, coded_vals);
return err;
}
grib_context_log(a->context, GRIB_LOG_DEBUG,
"grib_accessor_class_data_apply_bitmap: %s : creating %s, %d values",
__PRETTY_FUNCTION__,
a->name, n_vals);
for (i = 0; i < n_vals; i++) {
if (val[i] == 0) {
val[i] = missing_value;
}
else {
val[i] = coded_vals[j++];
if (j > coded_n_vals) {
grib_context_free(a->context, coded_vals);
grib_context_log(a->context, GRIB_LOG_ERROR,
"grib_accessor_class_data_apply_bitmap [%s]:"
" %s : number of coded values does not match bitmap %ld %ld",
a->name, __PRETTY_FUNCTION__, coded_n_vals, n_vals);
return GRIB_ARRAY_TOO_SMALL;
}
}
}
*len = n_vals;
grib_context_free(a->context, coded_vals);
return err;
}

271
src/grib_accessor_class_data_ccsds_packing.cc
View File

@ -8,7 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_accessor_class_data_ccsds_packing.h"
/*
This is used by make_class.pl
@ -49,36 +49,12 @@ or edit "accessor.class" and rerun ./make_class.pl
*/
static int pack_double(grib_accessor*, const double* val, size_t* len);
static int unpack_double(grib_accessor*, double* val, size_t* len);
static int unpack_float(grib_accessor*, float* val, size_t* len);
static int value_count(grib_accessor*, long*);
static void init(grib_accessor*, const long, grib_arguments*);
static void init_class(grib_accessor_class*);
static int unpack_double_element(grib_accessor*, size_t i, double* val);
static int unpack_double_element_set(grib_accessor*, const size_t* index_array, size_t len, double* val_array);
typedef struct grib_accessor_data_ccsds_packing
{
grib_accessor att;
/* Members defined in gen */
/* Members defined in values */
int carg;
const char* seclen;
const char* offsetdata;
const char* offsetsection;
int dirty;
/* Members defined in data_ccsds_packing */
const char* number_of_values;
const char* reference_value;
const char* binary_scale_factor;
const char* decimal_scale_factor;
const char* bits_per_value;
const char* number_of_data_points;
const char* ccsds_flags;
const char* ccsds_block_size;
const char* ccsds_rsi;
} grib_accessor_data_ccsds_packing;
extern grib_accessor_class* grib_accessor_class_values;
static grib_accessor_class _grib_accessor_class_data_ccsds_packing = {
@ -104,8 +80,8 @@ static grib_accessor_class _grib_accessor_class_data_ccsds_packing = {
0, /* grib_unpack procedures long */
&pack_double, /* grib_pack procedures double */
0, /* grib_pack procedures float */
&unpack_double, /* grib_unpack procedures double */
&unpack_float, /* grib_unpack procedures float */
&GribAccessorDataCcsdsPacking<double>::unpack, /* grib_unpack procedures double */
&GribAccessorDataCcsdsPacking<float>::unpack, /* grib_unpack procedures float */
0, /* grib_pack procedures string */
0, /* grib_unpack procedures string */
0, /* grib_pack array procedures string */
@ -200,7 +176,7 @@ static int value_count(grib_accessor* a, long* count)
#include <libaec.h>
static const char* aec_get_error_message(int code)
const char* aec_get_error_message(int code)
{
if (code == AEC_MEM_ERROR) return "AEC_MEM_ERROR";
if (code == AEC_DATA_ERROR) return "AEC_DATA_ERROR";
@ -209,7 +185,8 @@ static const char* aec_get_error_message(int code)
if (code == AEC_OK) return "AEC_OK";
return "Unknown error code";
}
static void print_aec_stream_info(struct aec_stream* strm, const char* func)
void print_aec_stream_info(struct aec_stream* strm, const char* func)
{
fprintf(stderr, "ECCODES DEBUG CCSDS %s aec_stream.flags=%u\n", func, strm->flags);
fprintf(stderr, "ECCODES DEBUG CCSDS %s aec_stream.bits_per_sample=%u\n", func, strm->bits_per_sample);
@ -219,232 +196,6 @@ static void print_aec_stream_info(struct aec_stream* strm, const char* func)
fprintf(stderr, "ECCODES DEBUG CCSDS %s aec_stream.avail_in=%lu\n", func, strm->avail_in);
}
static int unpack_double(grib_accessor* a, double* val, size_t* len)
{
grib_accessor_data_ccsds_packing* self = (grib_accessor_data_ccsds_packing*)a;
grib_handle* hand = grib_handle_of_accessor(a);
int err = GRIB_SUCCESS, i = 0;
size_t buflen = 0;
struct aec_stream strm;
double bscale = 0;
double dscale = 0;
unsigned char* buf = NULL;
size_t n_vals = 0;
size_t size = 0;
unsigned char* decoded = NULL;
/*unsigned char* p = NULL;*/
long pos = 0;
long nn = 0;
long binary_scale_factor = 0;
long decimal_scale_factor = 0;
double reference_value = 0;
long bits_per_value = 0;
long bits8;
long ccsds_flags;
long ccsds_block_size;
long ccsds_rsi;
self->dirty = 0;
if ((err = grib_value_count(a, &nn)) != GRIB_SUCCESS)
return err;
n_vals = nn;
if ((err = grib_get_long_internal(hand, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(hand, self->reference_value, &reference_value)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS)
return err;
/* ECC-477: Don't call grib_get_long_internal to suppress error message being output */
if ((err = grib_get_long(hand, self->ccsds_flags, &ccsds_flags)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->ccsds_block_size, &ccsds_block_size)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->ccsds_rsi, &ccsds_rsi)) != GRIB_SUCCESS)
return err;
/* TODO: This should be called upstream */
if (*len < n_vals)
return GRIB_ARRAY_TOO_SMALL;
/* Special case */
if (bits_per_value == 0) {
for (i = 0; i < n_vals; i++)
val[i] = reference_value;
*len = n_vals;
return GRIB_SUCCESS;
}
bscale = grib_power(binary_scale_factor, 2);
dscale = grib_power(-decimal_scale_factor, 10);
buflen = grib_byte_count(a);
buf = (unsigned char*)hand->buffer->data;
buf += grib_byte_offset(a);
strm.flags = ccsds_flags;
strm.bits_per_sample = bits_per_value;
strm.block_size = ccsds_block_size;
strm.rsi = ccsds_rsi;
strm.next_in = buf;
strm.avail_in = buflen;
bits8 = ((bits_per_value + 7) / 8) * 8;
size = n_vals * ((bits_per_value + 7) / 8);
decoded = (unsigned char*)grib_context_buffer_malloc_clear(a->context, size);
if (!decoded) {
err = GRIB_OUT_OF_MEMORY;
goto cleanup;
}
strm.next_out = decoded;
strm.avail_out = size;
if (hand->context->debug) print_aec_stream_info(&strm, "unpack_double");
if ((err = aec_buffer_decode(&strm)) != AEC_OK) {
grib_context_log(a->context, GRIB_LOG_ERROR, "CCSDS unpack_double: aec_buffer_decode error %d (%s)\n",
err, aec_get_error_message(err));
err = GRIB_ENCODING_ERROR;
goto cleanup;
}
/* printf("bscale=%g dscale=%g reference_value=%g\n",bscale,dscale,reference_value); */
pos = 0;
#if 0
p = decoded;
for (i = 0; i < n_vals; i++) {
val[i] = (double)(((grib_decode_unsigned_long(p, &pos, bits8) * bscale) + reference_value) * dscale);
}
#endif
/* ECC-1427: Performance improvement */
grib_decode_double_array(decoded, &pos, bits8 , reference_value, bscale, dscale, n_vals, val);
*len = n_vals;
cleanup:
grib_context_buffer_free(a->context, decoded);
return err;
}
// TODO(maee): ECC-1467: Copied the 'double' version and reused by copy/paste!
// Should be merged with unpack_double and refactored! Most probably using C++ templates
static int unpack_float(grib_accessor* a, float* val, size_t* len)
{
grib_accessor_data_ccsds_packing* self = (grib_accessor_data_ccsds_packing*)a;
grib_handle* hand = grib_handle_of_accessor(a);
int err = GRIB_SUCCESS, i = 0;
size_t buflen = 0;
struct aec_stream strm;
double bscale = 0;
double dscale = 0;
unsigned char* buf = NULL;
size_t n_vals = 0;
size_t size = 0;
unsigned char* decoded = NULL;
/*unsigned char* p = NULL;*/
long pos = 0;
long nn = 0;
long binary_scale_factor = 0;
long decimal_scale_factor = 0;
double reference_value = 0;
long bits_per_value = 0;
long bits8;
long ccsds_flags;
long ccsds_block_size;
long ccsds_rsi;
self->dirty = 0;
if ((err = grib_value_count(a, &nn)) != GRIB_SUCCESS)
return err;
n_vals = nn;
if ((err = grib_get_long_internal(hand, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(hand, self->reference_value, &reference_value)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS)
return err;
/* ECC-477: Don't call grib_get_long_internal to suppress error message being output */
if ((err = grib_get_long(hand, self->ccsds_flags, &ccsds_flags)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->ccsds_block_size, &ccsds_block_size)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->ccsds_rsi, &ccsds_rsi)) != GRIB_SUCCESS)
return err;
/* TODO: This should be called upstream */
if (*len < n_vals)
return GRIB_ARRAY_TOO_SMALL;
/* Special case */
if (bits_per_value == 0) {
for (i = 0; i < n_vals; i++)
val[i] = reference_value;
*len = n_vals;
return GRIB_SUCCESS;
}
bscale = grib_power(binary_scale_factor, 2);
dscale = grib_power(-decimal_scale_factor, 10);
buflen = grib_byte_count(a);
buf = (unsigned char*)hand->buffer->data;
buf += grib_byte_offset(a);
strm.flags = ccsds_flags;
strm.bits_per_sample = bits_per_value;
strm.block_size = ccsds_block_size;
strm.rsi = ccsds_rsi;
strm.next_in = buf;
strm.avail_in = buflen;
bits8 = ((bits_per_value + 7) / 8) * 8;
size = n_vals * ((bits_per_value + 7) / 8);
decoded = (unsigned char*)grib_context_buffer_malloc_clear(a->context, size);
if (!decoded) {
err = GRIB_OUT_OF_MEMORY;
goto cleanup;
}
strm.next_out = decoded;
strm.avail_out = size;
if (hand->context->debug) print_aec_stream_info(&strm, "unpack_float");
if ((err = aec_buffer_decode(&strm)) != AEC_OK) {
grib_context_log(a->context, GRIB_LOG_ERROR, "CCSDS unpack_float: aec_buffer_decode error %d (%s)\n",
err, aec_get_error_message(err));
err = GRIB_ENCODING_ERROR;
goto cleanup;
}
pos = 0;
/* ECC-1427: Performance improvement */
grib_decode_float_array(decoded, &pos, bits8 , reference_value, bscale, dscale, n_vals, val);
*len = n_vals;
cleanup:
grib_context_buffer_free(a->context, decoded);
return err;
}
static int pack_double(grib_accessor* a, const double* val, size_t* len)
{
grib_accessor_data_ccsds_packing* self = (grib_accessor_data_ccsds_packing*)a;
@ -787,16 +538,6 @@ static void print_error_feature_not_enabled(grib_context* c)
"grib_accessor_data_ccsds_packing: CCSDS support not enabled. "
"Please rebuild with -DENABLE_AEC=ON (Adaptive Entropy Coding library)");
}
static int unpack_double(grib_accessor* a, double* val, size_t* len)
{
print_error_feature_not_enabled(a->context);
return GRIB_FUNCTIONALITY_NOT_ENABLED;
}
static int unpack_float(grib_accessor* a, float* val, size_t* len)
{
print_error_feature_not_enabled(a->context);
return GRIB_FUNCTIONALITY_NOT_ENABLED;
}
static int pack_double(grib_accessor* a, const double* val, size_t* len)
{
print_error_feature_not_enabled(a->context);

160
src/grib_accessor_class_data_ccsds_packing.h Normal file
View File

@ -0,0 +1,160 @@
#pragma once
#include "grib_api_internal_cpp.h"
#if defined(HAVE_LIBAEC) || defined(HAVE_AEC)
#include <libaec.h>
const char* aec_get_error_message(int code);
void print_aec_stream_info(struct aec_stream* strm, const char* func);
#endif
//static int unpack_double(grib_accessor*, double* val, size_t* len);
//static int unpack_float(grib_accessor*, float* val, size_t* len);
typedef struct grib_accessor_data_ccsds_packing
{
grib_accessor att;
/* Members defined in gen */
/* Members defined in values */
int carg;
const char* seclen;
const char* offsetdata;
const char* offsetsection;
int dirty;
/* Members defined in data_ccsds_packing */
const char* number_of_values;
const char* reference_value;
const char* binary_scale_factor;
const char* decimal_scale_factor;
const char* bits_per_value;
const char* number_of_data_points;
const char* ccsds_flags;
const char* ccsds_block_size;
const char* ccsds_rsi;
} grib_accessor_data_ccsds_packing;
template <typename T>
class GribAccessorDataCcsdsPacking {
public:
#if defined(HAVE_LIBAEC) || defined(HAVE_AEC)
static int unpack(grib_accessor* a, T* val, size_t* len)
{
grib_accessor_data_ccsds_packing* self = (grib_accessor_data_ccsds_packing*)a;
grib_handle* hand = grib_handle_of_accessor(a);
int err = GRIB_SUCCESS, i = 0;
size_t buflen = 0;
struct aec_stream strm;
double bscale = 0;
double dscale = 0;
unsigned char* buf = NULL;
size_t n_vals = 0;
size_t size = 0;
unsigned char* decoded = NULL;
/*unsigned char* p = NULL;*/
long pos = 0;
long nn = 0;
long binary_scale_factor = 0;
long decimal_scale_factor = 0;
double reference_value = 0;
long bits_per_value = 0;
long bits8;
long ccsds_flags;
long ccsds_block_size;
long ccsds_rsi;
self->dirty = 0;
if ((err = grib_value_count(a, &nn)) != GRIB_SUCCESS)
return err;
n_vals = nn;
if ((err = grib_get_long_internal(hand, self->bits_per_value, &bits_per_value)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_double_internal(hand, self->reference_value, &reference_value)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->binary_scale_factor, &binary_scale_factor)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->decimal_scale_factor, &decimal_scale_factor)) != GRIB_SUCCESS)
return err;
/* ECC-477: Don't call grib_get_long_internal to suppress error message being output */
if ((err = grib_get_long(hand, self->ccsds_flags, &ccsds_flags)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->ccsds_block_size, &ccsds_block_size)) != GRIB_SUCCESS)
return err;
if ((err = grib_get_long_internal(hand, self->ccsds_rsi, &ccsds_rsi)) != GRIB_SUCCESS)
return err;
/* TODO: This should be called upstream */
if (*len < n_vals)
return GRIB_ARRAY_TOO_SMALL;
/* Special case */
if (bits_per_value == 0) {
for (i = 0; i < n_vals; i++)
val[i] = reference_value;
*len = n_vals;
return GRIB_SUCCESS;
}
bscale = grib_power(binary_scale_factor, 2);
dscale = grib_power(-decimal_scale_factor, 10);
buflen = grib_byte_count(a);
buf = (unsigned char*)hand->buffer->data;
buf += grib_byte_offset(a);
strm.flags = ccsds_flags;
strm.bits_per_sample = bits_per_value;
strm.block_size = ccsds_block_size;
strm.rsi = ccsds_rsi;
strm.next_in = buf;
strm.avail_in = buflen;
bits8 = ((bits_per_value + 7) / 8) * 8;
size = n_vals * ((bits_per_value + 7) / 8);
decoded = (unsigned char*)grib_context_buffer_malloc_clear(a->context, size);
if (!decoded) {
err = GRIB_OUT_OF_MEMORY;
goto cleanup;
}
strm.next_out = decoded;
strm.avail_out = size;
if (hand->context->debug) print_aec_stream_info(&strm, "unpack_*");
if ((err = aec_buffer_decode(&strm)) != AEC_OK) {
grib_context_log(a->context, GRIB_LOG_ERROR, "CCSDS %s: aec_buffer_decode error %d (%s)\n",
__PRETTY_FUNCTION__, err, aec_get_error_message(err));
err = GRIB_ENCODING_ERROR;
goto cleanup;
}
pos = 0;
/* ECC-1427: Performance improvement */
//grib_decode_float_array(decoded, &pos, bits8 , reference_value, bscale, dscale, n_vals, val);
grib_decode_array<T>(decoded, &pos, bits8 , reference_value, bscale, dscale, n_vals, val);
*len = n_vals;
cleanup:
grib_context_buffer_free(a->context, decoded);
return err;
}
#else
static int unpack(grib_accessor* a, float* val, size_t* len)
{
print_error_feature_not_enabled(a->context);
return GRIB_FUNCTIONALITY_NOT_ENABLED;
}
#endif
};

480
src/grib_accessor_class_data_complex_packing.cc
View File

@ -8,7 +8,7 @@
* virtue of its status as an intergovernmental organisation nor does it submit to any jurisdiction.
*/
#include "grib_api_internal.h"
#include "grib_accessor_class_data_complex_packing.h"
#include "grib_optimize_decimal_factor.h"
#include <math.h>
/*
@ -48,45 +48,10 @@ or edit "accessor.class" and rerun ./make_class.pl
*/
static int pack_double(grib_accessor*, const double* val, size_t* len);
static int unpack_double(grib_accessor*, double* val, size_t* len);
static int unpack_float(grib_accessor*, float* val, size_t* len);
static int value_count(grib_accessor*, long*);
static void init(grib_accessor*, const long, grib_arguments*);
static void init_class(grib_accessor_class*);
typedef struct grib_accessor_data_complex_packing
{
grib_accessor att;
/* Members defined in gen */
/* Members defined in values */
int carg;
const char* seclen;
const char* offsetdata;
const char* offsetsection;
int dirty;
/* Members defined in data_simple_packing */
int edition;
const char* units_factor;
const char* units_bias;
const char* changing_precision;
const char* number_of_values;
const char* bits_per_value;
const char* reference_value;
const char* binary_scale_factor;
const char* decimal_scale_factor;
const char* optimize_scaling_factor;
/* Members defined in data_complex_packing */
const char* GRIBEX_sh_bug_present;
const char* ieee_floats;
const char* laplacianOperatorIsSet;
const char* laplacianOperator;
const char* sub_j;
const char* sub_k;
const char* sub_m;
const char* pen_j;
const char* pen_k;
const char* pen_m;
} grib_accessor_data_complex_packing;
extern grib_accessor_class* grib_accessor_class_data_simple_packing;
@ -113,8 +78,8 @@ static grib_accessor_class _grib_accessor_class_data_complex_packing = {
0, /* grib_unpack procedures long */
&pack_double, /* grib_pack procedures double */
0, /* grib_pack procedures float */
&unpack_double, /* grib_unpack procedures double */
&unpack_float, /* grib_unpack procedures float */
&GribAccessorClassDataComplexPacking<double>::unpack, /* grib_unpack procedures double */
&GribAccessorClassDataComplexPacking<float>::unpack, /* grib_unpack procedures float */
0, /* grib_pack procedures string */
0, /* grib_unpack procedures string */
0, /* grib_pack array procedures string */
@ -181,8 +146,6 @@ static void init_class(grib_accessor_class* c)
/* END_CLASS_IMP */
typedef unsigned long (*encode_float_proc)(double);
typedef double (*decode_float_proc)(unsigned long);
static void init(grib_accessor* a, const long v, grib_arguments* args)
{
@ -234,443 +197,6 @@ static int value_count(grib_accessor* a, long* count)
return ret;
}
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);
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;
/* 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(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 = (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, 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 * (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];
/* 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++;
}
Assert(*len >= i);
*len = i;
grib_context_free(a->context, scals);
return ret;
}
// TODO(maee): ECC-1467
static int unpack_float(grib_accessor* a, float* 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;
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++;
}
Assert(*len >= i);
*len = i;
grib_context_free(a->context, scals);
return ret;
}
#define MAXVAL(a, b) a > b ? a : b
static double calculate_pfactor(grib_context* ctx, const double* spectralField, long fieldTruncation, long subsetTruncation)

269
src/grib_accessor_class_data_complex_packing.h Normal file
View File

@ -0,0 +1,269 @@
// ECC-1467
#pragma once
#include "grib_api_internal_cpp.h"
#include <type_traits>
typedef unsigned long (*encode_float_proc)(double);
typedef double (*decode_float_proc)(unsigned long);
typedef struct grib_accessor_data_complex_packing
{
grib_accessor att;
/* Members defined in gen */
/* Members defined in values */
int carg;
const char* seclen;
const char* offsetdata;
const char* offsetsection;
int dirty;
/* Members defined in data_simple_packing */
int edition;
const char* units_factor;
const char* units_bias;
const char* changing_precision;
const char* number_of_values;
const char* bits_per_value;
const char* reference_value;
const char* binary_scale_factor;
const char* decimal_scale_factor;
const char* optimize_scaling_factor;
/* Members defined in data_complex_packing */
const char* GRIBEX_sh_bug_present;
const char* ieee_floats;
const char* laplacianOperatorIsSet;
const char* laplacianOperator;
const char* sub_j;
const char* sub_k;
const char* sub_m;
const char* pen_j;
const char* pen_k;
const char* pen_m;
} grib_accessor_data_complex_packing;
template <typename T>
class GribAccessorClassDataComplexPacking {
public:
static int unpack(grib_accessor* a, T* val, size_t* len);
};
template <typename T>
int GribAccessorClassDataComplexPacking<T>::unpack(grib_accessor* a, T* 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;
T* scals = NULL;
T*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;
T 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<T>(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 = (T*)grib_context_malloc(a->context, maxv * sizeof(T));
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 * (T)((grib_decode_unsigned_long(lres, &lpos, bits_per_value) * s) + reference_value) * scals[lup];
val[i++] = d * (T)((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++;
}
Assert(*len >= i);
*len = i;
grib_context_free(a->context, scals);
return ret;
}

4
src/grib_accessor_class_data_raw_packing.cc
View File

@ -11,7 +11,7 @@
* Enrico Fucile
****************************/
#include "grib_api_internal.h"
#include "grib_api_internal_cpp.h"
#define PRE_PROCESSING_NONE 0
#define PRE_PROCESSING_DIFFERENCE 1
@ -210,7 +210,7 @@ static int unpack_double(grib_accessor* a, double* val, size_t* len)
if (*len < nvals)
return GRIB_ARRAY_TOO_SMALL;
code = grib_ieee_decode_array(a->context, buf, nvals, bytes, val);
code = grib_ieee_decode_array<double>(a->context, buf, nvals, bytes, val);
*len = nvals;

8
src/grib_accessor_class_data_simple_packing.cc
View File

@ -12,7 +12,7 @@
* Enrico Fucile
*******************************/
#include "grib_api_internal.h"
#include "grib_api_internal_cpp.h"
#include "grib_optimize_decimal_factor.h"
#include <float.h>
@ -442,9 +442,9 @@ static int _unpack_real(grib_accessor* a, double* dval, float* fval, size_t* len
"unpack_double: calling outline function : bpv %d, rv : %g, sf : %d, dsf : %d ",
bits_per_value, reference_value, binary_scale_factor, decimal_scale_factor);
if(dval)
grib_decode_double_array(buf, &pos, bits_per_value, reference_value, s, d, n_vals, dval);
grib_decode_array<double>(buf, &pos, bits_per_value, reference_value, s, d, n_vals, dval);
if(fval)
grib_decode_float_array(buf, &pos, bits_per_value, reference_value, s, d, n_vals, fval);
grib_decode_array<float>(buf, &pos, bits_per_value, reference_value, s, d, n_vals, fval);
*len = (long)n_vals;
@ -586,7 +586,7 @@ static int _unpack_double(grib_accessor* a, double* val, size_t* len, unsigned c
grib_context_log(a->context, GRIB_LOG_DEBUG,
"unpack_double: calling outline function : bpv %d, rv : %g, sf : %d, dsf : %d ",
bits_per_value, reference_value, binary_scale_factor, decimal_scale_factor);
grib_decode_double_array(buf, &pos, bits_per_value, reference_value, s, d, n_vals, val);
grib_decode_array<double>(buf, &pos, bits_per_value, reference_value, s, d, n_vals, val);
*len = (long)n_vals;

68
src/grib_accessor_class_gen.cc
View File

@ -13,7 +13,7 @@
* Enrico Fucile *
* Shahram Najm *
***************************************************************************/
#include "grib_api_internal.h"
#include "grib_accessor_class_gen.h"
/*
This is used by make_class.pl
@ -56,10 +56,6 @@ static int pack_string(grib_accessor*, const char*, size_t* len);
static int pack_string_array(grib_accessor*, const char**, size_t* len);
static int pack_expression(grib_accessor*, grib_expression*);
static int unpack_bytes(grib_accessor*, unsigned char*, size_t* len);
static int unpack_double(grib_accessor*, double* val, size_t* len);
static int unpack_float(grib_accessor*, float* val, size_t* len);
static int unpack_long(grib_accessor*, long* val, size_t* len);
static int unpack_string(grib_accessor*, char*, size_t* len);
static int unpack_string_array(grib_accessor*, char**, size_t* len);
static size_t string_length(grib_accessor*);
static long byte_count(grib_accessor*);
@ -108,13 +104,13 @@ static grib_accessor_class _grib_accessor_class_gen = {
0, /* grib_pack procedures long */
&is_missing, /* grib_pack procedures long */
&pack_long, /* grib_pack procedures long */
&unpack_long, /* grib_unpack procedures long */
&GribAccessorClassGen<long>::unpack, /* grib_unpack procedures long */
&pack_double, /* grib_pack procedures double */
0, /* grib_pack procedures float */
&unpack_double, /* grib_unpack procedures double */
&unpack_float, /* grib_unpack procedures float */
&GribAccessorClassGen<double>::unpack, /* grib_unpack procedures double */
&GribAccessorClassGen<float>::unpack, /* grib_unpack procedures float */
&pack_string, /* grib_pack procedures string */
&unpack_string, /* grib_unpack procedures string */
&GribAccessorClassGen<char>::unpack, /* grib_unpack procedures string */
&pack_string_array, /* grib_pack array procedures string */
&unpack_string_array, /* grib_unpack array procedures string */
&pack_bytes, /* grib_pack procedures bytes */
@ -274,10 +270,12 @@ static int clear(grib_accessor* a)
return GRIB_SUCCESS;
}
static int unpack_long(grib_accessor* a, long* v, size_t* len)
// ECC-1467
template <>
int GribAccessorClassGen<long>::unpack(grib_accessor* a, long* v, size_t* len)
{
int type = GRIB_TYPE_UNDEFINED;
if (a->cclass->unpack_double && a->cclass->unpack_double != &unpack_double) {
if (a->cclass->unpack_double && a->cclass->unpack_double != &GribAccessorClassGen<double>::unpack) {
double val = 0.0;
size_t l = 1;
grib_unpack_double(a, &val, &l);
@ -289,7 +287,7 @@ static int unpack_long(grib_accessor* a, long* v, size_t* len)
return GRIB_SUCCESS;
}
if (a->cclass->unpack_string && a->cclass->unpack_string != &unpack_string) {
if (a->cclass->unpack_string && a->cclass->unpack_string != &GribAccessorClassGen<char>::unpack) {
char val[1024];
size_t l = sizeof(val);
char* last = NULL;
@ -310,48 +308,18 @@ static int unpack_long(grib_accessor* a, long* v, size_t* len)
return GRIB_NOT_IMPLEMENTED;
}
static int unpack_double(grib_accessor* a, double* v, size_t* len)
{
int type = GRIB_TYPE_UNDEFINED;
if (a->cclass->unpack_long && a->cclass->unpack_long != &unpack_long) {
long val = 0;
size_t l = 1;
grib_unpack_long(a, &val, &l);
*v = val;
grib_context_log(a->context, GRIB_LOG_DEBUG, "Casting long %s to double", a->name);
return GRIB_SUCCESS;
}
if (a->cclass->unpack_string && a->cclass->unpack_string != &unpack_string) {
char val[1024];
size_t l = sizeof(val);
char* last = NULL;
grib_unpack_string(a, val, &l);
*v = strtod(val, &last);
if (*last == 0) { /* conversion of string to double worked */
grib_context_log(a->context, GRIB_LOG_DEBUG, "Casting string %s to long", a->name);
return GRIB_SUCCESS;
}
}
grib_context_log(a->context, GRIB_LOG_ERROR, "Cannot unpack %s as double", a->name);
if (grib_get_native_type(grib_handle_of_accessor(a), a->name, &type) == GRIB_SUCCESS) {
grib_context_log(a->context, GRIB_LOG_ERROR, "Hint: Try unpacking as %s", grib_get_type_name(type));
}
return GRIB_NOT_IMPLEMENTED;
}
// TODO(maee): ECC-1467: Copied the 'double' version and reused by copy/paste!
static int unpack_float(grib_accessor* a, float* v, size_t* len)
// ECC-1467
template <>
int GribAccessorClassGen<float>::unpack(grib_accessor* a, float* v, size_t* len)
{
return GRIB_NOT_IMPLEMENTED;
}
static int unpack_string(grib_accessor* a, char* v, size_t* len)
// ECC-1467
template <>
int GribAccessorClassGen<char>::unpack(grib_accessor* a, char* v, size_t* len)
{
if (a->cclass->unpack_double && a->cclass->unpack_double != &unpack_double) {
if (a->cclass->unpack_double && a->cclass->unpack_double != &GribAccessorClassGen<double>::unpack) {
double val = 0.0;
size_t l = 1;
grib_unpack_double(a, &val, &l);
@ -361,7 +329,7 @@ static int unpack_string(grib_accessor* a, char* v, size_t* len)
return GRIB_SUCCESS;
}
if (a->cclass->unpack_long && a->cclass->unpack_long != &unpack_long) {
if (a->cclass->unpack_long && a->cclass->unpack_long != &GribAccessorClassGen<long>::unpack) {
long val = 0;
size_t l = 1;
grib_unpack_long(a, &val, &l);

53
src/grib_accessor_class_gen.h Normal file
View File

@ -0,0 +1,53 @@
// ECC-1467
#pragma once
#include "grib_api_internal_cpp.h"
#include <type_traits>
#include <typeinfo>
template <typename T>
class GribAccessorClassGen {
public:
static int unpack(grib_accessor* a, T* v, size_t* len);
};
template <> int GribAccessorClassGen<long>::unpack(grib_accessor* a, long* v, size_t* len);
template <> int GribAccessorClassGen<char>::unpack(grib_accessor* a, char* v, size_t* len);
template <> int GribAccessorClassGen<float>::unpack(grib_accessor* a, float* v, size_t* len);
template <typename T>
int GribAccessorClassGen<T>::unpack(grib_accessor* a, T* v, size_t* len)
{
static_assert(std::is_floating_point<T>::value, "Requires floating point numbers");
int type = GRIB_TYPE_UNDEFINED;
if (a->cclass->unpack_long && a->cclass->unpack_long != &GribAccessorClassGen<long>::unpack) {
long val = 0;
size_t l = 1;
grib_unpack_long(a, &val, &l);
*v = val;
grib_context_log(a->context, GRIB_LOG_DEBUG, "Casting long %s to %s", a->name, typeid(T).name());
return GRIB_SUCCESS;
}
if (a->cclass->unpack_string && a->cclass->unpack_string != &GribAccessorClassGen<char>::unpack) {
char val[1024];
size_t l = sizeof(val);
char* last = NULL;
grib_unpack_string(a, val, &l);
*v = strtod(val, &last);
if (*last == 0) { /* conversion of string to double worked */
grib_context_log(a->context, GRIB_LOG_DEBUG, "Casting string %s to long", a->name);
return GRIB_SUCCESS;
}
}
grib_context_log(a->context, GRIB_LOG_ERROR, "Cannot unpack %s as %s", a->name, typeid(T).name());
if (grib_get_native_type(grib_handle_of_accessor(a), a->name, &type) == GRIB_SUCCESS) {
grib_context_log(a->context, GRIB_LOG_ERROR, "Hint: Try unpacking as %s", grib_get_type_name(type));
}
return GRIB_NOT_IMPLEMENTED;
}

8
src/grib_api_internal.h
View File

@ -16,7 +16,6 @@
#ifndef grib_api_internal_H
#define grib_api_internal_H
#ifdef __cplusplus
extern "C" {
#endif
@ -253,6 +252,8 @@ extern int pthread_mutexattr_settype(pthread_mutexattr_t* attr, int type);
#define GRIB_COMPARE_NAMES (1 << 0)
#define GRIB_COMPARE_TYPES (1 << 1)
extern const int max_nbits;
typedef struct grib_expression grib_expression;
typedef struct grib_arguments grib_arguments;
@ -1606,11 +1607,14 @@ typedef struct j2k_encode_helper
#include "eccodes_prototypes.h"
#ifdef __cplusplus
}
#endif
#endif
/* This part is automatically generated by ./errors.pl, do not edit */
#ifndef grib_errors_internal_H
#define grib_errors_internal_H

5
src/grib_api_internal_cpp.h Normal file
View File

@ -0,0 +1,5 @@
#pragma once
#include "grib_value.h"
#include "grib_bits_any_endian_simple.h"
#include "grib_ieeefloat.h"

2
src/grib_bits.cc
View File

@ -48,7 +48,7 @@ static const unsigned long dmasks[] = {
0x00,
};
static const int max_nbits = sizeof(unsigned long) * 8;
extern const int max_nbits = sizeof(unsigned long) * 8;
static const int max_nbits_size_t = sizeof(size_t) * 8;
unsigned long grib_decode_unsigned_byte_long(const unsigned char* p, long o, int l)

251
src/grib_bits_any_endian_simple.cc
View File

@ -13,10 +13,7 @@
* *
***************************************************************************/
/* A mask with x least-significant bits set, possibly 0 or >=32 */
/* -1UL is 1111111... in every bit in binary representation */
#define BIT_MASK1(x) \
(((x) >= max_nbits) ? (unsigned long)-1UL : (1UL << (x)) - 1)
#include "grib_bits_any_endian_simple.h"
/**
* decode an array of n_vals values from a octet-stream
@ -79,178 +76,90 @@ int grib_decode_long_array(const unsigned char* p, long* bitp, long bitsPerValue
* @param n_vals number of values to decode
* @param val output, values encoded as 32/64bit numbers
*/
int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerValue,
double reference_value, double s, double d,
size_t n_vals, double* val)
{
long i = 0;
unsigned long lvalue = 0;
double x;
//int grib_decode_double_array(const unsigned char* p, long* bitp, long bitsPerValue,
// double reference_value, double s, double d,
// size_t n_vals, double* val)
//{
// long i = 0;
// unsigned long lvalue = 0;
// double x;
#if 0
/* slow reference code */
int j=0;
for(i=0;i < n_vals;i++) {
lvalue=0;
for(j=0; j< bitsPerValue;j++){
lvalue <<= 1;
if(grib_get_bit( p, *bitp)) lvalue += 1;
*bitp += 1;
}
x=((lvalue*s)+reference_value)*d;
val[i] = (double)x;
}
#endif
if (bitsPerValue % 8 == 0) {
/* See ECC-386 */
int bc;
int l = bitsPerValue / 8;
size_t o = 0;
//#if 0
// [> slow reference code <]
// int j=0;
// for(i=0;i < n_vals;i++) {
// lvalue=0;
// for(j=0; j< bitsPerValue;j++){
// lvalue <<= 1;
// if(grib_get_bit( p, *bitp)) lvalue += 1;
// *bitp += 1;
// }
// x=((lvalue*s)+reference_value)*d;
// val[i] = (double)x;
// }
//#endif
// if (bitsPerValue % 8 == 0) {
// [> See ECC-386 <]
// int bc;
// int l = bitsPerValue / 8;
// size_t o = 0;
for (i = 0; i < n_vals; i++) {
lvalue = 0;
lvalue <<= 8;
lvalue |= p[o++];
// for (i = 0; i < n_vals; i++) {
// lvalue = 0;
// lvalue <<= 8;
// lvalue |= p[o++];
for (bc = 1; bc < l; bc++) {
lvalue <<= 8;
lvalue |= p[o++];
}
x = ((lvalue * s) + reference_value) * d;
val[i] = (double)x;
/* *bitp += bitsPerValue * n_vals; */
}
}
else {
unsigned long mask = BIT_MASK1(bitsPerValue);
// for (bc = 1; bc < l; bc++) {
// lvalue <<= 8;
// lvalue |= p[o++];
// }
// x = ((lvalue * s) + reference_value) * d;
// val[i] = (double)x;
// [> *bitp += bitsPerValue * n_vals; <]
// }
// }
// else {
// unsigned long mask = BIT_MASK1(bitsPerValue);
/* pi: position of bitp in p[]. >>3 == /8 */
long pi = *bitp / 8;
/* some bits might of the current byte at pi might be used */
/* by the previous number usefulBitsInByte gives remaining unused bits */
/* number of useful bits in current byte */
int usefulBitsInByte = 8 - (*bitp & 7);
for (i = 0; i < n_vals; i++) {
/* value read as long */
long bitsToRead = 0;
lvalue = 0;
bitsToRead = bitsPerValue;
/* read one byte after the other to lvalue until >= bitsPerValue are read */
while (bitsToRead > 0) {
lvalue <<= 8;
lvalue += p[pi];
pi++;
bitsToRead -= usefulBitsInByte;
usefulBitsInByte = 8;
}
*bitp += bitsPerValue;
/* bitsToRead is now <= 0, remove the last bits */
lvalue >>= -1 * bitsToRead;
/* set leading bits to 0 - removing bits used for previous number */
lvalue &= mask;
// [> pi: position of bitp in p[]. >>3 == /8 <]
// long pi = *bitp / 8;
// [> some bits might of the current byte at pi might be used <]
// [> by the previous number usefulBitsInByte gives remaining unused bits <]
// [> number of useful bits in current byte <]
// int usefulBitsInByte = 8 - (*bitp & 7);
// for (i = 0; i < n_vals; i++) {
// [> value read as long <]
// long bitsToRead = 0;
// lvalue = 0;
// bitsToRead = bitsPerValue;
// [> read one byte after the other to lvalue until >= bitsPerValue are read <]
// while (bitsToRead > 0) {
// lvalue <<= 8;
// lvalue += p[pi];
// pi++;
// bitsToRead -= usefulBitsInByte;
// usefulBitsInByte = 8;
// }
// *bitp += bitsPerValue;
// [> bitsToRead is now <= 0, remove the last bits <]
// lvalue >>= -1 * bitsToRead;
// [> set leading bits to 0 - removing bits used for previous number <]
// lvalue &= mask;
usefulBitsInByte = -1 * bitsToRead; /* prepare for next round */
if (usefulBitsInByte > 0) {
pi--; /* reread the current byte */
}
else {
usefulBitsInByte = 8; /* start with next full byte */
}
/* scaling and move value to output */
x = ((lvalue * s) + reference_value) * d;
val[i] = (double)x;
}
}
return 0;
}
// TODO(maee): ECC-1467: Copied the 'double' version and reused by copy/paste!
// This and the grib_decode_double_array function
// should be merged and refactored! Most probably using C++ templates
int grib_decode_float_array(const unsigned char* p, long* bitp, long bitsPerValue,
double reference_value, double s, double d,
size_t n_vals, float* val)
{
long i = 0;
unsigned long lvalue = 0;
double x;
#if 0
/* slow reference code */
int j=0;
for(i=0;i < n_vals;i++) {
lvalue=0;
for(j=0; j< bitsPerValue;j++){
lvalue <<= 1;
if(grib_get_bit( p, *bitp)) lvalue += 1;
*bitp += 1;
}
x=((lvalue*s)+reference_value)*d;
val[i] = (double)x;
}
#endif
if (bitsPerValue % 8 == 0) {
/* See ECC-386 */
int bc;
int l = bitsPerValue / 8;
size_t o = 0;
for (i = 0; i < n_vals; i++) {
lvalue = 0;
lvalue <<= 8;
lvalue |= p[o++];
for (bc = 1; bc < l; bc++) {
lvalue <<= 8;
lvalue |= p[o++];
}
x = ((lvalue * s) + reference_value) * d;
val[i] = (float)x;
/* *bitp += bitsPerValue * n_vals; */
}
}
else {
unsigned long mask = BIT_MASK1(bitsPerValue);
/* pi: position of bitp in p[]. >>3 == /8 */
long pi = *bitp / 8;
/* some bits might of the current byte at pi might be used */
/* by the previous number usefulBitsInByte gives remaining unused bits */
/* number of useful bits in current byte */
int usefulBitsInByte = 8 - (*bitp & 7);
for (i = 0; i < n_vals; i++) {
/* value read as long */
long bitsToRead = 0;
lvalue = 0;
bitsToRead = bitsPerValue;
/* read one byte after the other to lvalue until >= bitsPerValue are read */
while (bitsToRead > 0) {
lvalue <<= 8;
lvalue += p[pi];
pi++;
bitsToRead -= usefulBitsInByte;
usefulBitsInByte = 8;
}
*bitp += bitsPerValue;
/* bitsToRead is now <= 0, remove the last bits */
lvalue >>= -1 * bitsToRead;
/* set leading bits to 0 - removing bits used for previous number */
lvalue &= mask;
usefulBitsInByte = -1 * bitsToRead; /* prepare for next round */
if (usefulBitsInByte > 0) {
pi--; /* reread the current byte */
}
else {
usefulBitsInByte = 8; /* start with next full byte */
}
/* scaling and move value to output */
x = ((lvalue * s) + reference_value) * d;
val[i] = (float)x;
}
}
return 0;
}
// usefulBitsInByte = -1 * bitsToRead; [> prepare for next round <]
// if (usefulBitsInByte > 0) {
// pi--; [> reread the current byte <]
// }
// else {
// usefulBitsInByte = 8; [> start with next full byte <]
// }
// [> scaling and move value to output <]
// x = ((lvalue * s) + reference_value) * d;
// val[i] = (double)x;
// }
// }
// return 0;
//}
int grib_decode_double_array_complex(const unsigned char* p, long* bitp, long nbits, double reference_value, double s, double* d, size_t size, double* val)
{

97
src/grib_bits_any_endian_simple.h Normal file
View File

@ -0,0 +1,97 @@
// ECC-1467
#pragma once
#include "grib_api_internal.h"
/* A mask with x least-significant bits set, possibly 0 or >=32 */
/* -1UL is 1111111... in every bit in binary representation */
#define BIT_MASK1(x) \
(((x) >= max_nbits) ? (unsigned long)-1UL : (1UL << (x)) - 1)
template <typename T>
int grib_decode_array(const unsigned char* p, long* bitp, long bitsPerValue,
double reference_value, double s, double d,
size_t n_vals, T* val)
{
size_t i = 0;
unsigned long lvalue = 0;
double x;
#if 0
/* slow reference code */
int j=0;
for(i=0; i < n_vals; i++) {
lvalue=0;
for(j=0; j< bitsPerValue;j++){
lvalue <<= 1;
if(grib_get_bit( p, *bitp)) lvalue += 1;
*bitp += 1;
}
x=((lvalue*s)+reference_value)*d;
val[i] = (double)x;
}
#endif
if (bitsPerValue % 8 == 0) {
/* See ECC-386 */
int bc;
int l = bitsPerValue / 8;
size_t o = 0;
for (i = 0; i < n_vals; i++) {
lvalue = 0;
lvalue <<= 8;
lvalue |= p[o++];
for (bc = 1; bc < l; bc++) {
lvalue <<= 8;
lvalue |= p[o++];
}
x = ((lvalue * s) + reference_value) * d;
val[i] = (T)x;
/* *bitp += bitsPerValue * n_vals; */
}
}
else {
unsigned long mask = BIT_MASK1(bitsPerValue);
/* pi: position of bitp in p[]. >>3 == /8 */
long pi = *bitp / 8;
/* some bits might of the current byte at pi might be used */
/* by the previous number usefulBitsInByte gives remaining unused bits */
/* number of useful bits in current byte */
int usefulBitsInByte = 8 - (*bitp & 7);
for (i = 0; i < n_vals; i++) {
/* value read as long */
long bitsToRead = 0;
lvalue = 0;
bitsToRead = bitsPerValue;
/* read one byte after the other to lvalue until >= bitsPerValue are read */
while (bitsToRead > 0) {
lvalue <<= 8;
lvalue += p[pi];
pi++;
bitsToRead -= usefulBitsInByte;
usefulBitsInByte = 8;
}
*bitp += bitsPerValue;
/* bitsToRead is now <= 0, remove the last bits */
lvalue >>= -1 * bitsToRead;
/* set leading bits to 0 - removing bits used for previous number */
lvalue &= mask;
usefulBitsInByte = -1 * bitsToRead; /* prepare for next round */
if (usefulBitsInByte > 0) {
pi--; /* reread the current byte */
}
else {
usefulBitsInByte = 8; /* start with next full byte */
}
/* scaling and move value to output */
x = ((lvalue * s) + reference_value) * d;
val[i] = (T)x;
}
}
return 0;
}

8
src/grib_ieeefloat.cc
View File

@ -12,7 +12,7 @@
* Enrico Fucile - 06.01.2009 *
* *
***************************************************************************/
#include "grib_api_internal.h"
#include "grib_ieeefloat.h"
#if GRIB_PTHREADS
static pthread_once_t once = PTHREAD_ONCE_INIT;
@ -344,7 +344,8 @@ double grib_long_to_ieee64(unsigned long x)
return dval;
}
int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, double* val)
template <>
int grib_ieee_decode_array<double> (grib_context* c, unsigned char* buf, size_t nvals, int bytes, double* val)
{
int err = 0, i = 0, j = 0;
unsigned char s[8] = {0,};
@ -387,7 +388,8 @@ 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)
template <>
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,};

5
src/grib_ieeefloat.h Normal file
View File

@ -0,0 +1,5 @@
#pragma once
#include "grib_api_internal.h"
template <typename T> int grib_ieee_decode_array(grib_context* c, unsigned char* buf, size_t nvals, int bytes, T* val);

76
src/grib_value.cc
View File

@ -12,6 +12,7 @@
* Enrico Fucile *
***************************************************************************/
#include "grib_api_internal.h"
#include "grib_value.h"
#include <float.h>
/* Note: A fast cut-down version of strcmp which does NOT return -1 */
@ -1279,10 +1280,11 @@ const char* grib_get_accessor_class_name(grib_handle* h, const char* name)
return act ? act->cclass->name : NULL;
}
static int _grib_get_double_array_internal(const grib_handle* h, grib_accessor* a, double* val, size_t buffer_len, size_t* decoded_length)
template <>
int _grib_get_array_internal<double>(const grib_handle* h, grib_accessor* a, double* val, size_t buffer_len, size_t* decoded_length)
{
if (a) {
int err = _grib_get_double_array_internal(h, a->same, val, buffer_len, decoded_length);
int err = _grib_get_array_internal<double>(h, a->same, val, buffer_len, decoded_length);
if (err == GRIB_SUCCESS) {
size_t len = buffer_len - *decoded_length;
@ -1297,10 +1299,11 @@ static int _grib_get_double_array_internal(const grib_handle* h, grib_accessor*
}
}
int _grib_get_float_array_internal(const grib_handle* h, grib_accessor* a, float* val, size_t buffer_len, size_t* decoded_length)
template <>
int _grib_get_array_internal<float>(const grib_handle* h, grib_accessor* a, float* val, size_t buffer_len, size_t* decoded_length)
{
if (a) {
int err = _grib_get_float_array_internal(h, a->same, val, buffer_len, decoded_length);
int err = _grib_get_array_internal<float>(h, a->same, val, buffer_len, decoded_length);
if (err == GRIB_SUCCESS) {
size_t len = buffer_len - *decoded_length;
@ -1317,28 +1320,35 @@ int _grib_get_float_array_internal(const grib_handle* h, grib_accessor* a, float
int grib_get_double_array_internal(const grib_handle* h, const char* name, double* val, size_t* length)
{
int ret = grib_get_double_array(h, name, val, length);
if (ret != GRIB_SUCCESS)
grib_context_log(h->context, GRIB_LOG_ERROR,
"unable to get %s as double array (%s)",
name, grib_get_error_message(ret));
return ret;
return grib_get_array_internal<double>(h, name, val, length);
}
int grib_get_float_array_internal(const grib_handle* h, const char* name, float* val, size_t* length)
{
int ret = grib_get_float_array(h, name, val, length);
if (ret != GRIB_SUCCESS)
grib_context_log(h->context, GRIB_LOG_ERROR,
"unable to get %s as float array (%s)",
name, grib_get_error_message(ret));
return ret;
return grib_get_array_internal<float>(h, name, val, length);
}
int grib_get_double_array(const grib_handle* h, const char* name, double* val, size_t* length)
template <>
int grib_get_array<float>(const grib_handle* h, const char* name, float* val, size_t *length)
{
size_t len = *length;
grib_accessor* a = grib_find_accessor(h, name);
if (!a) return GRIB_NOT_FOUND;
//[> TODO: For now only GRIB supported... no BUFR keys <]
if (h->product_kind != PRODUCT_GRIB) {
grib_context_log(h->context, GRIB_LOG_ERROR, "grib_get_float_array only supported for GRIB");
return GRIB_NOT_IMPLEMENTED;
}
Assert(name[0]!='/');
Assert(name[0]!='#');
*length = 0;
return _grib_get_array_internal<float>(h,a,val,len,length);
}
template <>
int grib_get_array<double>(const grib_handle* h, const char* name, double* val, size_t* length)
{
size_t len = *length;
grib_accessor* a = NULL;
@ -1358,31 +1368,23 @@ int grib_get_double_array(const grib_handle* h, const char* name, double* val, s
if (!a)
return GRIB_NOT_FOUND;
if (name[0] == '#') {
return grib_unpack_double(a, val, length);
return grib_unpack_double(a, val, length); // TODO: change to template
}
else {
*length = 0;
return _grib_get_double_array_internal(h, a, val, len, length);
return _grib_get_array_internal<double>(h, a, val, len, length);
}
}
}
// TODO(maee): ECC-1467: Copied the 'double' version and reused by copy/paste!
int grib_get_double_array(const grib_handle* h, const char* name, double* val, size_t* length)
{
return grib_get_array<double>(h, name, val, length);
}
int grib_get_float_array(const grib_handle* h, const char* name, float* val, size_t *length)
{
size_t len = *length;
grib_accessor* a = grib_find_accessor(h, name);
if (!a) return GRIB_NOT_FOUND;
/* TODO: For now only GRIB supported... no BUFR keys */
if (h->product_kind != PRODUCT_GRIB) {
grib_context_log(h->context, GRIB_LOG_ERROR, "grib_get_float_array only supported for GRIB");
return GRIB_NOT_IMPLEMENTED;
}
Assert(name[0]!='/');
Assert(name[0]!='#');
*length = 0;
return _grib_get_float_array_internal(h,a,val,len,length);
return grib_get_array<float>(h, name, val, length);
}
int ecc__grib_get_string_length(grib_accessor* a, size_t* size)

24
src/grib_value.h Normal file
View File

@ -0,0 +1,24 @@
#pragma once
#include "grib_api_internal.h"
#include <typeinfo>
#include <type_traits>
template <typename T>
int grib_get_array(const grib_handle* h, const char* name, T* val, size_t* length);
template <typename T>
int _grib_get_array_internal(const grib_handle* h, grib_accessor* a, T* val, size_t buffer_len, size_t* decoded_length);
template <typename T>
int grib_get_array_internal(const grib_handle* h, const char* name, T* val, size_t* length)
{
int ret = grib_get_array<T>(h, name, val, length);
if (ret != GRIB_SUCCESS)
grib_context_log(h->context, GRIB_LOG_ERROR,
"unable to get %s as %s array (%s)",
name, typeid(T).name(), grib_get_error_message(ret));
return ret;
}