Calculator.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include <cmath>
 
#include <QRegExp>
#include <QStack>
#include <QVector>
 
#include "Calculator.h"
#include "InfixToPostfix.h"
#include "IException.h"
#include "SpecialPixel.h"
 
using namespace std;
 
namespace Isis {
  Calculator::Calculator() {
    p_valStack = NULL;
 
    p_valStack = new QStack< QVector<double> >();
  }
 
  Calculator::~Calculator() {
    if (p_valStack) {
      delete p_valStack;
      p_valStack = NULL;
    }
  }
 
  double NegateOperator(double a) {
    return -1 * a;
  }
 
 
  double MultiplyOperator(double a, double b) {
    return a * b;
  }
 
 
  double DivideOperator(double a, double b) {
    return a / b;
  }
 
 
  double AddOperator(double a, double b) {
    return a + b;
  }
 
 
  double SubtractOperator(double a, double b) {
    return a - b;
  }
 
 
  double GreaterThanOperator(double a, double b) {
    return a > b ? 1.0 : 0.0;
  }
 
 
  double LessThanOperator(double a, double b) {
    return a < b ? 1.0 : 0.0;
  }
 
 
  double EqualOperator(double a, double b) {
    return a == b ? 1.0 : 0.0;
  }
 
 
  double GreaterThanOrEqualOperator(double a, double b) {
    return a >= b ? 1.0 : 0.0;
  }
 
 
  double LessThanOrEqualOperator(double a, double b) {
    return a <= b ? 1.0 : 0.0;
  }
 
 
  double NotEqualOperator(double a, double b) {
    return a != b ? 1.0 : 0.0;
  }
 
 
  double CosecantOperator(double a) {
    return 1.0 / sin(a);
  }
 
 
  double SecantOperator(double a) {
    return 1.0 / cos(a);
  }
 
  double CotangentOperator(double a) {
    return 1.0 / tan(a);
  }
 
 
  int Round(double a) {
    return (a > 0) ? (int)(a + 0.5) : (int)(a - 0.5);
  }
 
 
  double BitwiseAndOperator(double a, double b) {
    return (double)(Round(a)&Round(b));
  }
 
 
  double BitwiseOrOperator(double a, double b) {
    return (double)(Round(a) | Round(b));
  }
 
 
  double ModulusOperator(double a, double b) {
    return (double)(Round(a) % Round(b));
  }
 
  double MaximumOperator(double a, double b) {
 
    if (std::isnan(a)) return (a);
    if (std::isnan(b)) return (b);
    return (a > b) ? a : b;
  }
 
  double MinimumOperator(double a, double b) {
    if (std::isnan(a)) return (a);
    if (std::isnan(b)) return (b);
    return (a < b) ? a : b;
  }
 
 
 
  void Calculator::Negative() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), NegateOperator);
    Push(result);
  }
 
 
  void Calculator::Multiply() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(), MultiplyOperator);
    Push(result);
  }
 
 
  void Calculator::Add() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(), AddOperator);
    Push(result);
  }
 
 
  void Calculator::Subtract() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(), SubtractOperator);
    Push(result);
  }
 
 
  void Calculator::Divide() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(), DivideOperator);
    Push(result);
  }
 
  void Calculator::Modulus() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(), ModulusOperator);
    Push(result);
  }
 
 
  void Calculator::Exponent() {
    QVector<double> exponent = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), exponent.begin(), exponent.end(), pow);
    Push(result);
  }
 
 
  void Calculator::SquareRoot() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), sqrt);
    Push(result);
  }
 
 
  void Calculator::AbsoluteValue() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), fabs);
    Push(result);
  }
 
 
  void Calculator::Log() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), log);
    Push(result);
  }
 
 
  void Calculator::Log10() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), log10);
    Push(result);
  }
 
 
  void Calculator::LeftShift() {
    QVector<double> y = Pop();
    if (y.size() != 1) {
      IString msg = "When trying to do a left shift calculation, a non-scalar "
                    "shift value was encountered. Shifting requires scalars.";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
    else {
      QVector<double> x = Pop();
 
      if ((int)y[0] > (int)x.size()) {
        IString msg = "When trying to do a left shift calculation, a shift "
                      "value greater than the data size was encountered. "
                      "Shifting by this value would erase all of the data.";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
      else {
        QVector<double> result;
        int shift = (int)y[0];
        result.resize(x.size());
 
        for (int i = 0; i < result.size(); i++) {
          if (i + shift < x.size() && i + shift >= 0)
            result[i] = x[i+shift];
          else
            result[i] = sqrt(-1.0); // create a NaN
        }
 
        Push(result);
      }
    }
  }
 
 
  void Calculator::RightShift() {
    QVector<double> y = Pop();
    if (y.size() != 1) {
      IString msg = "When trying to do a right shift calculation, a non-scalar "
                    "shift value was encountered. Shifting requires scalars.";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
    else {
      QVector<double> x = Pop();
 
      if ((int)y[0] > (int)x.size()) {
        IString msg = "When trying to do a right shift calculation, a shift "
                      "value greater than the data size was encountered. "
                      "Shifting by this value would erase all of the data.";
        throw IException(IException::Unknown, msg, _FILEINFO_);
      }
      else {
        QVector<double> result;
        int shift = (int)y[0];
        result.resize(x.size());
 
        for (int i = 0; i < (int)result.size(); i++) {
          if (i - shift < (int)x.size() && i - shift >= 0) {
            result[i] = x[i-shift];
          }
          else {
            result[i] = sqrt(-1.0); // create a NaN
          }
        }
 
        Push(result);
      }
    }
  }
 
  void Calculator::MinimumLine() {
    QVector<double> result = Pop();
 
    double minVal = result[0];
    for (int i = 0; i < result.size(); i++) {
      if (!IsSpecial(result[i])) {
        minVal = min(minVal, result[i]);
      }
    }
 
    result.clear();
    result.push_back(minVal);
    Push(result);
  }
 
 
  void Calculator::MaximumLine() {
    QVector<double> result = Pop();
 
    double maxVal = result[0];
    for (int i = 0; i < result.size(); i++) {
      if (!IsSpecial(result[i])) {
        maxVal = max(maxVal, result[i]);
      }
    }
 
    result.clear();
    result.push_back(maxVal);
    Push(result);
  }
 
 
  void Calculator::MinimumPixel() {
    QVector<double> x = Pop();
    QVector<double> y = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     MinimumOperator);
    Push(result);
  }
 
 
  void Calculator::MaximumPixel() {
    QVector<double> x = Pop();
    QVector<double> y = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     MaximumOperator);
    Push(result);
  }
 
 
  void Calculator::GreaterThan() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     GreaterThanOperator);
    Push(result);
  }
 
 
  void Calculator::LessThan() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     LessThanOperator);
    Push(result);
  }
 
 
  void Calculator::Equal() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     EqualOperator);
    Push(result);
  }
 
 
  void Calculator::GreaterThanOrEqual() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     GreaterThanOrEqualOperator);
    Push(result);
  }
 
 
  void Calculator::LessThanOrEqual() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     LessThanOrEqualOperator);
    Push(result);
  }
 
 
  void Calculator::NotEqual() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     NotEqualOperator);
    Push(result);
  }
 
 
  // Commented out because bitwise ops only work with integers instead of
  // doubles
 
  void Calculator::And() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     BitwiseAndOperator);
    Push(result);
  }
 
 
  void Calculator::Or() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(),
                     BitwiseOrOperator);
    Push(result);
  }
 
 
  void Calculator::Sine() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), sin);
    Push(result);
  }
 
 
  void Calculator::Cosine() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), cos);
    Push(result);
  }
 
 
  void Calculator::Tangent() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), tan);
    Push(result);
  }
 
 
  void Calculator::Cosecant() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), CosecantOperator);
    Push(result);
  }
 
 
  void Calculator::Secant() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), SecantOperator);
    Push(result);
  }
 
 
  void Calculator::Cotangent() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), CotangentOperator);
    Push(result);
  }
 
 
  void Calculator::Arcsine() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), asin);
    Push(result);
  }
 
 
  void Calculator::Arccosine() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), acos);
    Push(result);
  }
 
 
  void Calculator::Arctangent() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), atan);
    Push(result);
  }
 
 
  void Calculator::ArcsineH() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), asinh);
    Push(result);
  }
 
 
  void Calculator::ArccosineH() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), acosh);
    Push(result);
  }
 
 
  void Calculator::ArctangentH() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), atanh);
    Push(result);
  }
 
 
  void Calculator::Arctangent2() {
    QVector<double> y = Pop();
    QVector<double> x = Pop();
    QVector<double> result;
 
    PerformOperation(result, x.begin(), x.end(), y.begin(), y.end(), atan2);
    Push(result);
  }
 
 
  void Calculator::SineH() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), sinh);
    Push(result);
  }
 
 
  void Calculator::CosineH() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), cosh);
    Push(result);
  }
 
 
  void Calculator::TangentH() {
    QVector<double> result = Pop();
    PerformOperation(result, result.begin(), result.end(), tanh);
    Push(result);
  }
 
 
  // Stack methods
 
  int Calculator::StackSize() {
    return p_valStack->size();
  }
 
  void Calculator::Push(QVector<double> &vect) {
    p_valStack->push(vect);
  }
 
 
  void Calculator::Push(double scalar) {
    QVector<double> s;
    s.push_back(scalar);
    Push(s);
  }
 
 
  void Calculator::Push(Buffer &buff) {
    QVector<double> b(buff.size());
 
    for (int i = 0; i < buff.size(); i++) {
      // Test for special pixels and map them to valid values
      if (IsSpecial(buff[i])) {
        if (Isis::IsNullPixel(buff[i])) {
          //b[i] = NAN;
          b[i] = sqrt(-1.0);
        }
        else if (Isis::IsHrsPixel(buff[i])) {
          //b[i] = INFINITY;
          b[i] = DBL_MAX * 2;
        }
        else if (Isis::IsHisPixel(buff[i])) {
          //b[i] = INFINITY;
          b[i] = DBL_MAX * 2;
        }
        else if (Isis::IsLrsPixel(buff[i])) {
          //b[i] = -INFINITY;
          b[i] = -DBL_MAX * 2;
        }
        else if (Isis::IsLisPixel(buff[i])) {
          //b[i] = -INFINITY;
          b[i] = -DBL_MAX * 2;
        }
      }
      else
        b[i] = buff[i];
    }
 
    Push(b);
  }
 
 
  QVector<double> Calculator::Pop(bool keepSpecials) {
    QVector<double> top;
 
    if (p_valStack->empty()) {
      IString msg = "Math calculator stack is empty, cannot perform any "
                    "more operations.";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
 
    top = p_valStack->top();
 
    if (keepSpecials) {
      for (int i = 0; i < (int)top.size(); i++) {
        if (std::isnan(top[i])) {
          top[i] = Isis::Null;
        }
        // Test for +INFINITY
        else if (top[i] > DBL_MAX) {
          top[i] = Isis::Hrs;
        }
        // Test for -INFINITY)
        else if (top[i] < -DBL_MAX) {
          top[i] = Isis::Lrs;
        }
        else {
          // Do nothing
        }
      }
    }
 
    p_valStack->pop();
 
 
    return top;
  }
 
 
  void Calculator::PrintTop() {
    if (p_valStack->empty()) return;
 
    QString temp;
    temp += "[ ";
    QVector<double> top = p_valStack->top();
    for (int i = 0; i < (int)top.size(); i++) {
      temp += QString::number(top[i]);
      temp += " ";
    }
    temp += "]";
    // On some operating systems, -nan was being outputted. 
    // Because this method is only used as a cout in our tests, we do not 
    // care about the difference between nan and -nan; they are the same in this case.
    temp.replace(QRegExp("-nan"), "nan");
    std::cout<<temp<<std::endl;
  }
 
 
  bool Calculator::Empty() {
    return p_valStack->empty();
  }
 
 
  void Calculator::Clear() {
    while (!p_valStack->empty()) {
      p_valStack->pop();
    }
  }
 
 
  void Calculator::PerformOperation(QVector<double> &results,
                                    QVector<double>::iterator arg1Start,
                                    QVector<double>::iterator arg1End,
                                    double operation(double)) {
    results.resize(arg1End - arg1Start);
 
    for (int pos = 0; pos < results.size(); pos++) {
      results[pos] = operation(*arg1Start);
 
      arg1Start++;
    }
  }
 
 
  void Calculator::PerformOperation(QVector<double> &results,
                                    QVector<double>::iterator arg1Start,
                                    QVector<double>::iterator arg1End,
                                    QVector<double>::iterator arg2Start,
                                    QVector<double>::iterator arg2End,
                                    double operation(double, double)) {
    if (arg1End - arg1Start != 1 && arg2End - arg2Start != 1 &&
        arg1End - arg1Start != arg2End - arg2Start) {
      IString msg = "The stack based calculator cannot operate on vectors "
                    "of differing sizes.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    int iSize = max(arg1End - arg1Start, arg2End - arg2Start);
    results.resize(iSize);
 
    for (int pos = 0; pos < results.size(); pos++) {
      results[pos] = operation(*arg1Start, *arg2Start);
 
      if (arg1Start + 1 != arg1End) arg1Start++;
      if (arg2Start + 1 != arg2End) arg2Start++;
    }
  }
} // End of namespace Isis