logic.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2012 Stefan Buehler  <sbuehler@ltu.se>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
#include "logic.h"
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <stdexcept>
 
#include "sorting.h"
 
// For checking, if a Numeric equal zero we have to take into account the
// numerical precicion. If a value is smaller than *precision* it is
// taken to be 0.
#define precision 0.
 
 
bool is_bool(const Index& x) { return (x == 0 || x == 1); }
 
 
bool is_multiple(const Index& x, const Index& y) {
  ARTS_ASSERT(y != 0);
  return (0 == fmod(Numeric(x), Numeric(y)));
}
 
 
bool is_size(ConstVectorView x, const Index& n) { return (n == x.nelem()); }
 
 
bool is_size(ConstMatrixView x, const Index& r, const Index& c) {
  return (r == x.nrows() && c == x.ncols());
}
 
 
bool is_size(ConstTensor3View x,
             const Index& p,
             const Index& r,
             const Index& c) {
  return (p == x.npages() && r == x.nrows() && c == x.ncols());
}
 
 
bool is_size(ConstTensor4View x,
             const Index& b,
             const Index& p,
             const Index& r,
             const Index& c) {
  return (b == x.nbooks() && p == x.npages() && r == x.nrows() &&
          c == x.ncols());
}
 
 
bool is_size(ConstTensor5View x,
             const Index& s,
             const Index& b,
             const Index& p,
             const Index& r,
             const Index& c) {
  return (s == x.nshelves() && b == x.nbooks() && p == x.npages() &&
          r == x.nrows() && c == x.ncols());
}
 
 
bool is_size(ConstTensor6View x,
             const Index& v,
             const Index& s,
             const Index& b,
             const Index& p,
             const Index& r,
             const Index& c) {
  return (v == x.nvitrines() && s == x.nshelves() && b == x.nbooks() &&
          p == x.npages() && r == x.nrows() && c == x.ncols());
}
 
 
bool is_size(ConstTensor7View x,
             const Index& l,
             const Index& v,
             const Index& s,
             const Index& b,
             const Index& p,
             const Index& r,
             const Index& c) {
  return (l == x.nlibraries() && v == x.nvitrines() && s == x.nshelves() &&
          b == x.nbooks() && p == x.npages() && r == x.nrows() &&
          c == x.ncols());
}
 
 
bool is_sorted(ConstVectorView x) {
  if (x.nelem() > 1) {
    for (Index i = 1; i < x.nelem(); i++) {
      if (!(x[i] >= x[i - 1])) return false;
    }
  }
  return true;
}
 
 
bool is_increasing(ConstVectorView x) {
  if (x.nelem() > 1) {
    for (Index i = 1; i < x.nelem(); i++) {
      if (!(x[i] > x[i - 1])) return false;
    }
  }
  return true;
}
 
 
bool is_increasing(const ArrayOfIndex& x) {
  if (x.nelem() > 1) {
    for (Index i = 1; i < x.nelem(); i++) {
      if (x[i] <= x[i - 1]) return false;
    }
  }
  return true;
}
 
 
bool is_decreasing(ConstVectorView x) {
  if (x.nelem() > 1) {
    for (Index i = 1; i < x.nelem(); i++) {
      if (!(x[i] < x[i - 1])) return false;
    }
  }
  return true;
}
 
 
bool is_unique(const ArrayOfIndex& x) {
  // We simply compare the second element to the first,
  // the third to the first and second, and so on.
 
  for (Index i = 1; i < x.nelem(); ++i)
    for (Index s = 0; s < i; ++s)
      if (x[i] == x[s]) return false;
 
  return true;
}
 
 
bool is_singular(ConstMatrixView A) {
  ARTS_ASSERT(A.nrows() == A.ncols());
  Numeric temp = 0.;
 
  for (Index i = 0; i < A.nrows(); i++) {
    Numeric big = 0.;
    for (Index j = 0; j < A.nrows(); j++) {
      if ((temp = fabs(A(i, j))) > big) big = temp;
    }
    // Due to numerical precision the values can deviate from 0.0
    if (big < precision) {
      throw runtime_error("Matrix is singular.");
      return true;
    }
  }
  return false;
}
 
 
bool is_diagonal(ConstMatrixView A) {
  ARTS_ASSERT(A.nrows() == A.ncols());
 
  for (Index i = 1; i < A.ncols(); i++) {
    for (Index j = 0; j < i; j++) {
      if (fabs(A(i, j)) > precision || fabs(A(j, i)) > precision) return false;
    }
  }
  return true;
}
 
 
bool is_same_within_epsilon(const Numeric& a,
                            const Numeric& b,
                            const Numeric& epsilon) {
  return abs(a - b) <= epsilon * max(abs(a), abs(b));
}
 
 
bool is_lon_cyclic(ConstVectorView grid, const Numeric& epsilon) {
  return is_same_within_epsilon(
      grid[grid.nelem() - 1] - grid[0], 360., epsilon);
}
 
bool is_regularly_increasing_within_epsilon(ConstVectorView x,
                                            const Numeric epsilon) {
  if (x.size() < 2) return true;
  const Numeric dx0 = x[1] - x[0];
  if (dx0 <= 0) return false;
 
  for (Index i = 2; i < x.size(); i++) {
    if (Numeric dxi = x[i] - x[i - 1];
        dxi <= 0 or std::abs(dxi - dx0) >= epsilon) {
      return false;
    }
  }
  return true;
}