Commit 8cad1156 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Module Expr: added missing value support

parent 1c68b574
2010-07-22 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* Module Expr: added missing value support [request: Marco van Hulten]
2010-07-13 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* gridFromNCfile: check grid_dims
......
......@@ -7,6 +7,8 @@ Version 1.4.6 (?? September 2010):
* bandpass: Bandpass filtering
* lowpass: Lowpass filtering
* highpass: Highpass filtering
Changed operators:
* expr, exprf: added missing values support
Fixed bugs:
* seldate: open output file only when time steps found
......
......@@ -24,7 +24,7 @@
Exprf aexprf Append evaluated expressions from script file
*/
/*
Operatoren: +, -, *, \
Operatoren: +, -, *, \, ^
Functions: sqrt, exp, log, log10, sin, cos, tan, asin, acos, atan
Functions: min, max, avg, std, var
Constansts: M_PI, M_E
......
......@@ -58,7 +58,7 @@ double corr(double * restrict in0, double * restrict in1,
out = wsum0 ?
FDIV((sum01 * wsum0 - sum0 * sum1),
FROOT((sum00 * wsum0 - sum0 * sum0) *
FSQRT((sum00 * wsum0 - sum0 * sum0) *
(sum11 * wsum0 - sum1 * sum1))) : missval;
return (out);
......
......@@ -144,7 +144,7 @@ void *Math(void *argument)
break;
case SQRT:
for ( i = 0; i < gridsize; i++ )
array2[i] = DBL_IS_EQUAL(array1[i], missval1) ? missval1 : ROOT(array1[i]);
array2[i] = DBL_IS_EQUAL(array1[i], missval1) ? missval1 : SQRT(array1[i]);
break;
case EXP:
for ( i = 0; i < gridsize; i++ )
......
......@@ -208,7 +208,7 @@ void *Timstat2(void *argument)
temp5 = SUB(work[varID][levelID][3].ptr[i], DIV(temp3, nvals2));
temp6 = MUL(temp4, temp5);
work[varID][levelID][0].ptr[i] = DIV(temp1, ROOT(temp6));
work[varID][levelID][0].ptr[i] = DIV(temp1, SQRT(temp6));
/*
if ( work[varID][levelID][0].ptr[i] < -1)
work[varID][levelID][0].ptr[i] = -1;
......
......@@ -283,7 +283,7 @@ void *Timstat3(void *argument)
if ( !DBL_IS_EQUAL(temp0, missval1) && temp0 < 0 ) /* This is possible because */
temp0 = 0; /* of rounding errors */
stddev_estimator = ROOT(DIV(temp0, deg_of_freedom));
stddev_estimator = SQRT(DIV(temp0, deg_of_freedom));
mean_estimator = -rconst;
for ( j = 0; j < n_in; j++ )
{
......@@ -301,7 +301,7 @@ void *Timstat3(void *argument)
var_factor[j]), fnvals));
}
norm = ROOT(temp1);
norm = SQRT(temp1);
temp2 = DIV(DIV(mean_estimator, norm), stddev_estimator);
fractil = deg_of_freedom < 1 ? missval1 :
......
......@@ -54,7 +54,7 @@ void trms(field_t field1, field_t field2, double *dp, field_t *field3)
rsumw = ADD(rsumw, wp);
}
ravg = ROOT(DIV(rsum, rsumw));
ravg = SQRT(DIV(rsum, rsumw));
if ( DBL_IS_EQUAL(ravg, missval1) ) rnmiss++;
......
......@@ -2,6 +2,7 @@
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <errno.h>
#include "cdi.h"
#include "cdo.h"
......@@ -155,10 +156,9 @@ nodeType *expr_con_var(int oper, nodeType *p1, nodeType *p2)
case '^':
if ( nmiss > 0 )
{
UMISS(func, oper);
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->u.con.value, p2->data[i+k*ngp]);
p->data[i+k*ngp] = POW(p1->u.con.value, p2->data[i+k*ngp]);
}
else
{
......@@ -275,10 +275,9 @@ nodeType *expr_var_con(int oper, nodeType *p1, nodeType *p2)
case '^':
if ( nmiss > 0 )
{
UMISS(func, oper);
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->data[i+k*ngp], p2->u.con.value);
p->data[i+k*ngp] = POW(p1->data[i+k*ngp], p2->u.con.value);
}
else
{
......@@ -311,6 +310,13 @@ nodeType *expr_var_var(int oper, nodeType *p1, nodeType *p2)
int ngp, ngp1, ngp2, i;
int nlev, nlev1, nlev2, k;
int loff1, loff2;
int nmiss, nmiss1, nmiss2;
double missval1, missval2;
nmiss1 = p1->nmiss;
nmiss2 = p2->nmiss;
missval1 = p1->missval;
missval2 = p2->missval;
ngp1 = gridInqSize(p1->gridID);
ngp2 = gridInqSize(p2->gridID);
......@@ -366,30 +372,77 @@ nodeType *expr_var_var(int oper, nodeType *p1, nodeType *p2)
switch ( oper )
{
case '+':
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] + p2->data[i+loff2];
if ( nmiss1 > 0 || nmiss2 > 0 )
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = ADD(p1->data[i+loff1], p2->data[i+loff2]);
}
else
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] + p2->data[i+loff2];
}
break;
case '-':
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] - p2->data[i+loff2];
if ( nmiss1 > 0 || nmiss2 > 0 )
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = SUB(p1->data[i+loff1], p2->data[i+loff2]);
}
else
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] - p2->data[i+loff2];
}
break;
case '*':
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] * p2->data[i+loff2];
if ( nmiss1 > 0 || nmiss2 > 0 )
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = MUL(p1->data[i+loff1], p2->data[i+loff2]);
}
else
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] * p2->data[i+loff2];
}
break;
case '/':
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] / p2->data[i+loff2];
if ( nmiss1 > 0 || nmiss2 > 0 )
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = DIV(p1->data[i+loff1], p2->data[i+loff2]);
}
else
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = p1->data[i+loff1] / p2->data[i+loff2];
}
break;
case '^':
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->data[i+loff1], p2->data[i+loff2]);
if ( nmiss1 > 0 || nmiss2 > 0 )
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = POW(p1->data[i+loff1], p2->data[i+loff2]);
}
else
{
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = pow(p1->data[i+loff1], p2->data[i+loff2]);
}
break;
default:
cdoAbort("%s: operator %c unsupported!", func, oper);
}
}
nmiss = 0;
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
if ( DBL_IS_EQUAL(p->data[i+k*ngp], missval1) ) nmiss++;
p->nmiss = nmiss;
if ( p1->tmpvar ) free(p1->data);
if ( p2->tmpvar ) free(p2->data);
......@@ -491,10 +544,12 @@ nodeType *ex_fun_var(char *fun, nodeType *p1)
int nlev, k;
int gridID, zaxisID;
int funcID = -1;
int nmiss;
double missval;
gridID = p1->gridID;
zaxisID = p1->zaxisID;
nmiss = p1->nmiss;
missval = p1->missval;
ngp = gridInqSize(gridID);
......@@ -521,19 +576,34 @@ nodeType *ex_fun_var(char *fun, nodeType *p1)
if ( funcID == -1 )
cdoAbort("function %s not available!", fun);
if ( strcmp("log", fun) == 0 )
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = log(p1->data[i+k*ngp]);
{
errno = -1;
p->data[i+k*ngp] = DBL_IS_EQUAL(p1->data[i+k*ngp], missval) ? missval : fun_sym_tbl[funcID].func(p1->data[i+k*ngp]);
if ( errno == EDOM || errno == ERANGE ) p->data[i+k*ngp] = missval;
}
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = fun_sym_tbl[funcID].func(p1->data[i+k*ngp]);
{
errno = -1;
p->data[i+k*ngp] = fun_sym_tbl[funcID].func(p1->data[i+k*ngp]);
if ( errno == EDOM || errno == ERANGE ) p->data[i+k*ngp] = missval;
}
}
nmiss = 0;
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
if ( DBL_IS_EQUAL(p->data[i+k*ngp], missval) ) nmiss++;
p->nmiss = nmiss;
if ( p1->tmpvar ) free(p1->data);
return (p);
......@@ -567,12 +637,14 @@ nodeType *ex_uminus_var(nodeType *p1)
nodeType *p;
int ngp, i;
int nlev, k;
int nmiss;
int gridID, zaxisID;
double missval;
gridID = p1->gridID;
zaxisID = p1->zaxisID;
missval = p1->missval;
gridID = p1->gridID;
zaxisID = p1->zaxisID;
nmiss = p1->nmiss;
missval = p1->missval;
ngp = gridInqSize(gridID);
nlev = zaxisInqSize(zaxisID);
......@@ -588,10 +660,21 @@ nodeType *ex_uminus_var(nodeType *p1)
p->data = (double *) malloc(ngp*nlev*sizeof(double));
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = -(p1->data[i+k*ngp]);
if ( nmiss > 0 )
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = DBL_IS_EQUAL(p1->data[i+k*ngp], missval) ? missval : -(p1->data[i+k*ngp]);
}
else
{
for ( k = 0; k < nlev; k++ )
for ( i = 0; i < ngp; i++ )
p->data[i+k*ngp] = -(p1->data[i+k*ngp]);
}
p->nmiss = nmiss;
return (p);
}
......
......@@ -42,9 +42,14 @@ double _FDIV_(double x, double y, double missval1, double missval2)
return FDIV(x, y);
}
double _FROOT_(double x, double missval1)
double _FPOW_(double x, double y, double missval1, double missval2)
{
return FROOT(x);
return FPOW(x, y);
}
double _FSQRT_(double x, double missval1)
{
return FSQRT(x);
}
......@@ -327,7 +332,7 @@ void fldrms(field_t field, field_t field2, field_t *field3)
}
*/
ravg = ROOT(DIV(rsum, rsumw));
ravg = SQRT(DIV(rsum, rsumw));
if ( DBL_IS_EQUAL(ravg, missval1) ) rnmiss++;
......@@ -381,7 +386,7 @@ void varrms(field_t field, field_t field2, field_t *field3)
}
*/
ravg = ROOT(DIV(rsum, rsumw));
ravg = SQRT(DIV(rsum, rsumw));
if ( DBL_IS_EQUAL(ravg, missval1) ) rnmiss++;
......
......@@ -28,21 +28,24 @@
#define FSUB(x,y) (DBL_IS_EQUAL((x),missval1) || DBL_IS_EQUAL((y),missval2) ? missval1 : (x)-(y))
#define FMUL(x,y) (DBL_IS_EQUAL((x),0.)||IS_EQUAL((y),0.) ? 0 : DBL_IS_EQUAL((x),missval1) || DBL_IS_EQUAL((y),missval2) ? missval1 : (x)*(y))
#define FDIV(x,y) (DBL_IS_EQUAL((x),missval1) || DBL_IS_EQUAL((y),missval2) || DBL_IS_EQUAL((y),0.) ? missval1 : (x)/(y))
#define FROOT(x) (DBL_IS_EQUAL((x),missval1) || (x)<0 ? missval1 : sqrt(x))
#define FPOW(x,y) (DBL_IS_EQUAL((x),missval1) || DBL_IS_EQUAL((y),missval2) ? missval1 : pow((x),(y)))
#define FSQRT(x) (DBL_IS_EQUAL((x),missval1) || (x)<0 ? missval1 : sqrt(x))
double _FADD_(double x, double y, double missval1, double missval2);
double _FSUB_(double x, double y, double missval1, double missval2);
double _FMUL_(double x, double y, double missval1, double missval2);
double _FDIV_(double x, double y, double missval1, double missval2);
double _FROOT_(double x, double missval1);
double _FPOW_(double x, double y, double missval1, double missval2);
double _FSQRT_(double x, double missval1);
#define ADD(x,y) _FADD_(x, y, missval1, missval2)
#define SUB(x,y) _FSUB_(x, y, missval1, missval2)
#define MUL(x,y) _FMUL_(x, y, missval1, missval2)
#define DIV(x,y) _FDIV_(x, y, missval1, missval2)
#define ROOT(x) _FROOT_(x, missval1)
#define POW(x,y) _FPOW_(x, y, missval1, missval2)
#define SQRT(x) _FSQRT_(x, missval1)
typedef struct {
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment