Commit 705382fb authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Field refactoring.

parent e2bdbbbb
......@@ -73,15 +73,18 @@ vfldavgw(const Field &field)
}
static void
prevarsum(const double *restrict array, size_t len, size_t nmiss, double missval, double &rsum, double &rsumw, double &rsumq,
double &rsumwq)
prevarsum(const Field &field, double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
assert(array != nullptr);
const size_t len = field.gridsize;
const double missval = field.missval;
const auto &array = field.vec;
assert(array.size() >= len);
double xsum = 0, xsumw = 0;
double xsumq = 0, xsumwq = 0;
if (nmiss)
if (field.nmiss)
{
for (size_t i = 0; i < len; ++i)
if (!DBL_IS_EQUAL(array[i], missval))
......@@ -110,15 +113,18 @@ prevarsum(const double *restrict array, size_t len, size_t nmiss, double missval
}
static void
preskewsum(const double *restrict array, size_t len, size_t nmiss, const double mean, double missval, double *rsum3w,
double *rsum4w, double *rsum3diff, double *rsum2diff)
preskewsum(const Field &field, const double mean, double &rsum3w, double &rsum4w, double &rsum3diff, double &rsum2diff)
{
assert(array != nullptr);
const size_t len = field.gridsize;
const double missval = field.missval;
const auto &array = field.vec;
assert(array.size() >= len);
double xsum3w = 0, xsum3diff = 0;
double xsum4w = 0, xsum2diff = 0;
if (nmiss)
if (field.nmiss)
{
for (size_t i = 0; i < len; ++i)
if (!DBL_IS_EQUAL(array[i], missval))
......@@ -140,22 +146,25 @@ preskewsum(const double *restrict array, size_t len, size_t nmiss, const double
xsum4w = len;
}
*rsum3diff = xsum3diff;
*rsum2diff = xsum2diff;
*rsum3w = xsum3w;
*rsum4w = xsum4w;
rsum3diff = xsum3diff;
rsum2diff = xsum2diff;
rsum3w = xsum3w;
rsum4w = xsum4w;
}
static void
prekurtsum(const double *restrict array, size_t len, size_t nmiss, const double mean, double missval, double *rsum3w,
double *rsum4w, double *rsum2diff, double *rsum4diff)
prekurtsum(const Field &field, const double mean, double &rsum3w, double &rsum4w, double &rsum2diff, double &rsum4diff)
{
assert(array != nullptr);
const size_t len = field.gridsize;
const double missval = field.missval;
const auto &array = field.vec;
assert(array.size() >= len);
double xsum3w = 0, xsum4diff = 0;
double xsum4w = 0, xsum2diff = 0;
if (nmiss)
if (field.nmiss)
{
for (size_t i = 0; i < len; ++i)
if (!DBL_IS_EQUAL(array[i], missval))
......@@ -177,55 +186,19 @@ prekurtsum(const double *restrict array, size_t len, size_t nmiss, const double
xsum4w = len;
}
*rsum4diff = xsum4diff;
*rsum2diff = xsum2diff;
*rsum3w = xsum3w;
*rsum4w = xsum4w;
}
double
fldvar(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum, rsumw, rsumq, rsumwq;
prevarsum(field.ptr, len, 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;
return rvar;
rsum4diff = xsum4diff;
rsum2diff = xsum2diff;
rsum3w = xsum3w;
rsum4w = xsum4w;
}
double
vfldvar(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum, rsumw, rsumq, rsumwq;
prevarsum(field.vec.data(), len, 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;
return rvar;
}
double
fldvar1(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum, rsumw, rsumq, rsumwq;
prevarsum(field.ptr, len, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
prevarsum(field, rsum, rsumw, rsumq, rsumwq);
double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : field.missval;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
return rvar;
......@@ -234,14 +207,10 @@ fldvar1(const Field &field)
double
vfldvar1(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum, rsumw, rsumq, rsumwq;
prevarsum(field.vec.data(), len, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
prevarsum(field, rsum, rsumw, rsumq, rsumwq);
double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : field.missval;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
return rvar;
......@@ -250,18 +219,15 @@ vfldvar1(const Field &field)
double
vfldkurt(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum3w; /* 3rd moment variables */
double rsum4w; /* 4th moment variables */
double rsum2diff, rsum4diff;
double rsum, rsumw, rsumq, rsumwq;
prevarsum(field.vec.data(), len, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
prekurtsum(field.vec.data(), len, nmiss, (rsum / rsumw), missval, &rsum3w, &rsum4w, &rsum2diff, &rsum4diff);
prevarsum(field, rsum, rsumw, rsumq, rsumwq);
prekurtsum(field, (rsum / rsumw), rsum3w, rsum4w, rsum2diff, rsum4diff);
if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum2diff, 0.0)) return missval;
if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum2diff, 0.0)) return field.missval;
double rvar = ((rsum4diff / rsum3w) / pow((rsum2diff) / (rsum3w), 2)) - 3.0;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
......@@ -272,18 +238,15 @@ vfldkurt(const Field &field)
double
vfldskew(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum3w; /* 3rd moment variables */
double rsum4w; /* 4th moment variables */
double rsum3diff, rsum2diff;
double rsum, rsumw, rsumq, rsumwq;
prevarsum(field.vec.data(), len, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
preskewsum(field.vec.data(), len, nmiss, (rsum / rsumw), missval, &rsum3w, &rsum4w, &rsum3diff, &rsum2diff);
prevarsum(field, rsum, rsumw, rsumq, rsumwq);
preskewsum(field, (rsum / rsumw), rsum3w, rsum4w, rsum3diff, rsum2diff);
if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum3w, 1.0) || IS_EQUAL(rsum2diff, 0.0)) return missval;
if (IS_EQUAL(rsum3w, 0.0) || IS_EQUAL(rsum3w, 1.0) || IS_EQUAL(rsum2diff, 0.0)) return field.missval;
double rvar = (rsum3diff / rsum3w) / pow((rsum2diff) / (rsum3w - 1.0), 1.5);
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
......@@ -292,16 +255,20 @@ vfldskew(const Field &field)
}
static void
prevarsumw(const double *restrict array, const double *restrict w, size_t len, size_t nmiss, double missval, double &rsum,
double &rsumw, double &rsumq, double &rsumwq)
prevarsumw(const Field &field, double &rsum, double &rsumw, double &rsumq, double &rsumwq)
{
assert(array != nullptr);
assert(w != nullptr);
const size_t len = field.gridsize;
const double missval = field.missval;
const auto &array = field.vec;
const auto &w = field.weightv;
assert(array.size() >= len);
assert(w.size() >= len);
double xsum = 0, xsumw = 0;
double xsumq = 0, xsumwq = 0;
if (nmiss)
if (field.nmiss)
{
for (size_t i = 0; i < len; ++i)
if (!DBL_IS_EQUAL(array[i], missval) && !DBL_IS_EQUAL(w[i], missval))
......@@ -332,14 +299,10 @@ prevarsumw(const double *restrict array, const double *restrict w, size_t len, s
double
vfldvarw(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum, rsumw, rsumq, rsumwq;
prevarsumw(field.vec.data(), field.weightv.data(), len, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
prevarsumw(field, rsum, rsumw, rsumq, rsumwq);
double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : missval;
double rvar = IS_NOT_EQUAL(rsumw, 0) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw) : field.missval;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
return rvar;
......@@ -348,14 +311,10 @@ vfldvarw(const Field &field)
double
vfldvar1w(const Field &field)
{
const size_t nmiss = field.nmiss > 0;
const size_t len = field.gridsize;
const double missval = field.missval;
double rsum, rsumw, rsumq, rsumwq;
prevarsumw(field.vec.data(), field.weightv.data(), len, nmiss, missval, rsum, rsumw, rsumq, rsumwq);
prevarsumw(field, rsum, rsumw, rsumq, rsumwq);
double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : missval;
double rvar = (rsumw * rsumw > rsumwq) ? (rsumq * rsumw - rsum * rsum) / (rsumw * rsumw - rsumwq) : field.missval;
if (rvar < 0 && rvar > -1.e-5) rvar = 0;
return rvar;
......@@ -384,12 +343,6 @@ vfldstd(const Field &field)
return varToStd(vfldvar(field), field.missval);
}
double
fldstd1(const Field &field)
{
return varToStd(fldvar1(field), field.missval);
}
double
vfldstd1(const Field &field)
{
......@@ -421,7 +374,7 @@ vfldrms(const Field &field, const Field &field2, Field &field3)
const double *array2 = field2.vec.data();
const double missval1 = field.missval;
const double missval2 = field2.missval;
const std::vector<double> &w = field.weightv;
const auto &w = field.weightv;
double rsum = 0, rsumw = 0;
const size_t len = gridInqSize(grid1);
......
......@@ -179,7 +179,6 @@ void vfarmod(Field &field, double divisor);
void vfarcfuncplx(Field &field, const double rconstcplx[2], int function);
// field2.cc
void farfun(Field &field1, const Field &field2, int function);
void vfarfun(Field &field1, const Field &field2, int function);
void faradd(Field &field1, const Field &field2);
......@@ -201,28 +200,19 @@ void vfaratan2(Field &field1, const Field &field2);
void farsetmiss(Field &field1, const Field &field2);
void vfarsetmiss(Field &field1, const Field &field2);
void farsumq(Field &field1, const Field &field2);
void vfarsumq(Field &field1, const Field &field2);
void farsumw(Field &field1, const Field &field2, double w);
void vfarsumw(Field &field1, const Field &field2, double w);
void farsumqw(Field &field1, const Field &field2, double w);
void vfarsumqw(Field &field1, const Field &field2, double w);
void vfarsumtr(Field &field1, const Field &field2, double refval);
void vfarcpy(Field &field1, const Field &field2);
void vfarminidx(Field &field1, Field &field2, const Field &field3, int idx);
void vfarmaxidx(Field &field1, Field &field2, const Field &field3, int idx);
void farvar(Field &field1, const Field &field2, const Field &field3, int divisor);
void vfarvar(Field &field1, const Field &field2, const Field &field3, int divisor);
void farstd(Field &field1, const Field &field2, const Field &field3, int divisor);
void vfarstd(Field &field1, const Field &field2, const Field &field3, int divisor);
void farcvar(Field &field1, const Field &field2, int nsets, int divisor);
void vfarcvar(Field &field1, const Field &field2, int nsets, int divisor);
void farcstd(Field &field1, const Field &field2, int nsets, int divisor);
void vfarcstd(Field &field1, const Field &field2, int nsets, int divisor);
void farmoq(Field &field1, const Field &field2);
void vfarmoq(Field &field1, const Field &field2);
void farmoqw(Field &field1, const Field &field2, double w);
void vfarmoqw(Field &field1, const Field &field2, double w);
void vfarcount(Field &field1, const Field &field2);
......
This diff is collapsed.
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