"Fossies" - the Fresh Open Source Software Archive

Member "fimex-1.4.1/src/Units.cc" (30 Oct 2019, 11033 Bytes) of package /linux/privat/fimex-1.4.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "Units.cc" see the Fossies "Dox" file reference documentation.

    1 /*
    2  * Fimex
    3  *
    4  * (C) Copyright 2008, met.no
    5  *
    6  * Project Info:  https://wiki.met.no/fimex/start
    7  *
    8  * This library is free software; you can redistribute it and/or modify it
    9  * under the terms of the GNU Lesser General Public License as published by
   10  * the Free Software Foundation; either version 2.1 of the License, or
   11  * (at your option) any later version.
   12  *
   13  * This library is distributed in the hope that it will be useful, but
   14  * WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
   15  * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
   16  * License for more details.
   17  *
   18  * You should have received a copy of the GNU Lesser General Public
   19  * License along with this library; if not, write to the Free Software
   20  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301,
   21  * USA.
   22  */
   23 
   24 #include "fimex/Units.h"
   25 
   26 #include "fimex/Logger.h"
   27 #include "fimex/UnitsConverter.h"
   28 #include "fimex/UnitsException.h"
   29 
   30 #include <memory>
   31 
   32 #include "MutexLock.h"
   33 #include "fimex_config.h"
   34 #ifdef HAVE_UDUNITS2_H
   35 #include "udunits2.h"
   36 #include "converter.h"
   37 #else // UDUNITS1
   38 extern "C" {
   39 #include "udunits.h"
   40 }
   41 // add forgotten utIsInit fom udunits
   42 extern "C" int utIsInit();
   43 #endif // UDUNITS2
   44 
   45 #include <cmath>
   46 
   47 namespace MetNoFimex
   48 {
   49 #ifdef HAVE_UDUNITS2_H
   50 static ut_system* utSystem;
   51 #endif
   52 
   53 static OmpMutex unitsMutex;
   54 extern OmpMutex& getUnitsMutex()
   55 {
   56     return unitsMutex;
   57 }
   58 
   59 static Logger_p logger = getLogger("fimex.Units");
   60 
   61 void handleUdUnitError(int unitErrCode, const std::string& message)
   62 {
   63     switch (unitErrCode) {
   64 #ifdef HAVE_UDUNITS2_H
   65     case UT_SUCCESS: break;
   66     case UT_BAD_ARG:  throw UnitException("An argument violates the function's contract: " + message);
   67     case UT_EXISTS: throw UnitException("Unit, prefix, or identifier already exists: " + message);
   68     case UT_NO_UNIT: throw UnitException("No such unit exists: " + message);
   69     case UT_OS: throw UnitException("Operating-system error: " + message);
   70     case UT_NOT_SAME_SYSTEM: throw UnitException("The units belong to different unit-systems: " + message);
   71     case UT_MEANINGLESS: throw UnitException("The operation on the unit(s) is meaningless: " + message);
   72     case UT_NO_SECOND: throw UnitException("The unit-system doesn't have a unit named 'second': " + message);
   73     case UT_VISIT_ERROR: throw UnitException("An error occurred while visiting a unit: " + message);
   74     case UT_CANT_FORMAT: throw UnitException("A unit can't be formatted in the desired manner: " + message);
   75     case UT_SYNTAX: throw UnitException("string unit representation contains syntax error: " + message);
   76     case UT_UNKNOWN: throw UnitException("string unit representation contains unknown word: " + message);
   77     case UT_OPEN_ARG: throw UnitException("Can't open argument-specified unit database: " + message);
   78     case UT_OPEN_ENV: throw UnitException("Can't open environment-specified unit database: " + message);
   79     case UT_OPEN_DEFAULT: throw UnitException("Can't open installed, default, unit database: " + message);
   80     case UT_PARSE: throw UnitException("Error parsing unit specification: " + message);
   81 #else // udunits1
   82     case 0: break;
   83     case UT_EOF: throw UnitException("end-of-file encountered : " + message);
   84     case UT_ENOFILE: throw UnitException("no units-file : " + message);
   85     case UT_ESYNTAX: throw UnitException("syntax error : " + message);
   86     case UT_EUNKNOWN: throw UnitException("unknown specification : " + message);
   87     case UT_EIO: throw UnitException("I/O error : " + message);
   88     case UT_EINVALID: throw UnitException("invalid unit-structure : " + message);
   89     case UT_ENOINIT: throw UnitException("package not initialized : " + message);
   90     case UT_ECONVERT: throw UnitException("two units are not convertable : " + message);
   91     case UT_EALLOC: throw UnitException("memory allocation failure : " + message);
   92     case UT_ENOROOM: throw UnitException("insufficient room supplied : " + message);
   93     case UT_ENOTTIME: throw UnitException("not a unit of time : " + message);
   94     case UT_DUP: throw UnitException("duplicate unit : " + message);
   95 #endif // UDUNITS2
   96     default: throw UnitException("unknown error");
   97     }
   98 }
   99 
  100 class LinearUnitsConverter : public UnitsConverter{
  101     double dscale_;
  102     double doffset_;
  103     double fscale_;
  104     double foffset_;
  105 
  106 public:
  107     LinearUnitsConverter(double scale, double offset)
  108         : dscale_(scale)
  109         , doffset_(offset)
  110         , fscale_(scale)
  111         , foffset_(offset)
  112     {
  113     }
  114     ~LinearUnitsConverter() {}
  115     double convert(double from) override { return dscale_ * from + doffset_; }
  116     float convert(float from) override { return fscale_ * from + foffset_; }
  117     bool isLinear() override { return true; }
  118     void getScaleOffset(double& scale, double& offset) override
  119     {
  120         scale = dscale_;
  121         offset = doffset_;
  122     }
  123 };
  124 
  125 #ifdef HAVE_UDUNITS2_H
  126 class Ud2UnitsConverter : public UnitsConverter {
  127     cv_converter* conv_;
  128 public:
  129     Ud2UnitsConverter(cv_converter* conv) : conv_(conv) {}
  130     ~Ud2UnitsConverter() { cv_free(conv_); }
  131     double convert(double from) override
  132     {
  133         double retval;
  134 #pragma omp critical (cv_converter)
  135         {
  136             retval = cv_convert_double(conv_, from);
  137         }
  138         return retval;
  139     }
  140     float convert(float from) override
  141     {
  142         float retval;
  143 #pragma omp critical(cv_converter)
  144         {
  145             retval = cv_convert_float(conv_, from);
  146         }
  147         return retval;
  148     }
  149     bool isLinear() override
  150     {
  151         // check some points
  152         double offset = convert(0.0);
  153         if (!std::isfinite(offset)) return false;
  154         double slope = convert(1.0) - offset;
  155         if (!std::isfinite(slope)) return false;
  156 
  157         double val = 10.;
  158         double cval = convert(val);
  159         if ((!std::isfinite(cval)) || (std::fabs(cval - (val*slope+offset)) > 1e-5)) return false;
  160         val = 100.;
  161         cval = convert(val);
  162         if ((!std::isfinite(cval)) || (std::fabs(cval - (val*slope+offset)) > 1e-5)) return false;
  163         val = 1000.;
  164         cval = convert(val);
  165         if ((!std::isfinite(cval)) || (std::fabs(cval - (val*slope+offset)) > 1e-5)) return false;
  166         val = -1000.;
  167         cval = convert(val);
  168         if ((!std::isfinite(cval)) || (std::fabs(cval - (val*slope+offset)) > 1e-5)) return false;
  169 
  170         return true;
  171     }
  172     void getScaleOffset(double& scale, double& offset) override
  173     {
  174         if (! isLinear()) throw UnitException("cannot get scale and offset of non-linear function");
  175         offset = convert(0.0);
  176         scale = convert(1.) - offset;
  177         // make sure, scale is determined at a place where offset is no longer numerically superior
  178         if (scale != 0 && offset != 0) {
  179             double temp = -offset/scale;
  180             //std::cerr << temp << " offset:" << offset << " scale: " << scale <<  std::endl;
  181             double scale2 = (convert(temp)-offset)/temp;
  182             if (fabs(scale - scale2) < 1e-3) {
  183                 // should only increase precision, not nan/inf or other cases
  184                 scale = scale2;
  185             }
  186         }
  187     }
  188 };
  189 #endif
  190 
  191 Units::Units()
  192 {
  193     OmpScopedLock lock(unitsMutex);
  194 #ifdef HAVE_UDUNITS2_H
  195     if (utSystem == 0) {
  196         ut_set_error_message_handler(&ut_ignore);
  197         utSystem = ut_read_xml(0);
  198         handleUdUnitError(ut_get_status());
  199     }
  200 #else
  201     if (!utIsInit()) {
  202         handleUdUnitError(utInit(0));
  203     }
  204 #endif
  205 }
  206 
  207 Units::Units(const Units& u)
  208 {}
  209 
  210 Units& Units::operator=(const Units& rhs)
  211 {
  212     return *this; // no state! no increase/decrease to counter required
  213 }
  214 
  215 Units::~Units()
  216 {}
  217 
  218 bool Units::unload(bool force)
  219 {
  220     bool retVal = false;
  221     if (force) {
  222         OmpScopedLock lock(unitsMutex);
  223 #ifdef HAVE_UDUNITS2_H
  224         ut_free_system(utSystem);
  225 #else
  226         utTerm();
  227 #endif
  228         retVal = true;
  229     }
  230 
  231     return retVal;
  232 }
  233 
  234 
  235 void Units::convert(const std::string& from, const std::string& to, double& slope, double& offset)
  236 {
  237     LOG4FIMEX(logger, Logger::DEBUG, "convert from '" << from << "' to '" << to << "'");
  238     UnitsConverter_p conv = getConverter(from, to);
  239     conv->getScaleOffset(slope, offset);
  240 }
  241 
  242 UnitsConverter_p Units::getConverter(const std::string& from, const std::string& to)
  243 {
  244     LOG4FIMEX(logger, Logger::DEBUG, "getConverter from '" << from << "' to '" << to << "'");
  245     if (from == to) {
  246         return std::make_shared<LinearUnitsConverter>(1., 0.);
  247     }
  248     OmpScopedLock lock(unitsMutex);
  249 #ifdef HAVE_UDUNITS2_H
  250     std::shared_ptr<ut_unit> fromUnit(ut_parse(utSystem, from.c_str(), UT_UTF8), ut_free);
  251     handleUdUnitError(ut_get_status(), "'" + from + "'");
  252     std::shared_ptr<ut_unit> toUnit(ut_parse(utSystem, to.c_str(), UT_UTF8), ut_free);
  253     handleUdUnitError(ut_get_status(), "'" + to + "'");
  254     cv_converter* conv = ut_get_converter(fromUnit.get(), toUnit.get());
  255     handleUdUnitError(ut_get_status(), "'" + from + "' converted to '" +to + "'");
  256     return std::make_shared<Ud2UnitsConverter>(conv);
  257 #else
  258     double slope, offset;
  259     utUnit fromUnit, toUnit;
  260     handleUdUnitError(utScan(from.c_str(), &fromUnit), from);
  261     handleUdUnitError(utScan(to.c_str(), &toUnit), to);
  262     handleUdUnitError(utConvert(&fromUnit, &toUnit, &slope, &offset));
  263     return std::make_shared<LinearUnitsConverter>(slope, offset);
  264 #endif
  265 
  266 }
  267 
  268 bool Units::areConvertible(const std::string& unit1, const std::string& unit2) const
  269 {
  270     LOG4FIMEX(logger, Logger::DEBUG, "test convertibility of " << unit1 << " to " << unit2);
  271     int areConv = 0;
  272 #ifdef HAVE_UDUNITS2_H
  273     try {
  274         OmpScopedLock lock(unitsMutex);
  275         std::shared_ptr<ut_unit> fromUnit(ut_parse(utSystem, unit1.c_str(), UT_UTF8), ut_free);
  276         handleUdUnitError(ut_get_status(), "'" + unit1 + "'");
  277         std::shared_ptr<ut_unit> toUnit(ut_parse(utSystem, unit2.c_str(), UT_UTF8), ut_free);
  278         handleUdUnitError(ut_get_status(), "'" + unit2 + "'");
  279         areConv = ut_are_convertible(fromUnit.get(), toUnit.get());
  280     } catch (UnitException& ue) {
  281         LOG4FIMEX(logger, Logger::WARN, ue.what());
  282     }
  283 #else
  284     ScopedLock lock(unitsMutex);
  285     utUnit fromUnit, toUnit;
  286     double slope, offset;
  287     handleUdUnitError(utScan(unit1.c_str(), &fromUnit), unit1);
  288     handleUdUnitError(utScan(unit2.c_str(), &toUnit), unit2);
  289     int error = utConvert(&fromUnit, &toUnit, &slope, &offset);
  290     switch (error) {
  291     case 0: areConv = 1; break;
  292     case UT_ECONVERT: areConv = 0; break;
  293     default: handleUdUnitError(error);
  294     }
  295 #endif
  296 
  297     return areConv;
  298 }
  299 bool Units::isTime(const std::string& timeUnit) const
  300 {
  301 #ifdef HAVE_UDUNITS2_H
  302     return areConvertible(timeUnit, "seconds since 1970-01-01 00:00:00");
  303 #else
  304     bool isTime = false;
  305     class ScopedCritical lock(unitsMutex);
  306     utUnit unit;
  307     handleUdUnitError(utScan(timeUnit.c_str(), &unit), timeUnit);
  308     isTime = (utIsTime(&unit) != 0);
  309     return isTime;
  310 #endif
  311 }
  312 
  313 const void* Units::exposeInternals() const {
  314 #ifdef HAVE_UDUNITS2_H
  315     return utSystem;
  316 #else
  317     return 0;
  318 #endif
  319 }
  320 
  321 }