// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* // ** Copyright UCAR (c) 1990 - 2016 // ** University Corporation for Atmospheric Research (UCAR) // ** National Center for Atmospheric Research (NCAR) // ** Boulder, Colorado, USA // ** BSD licence applies - redistribution and use in source and binary // ** forms, with or without modification, are permitted provided that // ** the following conditions are met: // ** 1) If the software is modified to produce derivative works, // ** such modified software should be clearly marked, so as not // ** to confuse it with the version available from UCAR. // ** 2) Redistributions of source code must retain the above copyright // ** notice, this list of conditions and the following disclaimer. // ** 3) Redistributions in binary form must reproduce the above copyright // ** notice, this list of conditions and the following disclaimer in the // ** documentation and/or other materials provided with the distribution. // ** 4) Neither the name of UCAR nor the names of its contributors, // ** if any, may be used to endorse or promote products derived from // ** this software without specific prior written permission. // ** DISCLAIMER: THIS SOFTWARE IS PROVIDED "AS IS" AND WITHOUT ANY EXPRESS // ** OR IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED // ** WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. // *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=* /** * @file ThresholdBiasMapping.cc */ //------------------------------------------------------------------ #include #include #include #include #include #include #include #include #include //------------------------------------------------------------------ static int _hms(const time_t &t) { DateTime dt(t); int y, m, d; y = dt.getYear(); m = dt.getMonth(); d = dt.getDay(); DateTime dt0(y, m, d, 0, 0, 0); time_t t0 = dt0.utime(); int hms = t - t0; return hms; } //------------------------------------------------------------------ static int _bestGuessHourMinSec(const time_t >, int hourMinSec0, int deltaSec) { // want the gen that is >= one of the target gens, and < the next one. int hms = _hms(gt); // so what is the biggest hms in a day? int k=0; while (hourMinSec0 + k*deltaSec < 24*3600) { k++; } k--; if (hms < hourMinSec0 + deltaSec) { // go with smallest one return hourMinSec0; } else if (hms >= hourMinSec0 + k*deltaSec) { return k*deltaSec; } else { for (int i=1; i= hourMinSec0 + i*deltaSec && hms < hourMinSec0 + (i+1)*deltaSec) { return i*deltaSec; } } LOG(ERROR) << "Did not expect to be here"; return -1; } } //------------------------------------------------------------------ static bool _bestGuessResolution(const std::vector ×, int &hourMinSec0, int &deltaSec) { hourMinSec0 = -1; deltaSec = -1; if (times.empty()) { return false; } else if (times.size() == 1) { // not setup to handle this case LOG(ERROR) << "One gen time in SPDB, not able to determine gen resolution"; return false; } int numAtSame = 0; for (int i=0; i<(int)(times.size())-1; ++i) { int deltai = times[i+1] - times[i]; if (i == 0) { deltaSec = deltai; numAtSame = 1; } else { if (deltai != deltaSec) { if (deltai < deltaSec) { if (numAtSame > 1) { LOG(WARNING) << "Inconsistent SPDB time delta " << deltai << " " << deltaSec; } deltaSec = deltai; numAtSame = 1; } } else { numAtSame++; } } } // now that we have a delta, get an 'anchor' to get the hour/min/sec // values within each day. Use times[0] int hms = _hms(times[0]); hourMinSec0 = hms; while (hourMinSec0 >= 0) { hourMinSec0 -= deltaSec; } hourMinSec0 += deltaSec; // this is the earliest hour/min/sec in each day return true; } //------------------------------------------------------------------ ThresholdBiasMapping::ThresholdBiasMapping(void) { } //------------------------------------------------------------------ ThresholdBiasMapping::ThresholdBiasMapping(const std::string &path) { _url = path; } //------------------------------------------------------------------ ThresholdBiasMapping::~ThresholdBiasMapping() { } //------------------------------------------------------------------ void ThresholdBiasMapping::store(int lt, double bias, double w) { _biasMap[lt] = bias; _threshMap[lt] = w; } //------------------------------------------------------------------ void ThresholdBiasMapping::writeAndClear(const time_t >) { string xml = _toXml(); DsSpdb s; s.setPutMode(Spdb::putModeOver); s.clearPutChunks(); s.clearUrls(); s.addUrl(_url); MemBuf mem; mem.free(); mem.add(xml.c_str(), xml.size() + 1); if (s.put(SPDB_XML_ID, SPDB_XML_LABEL, 1, gt, gt, mem.getLen(), (void *)mem.getPtr())) { LOG(ERROR) << "problems writing out SPDB"; } } //------------------------------------------------------------------ bool ThresholdBiasMapping::readFirstBefore(const time_t >, int maxSecondsBack) { DsSpdb s; _threshMap.clear(); _biasMap.clear(); if (s.getFirstBefore(_url, gt-1, maxSecondsBack) != 0) { LOG(DEBUG) << "No SPDB data found in data base " << _url << " before " << DateTime::strn(gt) << " , within " << maxSecondsBack << " seconds"; return false; } else { return _load(s); } } //------------------------------------------------------------------ bool ThresholdBiasMapping::readFirstAfter(const time_t >, int maxSecondsAhead) { DsSpdb s; _threshMap.clear(); _biasMap.clear(); if (s.getFirstAfter(_url, gt, maxSecondsAhead) != 0) { LOG(DEBUG) << "No SPDB data found in data base " << _url << " after " << DateTime::strn(gt) << " , within " << maxSecondsAhead << " seconds"; return false; } else { return _load(s); } } //------------------------------------------------------------------ bool ThresholdBiasMapping::readNearest(const time_t >, int maxSecondsBack, int maxSecondsAhead) { // want gen time in nearest day that is >= gt and < next gen time in // that nearest day DsSpdb S; if (S.compileTimeList(_url, gt - maxSecondsBack, gt + maxSecondsAhead) != 0) { LOG(ERROR) << "No SPDB data in data base " << _url << " from " << DateTime::strn(gt - maxSecondsBack) << " to " << DateTime::strn(gt + maxSecondsAhead); return false; } vector times = S.getTimeList(); if (times.empty()) { LOG(ERROR) << "No SPDB data in data base " << _url << " from " << DateTime::strn(gt - maxSecondsBack) << " to " << DateTime::strn(gt + maxSecondsAhead); return false; } bool first = true; time_t best; if (times.size() == 1) { // this single time is it. first = false; best = times[0]; } else { // size > 1 // try to figure out the pattern of gen times (resolution, first gen in day) int hourMinSec0, deltaSec; if (!_bestGuessResolution(times, hourMinSec0, deltaSec)) { return false; } // now figure out which gen to search for int hms= _bestGuessHourMinSec(gt, hourMinSec0, deltaSec); if (hms < 0) { return false; } // now find closest day to day of gt with this hour/min/sec double delta; for (size_t i=0; i::const_iterator i; for (i=_threshMap.begin(); i!=_threshMap.end(); ++i) { char buf[1000]; sprintf(buf, "setenv THRESH_%d %10.6lf\n", i->first, i->second); ret += buf; } return ret; } //------------------------------------------------------------------ bool ThresholdBiasMapping::getThreshold(int lt, double &thresh) const { map::const_iterator i; i = _threshMap.find(lt); if (i == _threshMap.end()) // if (_threshMap.find(lt) == _threshMap.end() ) { LOG(ERROR) << "Lead time " << lt << " not in map"; return false; } else { thresh = i->second; // thresh = _threshMap.at(lt); return true; } } //------------------------------------------------------------------ bool ThresholdBiasMapping::getLeadtimeInterpolatedThreshold(int lt, double &thresh) const { if (_threshMap.empty()) { return false; } // pull out all the lead times map::const_iterator i; vector locLt; for (i=_threshMap.begin(); i!=_threshMap.end(); ++i) { if (lt == i->first) { thresh = i->second; return true; } locLt.push_back(i->first); } // now see where our lt fits in, since it was never an exact match // here we assume the leads are in ascending order, which is guaranteed // by the map class, so should be ok if (lt < locLt[0]) { LOG(DEBUG_VERBOSE) << "Lead time " << lt << " smaller than all leads in SPDB, use smallest lead times threshold"; return getThreshold(locLt[0], thresh); } if (lt > *(locLt.rbegin())) { LOG(DEBUG_VERBOSE) << "Lead time " << lt << " larger than all leads in SPDB, use largest lead times threshold"; return getThreshold(*(locLt.rbegin()), thresh); } for (int i=1; i(locLt.size()); ++i) { if (lt > locLt[i-1] && lt < locLt[i]) { // interpolate now double w0 = locLt[i] - lt; double w1 = lt - locLt[i-1]; double w = locLt[i] - locLt[i-1]; w0 /= w; w1 /= w; LOG(DEBUG_VERBOSE) << "Interpolating lead times " << locLt[i-1] << " " << locLt[i]; LOG(DEBUG_VERBOSE) << "Interpolating weights are " << w0 << " " << w1; double t0, t1; if (!getThreshold(locLt[i-1], t0)) { LOG(ERROR) << "Unexpected inability to get threshold at " << locLt[i-1]; return false; } if (!getThreshold(locLt[i], t1)) { LOG(ERROR) << "Unexpected inability to get threshold at " << locLt[i]; return false; } LOG(DEBUG_VERBOSE) << "Interpolating thresholds are " << t0 << " " << t1; thresh = (w0*t0 + w1*t1); return true; } } LOG(ERROR) << "Unexpected lead time interpolation failed, no window"; return false; } //------------------------------------------------------------------ bool ThresholdBiasMapping::getPartitionedThreshold(int lt, double &thresh) const { // produce the partition by pulling out all the lead times map::const_iterator i; vector locLt; for (i=_threshMap.begin(); i!=_threshMap.end(); ++i) { locLt.push_back(i->first); } int dt0 = 0; for (int j=0; j(locLt.size()); ++j) { if (j > 0) { int dt = locLt[j] - locLt[j-1]; if (j == 1) { dt0 = dt; } else { if (dt != dt0) { LOG(ERROR) << "Lead time partitioning uneven " << dt << " " << dt0; return false; } } } } // so here the partitions are [locLt[0]+i*dt0, locLt[0]+(i+1)*dt0), // i=0,1,...locLt.size()-1 // so want lt >= locLt[0] + i*dt0 // lt < locLt[0] + (i+1)*dt0 // I believe this works with truncation int j = (lt - locLt[0])/dt0; if (j >= static_cast(locLt.size())) { j = static_cast(locLt.size()) - 1; } if (j < 0) { j = 0; } int partitionLt = locLt[j]; LOG(DEBUG) << "Partition chosen = [" << partitionLt << "," << partitionLt+dt0 << ")"; return getThreshold(partitionLt, thresh); } //------------------------------------------------------------------ bool ThresholdBiasMapping::getBias(int lt, double &bias) const { map::const_iterator i; i = _biasMap.find(lt); // if (_biasMap.find(lt) == _biasMap.end() ) if (i == _biasMap.end()) { LOG(ERROR) << "Lead time " << lt << " not in map"; return false; } else { bias = i->second; // bias = _biasMap.at(lt); return true; } } //------------------------------------------------------------------ bool ThresholdBiasMapping::compare(const ThresholdBiasMapping &inp, double maxDiff, std::string &comparisons) const { bool status = true; std::map::const_iterator i, iInp; comparisons = ""; for (i = _threshMap.begin(); i != _threshMap.end(); ++i) { int lt = i->first; double thresh0 = i->second; iInp = inp._threshMap.find(lt); if (iInp == inp._threshMap.end()) // if (inp._threshMap.find(lt) == inp._threshMap.end() ) { LOG(ERROR) << "Lead time " << lt << " not in input map"; status = false; } else { // double thresh1 = inp._threshMap.at(lt); double thresh1 = iInp->second; char buf[1000]; sprintf(buf, "%6d %6.4lf %6.4lf diff:%6.4lf\n", lt, thresh0, thresh1, thresh1-thresh0); comparisons += buf; if (fabs(thresh1-thresh0) > maxDiff) { status = false; } } } for (i = inp._threshMap.begin(); i != inp._threshMap.end(); ++i) { int lt = i->first; if (_threshMap.find(lt) == _threshMap.end() ) { LOG(ERROR) << "Lead time " << lt << " not in local map"; status = false; } } return status; } //------------------------------------------------------------------ bool ThresholdBiasMapping::biasGood(double target, double tolerance) const { std::map::const_iterator i; for (i=_biasMap.begin(); i!=_biasMap.end(); ++i) { if (fabs(i->second - target) >= tolerance) { LOG(DEBUG) << "bias=" << i->second <<" at lead=" << i->first << " outside tolerance"; return false; } } LOG(DEBUG) << "All bias values within " << tolerance << " of " << target; return true; } //------------------------------------------------------------------ bool ThresholdBiasMapping::biasGood(double target, double tolerance, std::vector leadSeconds) const { if (leadSeconds.empty()) { return biasGood(target, tolerance); } std::map::const_iterator i; for (i=_biasMap.begin(); i!=_biasMap.end(); ++i) { if (find(leadSeconds.begin(), leadSeconds.end(), i->first) != leadSeconds.end()) { if (fabs(i->second - target) >= tolerance) { LOG(DEBUG) << "bias=" << i->second <<" at lead=" << i->first << " outside tolerance"; return false; } } } LOG(DEBUG) << "All bias values within " << tolerance << " of " << target; return true; } //------------------------------------------------------------------ void ThresholdBiasMapping::copyData(const ThresholdBiasMapping &inp) { _threshMap = inp._threshMap; _biasMap = inp._biasMap; } //------------------------------------------------------------------ std::string ThresholdBiasMapping::asciiTable(void) const { std::string ret = ""; // this assumes the same lead times for both mappings std::map::const_iterator it, ib; for (it=_threshMap.begin(), ib=_biasMap.begin(); it!=_threshMap.end() && ib != _biasMap.end(); ++it, ++ib) { int lt = it->first; int lt2 = ib->first; if (lt != lt2) { LOG(ERROR) << "Lead times do not match " << lt << " " << lt2; } else { double thresh = it->second; double bias = ib->second; char buf[1000]; sprintf(buf, "%6d %3.2lf %3.2lf\n", lt, thresh, bias); ret += buf; } } return ret; } //------------------------------------------------------------------ std::string ThresholdBiasMapping::asciiTable(const time_t &t) const { std::string stime = DateTime::strs(t); std::string ret = ""; // this assumes the same lead times for both mappings std::map::const_iterator it, ib; for (it=_threshMap.begin(), ib=_biasMap.begin(); it!=_threshMap.end() && ib != _biasMap.end(); ++it, ++ib) { int lt = it->first; int lt2 = ib->first; if (lt != lt2) { LOG(ERROR) << "Lead times do not match " << lt << " " << lt2; } else { double thresh = it->second; double bias = ib->second; char buf[1000]; sprintf(buf, "%s %6d %3.2lf %3.2lf\n", stime.c_str(), lt, thresh, bias); ret += buf; } } return ret; } //------------------------------------------------------------------ std::string ThresholdBiasMapping::_toXml(void) const { map::const_iterator i; string s = TaXml::writeStartTag("Thresholds", 0); for (i=_threshMap.begin(); i!=_threshMap.end(); ++i) { int lt = i->first; double t = i->second; s += _appendXml("ThreshMap", "Thresh", lt, t); } s += TaXml::writeEndTag("Thresholds", 0); s += TaXml::writeStartTag("WeightedBias", 0); for (i=_biasMap.begin(); i!=_biasMap.end(); ++i) { int lt = i->first; double t = i->second; s += _appendXml("BiasMap", "Bias", lt, t); } s += TaXml::writeEndTag("WeightedBias", 0); return s; } //------------------------------------------------------------------ std::string ThresholdBiasMapping::_appendXml(const std::string &pairKey, const std::string &valueKey, int lt, double w) const { string s; s = TaXml::writeStartTag(pairKey, 0); s = s.substr(0, s.size()-1); s += TaXml::writeInt("Lead", 0, lt); s = s.substr(0, s.size()-1); s += TaXml::writeDouble(valueKey, 0, w, "%.10lf"); s = s.substr(0, s.size()-1); s += TaXml::writeEndTag(pairKey, 0); return s; } //------------------------------------------------------------------ bool ThresholdBiasMapping::_fromXml(const std::string &xml) { _threshMap.clear(); _biasMap.clear(); string s; if (TaXml::readString(xml, "Thresholds", s)) { LOG(ERROR) << "No XML with key Thresholds"; return false; } vector m; if (TaXml::readStringArray(s, "ThreshMap", m)) { LOG(ERROR) << "No XML array with key ThreshMap"; return false; } for (size_t i=0; i &chunks = s.getChunks(); int nChunks = static_cast(chunks.size()); if (nChunks <= 0) { //LOG(WARNING) << "No chunks"; return false; } if (nChunks > 1) { LOG(WARNING) << "Too many chunks " << nChunks << " expected 1"; return false; } vector lchunks; for (int i=0; i