Commit 68c3366c authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Use VarList.

parent b1655787
Pipeline #2906 passed with stages
in 15 minutes and 13 seconds
......@@ -94,9 +94,7 @@ Arithc(void *process)
if (operatorArgc() > 2) cdoAbort("Too many arguments!");
const auto rconst = parameter2double(cdoOperatorArgv(0));
double rconstcplx[2];
rconstcplx[0] = rconst;
rconstcplx[1] = 0;
double rconstcplx[2] = { rconst, 0 };
if (operatorArgc() == 2) rconstcplx[1] = parameter2double(cdoOperatorArgv(1));
const auto streamID1 = cdoOpenRead(0);
......@@ -118,6 +116,9 @@ Arithc(void *process)
if (nwpv == 2 && !opercplx) cdoAbort("Fields with complex numbers are not supported by this operator!");
const auto gridsizemax = nwpv * vlistGridsizeMax(vlistID1);
VarList varList1;
varListInit(varList1, vlistID1);
Field field;
field.resize(gridsizemax);
......@@ -134,11 +135,9 @@ Arithc(void *process)
if (vars[varID])
{
nwpv = vlistInqNWPV(vlistID1, varID);
const auto gridID = vlistInqVarGrid(vlistID1, varID);
field.grid = gridID;
field.size = nwpv * gridInqSize(gridID);
field.missval = vlistInqVarMissval(vlistID1, varID);
field.grid = varList1[varID].gridID;
field.size = varList1[varID].nwpv * varList1[varID].gridsize;
field.missval = varList1[varID].missval;
if (nwpv == 2)
vfarcfuncplx(field, rconstcplx, operfunc);
......
......@@ -59,6 +59,9 @@ Arithlat(void *process)
const auto streamID2 = cdoOpenWrite(1);
cdoDefVlist(streamID2, vlistID2);
VarList varList1;
varListInit(varList1, vlistID1);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Varray<double> array(gridsizemax);
......@@ -74,8 +77,8 @@ Arithlat(void *process)
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, &array[0], &nmiss);
auto gridID = vlistInqVarGrid(vlistID1, varID);
const auto gridsize = gridInqSize(gridID);
auto gridID = varList1[varID].gridID;
const auto gridsize = varList1[varID].gridsize;
if (gridID != gridID0)
{
......
......@@ -62,6 +62,9 @@ Complextorect(void *process)
cdoDefVlist(streamID2, vlistID2);
cdoDefVlist(streamID3, vlistID3);
VarList varList1;
varListInit(varList1, vlistID1);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Varray<double> array1(2 * gridsizemax);
Varray<double> array2(gridsizemax), array3(gridsizemax);
......@@ -81,11 +84,11 @@ Complextorect(void *process)
cdoDefRecord(streamID2, varID, levelID);
cdoDefRecord(streamID3, varID, levelID);
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto gridsize = varList1[varID].gridsize;
cdoReadRecord(streamID1, array1.data(), &nmiss);
const auto missval1 = vlistInqVarMissval(vlistID1, varID);
const auto missval1 = varList1[varID].missval;
const auto missval2 = missval1;
if (operatorID == COMPLEXTORECT)
......
......@@ -60,7 +60,7 @@ enum ROC_ENUM_TYPE
};
static double
roc_curve_integrate(const double **roc, const int n)
roc_curve_integrate(const Varray2D<double> &roc, const int n)
{
double area = 0;
......@@ -93,12 +93,8 @@ Ensstat3(void *process)
int gridID, gridID2;
int have_miss = 0;
CdoStreamID streamID = 0, streamID2 = 0;
int **array2 = nullptr;
int **ctg_tab = nullptr, *hist = nullptr; // contingency table and histogram
double missval;
double *dat; // pointer to ensemble data for ROC
double *uThresh = nullptr, *lThresh = nullptr; // thresholds for histograms
double **roc = nullptr; // receiver operating characteristics table
double *dat; // pointer to ensemble data for ROC
double val;
int ival;
......@@ -106,7 +102,7 @@ Ensstat3(void *process)
{
CdoStreamID streamID;
int vlistID;
double *array;
Varray<double> array;
};
cdoInitialize(process);
......@@ -133,12 +129,12 @@ Ensstat3(void *process)
if (Options::cdoVerbose) cdoPrint("Ensemble over %d files.", nfiles);
const char *ofilename = cdoGetStreamName(nfiles);
auto ofilename = cdoGetStreamName(nfiles);
if (!Options::cdoOverwriteMode && fileExists(ofilename) && !userFileOverwrite(ofilename))
cdoAbort("Outputfile %s already exists!", ofilename);
ens_file_t *ef = (ens_file_t *) Malloc(nfiles * sizeof(ens_file_t));
std::vector<ens_file_t> ef(nfiles);
/* *************************************************** */
/* should each thread be allocating memory locally???? */
......@@ -167,13 +163,15 @@ Ensstat3(void *process)
const auto vlistID1 = ef[0].vlistID;
const auto vlistID2 = vlistCreate();
const auto nvars = vlistNvars(vlistID1);
int *varID2 = (int *) Malloc(nvars * sizeof(int));
std::vector<int> varIDs2(nvars);
double *levs = (double *) Calloc(nfiles, sizeof(double));
const auto zaxisID2 = zaxisCreate(ZAXIS_GENERIC, nfiles);
for (int i = 0; i < nfiles; i++) levs[i] = i;
zaxisDefLevels(zaxisID2, levs);
zaxisDefName(zaxisID2, "histogram_binID");
{
Varray<double> levs(nfiles, 0);
for (int i = 0; i < nfiles; i++) levs[i] = i;
zaxisDefLevels(zaxisID2, levs.data());
zaxisDefName(zaxisID2, "histogram_binID");
}
int time_mode = datafunc == TIME ? TIME_VARYING : TIME_CONSTANT;
......@@ -195,7 +193,7 @@ Ensstat3(void *process)
else // datafunc == SPACE
gridID2 = vlistInqVarGrid(vlistID1, varID);
varID2[varID] = vlistDefVar(vlistID2, gridID2, zaxisID2, time_mode);
varIDs2[varID] = vlistDefVar(vlistID2, gridID2, zaxisID2, time_mode);
}
const int taxisID1 = vlistInqTaxis(vlistID1);
......@@ -221,36 +219,41 @@ Ensstat3(void *process)
auto gridsize = vlistGridsizeMax(vlistID1);
for (int fileID = 0; fileID < nfiles; fileID++) ef[fileID].array = (double *) Malloc(gridsize * sizeof(double));
for (int fileID = 0; fileID < nfiles; fileID++) ef[fileID].array.resize(gridsize);
std::vector<std::vector<int>> array2(nfiles + 1);
if (operfunc == func_rank && datafunc == SPACE)
{ // need to memorize data for entire grid before writing
array2 = (int **) Malloc((nfiles + 1) * sizeof(int *));
for (binID = 0; binID < nfiles; binID++) array2[binID] = (int *) Calloc(gridsize, sizeof(int));
for (binID = 0; binID < nfiles; binID++) array2[binID].resize(gridsize, 0);
}
else if (operfunc == func_rank)
{ // can process data separately for each timestep and only need to cumulate values over the grid
array2 = (int **) Malloc((nfiles + 1) * sizeof(int *));
for (binID = 0; binID < nfiles; binID++) array2[binID] = (int *) Calloc(1, sizeof(int));
for (binID = 0; binID < nfiles; binID++) array2[binID].resize(1, 0);
}
std::vector<int> hist;
Varray2D<double> roc; // receiver operating characteristics table
Varray2D<int> ctg_tab; // contingency table and histogram
Varray<double> uThresh, lThresh; // thresholds for histograms
if (operfunc == func_roc)
{
ctg_tab = (int **) Malloc((nbins + 1) * sizeof(int *));
hist = (int *) Malloc(nbins * sizeof(int));
uThresh = (double *) Malloc(nbins * sizeof(double));
lThresh = (double *) Malloc(nbins * sizeof(double));
roc = (double **) Malloc((nbins + 1) * sizeof(double *));
hist.resize(nbins);
roc.resize(nbins + 1);
ctg_tab.resize(nbins + 1);
for (int i = 0; i <= nbins; i++)
{
roc[i].resize(2, 0);
ctg_tab[i].resize(4, 0);
}
uThresh.resize(nbins);
lThresh.resize(nbins);
for (int i = 0; i < nbins; i++)
{
ctg_tab[i] = (int *) Calloc(4, sizeof(int));
roc[i] = (double *) Calloc(2, sizeof(double));
uThresh[i] = ((double) i + 1) / nbins;
lThresh[i] = (double) i / nbins;
}
ctg_tab[nbins] = (int *) Calloc(4, sizeof(int));
roc[nbins] = (double *) Calloc(2, sizeof(double));
}
int tsID = 0;
......@@ -288,7 +291,7 @@ Ensstat3(void *process)
{
streamID = ef[fileID].streamID;
cdoInqRecord(streamID, &varID, &levelID);
cdoReadRecord(streamID, ef[fileID].array, &nmiss);
cdoReadRecord(streamID, ef[fileID].array.data(), &nmiss);
}
gridID = vlistInqVarGrid(vlistID1, varID);
......@@ -403,7 +406,7 @@ Ensstat3(void *process)
{
val = (double) array2[binID][0];
// fprintf(stderr,"%i ",(int)val);
cdoDefRecord(streamID2, varID2[varID], binID);
cdoDefRecord(streamID2, varIDs2[varID], binID);
cdoWriteRecord(streamID2, &val, nmiss);
}
// fprintf(stderr,"\n");
......@@ -429,7 +432,7 @@ Ensstat3(void *process)
ctg_tab[binID][0], ctg_tab[binID][1], ctg_tab[binID][2], ctg_tab[binID][3], tpr, fpr);
}
fprintf(stderr, "nbins %10i\n", nbins);
fprintf(stderr, "#ROC CurveArea: %10.6f\n", roc_curve_integrate((const double **) roc, nbins));
fprintf(stderr, "#ROC CurveArea: %10.6f\n", roc_curve_integrate(roc, nbins));
} // if ( DEBUG_ROC )
} // else if (operfunc == func_roc )
} // for ( recID=0; recID<nrecs; recID++ )
......@@ -446,7 +449,7 @@ Ensstat3(void *process)
{
for (int i = 0; i < osize; i++) tmpdoub[i] = (double) array2[binID][i];
cdoDefRecord(streamID2, varID2[varID], binID);
cdoDefRecord(streamID2, varIDs2[varID], binID);
cdoWriteRecord(streamID2, tmpdoub.data(), nmiss);
}
}
......@@ -471,28 +474,13 @@ Ensstat3(void *process)
ctg_tab[i][1], ctg_tab[i][2], ctg_tab[i][3], sum, tpr, fpr);
}
fprintf(stdout, "#ROC CurveArea: %10.6f\n", roc_curve_integrate((const double **) roc, nbins));
fprintf(stdout, "#ROC CurveArea: %10.6f\n", roc_curve_integrate(roc, nbins));
}
for (int fileID = 0; fileID < nfiles; fileID++)
{
streamID = ef[fileID].streamID;
cdoStreamClose(streamID);
}
for (int fileID = 0; fileID < nfiles; fileID++) cdoStreamClose(ef[fileID].streamID);
if (operfunc != func_roc) cdoStreamClose(streamID2);
for (int fileID = 0; fileID < nfiles; fileID++)
if (ef[fileID].array) Free(ef[fileID].array);
if (ef) Free(ef);
if (array2)
{
for (binID = 0; binID < nfiles; binID++) Free(array2[binID]);
Free(array2);
}
cdoFinish();
return nullptr;
......
......@@ -235,6 +235,9 @@ Ensval(void *process)
if (Options::cdoVerbose) cdoPrint(" sum_weights %10.6f", sum_weights);
VarList varList1;
varListInit(varList1, vlistID1);
int tsID = 0;
do
{
......@@ -263,8 +266,8 @@ Ensval(void *process)
if (fileID == 0)
{
gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
missval = vlistInqVarMissval(vlistID1, varID);
gridsize = varList1[varID].gridsize;
missval = varList1[varID].missval;
weights.resize(gridsize);
}
......@@ -362,8 +365,7 @@ Ensval(void *process)
}
}
// brs_o[i] - number of times that the obs is between
// Ensemble i-1 and i
// brs_o[i] - number of times that the obs is between Ensemble i-1 and i
if (1)
{
if (x[0] > xa)
......@@ -386,8 +388,7 @@ Ensval(void *process)
if (operfunc == CRPS)
{
// First Bin
double p = 0.;
double g = 0.;
double p = 0., g = 0.;
double o = heavyside0 / gridsize;
if (o > 0.) g = beta[0] / o;
crps_reli = g * (o - p) * (o - p);
......@@ -435,9 +436,7 @@ Ensval(void *process)
brs_resol = 0;
brs_uncty = 0;
double gsum = 0;
double obar = 0;
double osum = 0;
double gsum = 0, obar = 0, osum = 0;
for (k = 0; k <= nens; k++)
{
obar += brs_g[k] * brs_o[k];
......@@ -451,9 +450,7 @@ Ensval(void *process)
cdoAbort("This is likely due to missing values");
}
double o = 0;
double p = 0;
double g = 0;
double o = 0, p = 0, g = 0;
brs_uncty = obar * (1 - obar);
for (k = 0; k <= nens; k++)
......
......@@ -499,7 +499,6 @@ Fillmiss(void *process)
/* Argument handling */
{
const auto oargc = operatorArgc();
if (oargc == 1)
{
nfill = parameter2int(cdoOperatorArgv(0));
......
......@@ -194,6 +194,9 @@ Fldstat(void *process)
}
}
VarList varList1;
varListInit(varList1, vlistID1);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
......@@ -209,8 +212,8 @@ Fldstat(void *process)
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, field.vec_d.data(), &field.nmiss);
field.grid = vlistInqVarGrid(vlistID1, varID);
field.size = gridInqSize(field.grid);
field.grid = varList1[varID].gridID;
field.size = varList1[varID].gridsize;
if (needWeights && field.grid != lastgrid)
{
......@@ -236,7 +239,7 @@ Fldstat(void *process)
}
}
field.missval = vlistInqVarMissval(vlistID1, varID);
field.missval = varList1[varID].missval;
auto sglval = (operfunc == func_pctl) ? vfldpctl(field, pn) : vfldfun(field, operfunc);
if (Options::cdoVerbose && (operfunc == func_min || operfunc == func_max))
......
......@@ -80,6 +80,10 @@ Intyear(void *process)
cdoDefVlist(streamIDs[iy], vlistID3);
}
VarList varList1, varList2;
varListInit(varList1, vlistID1);
varListInit(varList2, vlistID2);
int tsID = 0;
while (true)
{
......@@ -113,7 +117,7 @@ Intyear(void *process)
cdoReadRecord(streamID1, array1.data(), &nmiss1);
cdoReadRecord(streamID2, array2.data(), &nmiss2);
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto gridsize = varList1[varID].gridsize;
for (int iy = 0; iy < nyears; iy++)
{
......@@ -124,8 +128,8 @@ Intyear(void *process)
if (nmiss1 || nmiss2)
{
const auto missval1 = vlistInqVarMissval(vlistID1, varID);
const auto missval2 = vlistInqVarMissval(vlistID2, varID);
const auto missval1 = varList1[varID].missval;
const auto missval2 = varList2[varID].missval;
for (size_t i = 0; i < gridsize; i++)
{
......
......@@ -267,6 +267,9 @@ Math(void *process)
const auto streamID2 = cdoOpenWrite(1);
cdoDefVlist(streamID2, vlistID2);
VarList varList1;
varListInit(varList1, vlistID1);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
......@@ -279,9 +282,9 @@ Math(void *process)
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, &array1[0], &nmiss);
const auto missval1 = vlistInqVarMissval(vlistID1, varID);
const auto n = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto number = vlistInqVarNumber(vlistID1, varID);
const auto missval1 = varList1[varID].missval;
const auto n = varList1[varID].gridsize;
const auto number = varList1[varID].nwpv;
if (number == CDI_REAL)
{
......
......@@ -76,6 +76,9 @@ Ninfo(void *process)
const auto ntsteps = vlistNtsteps(vlistID);
const auto ngrids = vlistNgrids(vlistID);
VarList varList;
varListInit(varList, vlistID);
switch (operfunc)
{
case NYEAR:
......@@ -152,13 +155,13 @@ Ninfo(void *process)
case NLEVEL:
for (varID = 0; varID < nvars; varID++)
{
fprintf(stdout, "%d\n", zaxisInqSize(vlistInqVarZaxis(vlistID, varID)));
fprintf(stdout, "%d\n", varList[varID].nlevels);
}
break;
case NGRIDPOINTS:
for (varID = 0; varID < nvars; varID++)
{
fprintf(stdout, "%zu\n", gridInqSize(vlistInqVarGrid(vlistID, varID)));
fprintf(stdout, "%zu\n", varList[varID].gridsize);
}
break;
case NGRIDS: fprintf(stdout, "%d\n", ngrids); break;
......
......@@ -349,7 +349,7 @@ Output(void *process)
const auto zaxisID = varList[varID].zaxisID;
const auto datatype = varList[varID].datatype;
const int dig = (datatype == CDI_DATATYPE_FLT64) ? Options::CDO_dbl_digits : Options::CDO_flt_digits;
gridsize = varList[varID].nwpv * gridInqSize(gridID);
gridsize = varList[varID].nwpv * varList[varID].gridsize;
auto nlon = gridInqXsize(gridID);
auto nlat = gridInqYsize(gridID);
const auto level = cdoZaxisInqLevel(zaxisID, levelID);
......
......@@ -54,9 +54,10 @@ Recttocomplex(void *process)
cdoDefVlist(streamID3, vlistID3);
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Varray<double> array1(gridsizemax);
Varray<double> array2(gridsizemax);
Varray<double> array3(2 * gridsizemax);
Varray<double> array1(gridsizemax), array2(gridsizemax), array3(2 * gridsizemax);
VarList varList1;
varListInit(varList1, vlistID1);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
......@@ -71,11 +72,10 @@ Recttocomplex(void *process)
cdoInqRecord(streamID2, &varID, &levelID);
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
cdoReadRecord(streamID1, array1.data(), &nmiss);
cdoReadRecord(streamID2, array2.data(), &nmiss);
const auto gridsize = varList1[varID].gridsize;
for (size_t i = 0; i < gridsize; ++i)
{
array3[2 * i] = array1[i];
......
......@@ -91,6 +91,9 @@ Replacevalues(void *process)
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Varray<double> array(gridsizemax);
VarList varList1;
varListInit(varList1, vlistID1);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
......@@ -102,8 +105,8 @@ Replacevalues(void *process)
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, array.data(), &nmiss);
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto missval = vlistInqVarMissval(vlistID1, varID);
const auto gridsize = varList1[varID].gridsize;
const auto missval = varList1[varID].missval;
if (operatorID == SETVALS)
{
......
......@@ -71,6 +71,9 @@ Runpctl(void *process)
for (int its = 0; its < ndates; its++) fieldsFromVlist(vlistID1, vars1[its], FIELD_VEC);
VarList varList1;
varListInit(varList1, vlistID1);
int tsID;
for (tsID = 0; tsID < ndates; tsID++)
{
......@@ -101,11 +104,11 @@ Runpctl(void *process)
{
for (int varID = 0; varID < nvars; varID++)
{
if (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT) continue;
if (varList1[varID].timetype == TIME_CONSTANT) continue;
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
const auto missval = vlistInqVarMissval(vlistID1, varID);
const auto gridsize = varList1[varID].gridsize;
const auto nlevels = varList1[varID].nlevels;
const auto missval = varList1[varID].missval;
for (int levelID = 0; levelID < nlevels; levelID++)
{
......
......@@ -91,11 +91,12 @@ Seaspctl(void *process)
fieldsFromVlist(vlistID1, vars1, FIELD_VEC);
HistogramSet hset(nvars);
VarList varList1;
varListInit(varList1, vlistID1);
for (int varID = 0; varID < nvars; varID++)
{
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
hset.createVarLevels(varID, nlevels, gridsize);
hset.createVarLevels(varID, varList1[varID].nlevels, varList1[varID].gridsize);
}
int tsID = 0;
......@@ -181,9 +182,9 @@ Seaspctl(void *process)
for (int varID = 0; varID < nvars; varID++)
{
if (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT) continue;
if (varList1[varID].timetype == TIME_CONSTANT) continue;
const auto nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
const auto nlevels = varList1[varID].nlevels;
for (int levelID = 0; levelID < nlevels; levelID++) hset.getVarLevelPercentiles(vars1[varID][levelID], varID, levelID, pn, FIELD_VEC);
}
......
......@@ -120,6 +120,9 @@ Setmiss(void *process)
const auto gridsizemax = vlistGridsizeMax(vlistID1);
Varray<double> array(gridsizemax);
VarList varList1;
varListInit(varList1, vlistID1);
int tsID = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID)))
{
......@@ -131,8 +134,8 @@ Setmiss(void *process)
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, array.data(), &nmiss);
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto missval = vlistInqVarMissval(vlistID1, varID);
const auto gridsize = varList1[varID].gridsize;
const auto missval = varList1[varID].missval;
if (operatorID == SETMISSVAL)
{
......
......@@ -361,6 +361,9 @@ Setpartab(void *process)
if (vlistNumber(vlistID1) != CDI_REAL) gridsizemax *= 2;
Varray<double> array(gridsizemax);
VarList varList2;
varListInit(varList2, vlistID2);
int tsID1 = 0;
while ((nrecs = cdoStreamInqTimestep(streamID1, tsID1)))
{
......@@ -393,9 +396,8 @@ Setpartab(void *process)