Commit cc5b137e authored by Uwe Schulzweida's avatar Uwe Schulzweida

Added operator timrms (status: experimental)

parent 5be622f9
Pipeline #5033 passed with stages
in 17 minutes and 43 seconds
......@@ -29,8 +29,8 @@
// correlation in time
static void
correlationInit(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals,
double *work0, double *work1, double *work2, double *work3, double *work4)
correlation_init(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals,
double *work0, double *work1, double *work2, double *work3, double *work4)
{
for (size_t i = 0; i < gridsize; ++i)
{
......@@ -85,8 +85,8 @@ correlation(size_t gridsize, double missval1, double missval2, size_t *nofvals,
// covariance in time
static void
covarianceInit(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals,
double *work0, double *work1, double *work2)
covariance_init(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals,
double *work0, double *work1, double *work2)
{
for (size_t i = 0; i < gridsize; ++i)
{
......@@ -128,6 +128,47 @@ covariance(size_t gridsize, double missval1, double missval2, size_t *nofvals, d
return nmiss;
}
// rms in time
static void
rms_init(size_t gridsize, const double *array1, const double *array2, double missval1, double missval2, size_t *nofvals,
double *work0)
{
for (size_t i = 0; i < gridsize; ++i)
{
if ((!DBL_IS_EQUAL(array1[i], missval1)) && (!DBL_IS_EQUAL(array2[i], missval2)))
{
work0[i] += (array2[i] - array1[i]) * (array2[i] - array1[i]);
nofvals[i]++;
}
}
}
static size_t
rms(size_t gridsize, double missval, size_t *nofvals, double *work0)
{
size_t nmiss = 0;
for (size_t i = 0; i < gridsize; ++i)
{
double rms;
const auto nvals = nofvals[i];
if (nvals > 0)
{
rms = std::sqrt(work0[i]);
}
else
{
nmiss++;
rms = missval;
}
work0[i] = rms;
}
return nmiss;
}
void *
Timstat2(void *process)
{
......@@ -147,6 +188,7 @@ Timstat2(void *process)
const auto operatorID = cdoOperatorID();
const auto operfunc = cdoOperatorF1(operatorID);
const auto nwork = cdoOperatorF2(operatorID);
const auto timeIsConst = (operfunc == func_rms);
operatorCheckArgc(0);
......@@ -172,6 +214,10 @@ Timstat2(void *process)
// const auto taxisID2 = vlistInqTaxis(vlistID2);
const auto taxisID3 = taxisDuplicate(taxisID1);
if (timeIsConst)
for (varID = 0; varID < nvars; ++varID)
vlistDefVarTimetype(vlistID3, varID, TIME_CONSTANT);
vlistDefTaxis(vlistID3, taxisID3);
const auto streamID3 = cdoOpenWrite(2);
cdoDefVlist(streamID3, vlistID3);
......@@ -230,13 +276,17 @@ Timstat2(void *process)
if (operfunc == func_cor)
{
correlationInit(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(),
rwork[0].data(), rwork[1].data(), rwork[2].data(), rwork[3].data(), rwork[4].data());
correlation_init(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(),
rwork[0].data(), rwork[1].data(), rwork[2].data(), rwork[3].data(), rwork[4].data());
}
else if (operfunc == func_covar)
{
covarianceInit(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(),
rwork[0].data(), rwork[1].data(), rwork[2].data());
covariance_init(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(),
rwork[0].data(), rwork[1].data(), rwork[2].data());
}
else if (operfunc == func_rms)
{
rms_init(gridsize, array1.data(), array2.data(), missval1, missval2, rnofvals.data(), rwork[0].data());
}
}
......@@ -270,6 +320,10 @@ Timstat2(void *process)
nmiss = covariance(gridsize, missval1, missval2, rnofvals.data(), rwork[0].data(),
rwork[1].data(), rwork[2].data());
}
else if (operfunc == func_rms)
{
nmiss = rms(gridsize, missval1, rnofvals.data(), rwork[0].data());
}
cdoDefRecord(streamID3, varID, levelID);
cdoWriteRecord(streamID3, rwork[0].data(), nmiss);
......
......@@ -1598,7 +1598,8 @@ void init_modules()
{"Daystat" ,module_Daystat },
{"Hourstat" ,module_Hourstat },
{"Timcor" ,module_Timcor },
{"Timscorvar" ,module_Timscorvar },
{"Timcorvar" ,module_Timcorvar },
{"Timrms" ,module_Timrms },
{"Timstat3" ,module_Timstat3 },
{"Tinfo" ,module_Tinfo },
{"Tocomplex" ,module_Tocomplex },
......
......@@ -206,7 +206,7 @@ static const module_t module_Monstat = {Timstat , MonstatHelp
static const module_t module_Daystat = {Timstat , DaystatHelp , DaystatOperators , EXPOSED , CDI_BOTH , 1 , 1 , NoRestriction };
static const module_t module_Hourstat = {Timstat , HourstatHelp , HourstatOperators , EXPOSED , CDI_BOTH , 1 , 1 , NoRestriction };
static const module_t module_Timcor = {Timstat2 , TimcorHelp , TimcorOperators , EXPOSED , CDI_REAL , 2 , 1 , NoRestriction };
static const module_t module_Timscorvar = {Timstat2 , TimcovarHelp , TimcovarOperators , EXPOSED , CDI_REAL , 2 , 1 , NoRestriction };
static const module_t module_Timcorvar = {Timstat2 , TimcovarHelp , TimcovarOperators , EXPOSED , CDI_REAL , 2 , 1 , NoRestriction };
static const module_t module_Timrms = {Timstat2 , {} , TimrmsOperators , EXPOSED , CDI_REAL , 2 , 1 , NoRestriction };
static const module_t module_Timstat3 = {Timstat3 , {} , Timstat3Operators , EXPOSED , CDI_REAL , 2 , 1 , NoRestriction };
static const module_t module_Tinfo = {Tinfo , {} , TinfoOperators , EXPOSED , CDI_BOTH , 1 , 0 , NoRestriction };
......
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