Commit 0300b798 authored by Fabian Wachsmann's avatar Fabian Wachsmann
Browse files

Try to calc cei by bootstrapping histo

parent f885e7e1
......@@ -328,7 +328,6 @@ AC_CONFIG_FILES([
test/Timselpctl.test
test/Timstat.test
test/Timstat2.test
test/UserInput.test
test/Varsstat.test
test/Vertint.test
test/Vertstat.test
......
......@@ -103,6 +103,7 @@ EcaEtccdi(void *process)
const auto taxisID4 = taxisDuplicate(taxisID1);
if (taxisHasBounds(taxisID4)) taxisDeleteBounds(taxisID4);
vlistDefTaxis(vlistID4, taxisID4);
int dpy = 0;
const auto streamID4 = cdoOpenWrite(3);
cdoDefVlist(streamID4, vlistID4);
......@@ -129,6 +130,7 @@ EcaEtccdi(void *process)
fieldsFromVlist(vlistID1, vars1[its], FIELD_VEC);
int tsID = 0;
int lastDay = 0;
while ((nrecs = cdoStreamInqTimestep(streamID2, tsID)))
{
if (nrecs != cdoStreamInqTimestep(streamID3, tsID))
......@@ -166,6 +168,7 @@ EcaEtccdi(void *process)
auto nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID2, varID));
hsets[dayoy].createVarLevels(varID, nlevels, gridsize);
}
lastDay = dayoy;
}
for (int recID = 0; recID < nrecs; recID++)
......@@ -190,8 +193,26 @@ EcaEtccdi(void *process)
tsID++;
}
for ( dayoy = 1; dayoy < MaxDays; dayoy++ )
for ( int dayoy = 1; dayoy < MaxDays; dayoy++)
{
if (wdaysSrc[dayoy])
dpy++;
}
int percyear[(MaxDays-1)*sumboot+1][ndates];
/* windowDays()
Saves for each day of year in the first dimension
and for each window day in the second dimension.
the day of year that needs to be taken into account for index=0 in the third dimension
and the number of window-day-exceeds-maxDays-or-1 for index=1 in the third dimension.
void
windowDays(int *wdays, int *wdaysSrc, int MaxDays, int ndates)
{ */
#ifdef _OPENMP
#pragma omp parallel shared(percyear, wdays, ndates, wdaysSrc, cei)
{
#pragma omp for schedule(dynamic)
#endif
for ( int dayoy = 1; dayoy < MaxDays; dayoy++ )
{
wdays[dayoy][0][1] = 0;
int wdayoy = dayoy;
......@@ -234,30 +255,38 @@ EcaEtccdi(void *process)
}
}
int daysPerYear = 0;
int lastDay = 0;
int percyear[(MaxDays-1)*sumboot+1][ndates];
/* windowYears()
Saves for each day of year in the first dimension
and for each window day of year in the second dimension
the year of the window day of year
void
windowYears(int *windowYears, int sumboot, int *wdays, int *wdaysSrc, int MaxDays, int ndates, int ytoadd)
{ */
#ifdef _OPENMP
#pragma omp for schedule(dynamic)
#endif
for ( int ytoadd = 0; ytoadd < sumboot; ytoadd++)
{
fieldsFromVlist(vlistID1, cei[ytoadd], FIELD_VEC);
for ( dayoy = 1; dayoy < MaxDays; dayoy++)
for ( int dayoy = 1; dayoy < MaxDays; dayoy++)
{
if (wdaysSrc[dayoy])
{
if ( ytoadd == 0 ) { daysPerYear++; lastDay = dayoy; }
int temp = dayoy+ytoadd*(MaxDays-1);
int dayInTimeseries = dayoy+ytoadd*(MaxDays-1);
for ( int ano = 0; ano < ndates; ano++)
{
percyear[temp][ano] = startboot + ytoadd + wdays[dayoy][ano][1];
if ( percyear[temp][ano] < startboot )
percyear[temp][ano] = endboot;
else if ( percyear[temp][ano] > endboot )
percyear[temp][ano] = startboot;
percyear[dayInTimeseries][ano] = startboot + ytoadd + wdays[dayoy][ano][1];
if ( percyear[dayInTimeseries][ano] < startboot )
percyear[dayInTimeseries][ano] = endboot;
else if ( percyear[dayInTimeseries][ano] > endboot )
percyear[dayInTimeseries][ano] = startboot;
}
}
}
}
#ifdef _OPENMP
}
#endif
tsID = 0;
while ( (nrecs = cdoStreamInqTimestep(streamID1, tsID++) ) )
{
......@@ -301,100 +330,84 @@ EcaEtccdi(void *process)
int bootsyear = 0;
int subyear = 0;
int otsID = 0;
for ( bootsyear = startboot; bootsyear < endboot+1; bootsyear++ )
{
for ( dayoy = 1; dayoy < MaxDays; dayoy++)
{
if (wdaysSrc[dayoy])
{
for (varID = 0; varID < nvars; varID++)
{
if (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT) continue;
nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevels; levelID++)
{
if ( bootsyear != startboot )
#ifdef _OPENMP
#pragma omp parallel for shared(sumboot, wdaysSrc, ndates, vars1, percyear, startboot, varID, levelID, hsets, wdays, varsPtemp, pn, vars2, cei) schedule(dynamic)
#endif
for ( int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
{
if (wdaysSrc[loopdoy])
{
/*Recreate orginal time series */
Field &source1 = varsTemp[dayoy][varID][levelID];
Field &target1 = vars1[dayoy+(bootsyear-1-startboot)*(MaxDays-1)][varID][levelID];
vfarcpy(target1, source1);
target1.nmiss = source1.nmiss;
}
/*Save bootstrapping year values in varsTemp */
if (vars1[dayoy+(bootsyear-startboot)*(MaxDays-1)][varID][levelID].vec.data())
for ( int ytoadd = 0; ytoadd < sumboot; ytoadd++)
{
Field &source = vars1[dayoy+(bootsyear-startboot)*(MaxDays-1)][varID][levelID];
Field &target = varsTemp[dayoy][varID][levelID];
vfarcpy(target, source);
target.nmiss = vars1[dayoy+(bootsyear-startboot)*(MaxDays-1)][varID][levelID].nmiss;
for ( int ano = 0; ano < ndates; ano++)
{
if (vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID].vec.data())
{
Field &source = vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID];
HistogramSet &hset = hsets[loopdoy];
hset.addSubVarLevelValues(varID, levelID, source, true);
/* cdoPrint("There are values out of bounds for day '%d', windowday: '%d', year '%d'\n", loopdoy, wdays[loopdoy][ano][0], percyear[loopdoy+ytoadd*(MaxDays-1)][ano]); */
}
else
cdoAbort("Save: Do not have data for: %d %d\n", dayoy, bootsyear);
cdoWarning("Addvarlevel: Do not have data for: %d %d\n", wdays[loopdoy][ano][0], percyear[loopdoy+ytoadd*(MaxDays-1)][ano]);
}
}
}
}
if (Options::cdoVerbose) cdoPrint("Bootsyear: %d\n", bootsyear);
for ( subyear = startboot; subyear < endboot+1; subyear++ )
{
if (Options::cdoVerbose) cdoPrint("Subyear: %d\n", subyear);
if ( subyear != bootsyear )
{
for ( dayoy = 1; dayoy < MaxDays; dayoy++)
{
if (wdaysSrc[dayoy])
}
}
for ( bootsyear = startboot; bootsyear < startboot+2; bootsyear++ )
{
if (Options::cdoVerbose) cdoPrint("Bootsyear: %d\n", bootsyear);
for (varID = 0; varID < nvars; varID++)
{
if (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT) continue;
nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevels; levelID++)
{
#ifdef _OPENMP
#pragma omp parallel sections
#pragma omp parallel for shared(sumboot, wdaysSrc, ndates, vars1, percyear, startboot, varID, levelID, hsets, wdays, varsPtemp, pn, vars2, cei) schedule(dynamic)
#endif
for ( int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
{
#ifdef _OPENMP
#pragma omp section
#endif
/*Reset histograms to begin with a new substitution year */
if (wdaysSrc[loopdoy])
{
hsets[dayoy].Reset(varID, levelID);
}
#ifdef _OPENMP
#pragma omp section
#endif
/*Copy one year in vars1 */
for ( int ytoadd = 0; ytoadd < sumboot; ytoadd++)
{
if (vars1[dayoy+(bootsyear-startboot)*(MaxDays-1)][varID][levelID].vec.data())
for ( int ano = 0; ano < ndates; ano++)
{
Field &source = vars1[dayoy+(subyear-startboot)*(MaxDays-1)][varID][levelID];
Field &target = vars1[dayoy+(bootsyear-startboot)*(MaxDays-1)][varID][levelID];
vfarcpy(target, source);
target.nmiss = vars1[dayoy+(subyear-startboot)*(MaxDays-1)][varID][levelID].nmiss;
}
else
cdoAbort("Copy bootyear: Do not have data for: %d %d\n", dayoy, bootsyear);
if ( percyear[loopdoy+ytoadd*(MaxDays-1)][ano] == bootsyear )
{
Field &source = vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID];
HistogramSet &hset = hsets[loopdoy];
hset.addSubVarLevelValues(varID, levelID, source, false);
}
/*percyear cannot be smaller than startboot */
else if ( percyear[loopdoy+ytoadd*(MaxDays-1)][ano] == bootsyear -1 )
{
Field &source = vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID];
HistogramSet &hset = hsets[loopdoy];
hset.addSubVarLevelValues(varID, levelID, source, true);
}
}
}
}
}
/* printf("Haben die Jahre ausgetauscht.\n");*/
/*Add values to histogram */
for (varID = 0; varID < nvars; varID++)
for ( subyear = startboot; subyear < endboot+1; subyear++ )
{
if (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT) continue;
nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevels; levelID++)
if (Options::cdoVerbose && varID == 0 && levelID == 0) cdoPrint("Subyear: %d\n", subyear);
if ( subyear != bootsyear )
{
#ifdef _OPENMP
#pragma omp parallel for shared(sumboot, wdaysSrc, ndates, vars1, percyear, startboot, varID, levelID, hsets, wdays, varsPtemp, varsTemp, pn, vars2, cei) schedule(dynamic)
#pragma omp parallel for shared(sumboot, wdaysSrc, ndates, vars1, percyear, startboot, varID, levelID, hsets, wdays, varsPtemp, pn, vars2, cei) schedule(dynamic)
#endif
for ( int loopdoy = 1; loopdoy < MaxDays; loopdoy++)
{
......@@ -404,31 +417,27 @@ EcaEtccdi(void *process)
{
for ( int ano = 0; ano < ndates; ano++)
{
/*
#ifdef _OPENMP
int tid = pthread_self();
printf("Schau: %d %d %d %d %d\n", loopdoy, wdays[loopdoy][ano][0], percyear[loopdoy+ytoadd*(MaxDays-1)][ano], wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1), ano);
#endif
*/
if ( percyear[loopdoy+ytoadd*(MaxDays-1)][ano] == subyear )
{
if (vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID].vec.data())
{
Field &source = vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID];
HistogramSet &hset = hsets[loopdoy];
hset.addVarLevelValues(varID, levelID, source);
hset.addSubVarLevelValues(varID, levelID, source, true);
/* cdoPrint("There are values out of bounds for day '%d', windowday: '%d', year '%d'\n", loopdoy, wdays[loopdoy][ano][0], percyear[loopdoy+ytoadd*(MaxDays-1)][ano]); */
}
else
cdoWarning("Addvarlevel: Do not have data for: %d %d\n", wdays[loopdoy][ano][0], percyear[loopdoy+ytoadd*(MaxDays-1)][ano]);
}
}
}
/* printf("Haben es zum temp array addiert.\n");*/
/*** Calculate percentile ***/
Field &pctls = varsPtemp[loopdoy][varID][levelID];
hsets[loopdoy].getVarLevelPercentiles(pctls, varID, levelID, pn);
/*** Compare data with percentile ***/
Field &source = varsTemp[loopdoy][varID][levelID];
Field &source = vars1[dayoy+(bootsyear-startboot)*(MaxDays-1)][varID][levelID];
Field &toCompare = vars2[loopdoy][varID][levelID];
vfarcpy(toCompare, source);
vfarselge(toCompare, pctls);
......@@ -440,23 +449,32 @@ EcaEtccdi(void *process)
#pragma omp critical
#endif
vfaradd(cei[bootsyear-startboot][varID][levelID], vars2[loopdoy][varID][levelID]);
for ( int ytoadd = 0; ytoadd < sumboot; ytoadd++)
{
for ( int ano = 0; ano < ndates; ano++)
{
if ( percyear[loopdoy+ytoadd*(MaxDays-1)][ano] == subyear )
{
if (vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID].vec.data())
{
Field &source = vars1[wdays[loopdoy][ano][0]+(percyear[loopdoy+ytoadd*(MaxDays-1)][ano]-startboot)*(MaxDays-1)][varID][levelID];
HistogramSet &hset = hsets[loopdoy];
hset.addSubVarLevelValues(varID, levelID, source, false);
/* cdoPrint("There are values out of bounds for day '%d', windowday: '%d', year '%d'\n", loopdoy, wdays[loopdoy][ano][0], percyear[loopdoy+ytoadd*(MaxDays-1)][ano]); */
}
else
cdoWarning("Addvarlevel: Do not have data for: %d %d\n", wdays[loopdoy][ano][0], percyear[loopdoy+ytoadd*(MaxDays-1)][ano]);
}
}
}
/* printf("Haben ein climaindex fuer ein Jahr berechnet.\n");*/
}
}
}
}
for (varID = 0; varID < nvars; varID++)
{
if (vlistInqVarTimetype(vlistID1, varID) == TIME_CONSTANT) continue;
nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevels; levelID++)
{
if (cei[bootsyear-startboot][varID][levelID].vec.data())
{
/* Divide vars2 to receive average */
vfarcdiv(cei[bootsyear-startboot][varID][levelID], (double) ((sumboot-1)*daysPerYear));
vfarcdiv(cei[bootsyear-startboot][varID][levelID], (double) ((sumboot-1)*dpy/100.0));
}
}
}
......
......@@ -170,7 +170,7 @@ Seaspctl(void *process)
cdoReadRecord(streamID1, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
hset.addVarLevelValues(varID, levelID, vars1[varID][levelID]);
hset.addSubVarLevelValues(varID, levelID, vars1[varID][levelID], true);
}
nsets++;
......
......@@ -158,7 +158,7 @@ timpctl(int operatorID)
cdoReadRecord(streamID1, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
hset.addVarLevelValues(varID, levelID, vars1[varID][levelID]);
hset.addSubVarLevelValues(varID, levelID, vars1[varID][levelID], true);
}
nsets++;
......
......@@ -186,7 +186,7 @@ Timselpctl(void *process)
cdoReadRecord(streamID1, vars1[varID][levelID].vec.data(), &nmiss);
vars1[varID][levelID].nmiss = nmiss;
hset.addVarLevelValues(varID, levelID, vars1[varID][levelID]);
hset.addSubVarLevelValues(varID, levelID, vars1[varID][levelID], true);
}
tsID++;
......
......@@ -171,7 +171,7 @@ Ydaypctl(void *process)
cdoReadRecord(streamID1, vars1[dayoy][varID][levelID].vec.data(), &nmiss);
vars1[dayoy][varID][levelID].nmiss = nmiss;
hsets[dayoy].addVarLevelValues(varID, levelID, vars1[dayoy][varID][levelID]);
hsets[dayoy].addSubVarLevelValues(varID, levelID, vars1[dayoy][varID][levelID], true);
}
nsets[dayoy]++;
......
......@@ -198,7 +198,7 @@ Ydrunpctl(void *process)
const auto nlevels = zaxisInqSize(vlistInqVarZaxis(vlistID1, varID));
for (levelID = 0; levelID < nlevels; levelID++)
for (inp = 0; inp < ndates; inp++) hsets[dayoy].addVarLevelValues(varID, levelID, vars1[inp][varID][levelID]);
for (inp = 0; inp < ndates; inp++) hsets[dayoy].addSubVarLevelValues(varID, levelID, vars1[inp][varID][levelID], true);
}
datetime[ndates] = datetime[0];
......
......@@ -171,7 +171,7 @@ Ymonpctl(void *process)
cdoReadRecord(streamID1, vars1[month][varID][levelID].vec.data(), &nmiss);
vars1[month][varID][levelID].nmiss = nmiss;
hsets[month].addVarLevelValues(varID, levelID, vars1[month][varID][levelID]);
hsets[month].addSubVarLevelValues(varID, levelID, vars1[month][varID][levelID], true);
}
nsets[month]++;
......
......@@ -177,7 +177,7 @@ Yseaspctl(void *process)
cdoReadRecord(streamID1, vars1[seas][varID][levelID].vec.data(), &nmiss);
vars1[seas][varID][levelID].nmiss = nmiss;
hsets[seas].addVarLevelValues(varID, levelID, vars1[seas][varID][levelID]);
hsets[seas].addSubVarLevelValues(varID, levelID, vars1[seas][varID][levelID], true);
}
nsets[seas]++;
......
......@@ -62,6 +62,15 @@ histBinValue(Histogram &hist, double value)
if (bin >= 0 && bin < hist.nbins) INT_PTR(hist.ptr)[bin]++;
}
static void
histBinSubValue(Histogram &hist, double value)
{
assert(hist.step > 0);
const int bin = cdo::min((int) ((value - hist.min) / hist.step), hist.nbins - 1);
if (bin >= 0 && bin < hist.nbins && INT_PTR(hist.ptr)[bin] > 0) INT_PTR(hist.ptr)[bin]--;
}
static void
histBin(Histogram &hist)
{
......@@ -124,6 +133,50 @@ histAddValue(Histogram &hist, double value)
return 0;
}
static void
histRemoveValue(Histogram &hist, double value)
{
int i = 0;
for ( i = 0; i < hist.nsamp; i++ )
{
if ( IS_EQUAL(DBL_PTR(hist.ptr)[i], value) )
{
if ( i != hist.nsamp-1 )
DBL_PTR(hist.ptr)[i] = DBL_PTR(hist.ptr)[hist.nsamp-1];
break;
}
}
if ( i == hist.nsamp )
cdoWarning("'%f' not found in histogram!", value);
else
hist.nsamp--;
}
static int
histSubValue(Histogram &hist, double value)
{
assert(hist.nbins > 0);
// 2011-08-01 Uwe Schulzweida: added check for rounding errors
if (value < hist.min && (hist.min - value) < 1e5) value = hist.min;
if (value > hist.max && (value - hist.max) < 1e5) value = hist.max;
if (IS_EQUAL(hist.min, hist.max)) return 0;
if (value < hist.min || value > hist.max) return 1;
if (hist.nsamp < DBL_CAPACITY(hist.nbins))
{
histRemoveValue(hist, value);
}
else if (hist.nsamp > DBL_CAPACITY(hist.nbins))
{
histBinSubValue(hist, value);
hist.nsamp--;
}
return 0;
}
static double
histGetPercentile(const Histogram &hist, double p)
{
......@@ -220,7 +273,7 @@ HistogramSet::defVarLevelBounds(int varID, int levelID, const Field &field1, con
}
void
HistogramSet::addVarLevelValues(int varID, int levelID, const Field &field)
HistogramSet::addSubVarLevelValues(int varID, int levelID, const Field &field, bool lisAdd)
{
const auto &array = field.vec;
assert(!array.empty());
......@@ -237,6 +290,8 @@ HistogramSet::addVarLevelValues(int varID, int levelID, const Field &field)
int nign = 0;
if ( lisAdd )
{
if (field.nmiss)
{
for (size_t i = 0; i < nhists; i++)
......@@ -246,10 +301,22 @@ HistogramSet::addVarLevelValues(int varID, int levelID, const Field &field)
{
for (size_t i = 0; i < nhists; i++) nign += histAddValue(hists[i], array[i]);
}
}
else
{
if (field.nmiss)
{
for (size_t i = 0; i < nhists; i++)
if (!DBL_IS_EQUAL(array[i], field.missval)) nign += histSubValue(hists[i], array[i]);
}
else
{
for (size_t i = 0; i < nhists; i++) nign += histSubValue(hists[i], array[i]);
}
}
if (nign) cdoWarning("%d out of %d grid values are out of bounds and have been ignored (%s)", nign, nhists, __func__);
}
void
HistogramSet::Reset(int varID, int levelID)
{
......
......@@ -79,7 +79,7 @@ public:
void createVarLevels(int varID, int nlevels, size_t nhists);
void defVarLevelBounds(int varID, int levelID, const Field &field1, const Field &field2);
void addVarLevelValues(int varID, int levelID, const Field &field);
void addSubVarLevelValues(int varID, int levelID, const Field &field, bool lisAdd);
void getVarLevelPercentiles(Field &field, int varID, int levelID, double p);
void Reset(int varID, int levelID);
};
......
......@@ -539,7 +539,7 @@ static void
printZaxisReferenceInfo(const int zaxisID)
{
int number = 0;
cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_NUMBEROFVGRIDUSED, &number);
/* cdiInqKeyInt(zaxisID, CDI_GLOBAL, CDI_KEY_NUMBEROFVGRIDUSED, &number)*/
if (number > 0)
{
fprintf(stdout, "%33s : ", "zaxis");
......
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