Commit 6e5ce653 authored by Uwe Schulzweida's avatar Uwe Schulzweida
Browse files

Replaced Malloc by std::vector.

parent 5400877a
......@@ -47,14 +47,14 @@
// NO MISSING VALUE SUPPORT ADDED SO FAR
static void
scale_eigvec_grid(double *restrict out, int tsID, size_t npack, const size_t *restrict pack, const double *restrict weight,
scale_eigvec_grid(Varray<double> &out, int tsID, size_t npack, const size_t *restrict pack, const double *restrict weight,
double **covar, double sum_w)
{
for (size_t i = 0; i < npack; ++i) out[pack[i]] = covar[tsID][i] / std::sqrt(weight[pack[i]] / sum_w);
}
static void
scale_eigvec_time(double *restrict out, int tsID, int nts, size_t npack, const size_t *restrict pack, const double *restrict weight,
scale_eigvec_time(Varray<double> &out, int tsID, int nts, size_t npack, const size_t *restrict pack, const double *restrict weight,
double **covar, double **data, double missval, double sum_w)
{
#ifdef _OPENMP
......@@ -258,7 +258,7 @@ EOFs(void *process)
// allocation of temporary fields and output structures
size_t npack = SIZE_MAX;
std::vector<size_t> pack(gridsize);
double *in = (double *) Malloc(gridsize * sizeof(double));
Varray<double> in(gridsize);
std::vector<std::vector<eofdata_t>> eofdata(nvars);
for (varID = 0; varID < nvars; ++varID)
......@@ -296,7 +296,7 @@ EOFs(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID1, &varID, &levelID);
cdoReadRecord(streamID1, in, &nmiss);
cdoReadRecord(streamID1, in.data(), &nmiss);
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto missval = vlistInqVarMissval(vlistID1, varID);
......@@ -407,7 +407,7 @@ EOFs(void *process)
int vtime = 0;
auto juldate = julianDateEncode(calendar, vdate, vtime);
double *out = in;
Varray<double> &out = in;
double *eig_val = nullptr;
int nts_out = nts;
......@@ -463,7 +463,6 @@ EOFs(void *process)
for (size_t ipack = 0; ipack < npack; ++ipack)
{
size_t j;
size_t i = pack[ipack];
for (size_t jpack = 0; jpack < npack; ++jpack)
{
......@@ -473,7 +472,7 @@ EOFs(void *process)
}
else
{
j = pack[jpack];
const auto j = pack[jpack];
covar[ipack][jpack] = covar[ipack][jpack] * // covariance
std::sqrt(weight[i]) * std::sqrt(weight[j]) / sum_w / // weights
nts; // number of data contributing
......@@ -544,9 +543,9 @@ EOFs(void *process)
else if (time_space)
scale_eigvec_time(out, tsID, nts, npack, pack.data(), weight.data(), covar, data, missval, sum_w);
nmiss = arrayNumMV(gridsize, out, missval);
nmiss = varrayNumMV(gridsize, out, missval);
cdoDefRecord(streamID3, varID, levelID);
cdoWriteRecord(streamID3, out, nmiss);
cdoWriteRecord(streamID3, out.data(), nmiss);
} // loop n_eig
nmiss = (DBL_IS_EQUAL(eig_val[tsID], missval)) ? 1 : 0;
......@@ -575,8 +574,6 @@ EOFs(void *process)
}
}
Free(in);
cdoStreamClose(streamID3);
cdoStreamClose(streamID2);
cdoStreamClose(streamID1);
......
......@@ -47,42 +47,33 @@ enum TDATA_TYPE
#define DEBUG_ROC 0
enum CONTINGENCY_TYPE
{
TP,
FP,
FN,
TN
TP, // TP - True positive ( event forecast and occured) HIT
FP, // FP - False positive ( event forecast and not occured) false ALARM
FN, // FN - False negative ( event not forecast and ocurred) MISSED
TN // TN - True negative ( event not forecast and not ocurred) CORRECT REJECTION
};
/* TP - True positive ( event forecast and occured) HIT */
/* FP - False positive ( event forecast and not occured) false ALARM */
/* TN - True negative ( event not forecast and not ocurred) CORRECT REJECTION
*/
/* FN - False negative ( event not forecast and ocurred) MISSED */
enum ROC_ENUM_TYPE
{
TPR,
FPR
TPR, // TPR = True Positive Rate = TP / ( TP + FN )
FPR // FNR = False Negtive Rate = FN / ( FP + TN )
};
/* TPR = True Positive Rate = TP / ( TP + FN ) */
/* FNR = False Negtive Rate = FN / ( FP + TN ) */
static double
roc_curve_integrate(const double **roc, const int n)
{
double y1, y0, x1, x0, dx, dy;
double step_area;
double area = 0;
for (int i = 1; i <= n; ++i)
{
x1 = roc[i][FPR];
x0 = roc[i - 1][FPR];
y1 = roc[i][TPR];
y0 = roc[i - 1][TPR];
dx = x1 - x0;
dy = y1 - y0;
step_area = -0.5 * dx * dy - dx * y0;
const auto x1 = roc[i][FPR];
const auto x0 = roc[i - 1][FPR];
const auto y1 = roc[i][TPR];
const auto y0 = roc[i - 1][TPR];
const auto dx = x1 - x0;
const auto dy = y1 - y0;
const auto step_area = -0.5 * dx * dy - dx * y0;
area += step_area;
}
......@@ -335,11 +326,11 @@ Ensstat3(void *process)
switch (operfunc)
{
case (func_rank):
/* ****************/
/* RANK HISTOGRAM */
/* ************** */
// for ( j=0; j<nfiles; j++ )
// fprintf(stderr,"%5.2g ",field[ompthID].vec[j]);
/* ****************/
/* RANK HISTOGRAM */
/* ************** */
// for ( j=0; j<nfiles; j++ )
// fprintf(stderr,"%5.2g ",field[ompthID].vec[j]);
#ifdef _OPENMP
#pragma omp critical
#endif
......
......@@ -179,8 +179,8 @@ const char *COLOUR = nullptr, *COLOUR_MIN = nullptr, *COLOUR_MAX = nullptr, *STY
*COLOUR_TRIAD = nullptr, *PROJECTION = nullptr;
static void
magplot(const char *plotfile, int operatorID, const char *varname, const char *units, long nlon, long nlat, double *grid_center_lon,
double *grid_center_lat, double *array, int nparam, const std::vector<std::string> &params, char *datetime, bool lregular)
magplot(const char *plotfile, int operatorID, const char *varname, const char *units, long nlon, long nlat, Varray<double> &grid_center_lon,
Varray<double> &grid_center_lat, Varray<double> &array, int nparam, const std::vector<std::string> &params, char *datetime, bool lregular)
{
long i;
......@@ -253,7 +253,7 @@ magplot(const char *plotfile, int operatorID, const char *varname, const char *u
/* Set the input data arrays to magics++ */
mag_set2r("input_field", array, nlon, nlat);
mag_set2r("input_field", array.data(), nlon, nlat);
if (lregular)
{
......@@ -270,8 +270,8 @@ magplot(const char *plotfile, int operatorID, const char *varname, const char *u
{
mag_setc("input_field_organization", "NONREGULAR");
mag_set2r("input_field_latitudes", grid_center_lat, nlon, nlat);
mag_set2r("input_field_longitudes", grid_center_lon, nlon, nlat);
mag_set2r("input_field_latitudes", grid_center_lat.data(), nlon, nlat);
mag_set2r("input_field_longitudes", grid_center_lon.data(), nlon, nlat);
}
/* magics_template_parser( magics_node ); */
......@@ -1144,18 +1144,18 @@ Magplot(void *process)
int nlat = gridInqYsize(gridID);
// int nlev = zaxisInqSize(zaxisID);
double *array = (double *) Malloc(gridsize * sizeof(double));
double *grid_center_lon = (double *) Malloc(gridsize * sizeof(double));
double *grid_center_lat = (double *) Malloc(gridsize * sizeof(double));
Varray<double> array(gridsize);
Varray<double> grid_center_lon(gridsize);
Varray<double> grid_center_lat(gridsize);
gridInqXvals(gridID, grid_center_lon);
gridInqYvals(gridID, grid_center_lat);
gridInqXvals(gridID, grid_center_lon.data());
gridInqYvals(gridID, grid_center_lat.data());
/* Convert lat/lon units if required */
gridInqXunits(gridID, units);
grid_to_degree(units, gridsize, grid_center_lon, "grid center lon");
grid_to_degree(units, gridsize, grid_center_lon.data(), "grid center lon");
gridInqYunits(gridID, units);
grid_to_degree(units, gridsize, grid_center_lat, "grid center lat");
grid_to_degree(units, gridsize, grid_center_lat.data(), "grid center lat");
int tsID = 0;
......@@ -1216,9 +1216,9 @@ Magplot(void *process)
for (int recID = 0; recID < nrecs; recID++)
{
cdoInqRecord(streamID, &varID, &levelID);
cdoReadRecord(streamID, array, &nmiss);
cdoReadRecord(streamID, array.data(), &nmiss);
if (nmiss) cdoSetNAN(vlistInqVarMissval(vlistID, varID), gridsize, array);
if (nmiss) cdoSetNAN(vlistInqVarMissval(vlistID, varID), gridsize, array.data());
vlistInqVarName(vlistID, varID, varname);
vlistInqVarUnits(vlistID, varID, units);
......@@ -1270,10 +1270,6 @@ Magplot(void *process)
cdoStreamClose(streamID);
if (array) Free(array);
if (grid_center_lon) Free(grid_center_lon);
if (grid_center_lat) Free(grid_center_lat);
/* quit_XMLtemplate_parser( ); */
#else
......
......@@ -49,8 +49,8 @@ extern char *DEVICE_TABLE;
extern int DEVICE_COUNT;
static void
magvector(const char *plotfile, int operatorID, const char *varname, long nlon, long nlat, double *grid_center_lon,
double *grid_center_lat, double *uarray, double *varray, int nparam, std::vector<std::string> &params, char *datetime)
magvector(const char *plotfile, int operatorID, const char *varname, long nlon, long nlat, Varray<double> &grid_center_lon,
Varray<double> &grid_center_lat, Varray<double> &uarray, Varray<double> &varray, int nparam, std::vector<std::string> &params, char *datetime)
{
long i;
......@@ -64,13 +64,13 @@ magvector(const char *plotfile, int operatorID, const char *varname, long nlon,
(void) varname;
if (uarray == nullptr && varray == nullptr)
if (uarray.empty() && varray.empty())
{
fprintf(stderr, " No Velocity Components in input file, cannot creaate Vector PLOT!\n");
return;
}
if (uarray == nullptr || varray == nullptr)
if (uarray.empty() || varray.empty())
{
fprintf(stderr, " Found only one Velocity Component in input file, cannot create Vector PLOT!\n");
return;
......@@ -152,8 +152,8 @@ magvector(const char *plotfile, int operatorID, const char *varname, long nlon,
mag_setr("input_field_initial_longitude", grid_center_lon[0]);
mag_setr("input_field_longitude_step", dlon);
mag_set2r("input_wind_u_component", uarray, nlon, nlat);
mag_set2r("input_wind_v_component", varray, nlon, nlat);
mag_set2r("input_wind_u_component", uarray.data(), nlon, nlat);
mag_set2r("input_wind_v_component", varray.data(), nlon, nlat);
mag_seti("map_label_latitude_frequency", 2);
mag_seti("map_label_longitude_frequency", 2);
......@@ -312,10 +312,10 @@ Magvector(void *process)
int nparam = operatorArgc();
std::vector<std::string> &pnames = cdoGetOperArgv();
int VECTOR = cdoOperatorAdd("vector", 0, 0, nullptr);
int STREAM = cdoOperatorAdd("stream", 0, 0, nullptr);
VECTOR = cdoOperatorAdd("vector", 0, 0, nullptr);
STREAM = cdoOperatorAdd("stream", 0, 0, nullptr);
int operatorID = cdoOperatorID();
const auto operatorID = cdoOperatorID();
if (nparam)
{
......@@ -347,19 +347,19 @@ Magvector(void *process)
int nlat = gridInqYsize(gridID);
// int nlev = zaxisInqSize(zaxisID);
double *uarray = (double *) Malloc(gridsize * sizeof(double));
double *varray = (double *) Malloc(gridsize * sizeof(double));
double *grid_center_lat = (double *) Malloc(gridsize * sizeof(double));
double *grid_center_lon = (double *) Malloc(gridsize * sizeof(double));
Varray<double> uarray(gridsize);
Varray<double> varray(gridsize);
Varray<double> grid_center_lat(gridsize);
Varray<double> grid_center_lon(gridsize);
gridInqYvals(gridID, grid_center_lat);
gridInqXvals(gridID, grid_center_lon);
gridInqYvals(gridID, grid_center_lat.data());
gridInqXvals(gridID, grid_center_lon.data());
/* Convert lat/lon units if required */
gridInqXunits(gridID, units);
grid_to_degree(units, gridsize, grid_center_lon, "grid center lon");
grid_to_degree(units, gridsize, grid_center_lon.data(), "grid center lon");
gridInqYunits(gridID, units);
grid_to_degree(units, gridsize, grid_center_lat, "grid center lat");
grid_to_degree(units, gridsize, grid_center_lat.data(), "grid center lat");
int tsID = 0;
......@@ -405,15 +405,15 @@ Magvector(void *process)
if (cstrIsEqual(varname, "var131") || cstrIsEqual(varname, "u")) // U Velocity as per GRIB is var131, as per NC 'u'
{
if (DBG) fprintf(stderr, "Found U VEL in Varname %s\n", varname);
cdoReadRecord(streamID, uarray, &nmiss);
if (nmiss) cdoSetNAN(vlistInqVarMissval(vlistID, varID), gridsize, uarray);
cdoReadRecord(streamID, uarray.data(), &nmiss);
if (nmiss) cdoSetNAN(vlistInqVarMissval(vlistID, varID), gridsize, uarray.data());
found++;
}
if (cstrIsEqual(varname, "var132") || cstrIsEqual(varname, "v")) // V Velocity as per GRIB is var132, as per NC 'v'
{
if (DBG) fprintf(stderr, "Found V VEL in Varname %s\n", varname);
cdoReadRecord(streamID, varray, &nmiss);
if (nmiss) cdoSetNAN(vlistInqVarMissval(vlistID, varID), gridsize, varray);
cdoReadRecord(streamID, varray.data(), &nmiss);
if (nmiss) cdoSetNAN(vlistInqVarMissval(vlistID, varID), gridsize, varray.data());
found++;
}
if (found == 2) break;
......@@ -459,11 +459,6 @@ Magvector(void *process)
cdoStreamClose(streamID);
if (uarray) Free(uarray);
if (varray) Free(varray);
if (grid_center_lon) Free(grid_center_lon);
if (grid_center_lat) Free(grid_center_lat);
/* quit_XMLtemplate_parser( ); */
quit_MAGICS();
......
......@@ -31,7 +31,7 @@
static void
rot_uv_back(int gridID, double *us, double *vs)
rot_uv_back(int gridID, Varray<double> &us, Varray<double> &vs)
{
double xpole = 0, ypole = 0, angle = 0;
if (gridInqType(gridID) == GRID_PROJECTION && gridInqProjType(gridID) == CDI_PROJ_RLL)
......@@ -74,13 +74,11 @@ void *
Rotuv(void *process)
{
int varID, levelID;
int varID1, varID2;
int nrecs;
int chcodes[MAXARG];
const char *chvars[MAXARG];
char varname[CDI_MAX_NAME];
char varname2[CDI_MAX_NAME];
double *single, *usvar = nullptr, *vsvar = nullptr;
cdoInitialize(process);
......@@ -118,8 +116,8 @@ Rotuv(void *process)
const auto maxrecs = vlistNrecs(vlistID1);
std::vector<RecordInfo> recList(maxrecs);
size_t **varnmiss = (size_t **) Malloc(nvars * sizeof(size_t *));
double **vardata = (double **) Malloc(nvars * sizeof(double *));
std::vector<std::vector<size_t>> varnmiss(nvars);
Varray3D<double> vardata(nvars);
bool lfound[MAXARG];
for (int i = 0; i < nch; i++) lfound[i] = false;
......@@ -155,8 +153,9 @@ Rotuv(void *process)
const auto gridsize = gridInqSize(gridID);
const auto nlevel = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
varnmiss[varID] = (size_t *) Malloc(nlevel * sizeof(size_t));
vardata[varID] = (double *) Malloc(gridsize * nlevel * sizeof(double));
varnmiss[varID].resize(nlevel);
vardata[varID].resize(nlevel);
for (levelID = 0; levelID < nlevel; levelID++) vardata[varID][levelID].resize(gridsize);
}
const auto taxisID1 = vlistInqTaxis(vlistID1);
......@@ -180,10 +179,7 @@ Rotuv(void *process)
recList[recID].levelID = levelID;
recList[recID].lconst = vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT;
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto offset = gridsize * levelID;
single = vardata[varID] + offset;
cdoReadRecord(streamID1, single, &varnmiss[varID][levelID]);
cdoReadRecord(streamID1, vardata[varID][levelID].data(), &varnmiss[varID][levelID]);
if (varnmiss[varID][levelID]) cdoAbort("Missing values unsupported for this operator!");
}
......@@ -205,8 +201,7 @@ Rotuv(void *process)
if (varID == nvars) cdoAbort("u-wind not found!");
usvar = vardata[varID];
varID1 = varID;
const auto usvarID = varID;
for (varID = 0; varID < nvars; varID++)
{
......@@ -224,46 +219,36 @@ Rotuv(void *process)
if (varID == nvars) cdoAbort("v-wind not found!");
vsvar = vardata[varID];
varID2 = varID;
const auto vsvarID = varID;
if (Options::cdoVerbose)
{
if (lvar)
{
vlistInqVarName(vlistID2, varID1, varname);
vlistInqVarName(vlistID2, varID2, varname2);
vlistInqVarName(vlistID2, usvarID, varname);
vlistInqVarName(vlistID2, vsvarID, varname2);
cdoPrint("Using var %s [%s](u) and var %s [%s](v)", varname, chvars[i], varname2, chvars[i + 1]);
}
else
cdoPrint("Using code %d [%d](u) and code %d [%d](v)", vlistInqVarCode(vlistID1, varID1), chcodes[i],
vlistInqVarCode(vlistID1, varID2), chcodes[i + 1]);
cdoPrint("Using code %d [%d](u) and code %d [%d](v)", vlistInqVarCode(vlistID1, usvarID), chcodes[i],
vlistInqVarCode(vlistID1, vsvarID), chcodes[i + 1]);
}
const auto gridID = vlistInqVarGrid(vlistID1, varID);
const auto gridsize = gridInqSize(gridID);
const auto nlevel1 = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID1));
const auto nlevel2 = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID2));
const auto nlevel1 = zaxisInqSize(vlistInqVarZaxis(vlistID1, usvarID));
const auto nlevel2 = zaxisInqSize(vlistInqVarZaxis(vlistID1, vsvarID));
if (nlevel1 != nlevel2) cdoAbort("u-wind and v-wind have different number of levels!");
for (levelID = 0; levelID < nlevel1; levelID++)
{
auto offset = gridsize * levelID;
rot_uv_back(gridID, usvar + offset, vsvar + offset);
}
rot_uv_back(gridID, vardata[usvarID][levelID], vardata[vsvarID][levelID]);
}
for (int recID = 0; recID < nrecs; recID++)
{
varID = recList[recID].varID;
levelID = recList[recID].levelID;
const auto gridsize = gridInqSize(vlistInqVarGrid(vlistID1, varID));
const auto offset = gridsize * levelID;
single = vardata[varID] + offset;
cdoDefRecord(streamID2, varID, levelID);
cdoWriteRecord(streamID2, single, varnmiss[varID][levelID]);
cdoWriteRecord(streamID2, vardata[varID][levelID].data(), varnmiss[varID][levelID]);
}
tsID++;
......@@ -272,12 +257,6 @@ Rotuv(void *process)
cdoStreamClose(streamID2);
cdoStreamClose(streamID1);
for (varID = 0; varID < nvars; varID++)
{
Free(varnmiss[varID]);
Free(vardata[varID]);
}
cdoFinish();
return nullptr;
......
......@@ -35,8 +35,8 @@ Setrcaname(void *process)
int scode, sltype, slevel;
int zaxisID, ltype, code, nlev;
int level;
size_t gridsize, nmiss;
double *array = nullptr;
size_t nmiss;
Varray<double> array;
cdoInitialize(process);
......@@ -118,8 +118,8 @@ Setrcaname(void *process)
if (!lcopy)
{
gridsize = vlistGridsizeMax(vlistID1);
array = (double *) Malloc(gridsize * sizeof(double));
const auto gridsizemax = vlistGridsizeMax(vlistID1);
array.resize(gridsizemax);
}
int tsID = 0;
......@@ -139,8 +139,8 @@ Setrcaname(void *process)
}
else
{
cdoReadRecord(streamID1, array, &nmiss);
cdoWriteRecord(streamID2, array, nmiss);
cdoReadRecord(streamID1, array.data(), &nmiss);
cdoWriteRecord(streamID2, array.data(), nmiss);
}
}
......@@ -152,9 +152,6 @@ Setrcaname(void *process)
vlistDestroy(vlistID2);
if (!lcopy)
if (array) Free(array);
cdoFinish();
return nullptr;
......
......@@ -188,6 +188,7 @@ Showattribute(void *process)
}
}
}
cdoStreamClose(streamID);
cdoFinish();
......
This diff is collapsed.
......@@ -532,7 +532,6 @@ cdoDefCompType(CdoStreamID p_streamID, int p_cdi_compression_type)
streamDefCompType(p_streamID->m_fileID, p_cdi_compression_type);
}
/* clang-format off */
bool
unchangedRecord()
{
......@@ -574,4 +573,3 @@ cdoSetNAN(double missval, size_t gridsize, double *array)
if (DBL_IS_EQUAL(array[i], missval)) array[i] = newmissval;
}
}
/* 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