Commit 5294f739 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Added new operator ensstd1 and ensvar1

parent 936b4b4c
......@@ -3,6 +3,11 @@
* using CDI library version 1.6.0
* Version 1.6.0 released
2013-01-23 Uwe Schulzweida <Uwe.Schulzweida@zmaw.de>
* New operator: ensstd1 - Ensemble standard deviation [Divisor is (n-1)]
* New operator: ensvar1 - Ensemble variance [Divisor is (n-1)]
2013-01-22 Helmut Haak <Helmut.Haak@zmaw.de>
* New operator: adisit - Potential temperature to in-situ temperature
......
......@@ -216,8 +216,10 @@ Operator catalog:
Ensstat enssum Ensemble sum
Ensstat ensmean Ensemble mean
Ensstat ensavg Ensemble average
Ensstat ensvar Ensemble variance
Ensstat ensstd Ensemble standard deviation
Ensstat ensstd1 Ensemble standard deviation
Ensstat ensvar Ensemble variance
Ensstat ensvar1 Ensemble variance
Ensstat enspctl Ensemble percentiles
Ensstat2 ensrkhistspace Ranked Histogram averaged over time
Ensstat2 ensrkhisttime Ranked Histogram averaged over space
......
......@@ -288,8 +288,8 @@ while (<MOFILE>) {
print TRCARD "sum & {\\bf sum} \\\\ \n";
print TRCARD "mean & {\\bf mean} \\\\ \n";
print TRCARD "average & {\\bf avg} \\\\ \n";
print TRCARD "variance & {\\bf var} \\\\ \n";
print TRCARD "standard deviation & {\\bf std} \\\\ \n";
print TRCARD "variance & {\\bf var, var1} \\\\ \n";
print TRCARD "standard deviation & {\\bf std, std1} \\\\ \n";
print TRCARD "\\hline\n";
print TRCARD "\\end{tabular*}\n";
print TRCARD "\\vspace{2mm}\n\n";
......
......@@ -5,7 +5,7 @@
@Section = Statistical values
@Class = Statistic
@Arguments = ifiles ofile
@Operators = ensmin ensmax enssum ensmean ensavg ensvar ensstd enspctl
@Operators = ensmin ensmax enssum ensmean ensavg ensstd ensstd1 ensvar ensvar1 enspctl
@BeginDescription
This module computes statistical values over an ensemble of input files.
......@@ -97,10 +97,48 @@ o(t,x) = \mbox{\bf avg}\{i_1(t,x), i_2(t,x), \cdots, i_n(t,x)\}
@EndOperator
@BeginOperator_ensstd
@Title = Ensemble standard deviation
@BeginDescription
Divisor is n.
@IfMan
o(t,x) = std{i1(t,x), i2(t,x), ..., in(t,x)}
@EndifMan
@IfDoc
@BeginMath
o(t,x) = \mbox{\bf std}\{i_1(t,x), i_2(t,x), \cdots, i_n(t,x)\}
@EndMath
@EndifDoc
@EndDescription
@EndOperator
@BeginOperator_ensstd1
@Title = Ensemble standard deviation
@BeginDescription
Divisor is (n-1).
@IfMan
o(t,x) = std1{i1(t,x), i2(t,x), ..., in(t,x)}
@EndifMan
@IfDoc
@BeginMath
o(t,x) = \mbox{\bf std1}\{i_1(t,x), i_2(t,x), \cdots, i_n(t,x)\}
@EndMath
@EndifDoc
@EndDescription
@EndOperator
@BeginOperator_ensvar
@Title = Ensemble variance
@BeginDescription
Divisor is n.
@IfMan
o(t,x) = var{i1(t,x), i2(t,x), ..., in(t,x)}
@EndifMan
......@@ -113,16 +151,18 @@ o(t,x) = \mbox{\bf var}\{i_1(t,x), i_2(t,x), \cdots, i_n(t,x)\}
@EndOperator
@BeginOperator_ensstd
@Title = Ensemble standard deviation
@BeginOperator_ensvar1
@Title = Ensemble variance
@BeginDescription
Divisor is (n-1).
@IfMan
o(t,x) = std{i1(t,x), i2(t,x), ..., in(t,x)}
o(t,x) = var1{i1(t,x), i2(t,x), ..., in(t,x)}
@EndifMan
@IfDoc
@BeginMath
o(t,x) = \mbox{\bf std}\{i_1(t,x), i_2(t,x), \cdots, i_n(t,x)\}
o(t,x) = \mbox{\bf var1}\{i_1(t,x), i_2(t,x), \cdots, i_n(t,x)\}
@EndMath
@EndifDoc
@EndDescription
......
......@@ -42,6 +42,12 @@ In this section the abbreviations as in the following table are used:
n^{-1} \sum_{i=1}^{n} (x_i - \overline{x})^2
\\
\begin{array}{l}
\makebox[3cm][l]{{\bf var1}} \\
\end{array}
& &
(n-1)^{-1} \sum_{i=1}^{n} (x_i - \overline{x})^2
\\
\begin{array}{l}
\makebox[3cm][l]{{\bf var} weighted by} \\
\{w_i, i=1, ..., n\} \\
\end{array}
......@@ -57,6 +63,12 @@ In this section the abbreviations as in the following table are used:
\sqrt{ n^{-1} \sum_{i=1}^{n} (x_i - \overline{x})^2 }
\\
\begin{array}{l}
\makebox[3cm][l]{{\bf std1}} \\
\end{array}
& &
\sqrt{ (n-1)^{-1} \sum_{i=1}^{n} (x_i - \overline{x})^2 }
\\
\begin{array}{l}
\makebox[3cm][l]{{\bf std} weighted by} \\
\{w_i, i=1, ..., n\} \\
\end{array}
......
......@@ -24,7 +24,9 @@
Ensstat ensmean Ensemble mean
Ensstat ensavg Ensemble average
Ensstat ensstd Ensemble standard deviation
Ensstat ensstd1 Ensemble standard deviation
Ensstat ensvar Ensemble variance
Ensstat ensvar1 Ensemble variance
Ensstat enspctl Ensemble percentiles
Ensstat enscrps Ensemble cumulative ranked probability score
......@@ -70,9 +72,7 @@ void *Ensstat(void *argument)
double *array;
} ens_file_t;
ens_file_t *ef = NULL;
/* RQ */
int pn = 0;
/* QR */
cdoInitialize(argument);
......@@ -81,21 +81,17 @@ void *Ensstat(void *argument)
cdoOperatorAdd("enssum", func_sum, 0, NULL);
cdoOperatorAdd("ensmean", func_mean, 0, NULL);
cdoOperatorAdd("ensavg", func_avg, 0, NULL);
cdoOperatorAdd("ensvar", func_var, 0, NULL);
cdoOperatorAdd("ensstd", func_std, 0, NULL);
/* RQ */
cdoOperatorAdd("ensstd", func_std1, 0, NULL);
cdoOperatorAdd("ensvar", func_var, 0, NULL);
cdoOperatorAdd("ensvar1", func_var1, 0, NULL);
cdoOperatorAdd("enspctl", func_pctl, 0, NULL);
/* QR */
/* >>> Cedrick Ansorge 18.10.2010 */
cdoOperatorAdd("enscrps", func_crps, 0, NULL);
cdoOperatorAdd("ensbrs", func_brs, 0, NULL);
/* <<< Cedrick Ansorge 18.10.2010 */
operatorID = cdoOperatorID();
operfunc = cdoOperatorF1(operatorID);
/* RQ */
if ( operfunc == func_pctl )
{
operatorInputArg("percentile number");
......@@ -104,7 +100,6 @@ void *Ensstat(void *argument)
if ( pn < 1 || pn > 99 )
cdoAbort("Illegal argument: percentile number %d is not in the range 1..99!", pn);
}
/* QR */
nfiles = cdoStreamCnt() - 1;
......@@ -222,12 +217,10 @@ void *Ensstat(void *argument)
field[ompthID].nmiss++;
}
/* RQ */
if ( operfunc == func_pctl )
array2[i] = fldpctl(field[ompthID], pn);
else
array2[i] = fldfun(field[ompthID], operfunc);
/* QR */
if ( DBL_IS_EQUAL(array2[i], field[ompthID].missval) )
{
......
......@@ -67,7 +67,9 @@ double fldfun(field_t field, int function)
else if ( function == func_mean ) rval = fldmean(field);
else if ( function == func_avg ) rval = fldavg(field);
else if ( function == func_std ) rval = fldstd(field);
else if ( function == func_std1 ) rval = fldstd1(field);
else if ( function == func_var ) rval = fldvar(field);
else if ( function == func_var1 ) rval = fldvar1(field);
else if ( function == func_crps ) rval = fldcrps(field);
else if ( function == func_brs ) rval = fldbrs(field);
......@@ -326,39 +328,49 @@ double fldavg(field_t field)
return (ravg);
}
double fldvar(field_t field)
{
long i;
long len = field.size;
int nmiss = field.nmiss;
double missval = field.missval;
double *array = field.ptr;
double *w = field.weight;
double rsum = 0, rsumw = 0, rvar = 0;
double rsumq = 0, rsumwq = 0;
static
void prevarsum(const double *restrict array, const double *restrict w, long len, int nmiss,
double missval, double *rsum, double *rsumw, double *rsumq, double *rsumwq)
{
long i;
*rsum = *rsumw = 0;
*rsumq = *rsumwq = 0;
if ( nmiss > 0 )
{
for ( i = 0; i < len; i++ )
if ( !DBL_IS_EQUAL(array[i], missval) && !DBL_IS_EQUAL(w[i], missval) )
{
rsum += w[i] * array[i];
rsumq += w[i] * array[i] * array[i];
rsumw += w[i];
rsumwq += w[i] * w[i];
*rsum += w[i] * array[i];
*rsumq += w[i] * array[i] * array[i];
*rsumw += w[i];
*rsumwq += w[i] * w[i];
}
}
else
{
for ( i = 0; i < len; i++ )
{
rsum += w[i] * array[i];
rsumq += w[i] * array[i] * array[i];
rsumw += w[i];
rsumwq += w[i] * w[i];
*rsum += w[i] * array[i];
*rsumq += w[i] * array[i] * array[i];
*rsumw += w[i];
*rsumwq += w[i] * w[i];
}
}
}
double fldvar(field_t field)
{
long len = field.size;
int nmiss = field.nmiss;
double missval = field.missval;
double *array = field.ptr;
double *w = field.weight;
double rsum, rsumw, rvar;
double rsumq, rsumwq;
prevarsum(array, w, len, nmiss, missval, &rsum, &rsumw, &rsumq, &rsumwq);
rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq*rsumw - rsum*rsum) / (rsumw*rsumw) : missval;
if ( rvar < 0 && rvar > -1.e-5 ) rvar = 0;
......@@ -367,6 +379,25 @@ double fldvar(field_t field)
}
double fldvar1(field_t field)
{
long len = field.size;
int nmiss = field.nmiss;
double missval = field.missval;
double *array = field.ptr;
double *w = field.weight;
double rsum, rsumw, rvar;
double rsumq, rsumwq;
prevarsum(array, w, len, nmiss, missval, &rsum, &rsumw, &rsumq, &rsumwq);
rvar = (rsumw*rsumw > rsumwq) ? (rsumq*rsumw - rsum*rsum) / (rsumw*rsumw - rsumwq) : missval;
if ( rvar < 0 && rvar > -1.e-5 ) rvar = 0;
return (rvar);
}
double fldstd(field_t field)
{
double missval = field.missval;
......@@ -387,6 +418,26 @@ double fldstd(field_t field)
}
double fldstd1(field_t field)
{
double missval = field.missval;
double rvar, rstd;
rvar = fldvar1(field);
if ( DBL_IS_EQUAL(rvar, missval) || rvar < 0 )
{
rstd = missval;
}
else
{
rstd = IS_NOT_EQUAL(rvar, 0) ? sqrt(rvar) : 0;
}
return (rstd);
}
void fldrms(field_t field, field_t field2, field_t *field3)
{
long i, len;
......
......@@ -75,7 +75,9 @@ double fldsum(field_t field);
double fldavg(field_t field);
double fldmean(field_t field);
double fldstd(field_t field);
double fldstd1(field_t field);
double fldvar(field_t field);
double fldvar1(field_t field);
/* RQ */
double fldpctl(field_t field, int k);
/* QR */
......
......@@ -11,12 +11,15 @@
#define func_sum 13
#define func_avg 14
#define func_mean 15
#define func_var 16
#define func_std 17
#define func_pctl 18
#define func_std1 18
#define func_var 19
#define func_var1 20
#define func_cor 20
#define func_covar 21
#define func_pctl 21
#define func_cor 22
#define func_covar 23
#define func_crps 30
#define func_brs 31
......
......@@ -299,7 +299,7 @@ void *Maggraph(void *argument);
"export_e5ml", "export_e5res"}
#define EnlargeOperators {"enlarge"}
#define EnlargegridOperators {"enlargegrid"}
#define EnsstatOperators {"ensmin", "ensmax", "enssum", "ensmean", "ensavg", "ensvar", "ensstd", "enspctl"}
#define EnsstatOperators {"ensmin", "ensmax", "enssum", "ensmean", "ensavg", "ensvar", "ensvar1", "ensstd", "ensstd1", "enspctl"}
#define Ensstat3Operators {"ensrkhist_space","ensrkhist_time","ensroc"}
#define EnsvalOperators {"enscrps","ensbrs"}
#define EofcoeffOperators {"eofcoeff"}
......
......@@ -1448,8 +1448,8 @@ static char *ConsecstatHelp[] = {
static char *EnsstatHelp[] = {
"NAME",
" ensmin, ensmax, enssum, ensmean, ensavg, ensvar, ensstd, enspctl - ",
" Statistical values over an ensemble",
" ensmin, ensmax, enssum, ensmean, ensavg, ensstd, ensstd1, ensvar, ensvar1, ",
" enspctl - Statistical values over an ensemble",
"",
"SYNOPSIS",
" ensmin ifiles ofile",
......@@ -1457,8 +1457,10 @@ static char *EnsstatHelp[] = {
" enssum ifiles ofile",
" ensmean ifiles ofile",
" ensavg ifiles ofile",
" ensvar ifiles ofile",
" ensstd ifiles ofile",
" ensstd1 ifiles ofile",
" ensvar ifiles ofile",
" ensvar1 ifiles ofile",
" enspctl,p ifiles ofile",
"",
"DESCRIPTION",
......@@ -1480,10 +1482,22 @@ static char *EnsstatHelp[] = {
" o(t,x) = mean{i1(t,x), i2(t,x), ..., in(t,x)}",
" ensavg Ensemble average",
" o(t,x) = avg{i1(t,x), i2(t,x), ..., in(t,x)}",
" ensvar Ensemble variance",
" o(t,x) = var{i1(t,x), i2(t,x), ..., in(t,x)}",
" ensstd Ensemble standard deviation",
" Divisor is n.",
" ",
" o(t,x) = std{i1(t,x), i2(t,x), ..., in(t,x)}",
" ensstd1 Ensemble standard deviation",
" Divisor is (n-1).",
" ",
" o(t,x) = std1{i1(t,x), i2(t,x), ..., in(t,x)}",
" ensvar Ensemble variance",
" Divisor is n.",
" ",
" o(t,x) = var{i1(t,x), i2(t,x), ..., in(t,x)}",
" ensvar1 Ensemble variance",
" Divisor is (n-1).",
" ",
" o(t,x) = var1{i1(t,x), i2(t,x), ..., in(t,x)}",
" enspctl Ensemble percentiles",
" o(t,x) = pth percentile {i1(t,x), i2(t,x), ..., in(t,x)}",
"",
......
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