ICalculator.hpp

#pragma once
 
#include "DM/IGeometry.hpp"
 
DM_NAMESPACE_BEGIN
 
/// \brief Generic calculator for evaluating formulas or manipulating objects based on the OPALS \ref ref_filter_generic "generic filter syntax"
class DM_API ICalculator : public ObjectBase
{
public:
 
  /// Choose the supported input data types.
  enum struct ReadAccess
  {
    none,                                    ///< support querying no data types at all. Hence, only constants and functions: pi atan(.23)
    coordinates,                             ///< support querying coordinates: x y z
    attributes   = coordinates << 1,         ///< support querying attributes: sigmaX _customAttr
    vectors      = coordinates | attributes, ///< shortcut: support querying coordinates or attributes
    rasters      = attributes << 1,          ///< support querying input raster data: r[0]
    neighbors    = rasters    << 1,          ///< support querying neighbors: n[0]
    full         = (neighbors << 1) - 1      ///< support querying all input data types
  };
 
  /// Choose which types of data may be assigned to.
  enum struct WriteAccess
  {
    none,                              ///< support no assignments at all
    coordinates,                       ///< support assignment to coordinates
    attributes  = coordinates << 1,    ///< support assignment to attributes
    rasters     = attributes  << 1,    ///< support assignment to output raster nodes: [0] [red] [](float32)
    full        = (rasters  << 1) - 1  ///< support assignment to all output data types
  };
 
  /// The result type of the calculator.
  /** If custom attributes are to be queried, the result type is unknown beforehand.*/
  enum struct DataTypeSuperset
  {
    arithmetic, string, unknown
  };
 
  /// Copy of GDALDataType, but without complex types.
  enum struct RasterDataType
  {
    unknown, ///< Unknown or unspecified type
    uint8,   ///< 8 bit unsigned integer
    uint16 , ///< 16 bit unsigned integer
    int16  , ///< 16 bit signed integer
    uint32 , ///< 32 bit unsigned integer
    int32  , ///< 32 bit signed integer
    float32, ///< 32 bit floating point
    float64  ///< 64 bit floating point
  };
 
  // Round \op value to \p type. Sets \p value to std::numeric_limits<double>::quiet_NaN() if it is outside the domain of \p type.
  static void convertOrNaN(double &value, RasterDataType type);
 
  static ICalculator* New();
 
  /// Construct a calculator.
  /// \param text Text that defines the calculator.
  /// \param readAccess  Choose which types of data the calculator may query i.e. read from.
  ///                    If supported, make sure to provide inRaster and/or neighbors in calls to resultDouble()!
  ///                    Throws if \p text would read from an unsupported type.
  /// \param writeAccess Choose which types of data the calculator may write to i.e. assign to.
  ///                    If supported, make sure to provide outRaster in calls to resultDouble()!
  ///                    Throws if \p text would assign to an unsupported type.
  static ICalculator* New(const char *text, ReadAccess readAccess, WriteAccess writeAccess);
 
  /// Deep-copy. Use for multi-threading!
  virtual ICalculator* clone() const = 0;
 
  virtual ICalculator* operator*(double factor) const = 0;
  virtual ICalculator* operator*(const ICalculator& that) const = 0;
 
  /// \name Logging of evaluation errors
  /// Set a functor to be called upon evaluation errors, e.g. division by zero, conversion failure to target type.
  /// \p msg describes the problem, and it may contain contextual information. Hence, similar problems may have different messages.
  ///@{
 
  /// Abstract functor to derive from.
  struct EvaluationErrorCallback
  {
    virtual void operator()(const char *msg) = 0;
  };
 
  /// Make this use your functor instance as callback. Its memory remains managed by you.
  virtual void evaluationErrorCallback(EvaluationErrorCallback *cb = nullptr) = 0;
 
  ///@}
 
  /// Get notified if an evaluated geometry has been changed.
  /** After each call of resultDouble, this flag will be true only if this has changed the coordinates and/or attributes of \p geom .
      An attribute is considered as changed if:
      - a new AddInfo has been assigned, or
      - its layout has been changed (column of same name but different type) or extended (additional column), or
      - an attribute that was null has been set to non-null, or
      - a non-null attribute has been set null, or
      - the value of a non-null attribute has been changed.
      The memory of \p flag remains managed by you.*/
  virtual void geometryChanged(bool *flag = nullptr) = 0;
 
  /// Initialize internal data. resultDouble calls this implicitly / lazily if some non-const method has been called since the last call to this.
  /** This does various checks and may throw.
      Call this explicitly to generate these errors beforehand. */
  virtual void init() const = 0;
 
  /// \name Input rasters
  /// Calculators process zero or more input raster datasets with zero or more bands.
  /// Each band has an invalid value, and it may have a name.
  /// Named bands may be accessed in the \p text by name instead of by index.
  /// While the number of neighbors may change for each call of resultDouble(),
  /// the number of input raster datasets, their band counts, and each band's name and nodata value typically do not change
  /// for many consecutive calls of resultDouble(). Hence, this information must be specified beforehand, if input rasters shall be used.
  ///@{
    
  /// \brief Tells the calculator for how many datasets the argument \p inRaster of resultDouble() will subsequently provide data,
  /// if the respective datasets are used (see inputDatasetSubset()).
  /** Datasets by default have a single band.
      For each dataset, \p inRaster expects an element for each band, sorted first by dataset and second by band index, in ascending order.
      Hence, with
     
          CalculatorHandle calc = ICalculator::New(...);
          calc->inputDatasetCount(2);
     
      resultDouble() will expect \p inRaster to point to 2 elements in total:
      one for the first dataset, followed by another element for the second dataset (since both datasets are single band).*/
  virtual void inputDatasetCount(unsigned nDatasets) = 0;
  /// Tells the calculator that the \p iDataset 'th dataset has \p nBands instead of 1.
  /** Hence, with
     
          CalculatorHandle calc = ICalculator::New(...);
          calc->inputDatasetCount(3);
          calc->inputBandCount(0, 2);
          calc->inputBandCount(2, 4);
     
      resultDouble() will expect \p inRaster to point to 7 elements, where 
      - the first 2 elements are associated with the 2 bands of the first dataset,
      - the third element is associated with the only band of the second dataset (because its band count has been left to its default 1), and
      - the trailing 4 elements are associated with the four bands of the third dataset.*/
  virtual void inputBandCount(unsigned iDataset, unsigned nBands) = 0;
  /// Associates a band index with a band name, so the band can be accessed by name in \p text.
  virtual void inputBandName(unsigned iDataset, unsigned iBand, char const *bandName) = 0;
  /// Associates the given value as invalid for the given input raster band. Default: NaN.
  virtual void inputBandNoData(unsigned iDataset, unsigned iBand, double noData) = 0;
  /// Temporarily restricts the datasets for which \p inRaster must provide elements to those specified in \p iDatasets.
  /** Input raster indices (r[<idx>]) and index ranges (mean(r[<iBegin> : <iEnd>])) then index into this reduced list of datasets. 
      \p iDatasets must point to \p nDatasets elements.
      Input band counts, names, and nodata-values of temporarily de-activated datasets are preserved and will be re-used
      once a dataset gets re-activated by calling inputDatasetSubset another time.
      Hence, with
     
          CalculatorHandle calc = ICalculator::New(...);
          calc->inputDatasetCount(3);
          calc->inputBandCount(0, 2);
          calc->inputBandName(0, 0, "red");
          calc->inputBandCount(2, 4);
          std::vector<unsigned> subset = {1, 2};
          calc->inputDatasetSubset(subset.data(), subset.size());
     
      resultDouble() expects \p inRaster to point to 5 elements, where 
      - the first element is associated with the only band of the second dataset, and
      - the trailing 4 elements are associated with the four bands of the third dataset.
      
      Because the first dataset has been excluded / de-activated, r[0] now refers to the single-band (second) dataset.
      If later on,
     
          subset = {0, 1};
          calc->inputDatasetSubset(subset.data(), subset.size());
     
      gets executed, resultDouble() expects \p inRaster to point to 3 elements, where 
      - the first 2 elements are associated with the 2 bands of the first dataset, and
      - the third element is associated with the only band of the second dataset.
      
      Again, the first band of the first dataset is associated with the name "red".
      Because the first dataset has been activated again, r[0] refers to the two-band (first) dataset.*/
  virtual void inputDatasetSubset(unsigned const *iDatasets = 0, unsigned nDatasets = 0) = 0;
 
  ///@}
 
 
  /// \name Output raster
  /// \p text defines the number of output bands, their names and types.
  /// Use this information to size \p outRaster and create the output raster datasets.
  ///@{
    
  /// Returns the number of elements that \p outRaster must point to.
  /** Hence, returns 0 if this calculator does not write to any output rasters. */
  virtual unsigned outputBandCount() const = 0;
  /// Returns the name of the \p iBand'th output band if specified in \p text, or otherwise an empty string.
  virtual const char * outputBandName(unsigned iBand) const = 0;
  /// Returns the actual type assumed for the \p iBand'th output band. If not specified in \p text, returns the one passed to outputBandTypeDefault, or otherwise RasterDataType::float64.
  virtual RasterDataType outputBandType(unsigned iBand) const = 0;
  /// Sets a value to be written to the \p iBand'th element of \p outRaster if invalid (aka NO_DATA). Default: std::numeric_limits<output band type>::max().
  virtual void outputBandNoData(unsigned iBand, double noData) = 0;
  /// Clears the NO_DATA value for the \p iBand'th element of \p outRaster. Evaluation will throw upon encountering an invalid value!
  virtual void clearOutputBandNoData(unsigned iBand) = 0;
  /// Returns true if a NO_DATA value has been set for the \p iBand'th element of \p outRaster.
  virtual bool hasOutputBandNoData(unsigned iBand) const = 0;
  /// Returns the NO_DATA value for the \p iBand'th element of \p outRaster. Throws if none has been set.
  virtual double outputBandNoData(unsigned iBand) const = 0;
  /// Sets the default data type for output bands, to be used if \p text does not specify one. Default: RasterDataType::float64.
  virtual void outputBandTypeDefault(RasterDataType) = 0;
  /// Query the statistics of values assigned so far to the \p iBand'th output band.
  /** \p nInvalid is the number of times that the value to be assigned has been invalid. 
      \p nNoData is the number of times that the value to be assigned has been valid and could be converted to the band's data type, but happened to be equal to the band's NoData value.
      \p nOverflow is the number of times that the value to be assigned has been valid, but outside the domain (i.e. too small or large)
         of the band's data type or not numeric (i.e. ColumnType::string or ColumnType::cstr).
 
      nInvalid + nNoData + nOverflow is the total number of times that the NoData value has been assigned to \p outRaster.
  */
  virtual void outputBandStats(unsigned iBand, size_t& nInvalid, size_t& nNoData, size_t& nOverflow) const = 0;
 
  ///@}
 
 
  /// \name Attribute assignment statistics
  /// Query the statistics of attributes assigned so far.
  ///@{
 
  /// Query the number of attribute assignment statistics available.
  virtual size_t nAttributeAssignmentStats() const = 0;
  /// Query the statistics of values assigned so far to the \p iAttribute'th attribute.
  /** \p name the name of this attribute
      \p nInvalid is the number of times that the value to be assigned has been invalid.
      \p nOverflow is the number of times that the value to be assigned has been valid, but outside the domain of this attribute's data type,
                   i.e. too small or too large, or numeric vs. string.
      nInvalid + nOverflow is the total number of evaluated geometries for which this attribute is now inexistent
      i.e. the attribute is either not part of the layout, or it is null.
  */
  virtual void attributeAssignmentStats(unsigned iAttribute, char const*& name, size_t& nInvalid, size_t& nOverflow) const = 0;
 
  ///@}
 
 
  /// \name Evaluation
  /// Evaluate this with the given geometry, input rasters, neighbors, and output rasters.
  /// If this has been created without ReadAccess::rasters, then \p inRaster will not be accessed and may be NULL.
  /// Otherwise, preceding calls to methods in the group for input rasters
  /// define the number of elements that \p inRaster must point to, and the dataset and band that each of its elements is associated with.
  /// If this has been created without ReadAccess::neighbors, then \p neighbors will not be accessed and may be NULL.
  /// \p neighbors must point to \p nNeighbors elements.
  /// If this has been created without WriteAccess::rasters, then \p outRaster will not be accessed and may be NULL.
  /// Otherwise, \p outRaster must point to outputBandCount() elements.
  ///
  /// Returns the result of evaluating the last statement.
  /// Assignments always return 1.
  /// The following operations yield invalid values:
  ///    - Access to inexistent or null attributes and invalid coordinates (DBL_MAX).
  ///    - Operations on invalid values.
  ///    - Domain errors, e.g. 1/0, acos(2).
  ///    - Conversion of values outside the value range of target types, e.g. uint(-1.), uint8(256.)
  /// Returns NaN for invalid results.
  /// Hence,
  ///
  ///     x
  ///
  /// returns x if x is valid (i.e. not DBL_MAX, not NaN), NaN otherwise.
  ///
  ///     3 * sigmaX
  ///
  /// returns 3*sigmaX if sigmaX is valid (attribute is part of the layout and not NULL), NaN otherwise.
  ///
  ///     3 * sigmaX; 5
  ///
  /// returns 5.
  ///
  ///     _val = 3 * sigmaX; sigmaY
  ///
  /// Sets _val and returns sigmaY if sigmaY is valid, otherwise NaN.
  /// If init() has been called just before, then resultDouble only
  /// throws on input raster dataset / band index out-of-range, e.g.:
  ///    - r[1] while only a single input raster dataset has been specified, or
  ///    - r[0][1] while only a single band has been specified for the first input raster dataset.
  ///@{
    
  /// Evaluate this with a mutable geometry.
  /** The geometry coordinates and/or attributes may get changed
      only if this has been created with WriteAccess::coordinates and/or WriteAccess::attributes.*/
  virtual double resultDouble(IGeometry &geom,
                              double const *inRaster = 0,
                              IGeometry const *const *neighbors = 0, unsigned nNeighbors = 0,
                              double *outRaster = 0) const = 0;
  /// Evaluate this with a constant geometry.
  /** Throws if this has been created with WriteAccess::coordinates and/or WriteAccess::attributes.*/
  virtual double resultDouble(IGeometry const& geom,
                              double const* inRaster = 0,
                              IGeometry const* const* neighbors = 0, unsigned nNeighbors = 0,
                              double* outRaster = 0) const = 0;
 
  ///@}
 
 
  /// Returns true iff resultDouble surely does not change any arguments and its return value is always the same, independent of its arguments.
  /** Returns true for constants and combinations thereof:
      New("1")->isConstant() == true
      New("3 * atan2(0.2, pi)")->isConstant() == true
      New("2 ? 3 : x")->isConstant() == true
      New("x = 1")->isConstant() == false
      New("3 || 4")->isConstant() == true
      New("0 && 4")->isConstant() == true
 
      Note that the following is correct output, because the validity of x is unknown beforehand and hence, resultDouble may return 0 or NaN.
      New("0 * x")->isConstant() == false
 
      For operator&& and operator||, this returns certain false negatives, e.g.:
      New("0 && x")->isConstant() == false
      New("3 || x")->isConstant() == false */
  virtual bool isConstant() const = 0;
 
  /// returns true iff the calculator reads an input raster.
  virtual bool readsRasters()   const = 0;
  /// returns true iff this reads one or more input raster datasets specified by a single index or by a slice that defines start and/or stop indices.
  virtual bool readsRastersByIndexOrLimitedSlice() const = 0;
  /// returns true iff the calculator reads a neighbor.
  virtual bool readsNeighbors() const = 0;
 
  /// returns the number of (statements that are) assignments.
  virtual unsigned assignmentCount() const = 0;
  /// returns the number of statements. Empty statements do not count.
  virtual unsigned statementCount() const = 0;
 
  virtual DataTypeSuperset typeSuperset() const = 0;
 
  virtual void print(std::ostream&) const = 0;
 
  virtual const char *text() const = 0;
};
 
typedef Handle< ICalculator > CalculatorHandle;
 
constexpr ICalculator::ReadAccess  operator| (const ICalculator::ReadAccess  &left, const ICalculator::ReadAccess  &right)
{
  return static_cast<ICalculator::ReadAccess>(static_cast<int>(left) | static_cast<int>(right));
}
 
constexpr ICalculator::WriteAccess operator| (const ICalculator::WriteAccess &left, const ICalculator::WriteAccess &right)
{
  return static_cast<ICalculator::WriteAccess>(static_cast<int>(left) | static_cast<int>(right));
}
 
DM_NAMESPACE_END