ECC-551: Standard deviation is shown as NAN

This commit is contained in:
Shahram Najm 2017-10-05 11:38:00 +01:00
parent 3fc6ba2265
commit 41ab2f7fb0
1 changed files with 99 additions and 100 deletions

View File

@ -150,151 +150,150 @@ static void init_class(grib_accessor_class* c)
static void init(grib_accessor* a,const long l, grib_arguments* c)
{
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
int n = 0;
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
int n = 0;
self->values = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->J = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->K = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->M = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->JS = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
a->flags |= GRIB_ACCESSOR_FLAG_READ_ONLY;
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
self->values = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->J = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->K = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->M = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->JS = grib_arguments_get_name(grib_handle_of_accessor(a),c,n++);
self->number_of_elements=4;
self->v=(double*)grib_context_malloc(a->context,
sizeof(double)*self->number_of_elements);
a->flags |= GRIB_ACCESSOR_FLAG_READ_ONLY;
a->flags |= GRIB_ACCESSOR_FLAG_FUNCTION;
a->length=0;
a->dirty=1;
self->number_of_elements=4;
self->v=(double*)grib_context_malloc(a->context,
sizeof(double)*self->number_of_elements);
a->length=0;
a->dirty=1;
}
static int unpack_double (grib_accessor* a, double* val, size_t *len)
static int unpack_double (grib_accessor* a, double* val, size_t *len)
{
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
int ret = 0,i=0;
double* values;
size_t size=0;
long J,K,M,N;
double avg,enorm,sd;
grib_context* c=a->context;
grib_handle* h=grib_handle_of_accessor(a);
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
int ret = 0,i=0;
double* values;
size_t size=0;
long J,K,M,N;
double avg,enorm,sd;
grib_context* c=a->context;
grib_handle* h=grib_handle_of_accessor(a);
if (!a->dirty) return GRIB_SUCCESS;
if (!a->dirty) return GRIB_SUCCESS;
if ( (ret=grib_get_size(h,self->values,&size)) != GRIB_SUCCESS) return ret;
if ( (ret=grib_get_size(h,self->values,&size)) != GRIB_SUCCESS) return ret;
if((ret=grib_get_long(grib_handle_of_accessor(a),self->J,&J))
!= GRIB_SUCCESS) return ret;
if((ret=grib_get_long(grib_handle_of_accessor(a),self->J,&J))
!= GRIB_SUCCESS) return ret;
if((ret=grib_get_long(grib_handle_of_accessor(a),self->K,&K))
!= GRIB_SUCCESS) return ret;
if((ret=grib_get_long(grib_handle_of_accessor(a),self->K,&K))
!= GRIB_SUCCESS) return ret;
if((ret=grib_get_long(grib_handle_of_accessor(a),self->M,&M))
!= GRIB_SUCCESS) return ret;
if((ret=grib_get_long(grib_handle_of_accessor(a),self->M,&M))
!= GRIB_SUCCESS) return ret;
if (J != M || M != K) return GRIB_NOT_IMPLEMENTED;
if (J != M || M != K) return GRIB_NOT_IMPLEMENTED;
N=(M+1)*(M+2)/2;
N=(M+1)*(M+2)/2;
if (2*N != size) {
grib_context_log(a->context,GRIB_LOG_ERROR,
"wrong number of components for spherical harmonics %ld != %ld",2*N,size);
return GRIB_WRONG_ARRAY_SIZE;
}
values=(double*)grib_context_malloc(c,size*sizeof(double));
if (!values) return GRIB_OUT_OF_MEMORY;
if (2*N != size) {
grib_context_log(a->context,GRIB_LOG_ERROR,
"wrong number of components for spherical harmonics %ld != %ld",2*N,size);
return GRIB_WRONG_ARRAY_SIZE;
}
if((ret = grib_get_double_array_internal(h,self->values,values,&size))
!= GRIB_SUCCESS) {
values=(double*)grib_context_malloc(c,size*sizeof(double));
if (!values) return GRIB_OUT_OF_MEMORY;
if((ret = grib_get_double_array_internal(h,self->values,values,&size))
!= GRIB_SUCCESS) {
grib_context_free(c,values);
return ret;
}
}
avg=values[0];
enorm=0;
sd=0;
avg=values[0];
enorm=0;
sd=0;
for (i=2; i < 2*J ; i+=2)
sd += values[i]*values[i];
for (i= 2*J; i < size; i+=2)
sd += values[i]*values[i] - values[i+1]*values[i+1];
for (i=2; i < 2*J ; i+=2)
sd += values[i]*values[i];
enorm=sd+avg*avg;
for (i= 2*J; i < size; i+=2)
sd += 2*values[i]*values[i] + 2*values[i+1]*values[i+1];
sd=sqrt(sd);
enorm=sqrt(enorm);
enorm=sd+avg*avg;
a->dirty=0;
sd=sqrt(sd);
enorm=sqrt(enorm);
grib_context_free(c,values);
self->v[0]=avg;
self->v[1]=enorm;
self->v[2]=sd;
self->v[3]= sd == 0 ? 1 : 0;
a->dirty=0;
for (i=0;i<self->number_of_elements;i++)
val[i]=self->v[i];
grib_context_free(c,values);
return ret;
self->v[0]=avg;
self->v[1]=enorm;
self->v[2]=sd;
self->v[3]= sd == 0 ? 1 : 0;
for (i=0;i<self->number_of_elements;i++)
val[i]=self->v[i];
return ret;
}
static int value_count(grib_accessor* a,long* count)
{
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
*count= self->number_of_elements;
return 0;
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
*count= self->number_of_elements;
return 0;
}
static void destroy(grib_context* c,grib_accessor* a)
{
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
grib_context_free(c,self->v);
grib_accessor_statistics_spectral* self = (grib_accessor_statistics_spectral*)a;
grib_context_free(c,self->v);
}
static int compare(grib_accessor* a, grib_accessor* b)
{
int retval=GRIB_SUCCESS;
double *aval=0;
double *bval=0;
int retval=GRIB_SUCCESS;
double *aval=0;
double *bval=0;
size_t alen = 0;
size_t blen = 0;
int err=0;
long count=0;
size_t alen = 0;
size_t blen = 0;
int err=0;
long count=0;
err=grib_value_count(a,&count);
if (err) return err;
alen=count;
err=grib_value_count(a,&count);
if (err) return err;
alen=count;
err=grib_value_count(b,&count);
if (err) return err;
blen=count;
err=grib_value_count(b,&count);
if (err) return err;
blen=count;
if (alen != blen) return GRIB_COUNT_MISMATCH;
if (alen != blen) return GRIB_COUNT_MISMATCH;
aval=(double*)grib_context_malloc(a->context,alen*sizeof(double));
bval=(double*)grib_context_malloc(b->context,blen*sizeof(double));
aval=(double*)grib_context_malloc(a->context,alen*sizeof(double));
bval=(double*)grib_context_malloc(b->context,blen*sizeof(double));
b->dirty=1;
a->dirty=1;
b->dirty=1;
a->dirty=1;
grib_unpack_double(a,aval,&alen);
grib_unpack_double(b,bval,&blen);
grib_unpack_double(a,aval,&alen);
grib_unpack_double(b,bval,&blen);
retval = GRIB_SUCCESS;
while (alen != 0) {
if (*bval != *aval) retval = GRIB_DOUBLE_VALUE_MISMATCH;
alen--;
}
retval = GRIB_SUCCESS;
while (alen != 0) {
if (*bval != *aval) retval = GRIB_DOUBLE_VALUE_MISMATCH;
alen--;
}
grib_context_free(a->context,aval);
grib_context_free(b->context,bval);
grib_context_free(a->context,aval);
grib_context_free(b->context,bval);
return retval;
return retval;
}