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) +
           exp(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 <typename Type>
inline fims::Vector<Type> apply_ageing_error(
    const fims::Vector<Type> &proportion_at_true_age,
    const fims::Vector<Type> &ageing_error, size_t year = 0) {
  const size_t n_ages = proportion_at_true_age.size();
  const size_t matrix_size = n_ages * n_ages;
  const size_t matrix_offset =
      (ageing_error.size() == matrix_size) ? 0 : year * matrix_size;
 
  fims::Vector<Type> proportion_at_observed_age(n_ages, static_cast<Type>(0.0));
  for (size_t true_age = 0; true_age < n_ages; true_age++) {
    for (size_t observed_age = 0; observed_age < n_ages; observed_age++) {
      proportion_at_observed_age[observed_age] +=
          proportion_at_true_age[true_age] *
          ageing_error[matrix_offset + true_age * n_ages + observed_age];
    }
  }
  return proportion_at_observed_age;
}

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 */