Commit 13d51674 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Set constant variables to const.

parent 4e5016ee
......@@ -42,11 +42,7 @@ Timstat3(void *process)
int vlistID[NIN], vlistID2 = -1;
int64_t vdate = 0;
int vtime = 0;
int nlevs;
int is;
int gridID;
double missval, missval1, missval2;
double fractil_1, fractil_2, statistic;
Varray3D<int> iwork[NIWORK];
FieldVector2D fwork[NFWORK];
int reached_eof[NIN];
......@@ -86,7 +82,7 @@ Timstat3(void *process)
const auto vlistID3 = vlistDuplicate(vlistID[0]);
auto gridsize = vlistGridsizeMax(vlistID[0]);
const auto gridsizemax = vlistGridsizeMax(vlistID[0]);
const auto nvars = vlistNvars(vlistID[0]);
const auto maxrecs = vlistNrecs(vlistID[0]);
......@@ -102,23 +98,23 @@ Timstat3(void *process)
for (int i = 0; i < NIN; ++i) reached_eof[i] = 0;
Field in[NIN], out[NOUT];
for (int i = 0; i < NIN; ++i) in[i].resize(gridsize);
for (int i = 0; i < NOUT; ++i) out[i].resize(gridsize);
for (int i = 0; i < NIN; ++i) in[i].resize(gridsizemax);
for (int i = 0; i < NOUT; ++i) out[i].resize(gridsizemax);
for (int iw = 0; iw < NFWORK; ++iw) fwork[iw].resize(nvars);
for (int iw = 0; iw < NIWORK; ++iw) iwork[iw].resize(nvars);
for (int varID = 0; varID < nvars; ++varID)
{
gridID = vlistInqVarGrid(vlistID[0], varID);
gridsize = vlistGridsizeMax(vlistID[0]);
nlevs = zaxisInqSize(vlistInqVarZaxis(vlistID[0], varID));
missval = missval1 = vlistInqVarMissval(vlistID[0], varID);
const auto gridID = vlistInqVarGrid(vlistID[0], varID);
const auto gridsize = vlistGridsizeMax(vlistID[0]);
const auto nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID[0], varID));
const auto missval = vlistInqVarMissval(vlistID[0], varID);
for (int iw = 0; iw < NFWORK; ++iw) fwork[iw][varID].resize(nlevs);
for (int iw = 0; iw < NIWORK; ++iw) iwork[iw][varID].resize(nlevs);
for (int iw = 0; iw < NFWORK; ++iw) fwork[iw][varID].resize(nlevels);
for (int iw = 0; iw < NIWORK; ++iw) iwork[iw][varID].resize(nlevels);
for (int levelID = 0; levelID < nlevs; ++levelID)
for (int levelID = 0; levelID < nlevels; ++levelID)
{
for (int iw = 0; iw < NFWORK; ++iw)
{
......@@ -127,10 +123,7 @@ Timstat3(void *process)
fwork[iw][varID][levelID].resize(gridsize);
}
for (int iw = 0; iw < NIWORK; ++iw)
{
iwork[iw][varID][levelID].resize(gridsize, 0);
}
for (int iw = 0; iw < NIWORK; ++iw) iwork[iw][varID][levelID].resize(gridsize, 0);
}
}
......@@ -141,7 +134,7 @@ Timstat3(void *process)
{
if (reached_eof[is]) continue;
auto nrecs = cdoStreamInqTimestep(streamID[is], tsID);
const auto nrecs = cdoStreamInqTimestep(streamID[is], tsID);
if (nrecs == 0)
{
reached_eof[is] = 1;
......@@ -156,7 +149,7 @@ Timstat3(void *process)
int varID, levelID;
cdoInqRecord(streamID[is], &varID, &levelID);
gridsize = gridInqSize(vlistInqVarGrid(vlistID[is], varID));
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID[is], varID));
in[is].missval = vlistInqVarMissval(vlistID[is], varID);
......@@ -198,42 +191,41 @@ Timstat3(void *process)
const auto varID = recList[recID].varID;
const auto levelID = recList[recID].levelID;
missval1 = fwork[0][varID][levelID].missval;
missval2 = missval1;
const auto missval1 = fwork[0][varID][levelID].missval;
const auto missval2 = missval1;
if (operatorID == VARQUOT2TEST)
{
for (size_t i = 0; i < gridsize; ++i)
for (size_t i = 0; i < gridsizemax; ++i)
{
auto fnvals0 = iwork[0][varID][levelID][i];
auto fnvals1 = iwork[1][varID][levelID][i];
auto temp0 = DIVMN(MULMN(fwork[0][varID][levelID].vec_d[i], fwork[0][varID][levelID].vec_d[i]), fnvals0);
auto temp1 = DIVMN(MULMN(fwork[2][varID][levelID].vec_d[i], fwork[2][varID][levelID].vec_d[i]), fnvals1);
auto temp2 = SUBMN(fwork[1][varID][levelID].vec_d[i], temp0);
auto temp3 = SUBMN(fwork[3][varID][levelID].vec_d[i], temp1);
statistic = DIVMN(temp2, ADDMN(temp2, MULMN(rconst, temp3)));
if (fnvals0 <= 1 || fnvals1 <= 1)
fractil_1 = fractil_2 = missval1;
else
const auto fnvals0 = iwork[0][varID][levelID][i];
const auto fnvals1 = iwork[1][varID][levelID][i];
const auto temp0 = DIVMN(MULMN(fwork[0][varID][levelID].vec_d[i], fwork[0][varID][levelID].vec_d[i]), fnvals0);
const auto temp1 = DIVMN(MULMN(fwork[2][varID][levelID].vec_d[i], fwork[2][varID][levelID].vec_d[i]), fnvals1);
const auto temp2 = SUBMN(fwork[1][varID][levelID].vec_d[i], temp0);
const auto temp3 = SUBMN(fwork[3][varID][levelID].vec_d[i], temp1);
const auto statistic = DIVMN(temp2, ADDMN(temp2, MULMN(rconst, temp3)));
double fractil_1 = missval1, fractil_2 = missval1;
if (fnvals0 > 1 && fnvals1 > 1)
cdo::beta_distr_constants((fnvals0 - 1) / 2, (fnvals1 - 1) / 2, 1 - risk, &fractil_1, &fractil_2, __func__);
out[0].vec_d[i]
= DBL_IS_EQUAL(statistic, missval1) ? missval1 : statistic <= fractil_1 || statistic >= fractil_2 ? 1 : 0;
= DBL_IS_EQUAL(statistic, missval1) ? missval1 : (statistic <= fractil_1 || statistic >= fractil_2) ? 1.0 : 0.0;
}
}
else if (operatorID == MEANDIFF2TEST)
{
double fractil;
double mean_factor[NIN], var_factor[NIN];
mean_factor[0] = 1;
mean_factor[1] = -1;
var_factor[0] = var_factor[1] = 1;
mean_factor[0] = 1.0;
mean_factor[1] = -1.0;
var_factor[0] = var_factor[1] = 1.0;
for (size_t i = 0; i < gridsize; ++i)
for (size_t i = 0; i < gridsizemax; ++i)
{
double temp0 = 0;
double temp0 = 0.0;
double deg_of_freedom = -n_in;
for (int j = 0; j < n_in; j++)
{
......@@ -246,7 +238,7 @@ Timstat3(void *process)
if (!DBL_IS_EQUAL(temp0, missval1) && temp0 < 0) // This is possible because of rounding errors
temp0 = 0;
auto stddev_estimator = SQRTMN(DIVMN(temp0, deg_of_freedom));
const auto stddev_estimator = SQRTMN(DIVMN(temp0, deg_of_freedom));
auto mean_estimator = -rconst;
for (int j = 0; j < n_in; j++)
{
......@@ -254,7 +246,7 @@ Timstat3(void *process)
mean_estimator = ADDMN(mean_estimator, MULMN(mean_factor[j], DIVMN(fwork[2 * j][varID][levelID].vec_d[i], fnvals)));
}
double temp1 = 0;
double temp1 = 0.0;
for (int j = 0; j < n_in; j++)
{
const auto fnvals = iwork[j][varID][levelID][i];
......@@ -264,7 +256,7 @@ Timstat3(void *process)
const auto norm = SQRTMN(temp1);
const auto temp2 = DIVMN(DIVMN(mean_estimator, norm), stddev_estimator);
fractil = deg_of_freedom < 1 ? missval1 : cdo::student_t_inv(deg_of_freedom, 1 - risk / 2, __func__);
const auto fractil = (deg_of_freedom < 1) ? missval1 : cdo::student_t_inv(deg_of_freedom, 1 - risk / 2, __func__);
out[0].vec_d[i]
= DBL_IS_EQUAL(temp2, missval1) || DBL_IS_EQUAL(fractil, missval1) ? missval1 : std::fabs(temp2) >= fractil;
......
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