fims_math.hpp

Go to the documentation of this file.

 
#ifndef FIMS_MATH_HPP
#define FIMS_MATH_HPP
 
// note: this is modeling platform specific, must be controlled by
// preprocessing macros
#include <cmath>
#include <random>
#include <sstream>
 
#include "../interface/interface.hpp"
#include "fims_vector.hpp"
#include "def.hpp"
 
namespace fims_math {
#ifdef STD_LIB
 
template <class Type>
inline const Type exp(const Type &x) {
  return std::exp(x);
}
 
template <class Type>
inline const Type log(const Type &x) {
  return std::log(x);
}
 
template <class Type>
inline const Type cos(const Type &x) {
  return std::cos(x);
}
 
template <class Type>
inline const Type sqrt(const Type &x) {
  return std::sqrt(x);
}
 
template <class Type>
inline const Type pow(const Type &x, const Type &y) {
  return std::pow(x, y);
}
 
template <class Type>
inline const Type lgamma(const Type &x) {
  return std::lgamma(x);
}
 
inline double Value(const double &x) { return x; }
 
template <class Type>
inline const Type pnorm(const Type &x, const Type &mean, const Type &sd) {
  const double z = static_cast<double>(x - mean) /
                   (static_cast<double>(sd) * std::sqrt(2.0));
  return static_cast<Type>(0.5 * (1.0 + std::erf(z)));
}
 
#endif
 
#ifdef TMB_MODEL
 
inline double Value(const TMBAD_FIMS_TYPE &x) { return x.Value(); }
 
inline double Value(const double &x) { return x; }
 
// Add the following line to CMakeLists.txt to enable documentation of TMB_MODEL
// in doxygen or none of the following is rendered.
// set(DOXYGEN_PREDEFINED "TMB_MODEL=1" "ENABLE_TMB_CODE")
 
template <class Type>
inline const Type exp(const Type &x) {
  // use std::exp for double type, look for TMB version of exp if AD type
  using std::exp;
  return exp(x);
}
 
template <class Type>
inline const Type log(const Type &x) {
  // use std::log for double type, look for TMB version of log if AD type
  using std::log;
  return log(x);
}
 
template <class Type>
inline const Type cos(const Type &x) {
  // use std::cos for double type, look for TMB version of cos if AD type
  using std::cos;
  return cos(x);
}
 
template <class Type>
inline const Type sqrt(const Type &x) {
  // use std::std for double type, look for TMB version of std if AD type
  using std::sqrt;
  return sqrt(x);
}
 
template <class Type>
inline const Type pow(const Type &x, const Type &y) {
  // use std::pow for double type, look for TMB version of pow if AD type
  using std::pow;
  return pow(x, y);
}
 
template <class Type>
inline const Type lgamma(const Type &x) {
  // use std::lgamma for double type, look for TMB version of lgamma if AD type
  using std::lgamma;
  return lgamma(x);
}
 
template <class Type>
inline const Type pnorm(const Type &x, const Type &mean, const Type &sd) {
  // Use :: to call TMB's global pnorm; an unqualified call recurses into this
  // wrapper and can hang CreateTMBModel() when growth uses normal
  // probabilities.
  return ::pnorm(x, mean, sd);
}
 
#endif
 
template <class Type>
inline const Type logistic(const Type &inflection_point, const Type &slope,
                           const Type &x) {
  return static_cast<Type>(1.0) /
         (static_cast<Type>(1.0) +
          exp(Type(-1.0) * slope * (x - inflection_point)));
}
 
template <class Type>
inline const Type logit(const Type &a, const Type &b, const Type &x) {
  return -fims_math::log(b - x) + fims_math::log(x - a);
}
 
template <class Type>
inline const Type inv_logit(const Type &a, const Type &b, const Type &logit_x) {
  return a + (b - a) / (static_cast<Type>(1.0) + fims_math::exp(-logit_x));
}
 
template <class Type>
inline const Type double_logistic(const Type &inflection_point_asc,
                                  const Type &slope_asc,
                                  const Type &inflection_point_desc,
                                  const Type &slope_desc, const Type &x) {
  return (static_cast<Type>(1.0)) /
         (static_cast<Type>(1.0) +
          exp(Type(-1.0) * slope_asc * (x - inflection_point_asc))) *
         (static_cast<Type>(1.0) -
          (static_cast<Type>(1.0)) /
              (static_cast<Type>(1.0) +
               exp(Type(-1.0) * slope_desc * (x - inflection_point_desc))));
}
 
template <class Type>
const Type ad_fabs(const Type &x, Type C = 1e-5) {
  return sqrt((x * x) + C);
}
 
template <typename Type>
inline const Type ad_min(const Type &a, const Type &b, Type C = 1e-5) {
  return (a + b - fims_math::ad_fabs(a - b, C)) * static_cast<Type>(0.5);
}
 
template <typename Type>
inline const Type ad_max(const Type &a, const Type &b, Type C = 1e-5) {
  return (a + b + fims_math::ad_fabs(a - b, C)) * static_cast<Type>(.5);
}
 
template <class T>
T sum(const std::vector<T> &v) {
  T ret = 0.0;
  for (size_t i = 0; i < v.size(); i++) {
    ret += v[i];
  }
  return ret;
}
 
template <class T>
T sum(const fims::Vector<T> &v) {
  T ret = 0.0;
  for (size_t i = 0; i < v.size(); i++) {
    ret += v[i];
  }
  return ret;
}
 
}  // namespace fims_math
 
#endif /* FIMS_MATH_HPP */