Commit e6bfc6de authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

fieldzon: added temporary vector.

parent 66ecbe7c
......@@ -81,35 +81,32 @@ prevarsum(const Field &field, double &rsum, double &rsumw, double &rsumq, double
assert(array.size() >= len);
double xsum = 0, xsumw = 0;
double xsumq = 0, xsumwq = 0;
rsum = 0;
rsumq = 0;
rsumw = 0;
rsumwq = 0;
if (field.nmiss)
{
for (size_t i = 0; i < len; ++i)
if (!DBL_IS_EQUAL(array[i], missval))
{
xsum += array[i];
xsumq += array[i] * array[i];
xsumw += 1;
xsumwq += 1;
rsum += array[i];
rsumq += array[i] * array[i];
rsumw += 1;
rsumwq += 1;
}
}
else
{
for (size_t i = 0; i < len; ++i)
{
xsum += array[i];
xsumq += array[i] * array[i];
rsum += array[i];
rsumq += array[i] * array[i];
}
xsumw = len;
xsumwq = len;
rsumw = len;
rsumwq = len;
}
rsum = xsum;
rsumq = xsumq;
rsumw = xsumw;
rsumwq = xsumwq;
}
static void
......@@ -259,41 +256,38 @@ prevarsumw(const Field &field, double &rsum, double &rsumw, double &rsumq, doubl
{
const size_t len = field.gridsize;
const double missval = field.missval;
const auto &array = field.vec;
const auto &v = field.vec;
const auto &w = field.weightv;
assert(array.size() >= len);
assert(v.size() >= len);
assert(w.size() >= len);
double xsum = 0, xsumw = 0;
double xsumq = 0, xsumwq = 0;
rsum = 0;
rsumq = 0;
rsumw = 0;
rsumwq = 0;
if (field.nmiss)
{
for (size_t i = 0; i < len; ++i)
if (!DBL_IS_EQUAL(array[i], missval) && !DBL_IS_EQUAL(w[i], missval))
if (!DBL_IS_EQUAL(v[i], missval) && !DBL_IS_EQUAL(w[i], missval))
{
xsum += w[i] * array[i];
xsumq += w[i] * array[i] * array[i];
xsumw += w[i];
xsumwq += w[i] * w[i];
rsum += w[i] * v[i];
rsumq += w[i] * v[i] * v[i];
rsumw += w[i];
rsumwq += w[i] * w[i];
}
}
else
{
for (size_t i = 0; i < len; ++i)
{
xsum += w[i] * array[i];
xsumq += w[i] * array[i] * array[i];
xsumw += w[i];
xsumwq += w[i] * w[i];
rsum += w[i] * v[i];
rsumq += w[i] * v[i] * v[i];
rsumw += w[i];
rsumwq += w[i] * w[i];
}
}
rsum = xsum;
rsumq = xsumq;
rsumw = xsumw;
rsumwq = xsumwq;
}
double
......
......@@ -194,7 +194,7 @@ meravgw(const Field &field1, Field &field2)
}
static void
prevarsum_merw(const double *restrict array, const double *restrict w, size_t nx, size_t ny, size_t nmiss, double missval,
prevarsum_merw(const std::vector<double> &v, const std::vector<double> &w, size_t ny, size_t nmiss, double missval,
double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
rsum = 0;
......@@ -205,22 +205,22 @@ prevarsum_merw(const double *restrict array, const double *restrict w, size_t nx
if (nmiss)
{
for (size_t j = 0; j < ny; j++)
if (!DBL_IS_EQUAL(array[j * nx], missval) && !DBL_IS_EQUAL(w[j * nx], missval))
if (!DBL_IS_EQUAL(v[j], missval) && !DBL_IS_EQUAL(w[j], missval))
{
rsum += w[j * nx] * array[j * nx];
rsumq += w[j * nx] * array[j * nx] * array[j * nx];
rsumw += w[j * nx];
rsumwq += w[j * nx] * w[j * nx];
rsum += w[j] * v[j];
rsumq += w[j] * v[j] * v[j];
rsumw += w[j];
rsumwq += w[j] * w[j];
}
}
else
{
for (size_t j = 0; j < ny; j++)
{
rsum += w[j * nx] * array[j * nx];
rsumq += w[j * nx] * array[j * nx] * array[j * nx];
rsumw += w[j * nx];
rsumwq += w[j * nx] * w[j * nx];
rsum += w[j] * v[j];
rsumq += w[j] * v[j] * v[j];
rsumw += w[j];
rsumwq += w[j] * w[j];
}
}
}
......@@ -232,17 +232,20 @@ mervarw(const Field &field1, Field &field2)
int grid = field1.grid;
size_t nmiss = field1.nmiss;
double missval = field1.missval;
const double *array = field1.vec.data();
const double *w = field1.weightv.data();
double rsum = 0, rsumw = 0;
double rsumq = 0, rsumwq = 0;
const size_t nx = gridInqXsize(grid);
const size_t ny = gridInqYsize(grid);
std::vector<double> v(ny), w(ny);
for (size_t i = 0; i < nx; i++)
{
prevarsum_merw(array + i, w + i, nx, ny, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
for (size_t j = 0; j < ny; j++) v[j] = field1.vec[j * nx + i];
for (size_t j = 0; j < ny; j++) w[j] = field1.weightv[j * nx + i];
prevarsum_merw(v, w, ny, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : missval;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
......@@ -260,17 +263,20 @@ mervar1w(const Field &field1, Field &field2)
int grid = field1.grid;
size_t nmiss = field1.nmiss;
double missval = field1.missval;
const double *array = field1.vec.data();
const double *w = field1.weightv.data();
double rsum = 0, rsumw = 0;
double rsumq = 0, rsumwq = 0;
const size_t nx = gridInqXsize(grid);
const size_t ny = gridInqYsize(grid);
std::vector<double> v(ny), w(ny);
for (size_t i = 0; i < nx; i++)
{
prevarsum_merw(array + i, w + i, nx, ny, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
for (size_t j = 0; j < ny; j++) v[j] = field1.vec[j * nx + i];
for (size_t j = 0; j < ny; j++) w[j] = field1.weightv[j * nx + i];
prevarsum_merw(v, w, ny, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
......
......@@ -19,26 +19,6 @@
#include "cdo_int.h"
#include "percentiles.h"
void
zonfun(const Field &field1, Field &field2, int function)
{
// clang-format off
switch (function)
{
case func_min: zonmin(field1, field2); break;
case func_max: zonmax(field1, field2); break;
case func_range: zonrange(field1, field2); break;
case func_sum: zonsum(field1, field2); break;
case func_mean: zonmean(field1, field2); break;
case func_avg: zonavg(field1, field2); break;
case func_std: zonstd(field1, field2); break;
case func_std1: zonstd1(field1, field2); break;
case func_var: zonvar(field1, field2); break;
case func_var1: zonvar1(field1, field2); break;
default: cdoAbort("function %d not implemented!", function);
}
// clang-format on
}
void
zonmin(const Field &field1, Field &field2)
......@@ -49,16 +29,20 @@ zonmin(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (field1.nmiss)
{
rmin = arrayMinMV(nx, &field1.vec[j * nx], field1.missval);
rmin = arrayMinMV(nx, v.data(), field1.missval);
if (DBL_IS_EQUAL(rmin, field1.missval)) rnmiss++;
}
else
{
rmin = arrayMin(nx, &field1.vec[j * nx]);
rmin = arrayMin(nx, v.data());
}
field2.vec[j] = rmin;
......@@ -76,16 +60,20 @@ zonmax(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (field1.nmiss)
{
rmax = arrayMaxMV(nx, &field1.vec[j * nx], field1.missval);
rmax = arrayMaxMV(nx, v.data(), field1.missval);
if (DBL_IS_EQUAL(rmax, field1.missval)) rnmiss++;
}
else
{
rmax = arrayMax(nx, &field1.vec[j * nx]);
rmax = arrayMax(nx, v.data());
}
field2.vec[j] = rmax;
......@@ -103,16 +91,20 @@ zonrange(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (field1.nmiss)
{
range = arrayRangeMV(nx, &field1.vec[j * nx], field1.missval);
range = arrayRangeMV(nx, v.data(), field1.missval);
if (DBL_IS_EQUAL(range, field1.missval)) rnmiss++;
}
else
{
range = arrayRange(nx, &field1.vec[j * nx]);
range = arrayRange(nx, v.data());
}
field2.vec[j] = range;
......@@ -130,16 +122,20 @@ zonsum(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (field1.nmiss)
{
rsum = arraySumMV(nx, &field1.vec[j * nx], field1.missval);
rsum = arraySumMV(nx, v.data(), field1.missval);
if (DBL_IS_EQUAL(rsum, field1.missval)) rnmiss++;
}
else
{
rsum = arraySum(nx, &field1.vec[j * nx]);
rsum = arraySum(nx, v.data());
}
field2.vec[j] = rsum;
......@@ -157,15 +153,19 @@ zonmean(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (field1.nmiss)
{
rmean = arrayMeanMV(nx, &field1.vec[j * nx], field1.missval);
rmean = arrayMeanMV(nx, v.data(), field1.missval);
}
else
{
rmean = arrayMean(nx, &field1.vec[j * nx]);
rmean = arrayMean(nx, v.data());
}
if (DBL_IS_EQUAL(rmean, field1.missval)) rnmiss++;
......@@ -185,15 +185,19 @@ zonavg(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(field1.grid);
const size_t ny = gridInqYsize(field1.grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; ++j)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (field1.nmiss)
{
ravg = arrayAvgMV(nx, &field1.vec[j * nx], field1.missval);
ravg = arrayAvgMV(nx, v.data(), field1.missval);
}
else
{
ravg = arrayMean(nx, &field1.vec[j * nx]);
ravg = arrayMean(nx, v.data());
}
if (DBL_IS_EQUAL(ravg, field1.missval)) rnmiss++;
......@@ -205,7 +209,7 @@ zonavg(const Field &field1, Field &field2)
}
static void
prevarsum_zon(const double *restrict array, size_t nx, size_t nmiss, double missval, double &rsum, double &rsumw, double &rsumq,
prevarsum_zon(const std::vector<double> v, size_t nx, size_t nmiss, double missval, double &rsum, double &rsumw, double &rsumq,
double &rsumwq)
{
const double w = 1. / nx;
......@@ -218,10 +222,10 @@ prevarsum_zon(const double *restrict array, size_t nx, size_t nmiss, double miss
if (nmiss)
{
for (size_t i = 0; i < nx; i++)
if (!DBL_IS_EQUAL(array[i], missval))
if (!DBL_IS_EQUAL(v[i], missval))
{
rsum += w * array[i];
rsumq += w * array[i] * array[i];
rsum += w * v[i];
rsumq += w * v[i] * v[i];
rsumw += w;
rsumwq += w * w;
}
......@@ -230,8 +234,8 @@ prevarsum_zon(const double *restrict array, size_t nx, size_t nmiss, double miss
{
for (size_t i = 0; i < nx; i++)
{
rsum += w * array[i];
rsumq += w * array[i] * array[i];
rsum += w * v[i];
rsumq += w * v[i] * v[i];
rsumw += w;
rsumwq += w * w;
}
......@@ -251,9 +255,13 @@ zonvar(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(grid);
const size_t ny = gridInqYsize(grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; j++)
{
prevarsum_zon(&field1.vec[j * nx], nx, nmiss, missval1, rsum, rsumw, rsumq, rsumwq);
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
prevarsum_zon(v, nx, nmiss, missval1, rsum, rsumw, rsumq, rsumwq);
double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : missval1;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
......@@ -277,9 +285,13 @@ zonvar1(const Field &field1, Field &field2)
const size_t nx = gridInqXsize(grid);
const size_t ny = gridInqYsize(grid);
std::vector<double> v(nx);
for (size_t j = 0; j < ny; j++)
{
prevarsum_zon(&field1.vec[j * nx], nx, nmiss, missval1, rsum, rsumw, rsumq, rsumwq);
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
prevarsum_zon(v, nx, nmiss, missval1, rsum, rsumw, rsumq, rsumwq);
double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval1;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
......@@ -342,15 +354,19 @@ zonpctl(Field &field1, Field &field2, int p)
const size_t nx = gridInqXsize(grid);
const size_t ny = gridInqYsize(grid);
std::vector<double> v(nx);
if (field1.nmiss)
{
std::vector<double> array2(nx);
for (size_t j = 0; j < ny; j++)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
size_t l = 0;
for (size_t i = 0; i < nx; i++)
if (!DBL_IS_EQUAL(field1.vec[j * nx + i], missval)) array2[l++] = field1.vec[j * nx + i];
if (!DBL_IS_EQUAL(field1.vec[j * nx + i], missval)) array2[l++] = v[i];
if (l > 0)
{
......@@ -367,9 +383,11 @@ zonpctl(Field &field1, Field &field2, int p)
{
for (size_t j = 0; j < ny; j++)
{
for (size_t i = 0; i < nx; i++) v[i] = field1.vec[j * nx + i];
if (nx > 0)
{
field2.vec[j] = percentile(&field1.vec[j * nx], nx, p);
field2.vec[j] = percentile(v.data(), nx, p);
}
else
{
......@@ -381,3 +399,24 @@ zonpctl(Field &field1, Field &field2, int p)
field2.nmiss = rnmiss;
}
void
zonfun(const Field &field1, Field &field2, int function)
{
// clang-format off
switch (function)
{
case func_min: zonmin(field1, field2); break;
case func_max: zonmax(field1, field2); break;
case func_range: zonrange(field1, field2); break;
case func_sum: zonsum(field1, field2); break;
case func_mean: zonmean(field1, field2); break;
case func_avg: zonavg(field1, field2); break;
case func_std: zonstd(field1, field2); break;
case func_std1: zonstd1(field1, field2); break;
case func_var: zonvar(field1, field2); break;
case func_var1: zonvar1(field1, field2); break;
default: cdoAbort("function %d not implemented!", function);
}
// clang-format on
}
Markdown is supported
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