From b0efac2599e63f5cac3328a1cd6e6d7b6861fe55 Mon Sep 17 00:00:00 2001 From: kellijohnson-NOAA Date: Mon, 23 Sep 2024 19:57:15 -0700 Subject: [PATCH] feat(ParameterVector): All parameters are vectors to allow for time-varying parameters. Also, removed old log because it was conflicting on Windows. Tests are update as well as vignettes. --- CMakeLists.txt | 2 +- NAMESPACE | 4 + R/zzz.R | 12 +- inst/include/common/def.hpp | 27 +- inst/include/common/fims_math.hpp | 837 +++++++++++---- inst/include/common/fims_vector.hpp | 19 +- inst/include/common/information.hpp | 47 +- .../distributions/functors/lognormal_lpdf.hpp | 4 +- .../functors/multinomial_lpmf.hpp | 4 +- .../distributions/functors/normal_lpdf.hpp | 4 +- .../include/interface/rcpp/rcpp_interface.hpp | 23 +- .../rcpp/rcpp_objects/rcpp_fleet.hpp | 351 +++++-- .../rcpp/rcpp_objects/rcpp_interface_base.hpp | 171 ++-- .../rcpp/rcpp_objects/rcpp_maturity.hpp | 118 ++- .../rcpp/rcpp_objects/rcpp_population.hpp | 76 +- .../rcpp/rcpp_objects/rcpp_recruitment.hpp | 180 ++-- .../rcpp/rcpp_objects/rcpp_selectivity.hpp | 338 +++--- .../rcpp_objects/rcpp_tmb_distribution.hpp | 966 +++++++++++------- .../population_dynamics/fleet/fleet.hpp | 292 +++--- .../growth/functors/ewaa.hpp | 55 +- .../population/population.hpp | 4 +- tests/gtest/integration_test_population.cpp | 4 +- tests/gtest/test_population_Index.cpp | 2 +- ...tion_dynamics_fleet_initialize_prepare.cpp | 9 +- tests/gtest/test_population_test_fixture.hpp | 674 ++++++------ tests/integration/integration_class.hpp | 12 +- .../testthat/helper-integration-tests-setup.R | 52 +- tests/testthat/test-rcpp-fims.R | 6 +- tests/testthat/test-rcpp-get_fixed.R | 52 +- tests/testthat/test-rcpp-maturity-interface.R | 24 +- .../test-rcpp-recruitment-interface.R | 26 +- .../test-rcpp-selectivity-interface.R | 36 +- ...test-unit-rcpp-interface-variable-vector.R | 54 +- vignettes/fims-demo.Rmd | 52 +- vignettes/fims-path-maturity.Rmd | 12 +- 35 files changed, 2750 insertions(+), 1799 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 64d98dd3..008dea9e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -7,7 +7,7 @@ project(FIMS ) # CXX is the language name # GoogleTest requires at least C++11 -set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD 17) include(FetchContent) diff --git a/NAMESPACE b/NAMESPACE index 3c7c0bd8..105b76b7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -3,6 +3,8 @@ export(AgeComp) export(BevertonHoltRecruitment) export(CreateTMBModel) +export(SetFIMSFunctions) +export(Finalize) export(DoubleLogisticSelectivity) export(EWAAgrowth) export(FIMSFrame) @@ -16,7 +18,9 @@ export(Population) export(TMBDlnormDistribution) export(TMBDmultinomDistribution) export(TMBDnormDistribution) +export(ToJSON) export(clear) +export(get_log) export(clear_logs) export(get_fixed) export(get_random) diff --git a/R/zzz.R b/R/zzz.R index a3f6e85c..c03c6807 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -4,15 +4,15 @@ Rcpp::loadModule(module = "fims", what = TRUE) library.dynam.unload("FIMS", libpath) } -setMethod("[<-", signature(x = "Rcpp_ParameterVector", i = "numeric"), - function(x, i) { - (x$at(i)) - return(x) - }) + +setMethod("[<-", signature(x = "Rcpp_ParameterVector"), function(x, i, j, value) { + x$set(i - 1, value) # R uses 1-based indexing, C++ uses 0-based indexing + x # Return the modified object +}) setMethod("[", signature(x = "Rcpp_ParameterVector", i = "numeric"), function(x, i) { - return(x$at(i)) + return(x$get(i-1)) }) # setMethod("lapply", signature(X = "Rcpp_ParameterVector", FUN = "sum"), diff --git a/inst/include/common/def.hpp b/inst/include/common/def.hpp index 7d671c7a..d895aaf8 100644 --- a/inst/include/common/def.hpp +++ b/inst/include/common/def.hpp @@ -21,7 +21,7 @@ // comments used to assist in diagnosing model issues and tracking progress. // These files will only be created if a logs folder is added to the root model // directory. -std::ofstream FIMS_LOG("logs/fims.log"); /**< Generic log file */ +std::ofstream FIMS_LOG_OLD("logs/fims.log"); /**< Generic log file */ std::ofstream INFO_LOG("logs/info.log"); /**< Information.hpp log file */ std::ofstream ERROR_LOG("logs/error.log"); /**< Error tracking log file */ std::ofstream DATA_LOG("logs/data.log"); /**< Data input tracking log file */ @@ -50,31 +50,6 @@ std::ofstream DEBUG_LOG( namespace fims { -/** - * A static class for FIMS logging. - */ - -class fims_log { - public: - static std::map - FIMS_LOGS; /**< Map Log of files */ - /** - * Static getter for retrieving a specific log file. - */ - static std::ofstream& get(const std::string& l) { - typename std::map::iterator it; - it = fims_log::FIMS_LOGS.find(l); - if (it == fims_log::FIMS_LOGS.end()) { - std::ofstream& of = fims_log::FIMS_LOGS[l]; - of.open(l.c_str()); - } - - return fims_log::FIMS_LOGS[l]; - } -}; - -std::map fims_log::FIMS_LOGS; - } // namespace fims #endif /* TRAITS_HPP */ diff --git a/inst/include/common/fims_math.hpp b/inst/include/common/fims_math.hpp index dcdf42d0..aee57c45 100644 --- a/inst/include/common/fims_math.hpp +++ b/inst/include/common/fims_math.hpp @@ -18,208 +18,677 @@ // preprocessing macros //#include "def.hpp" #include +#include +#include #include "../interface/interface.hpp" +#include "fims_vector.hpp" namespace fims_math { #ifdef STD_LIB -/** - * @brief The exponential function. - * - * @param x value to exponentiate. Please use fims_math::exp(x) if x is - * an integer. - * @return the exponentiated value - */ -template -inline const Type exp(const Type &x) { - return std::exp(x); -} - -/** - * @brief The natural log function (base e) - * @param x the value to take the log of. Please use fims_math::log(x) - * if x is an integer. - * @return - */ -template -inline const Type log(const Type &x) { - return std::log(x); -} + + /** + * @brief The exponential function. + * + * @param x value to exponentiate. Please use fims_math::exp(x) if x is + * an integer. + * @return the exponentiated value + */ + template + inline const Type exp(const Type &x) { + return std::exp(x); + } + + /** + * @brief The natural log function (base e) + * @param x the value to take the log of. Please use fims_math::log(x) + * if x is an integer. + * @return + */ + template + inline const Type log(const Type &x) { + return std::log(x); + } + + /** + * @brief The natural log function (base e) + * @param x the value to take the log of. Please use fims_math::log(x) + * if x is an integer. + * @return + */ + template + inline const Type cos(const Type &x) { + return std::cos(x); + } + + template + inline const Type sqrt(const Type &x) { + return std::sqrt(x); + } + + template + inline const Type pow(const Type &x, const Type &y) { + return std::pow(x, y); + } #endif #ifdef TMB_MODEL -// #include - -/** - * @brief The exponential function. - * The code cannot be tested using the compilation flag - * -DTMB_MODEL through CMake and Google Test - * @param x value to exponentiate. Please use fims_math::exp(x) if x is - * an integer. - * @return the exponentiated value - */ -template -inline const Type exp(const Type &x) { - using ::exp; - return exp(x); -} - -template <> -inline const double exp(const double &x) { - return std::exp(x); -} - -/** - * @brief The natural log function (base e) - * The code cannot be tested using the compilation flag - * -DTMB_MODEL through CMake and Google Test. - * @param x the value to take the log of. Please use fims_math::log(x) - * if x is an integer. - * @return the log of the value - */ -template -inline const Type log(const Type &x) { - return log(x); -} + // #include + + /** + * @brief The exponential function. + * The code cannot be tested using the compilation flag + * -DTMB_MODEL through CMake and Google Test + * @param x value to exponentiate. Please use fims_math::exp(x) if x is + * an integer. + * @return the exponentiated value + */ + template + inline const Type exp(const Type &x) { + using ::exp; + return exp(x); + } + + template <> + inline const double exp(const double &x) { + return std::exp(x); + } -template <> -inline const double log(const double &x) { - return std::log(x); -} + /** + * @brief The natural log function (base e) + * The code cannot be tested using the compilation flag + * -DTMB_MODEL through CMake and Google Test. + * @param x the value to take the log of. Please use fims_math::log(x) + * if x is an integer. + * @return the log of the value + */ + template + inline const Type log(const Type &x) { + return log(x); + } + + template <> + inline const double log(const double &x) { + return std::log(x); + } + + template + inline const Type cos(const Type &x) { + return cos(x); + } + + template <> + inline const double cos(const double &x) { + return std::cos(x); + } + + template + inline const Type sqrt(const Type &x) { + return sqrt(x); + } + + template <> + inline const double sqrt(const double &x) { + return std::sqrt(x); + } + + template + inline const Type pow(const Type &x, const Type &y) { + return pow(x, y); + } + + template <> + inline const double pow(const double &x, const double &y) { + return std::pow(x, y); + } #endif -/** - * @brief The general logistic function - * - * \f$ \frac{1.0}{ 1.0 + exp(-1.0 * slope (x - inflection_point))} \f$ - * - * @param inflection_point the inflection point of the logistic function - * @param slope the slope of the logistic function - * @param x the index the logistic function should be evaluated at - * @return - */ -template -inline const Type logistic(const Type &inflection_point, const Type &slope, - const Type &x) { - return (1.0) / (1.0 + exp(-1.0 * slope * (x - inflection_point))); -} - -/** - * @brief A logit function for bounding of parameters - * - * \f$ -\mathrm{log}(b-x) + \mathrm{log}(x-a) \f$ - * @param a lower bound - * @param b upper bound - * @param x the parameter in bounded space - * @return the parameter in real space - * - */ -template -inline const Type logit(const Type &a, const Type &b, const Type &x) { - return -fims_math::log(b - x) + fims_math::log(x - a); -} + /** + * @brief The general logistic function + * + * \f$ \frac{1.0}{ 1.0 + exp(-1.0 * slope (x - inflection_point))} \f$ + * + * @param inflection_point the inflection point of the logistic function + * @param slope the slope of the logistic function + * @param x the index the logistic function should be evaluated at + * @return + */ + template + inline const Type logistic(const Type &inflection_point, const Type &slope, + const Type &x) { + return (1.0) / (1.0 + exp(-1.0 * slope * (x - inflection_point))); + } -/** - * @brief An inverse logit function for bounding of parameters - * - * \f$ a+\frac{b-a}{1+\mathrm{exp}(-\mathrm{logit}(x))}\f$ - * @param a lower bound - * @param b upper bound - * @param logit_x the parameter in real space - * @return the parameter in bounded space - * - */ -template -inline const Type inv_logit(const Type &a, const Type &b, const Type &logit_x) { - return a + (b - a) / (1.0 + fims_math::exp(-logit_x)); -} + /** + * @brief A logit function for bounding of parameters + * + * \f$ -\mathrm{log}(b-x) + \mathrm{log}(x-a) \f$ + * @param a lower bound + * @param b upper bound + * @param x the parameter in bounded space + * @return the parameter in real space + * + */ + template + inline const Type logit(const Type &a, const Type &b, const Type &x) { + return -fims_math::log(b - x) + fims_math::log(x - a); + } -/** - * @brief The general double logistic function - * - * \f$ \frac{1.0}{ 1.0 + exp(-1.0 * slope_{asc} (x - inflection_point_{asc}))} - * \left(1-\frac{1.0}{ 1.0 + exp(-1.0 * slope_{desc} (x - - * inflection_point_{desc}))} \right)\f$ - * - * @param inflection_point_asc the inflection point of the ascending limb of the - * double logistic function - * @param slope_asc the slope of the ascending limb of the double logistic - * function - * @param inflection_point_desc the inflection point of the descending limb of - * the double logistic function, where inflection_point_desc > - * inflection_point_asc - * @param slope_desc the slope of the descending limb of the double logistic - * function - * @param x the index the logistic function should be evaluated at - * @return - */ + /** + * @brief An inverse logit function for bounding of parameters + * + * \f$ a+\frac{b-a}{1+\mathrm{exp}(-\mathrm{logit}(x))}\f$ + * @param a lower bound + * @param b upper bound + * @param logit_x the parameter in real space + * @return the parameter in bounded space + * + */ + template + inline const Type inv_logit(const Type &a, const Type &b, const Type &logit_x) { + return a + (b - a) / (1.0 + fims_math::exp(-logit_x)); + } -template -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 (1.0) / (1.0 + exp(-1.0 * slope_asc * (x - inflection_point_asc))) * - (1.0 - - (1.0) / (1.0 + exp(-1.0 * slope_desc * (x - inflection_point_desc)))); -} - -/** - * - * Used when x could evaluate to zero, which will result in a NaN for - * derivative values. - * - * Evaluates: - * - * \f$ (x^2+C)^.5 \f$ - * - * @param x value to keep positive - * @param C default = 1e-5 - * @return - */ -template -const Type ad_fabs(const Type &x, Type C = 1e-5) { - return sqrt((x * x) + C); //, .5); -} + /** + * @brief The general double logistic function + * + * \f$ \frac{1.0}{ 1.0 + exp(-1.0 * slope_{asc} (x - inflection_point_{asc}))} + * \left(1-\frac{1.0}{ 1.0 + exp(-1.0 * slope_{desc} (x - + * inflection_point_{desc}))} \right)\f$ + * + * @param inflection_point_asc the inflection point of the ascending limb of the + * double logistic function + * @param slope_asc the slope of the ascending limb of the double logistic + * function + * @param inflection_point_desc the inflection point of the descending limb of + * the double logistic function, where inflection_point_desc > + * inflection_point_asc + * @param slope_desc the slope of the descending limb of the double logistic + * function + * @param x the index the logistic function should be evaluated at + * @return + */ -/** - * - * Returns the minimum between a and b in a continuous manner using: - * - * (a + b - fims_math::ad_fabs(a - b))*.5; - * Reference: \ref fims_math::ad_fabs() - * - * This is an approximation with minimal error. - * - * @param a - * @param b - * @param C default = 1e-5 - * @return - */ + template + 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 (1.0) / (1.0 + exp(-1.0 * slope_asc * (x - inflection_point_asc))) * + (1.0 - + (1.0) / (1.0 + exp(-1.0 * slope_desc * (x - inflection_point_desc)))); + } + + /** + * + * Used when x could evaluate to zero, which will result in a NaN for + * derivative values. + * + * Evaluates: + * + * \f$ (x^2+C)^.5 \f$ + * + * @param x value to keep positive + * @param C default = 1e-5 + * @return + */ + template + const Type ad_fabs(const Type &x, Type C = 1e-5) { + return sqrt((x * x) + C); //, .5); + } + + /** + * + * Returns the minimum between a and b in a continuous manner using: + * + * (a + b - fims_math::ad_fabs(a - b))*.5; + * Reference: \ref fims_math::ad_fabs() + * + * This is an approximation with minimal error. + * + * @param a + * @param b + * @param C default = 1e-5 + * @return + */ + + template + 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)) * 0.5; + } + + /** + * Returns the maximum between a and b in a continuous manner using: + * + * (a + b + fims_math::ad_fabs(a - b)) *.5; + * Reference: \ref fims_math::ad_fabs() + * This is an approximation with minimal error. + * + * @param a + * @param b + * @param C default = 1e-5 + * @return + */ + template + 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 (.5); + } + + /** + * \ingroup NormalDistribution + * + * @brief The probability density function for a standard normal distribution. + * + * \f$ + * f(x) = \frac{1}{\sqrt{2\pi\theta^2}}exp(-\frac{(x-\mu)^2}{2\theta^2}) + * \f$
+ * Where,
+ * \f$\ \mu \f$ is the mean.
+ * \f$\ \theta\f$ is the standard deviation.
+ * + * @image html http://upload.wikimedia.org/wikipedia/commons/7/74/Normal_Distribution_PDF.svg + * + *
Source: http://en.wikipedia.org/wiki/Normal_distribution
+ * + * + * + * @author Matthew Supernaw + * @date June 6, 2013 + * + * @param x + * @param mean + * @param sd + * @return relative likelihood for this random variable to have value x. + */ + template + const T dnorm(const T &x, const T &mean, const T &sd, bool ret_log = false) { + + // Check if the standard deviation is valid + if (sd <= static_cast (0)) { + throw std::invalid_argument("Standard deviation must be positive."); + } + + // Constants + const T inv_sqrt_2pi = 1.0 / std::sqrt(2.0 * M_PI); + + // Compute the exponent part + T z = (x - mean) / sd; + T exponent = -1.0*static_cast (0.5) * z * z; + + // Compute log density + if (ret_log) { + return exponent - log(sd) - log(inv_sqrt_2pi); + } + + // Otherwise, return the density + return inv_sqrt_2pi / sd * exp(exponent); + } + + /** + * @brief Normal distribution random sampling. + * + * @param mean The mean of the distribution \( \mu \). + * @param sd The standard deviation of the distribution \( \sigma \). + * @param gen Random number generator. + * @return T Random sample from the Normal distribution. + */ + template + T rnorm(T mean, T sd) { + static std::random_device rd; + static std::mt19937 gen(rd()); + static std::uniform_real_distribution dist(0.0, 1.0); + + T u1 = static_cast (dist(gen)); + T u2 = static_cast (dist(gen)); + + T z0 = fims_math::sqrt(-2.0 * log(u1)) * fims_math::cos(2.0 * M_PI * u2); + + return z0 * sd + mean; + } + + /** + * @ingroup LogNormal + * + * @brief The probability distribution function for a log-normal distribution. + * + * \f$ + * f(x) = \frac{1}{x\sqrt{2\pi}\theta}e^{-\frac{(ln x - \mu)^{2}}{2\theta^{2}}} + * \f$ + * + *
Where,
+ * \f$\mu\f$ is the log of the mean.
+ * \f$\theta\f$ is the log of the standard deviation.
+ * + * @image html http://upload.wikimedia.org/wikipedia/commons/thumb/8/80/Some_log-normal_distributions.svg/500px-Some_log-normal_distributions.svg.png + * + *
Source:http://en.wikipedia.org/wiki/Log-normal_distribution

+ * + * + * @param x + * @param meanLog + * @param sdLog + * @return + */ + template + const T dlnorm(const T &x, const T &meanLog, const T &sdLog, bool ret_log = false) { + if (x > T(0)) { + T ret = (T(1) / (x * sdLog * pow(T(2) * T(M_PI), 0.5))) * + exp(T(-1)*((std::log(x) - meanLog)* + (log(x) - meanLog)) / (T(2) * sdLog)); + + if (!ret_log) { + return ret; + } else { + return log(ret); + } + } + return T(0); + } + + /** + * @brief Computes the Gamma function using the Lanczos approximation. + * + * \f[ + * \Gamma(z) \approx \sqrt{2 \pi} \left( z + \frac{g + 0.5}{z} \right)^{z + 0.5} e^{-(z + g + 0.5)} A_g(z) + * \f] + * + * Where \( A_g(z) \) is a series with coefficients. + * + * @param z Input value. + * @return T Value of the Gamma function at z. + */ + template + T gamma(T z) { + if (z <= 0) { + throw std::invalid_argument("z must be greater than 0"); + } + + static const std::vector coef = { + 676.5203681218851, -1259.1392167224028, + 771.32342877765313, -176.61502916214059, + 12.507343278686905, -0.13857109526572012, + 9.9843695780195716e-6, 1.5056327351493116e-7 + }; + + static const T sqrt_two_pi = fims_math::sqrt(2 * M_PI); + + T x = 0.99999999999980993; + for (size_t i = 0; i < coef.size(); ++i) { + x += coef[i] / (z + i + 1); + } + + T t = z + coef.size() - 0.5; + return sqrt_two_pi * fims_math::pow(t, z + 0.5) * fims_math::exp(-t) * x; + } + + template + T lgamma(T x) {// x must be positive + + if (x <= 0.0) { + std::stringstream os; + os << "Invalid input argument " << x << ". Argument must be positive."; + throw std::invalid_argument(os.str()); + } + + if (x < 12.0) { + return fims_math::log(fims_math::ad_fabs(fims_math::gamma(x))); + } + + // Abramowitz and Stegun 6.1.41 + // Asymptotic series should be good to at least 11 or 12 figures + // For error analysis, see Whittiker and Watson + // A Course in Modern Analysis (1927), page 252 + + static const T c[8] = { + 1.0 / 12.0, + -1.0 / 360.0, + 1.0 / 1260.0, + -1.0 / 1680.0, + 1.0 / 1188.0, + -691.0 / 360360.0, + 1.0 / 156.0, + -3617.0 / 122400.0 + }; + T z = 1.0 / (x * x); + T sum = c[7]; + for (int i = 6; i >= 0; i--) { + sum *= z; + sum += c[i]; + } + T series = sum / x; + + + static const T halfLogTwoPi = 0.91893853320467274178032973640562; + T logGamma = (x - 0.5) * log(x) - x + halfLogTwoPi + series; + return logGamma; + } + + template + T gamma_p(T a, T x) { + if (x < 0 || a <= 0) { + return T(NAN); // Not a Number if input is invalid + } + + if (x == 0) { + return 0.0; // P(a, 0) is always 0 + } + + T result; + if (x < a + 1.0) { + // Series representation + T sum = 1.0 / a; + T term = 1.0 / a; + for (int n = 1; n < 100; ++n) { // Number of iterations for series expansion + term *= x / (a + n); + sum += term; + if (term < 1e-10) { + break; + } + } + result = sum * exp(-x + a * log(x) - fims_math::lgamma(a)); + } else { + // Use the continued fraction representation + T b = x + 1.0 - a; + T c = 1.0 / 1e-30; + T d = 1.0 / b; + T h = d; + for (int i = 1; i < 100; ++i) { + T an = -i * (i - a); + b += 2.0; + d = an * d + b; + if (d == 0) d = 1e-30; + c = b + an / c; + if (c == 0) c = 1e-30; + d = 1.0 / d; + T del = d * c; + h *= del; + if (fims_math::ad_fabs(del - 1.0) < 1e-10) { + break; + } + } + result = 1.0 - fims_math::exp(-x + a * fims_math::log(x) - fims_math::lgamma(a)) * h; + } + return result; + } + + template + T LogGammaLanczos(const T& x) { + // Log of Gamma from Lanczos with g=5, n=6/7 + // not in A & S + static const T coef[6] = {76.18009172947146, + -86.50532032941677, 24.01409824083091, + -1.231739572450155, 0.1208650973866179E-2, + -0.5395239384953E-5}; + T LogSqrtTwoPi = T(0.91893853320467274178); + T denom = x + T(1.0); + T y = x + T(5.5); + T series = T(1.000000000190015); + for (int i = 0; i < 6; ++i) { + series += coef[i] / denom; + denom += 1.0; + } + return (LogSqrtTwoPi + (x + 0.5) * fims_math::log(y) - + y + fims_math::log(series / x)); + } + + template + fims::Vector LogGammaLanczos(const fims::Vector& v) { + fims::Vector ret(v.size()); + for (int i = 0; i < v.size(); i++) { + ret[i] = LogGammaLanczos(v[i]); + } + return ret; + } + + template + T LogGammaSeries(const T& z) { + // A & S 6.1.41 (Stirling's approximation) + T x1 = (z - 0.5) * log(z); + T x3 = 0.5 * fims_math::log(2.0 * M_PI); + + T x4 = 1 / (12 * z); + T x5 = 1 / (360 * z * z * z); + T x6 = 1 / (1260 * z * z * z * z * z); + T x7 = 1 / (1680 * z * z * z * z * z * z * z); + // more terms possible + return x1 - z + x3 + x4 - x5 + x6 - x7; + } + + template + fims::Vector LogGammaSeries(const fims::Vector& v) { + fims::Vector ret(v.size()); + for (int i = 0; i < v.size(); i++) { + ret[i] = LogGammaSeries(v[i]); + } + return ret; + } + + template + std::vector lgamma(const std::vector& v) { + std::vector ret(v.size()); + for (int i = 0; i < v.size(); i++) { + ret[i] = lgamma(v[i]); + } + return ret; + } + + template + fims::Vector lgamma(const fims::Vector& v) { + fims::Vector ret(v.size()); + for (int i = 0; i < v.size(); i++) { + ret[i] = lgamma(v[i]); + } + return ret; + } + + template + T lgamma_ad(T z) { + if (z <= 0) { + throw std::invalid_argument("z must be greater than 0"); + } + + static const T g = 7; + static const std::vector coef = { + 0.99999999999980993, 676.5203681218851, + -1259.1392167224028, 771.32342877765313, + -176.61502916214059, 12.507343278686905, + -0.13857109526572012, 9.9843695780195716e-6, + 1.5056327351493116e-7 + }; + + T x = coef[0]; + for (size_t i = 1; i < coef.size(); ++i) { + x += coef[i] / (z + i - 1.0); + } + + T t = z + g - 0.5; + return (fims_math::log(2.0 * M_PI) / 2.0) + ((z - 0.5) * fims_math::log(t)) - t + fims_math::log(x); + } + + template + fims::Vector lgamma_ad(const fims::Vector& v) { + fims::Vector ret(v.size()); + for (int i = 0; i < v.size(); i++) { + ret[i] = lgamma_ad(v[i]); + } + return ret; + } + + template + T sum(const std::vector& v) { + T ret = 0.0; + for (int i = 0; i < v.size(); i++) { + ret += v[i]; + } + return ret; + } + + template + T sum(const fims::Vector& v) { + T ret = 0.0; + for (int i = 0; i < v.size(); i++) { + ret += v[i]; + } + return ret; + } + + /** + * Multinomial Probability Density function. p is internally normalized to sum 1. + * + * @brief + * + * @param x + * @param p + * @param ret_log + * @return + */ + template + T dmultinom(const fims::Vector& x, fims::Vector p, bool ret_log = false) { + if (x.size() != p.size()) { + throw std::invalid_argument("Size of x and p must be the same."); + } + + T sum_x = fims_math::sum(x); + T log_factorial_sum_x = fims_math::lgamma_ad(sum_x + 1.0); // log(sum(x)!) + + T log_factorial_x = 0; + T sum_x_log_p = 0; + + for (size_t i = 0; i < x.size(); ++i) { + if (p[i] <= 0 || p[i] >= 1) { + throw std::invalid_argument("Probabilities must be between 0 and 1, exclusive."); + } + log_factorial_x += fims_math::lgamma_ad(x[i] + 1); // log(x_i!) + sum_x_log_p += x[i] * fims_math::log(p[i]); // x_i * log(p_i) + } + T ret = log_factorial_sum_x - log_factorial_x + sum_x_log_p; + + if (ret_log) { + return ret; + } else { + return fims_math::exp(ret); + } + } + + template + fims::Vector rmultinom(const fims::Vector& prob) { + std::random_device rd; + std::mt19937 gen(rd()); + std::discrete_distribution<> dist(prob.begin(), prob.end()); + + size_t n = prob.size(); + fims::Vector result(prob.size(), 0); + for (int i = 0; i < n; ++i) { + int index = dist(gen); + ++result[index]; + } + + return result; + } -template -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)) * 0.5; -} -/** - * Returns the maximum between a and b in a continuous manner using: - * - * (a + b + fims_math::ad_fabs(a - b)) *.5; - * Reference: \ref fims_math::ad_fabs() - * This is an approximation with minimal error. - * - * @param a - * @param b - * @param C default = 1e-5 - * @return - */ -template -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(.5); -} -} // namespace fims_math +} // namespace fims_math #endif /* FIMS_MATH_HPP */ diff --git a/inst/include/common/fims_vector.hpp b/inst/include/common/fims_vector.hpp index 0892846d..3d619ff2 100644 --- a/inst/include/common/fims_vector.hpp +++ b/inst/include/common/fims_vector.hpp @@ -2,7 +2,7 @@ #define FIMS_VECTOR_HPP #include "../interface/interface.hpp" - +#include namespace fims { /** @@ -365,4 +365,21 @@ bool operator==(const fims::Vector& lhs, const fims::Vector& rhs) { } // namespace fims +template +std::ostream& operator <<(std::ostream& out, fims::Vector& v){ + out<<"["; + + if(v.size() == 0){ + out<<"]"; + return out; + } + for(size_t i = 0; i< v.size()-1; i++){ + out< > d = (*it).second; if(d->input_type == "prior"){ @@ -193,6 +194,7 @@ class Information { * Loop over distributions and set links to distribution x value if distribution is a random effects type. */ void setup_random_effects(){ + FIMS_INFO_LOG("segment"); for(density_components_iterator it = this->density_components.begin(); it!= this->density_components.end(); ++it){ std::shared_ptr > d = (*it).second; if(d->input_type == "random_effects"){ @@ -216,6 +218,7 @@ class Information { * Loop over distributions and set links to distribution expected value if distribution is a data type. */ void setup_data(){ + FIMS_INFO_LOG("segment"); for(density_components_iterator it = this->density_components.begin(); it!= this->density_components.end(); ++it){ std::shared_ptr > d = (*it).second; if(d->input_type == "data"){ @@ -256,12 +259,13 @@ class Information { << std::endl; INFO_LOG << "Initializing fleet objects for " << this->fleets.size() << " fleets." << std::endl; + FIMS_INFO_LOG("segment"); for (fleet_iterator it = this->fleets.begin(); it != this->fleets.end(); ++it) { std::shared_ptr > f = (*it).second; INFO_LOG << "Initializing fleet " << f->id << "." << std::endl; - + FIMS_INFO_LOG("segment"); f->Initialize(f->nyears, f->nages); INFO_LOG << "Checking for available fleet selectivity pattern." @@ -274,18 +278,20 @@ class Information { sel_id); // if find, set it, otherwise invalid INFO_LOG << "Input fleet selectivity pattern id = " << sel_id << "." << std::endl; - + FIMS_INFO_LOG("segment"); if (it != this->selectivity_models.end()) { f->selectivity = (*it).second; // elements in container held in pair // (first is id, second is object - // shared pointer to distribution) INFO_LOG << "Selectivity successfully set." << std::endl; + FIMS_INFO_LOG("segment"); } else { valid_model = false; + FIMS_INFO_LOG("segment"); ERROR_LOG << "Error: Expected selectivity pattern not defined for fleet " << f->id << ", selectivity pattern " << sel_id << std::endl; - exit(1); + //exit(1); } } else { @@ -293,39 +299,43 @@ class Information { ERROR_LOG << "Error: No selectivity pattern defined for fleet " << f->id << ". FIMS requires selectivity be defined for all fleets." << std::endl; - exit(1); + //exit(1); } // end set selectivity } INFO_LOG << "Expecting to import " << this->data_objects.size() << " data objects." << std::endl; + FIMS_INFO_LOG("segment"); for(density_components_iterator it = this->density_components.begin(); it!= this->density_components.end(); ++it){ std::shared_ptr > d = (*it).second; INFO_LOG << "Checking for available density components data objects." << std::endl; + FIMS_INFO_LOG("segment"); //set data objects if distribution is a data type if(d->input_type == "data"){ if(d->observed_data_id_m != -999){ uint32_t observed_data_id = static_cast(d->observed_data_id_m); data_iterator it = this->data_objects.find(observed_data_id); INFO_LOG << "Input data id = " << observed_data_id << "." << std::endl; - + FIMS_INFO_LOG("segment"); if (it != this->data_objects.end()) { d->observed_values = (*it).second; INFO_LOG << "Data for density component, " << d->id << " successfully set." << std::endl; DATA_LOG << "" << std::endl; + FIMS_INFO_LOG("segment"); } else { valid_model = false; ERROR_LOG << "Error: Expected data observations not defined for density component " << d->id << ", observed data " << observed_data_id << std::endl; - exit(1); + FIMS_INFO_LOG("segment"); + //exit(1); } } else { valid_model = false; ERROR_LOG << "Error: No data input for density " << d->id << std::endl; - exit(1); + //exit(1); } } // end set data @@ -335,6 +345,7 @@ class Information { << this->populations.size() << " populations." << std::endl; for (population_iterator it = this->populations.begin(); it != this->populations.end(); ++it) { + FIMS_INFO_LOG("segment"); std::shared_ptr > p = (*it).second; INFO_LOG << "Setting up links from population " << p->id @@ -349,6 +360,7 @@ class Information { // from population to fleets? // any shared member in p (population is pushed into fleets) p->fleets.push_back(f); + FIMS_INFO_LOG("segment"); INFO_LOG << f->id << " " << std::flush; } INFO_LOG << "]" << std::endl; @@ -359,6 +371,7 @@ class Information { INFO_LOG << "Checking for available recruitment function." << std::endl; // set recruitment if (p->recruitment_id != -999) { + FIMS_INFO_LOG("segment"); uint32_t recruitment_uint = static_cast(p->recruitment_id); recruitment_models_iterator it = this->recruitment_models.find(recruitment_uint); @@ -374,7 +387,7 @@ class Information { "population " << p->id << ", recruitment function " << recruitment_uint << std::endl; - exit(1); + //exit(1); } } else { @@ -384,7 +397,7 @@ class Information { << ". FIMS requires recruitment functions be defined for all " "populations." << std::endl; - exit(1); + //exit(1); } INFO_LOG << "Checking for available growth function." << std::endl; @@ -393,6 +406,7 @@ class Information { uint32_t growth_uint = static_cast(p->growth_id); growth_models_iterator it = this->growth_models.find( growth_uint); // growth_models is specified in information.hpp + FIMS_INFO_LOG("segment"); // and used in rcpp // at the head of information.hpp; are the // dimensions of ages defined in rcpp or where? @@ -404,10 +418,11 @@ class Information { INFO_LOG << "Growth function successfully set." << std::endl; } else { valid_model = false; + FIMS_INFO_LOG("segment"); ERROR_LOG << "Error: Expected growth function not defined for population " << p->id << ", growth function " << growth_uint << std::endl; - exit(1); + //exit(1); } } else { @@ -417,12 +432,13 @@ class Information { << ". FIMS requires growth functions be defined for all " "populations." << std::endl; - exit(1); + //exit(1); } INFO_LOG << "Checking for available maturity function." << std::endl; // set maturity if (p->maturity_id != -999) { + FIMS_INFO_LOG("segment"); uint32_t maturity_uint = static_cast(p->maturity_id); maturity_models_iterator it = this->maturity_models.find( maturity_uint); // >maturity_models is specified in @@ -436,9 +452,9 @@ class Information { ERROR_LOG << "Error: Expected maturity function not defined for population " << p->id << ", maturity function " << maturity_uint << std::endl; - exit(1); + //exit(1); } - + FIMS_INFO_LOG("segment"); } else { valid_model = false; ERROR_LOG << "Error: No maturity function defined for population " @@ -446,13 +462,14 @@ class Information { << ". FIMS requires maturity functions be defined for all " "populations." << std::endl; - exit(1); + //exit(1); } + FIMS_INFO_LOG("segment"); INFO_LOG << "Completed initialization for population " << p->id << "." << std::endl; } INFO_LOG << "Completed initialization of all populations." << std::endl; - + FIMS_INFO_LOG("segment"); //setup priors, random effect, and data density components setup_priors(); diff --git a/inst/include/distributions/functors/lognormal_lpdf.hpp b/inst/include/distributions/functors/lognormal_lpdf.hpp index 46fbce65..f845af39 100644 --- a/inst/include/distributions/functors/lognormal_lpdf.hpp +++ b/inst/include/distributions/functors/lognormal_lpdf.hpp @@ -88,13 +88,13 @@ namespace fims_distributions if(this->input_type == "data"){ if(this->observed_values->at(i) != this->observed_values->na_value){ // this->lpdf_vec[i] = this->keep[i] * -dnorm(log(this->observed_values->at(i)), logmu[i], logsd[i], true) - log(this->observed_values->->at(i)); - this->lpdf_vec[i] = dnorm(log(this->observed_values->at(i)), logmu[i], logsd[i], true) - log(this->observed_values->at(i)); + this->lpdf_vec[i] = fims_math::dnorm(log(this->observed_values->at(i)), logmu[i], logsd[i], true) - log(this->observed_values->at(i)); } else { this->lpdf_vec[i] = 0; MODEL_LOG << "lpdf_vec for obs " << i << " is: " << this->lpdf_vec[i] <lpdf_vec[i] = dnorm(log(this->x[i]), logmu[i], logsd[i], true); + this->lpdf_vec[i] = fims_math::dnorm(log(this->x[i]), logmu[i], logsd[i], true); } lpdf += this->lpdf_vec[i]; diff --git a/inst/include/distributions/functors/multinomial_lpmf.hpp b/inst/include/distributions/functors/multinomial_lpmf.hpp index 87b5035b..a380ec13 100644 --- a/inst/include/distributions/functors/multinomial_lpmf.hpp +++ b/inst/include/distributions/functors/multinomial_lpmf.hpp @@ -53,7 +53,7 @@ namespace fims_distributions dims[1] = this->observed_values->get_jmax(); Type lpdf = 0.0; /**< total log probability mass contribution of the distribution */ - this->lpdf_vec.resize(dims[0]*dims[1]); + this->lpdf_vec.resize(dims[0]);//*dims[1]); std::fill(this->lpdf_vec.begin(), this->lpdf_vec.end(), 0); @@ -88,7 +88,7 @@ namespace fims_distributions } } - this->lpdf_vec[i] = dmultinom((vector)x_vector, (vector)prob_vector, true); + this->lpdf_vec[i] = fims_math::dmultinom(x_vector, prob_vector, true); lpdf += this->lpdf_vec[i]; /* if (this->simulate_flag) diff --git a/inst/include/distributions/functors/normal_lpdf.hpp b/inst/include/distributions/functors/normal_lpdf.hpp index afa24df8..60aef4c2 100644 --- a/inst/include/distributions/functors/normal_lpdf.hpp +++ b/inst/include/distributions/functors/normal_lpdf.hpp @@ -83,7 +83,7 @@ struct NormalLPDF : public DensityComponentBase { if(this->input_type == "data"){ if(this->observed_values->at(i) != this->observed_values->na_value){ // this->lpdf_vec[i] = this->keep[i] * -dnorm(this->observed_values->at(i), mu[i], sd[i], true); - this->lpdf_vec[i] = dnorm(this->observed_values->at(i), mu[i], sd[i], true); + this->lpdf_vec[i] = fims_math::dnorm(this->observed_values->at(i), mu[i], sd[i], true); DISTRIBUTIONS_LOG << "obsered_values " << i << " is: " << this->observed_values->at(i) << std::endl; DISTRIBUTIONS_LOG << "mu " << i << " is: " << mu[i] << std::endl; DISTRIBUTIONS_LOG << "sd " << i << " is: " << sd[i] << std::endl; @@ -93,7 +93,7 @@ struct NormalLPDF : public DensityComponentBase { } } else { - this->lpdf_vec[i] = dnorm(this->x[i], mu[i], sd[i], true); + this->lpdf_vec[i] = fims_math::dnorm(this->x[i], mu[i], sd[i], true); DISTRIBUTIONS_LOG << "x " << i << " is: " << this->x[i] << std::endl; DISTRIBUTIONS_LOG << "mu " << i << " is: " << mu[i] << std::endl; DISTRIBUTIONS_LOG << "sd " << i << " is: " << sd[i] << std::endl; diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index ac287bf6..755899f2 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/inst/include/interface/rcpp/rcpp_interface.hpp @@ -34,6 +34,7 @@ bool FIMS_finalized = false; * @brief Create the TMB model object and add interface objects to it. */ bool CreateTMBModel() { + FIMS_INFO_LOG("adding FIMS objects to TMB"); for (size_t i = 0; i < FIMSRcppInterfaceBase::fims_interface_objects.size(); i++) { FIMSRcppInterfaceBase::fims_interface_objects[i]->add_to_fims_tmb(); @@ -143,11 +144,14 @@ std::string ToJSON() { ss << "\"objective_function_value\": " << FIMS_function_value << ",\n"; ss << "\"max_gradient_component\": " << FIMS_mgc_value << ",\n"; ss << "\"final_gradient\": ["; - for (size_t i = 0; i < FIMS_function_gradient.size() - 1; i++) { - ss << FIMS_function_gradient[i] << ", "; + if (FIMS_function_gradient.size() > 0) { + for (size_t i = 0; i < FIMS_function_gradient.size() - 1; i++) { + ss << FIMS_function_gradient[i] << ", "; + } + ss << FIMS_function_gradient[FIMS_function_gradient.size() - 1] << "],\n"; + } else { + ss << "],"; } - ss << FIMS_function_gradient[FIMS_function_gradient.size() - 1] << "],\n"; - size_t length = FIMSRcppInterfaceBase::fims_interface_objects.size(); for (size_t i = 0; i < length - 1; i++) { ss << FIMSRcppInterfaceBase::fims_interface_objects[i]->to_json() << ",\n"; @@ -218,10 +222,10 @@ void clear_info_log() { * Clears the contents of fims log file. */ void clear_fims_log() { - FIMS_LOG.flush(); + FIMS_LOG_OLD.flush(); std::ofstream CLEAR_LOG("logs/fims.log"); CLEAR_LOG.close(); - FIMS_LOG.seekp(0); + FIMS_LOG_OLD.seekp(0); } /** @@ -571,6 +575,7 @@ void log_error(std::string log_entry) { RCPP_EXPOSED_CLASS(Parameter) RCPP_EXPOSED_CLASS(ParameterVector) + RCPP_MODULE(fims) { Rcpp::function("CreateTMBModel", &CreateTMBModel); Rcpp::function("SetFIMSFunctions", &SetFIMSFunctions); @@ -610,6 +615,7 @@ RCPP_MODULE(fims) { .constructor() .constructor() .field("value", &Parameter::value_m, "numeric parameter value") + .field("value", &Parameter::estimated_value_m, "numeric estimated parameter value") .field("min", &Parameter::min_m, "minimum parameter value") .field("max", &Parameter::max_m, "maximum parameter value") .field("id", &Parameter::id_m, "unique id for parameter class") @@ -620,7 +626,10 @@ RCPP_MODULE(fims) { .constructor() .constructor() .constructor() - .field("data", &ParameterVector::storage_m, "list where each element is a Parameter class") + .method("get", &ParameterVector::get) + .method("set", &ParameterVector::set) + .method("show", &ParameterVector::show) + // .field("data", &ParameterVector::storage_m, "list where each element is a Parameter class") .method("at", &ParameterVector::at, "returns a Parameter at the indicated position given the index argument") .method("size", &ParameterVector::size, "returns the size of the Parameter Vector") .method("resize", &ParameterVector::resize, "resizes the Parameter Vector given the provided length argument") diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp index 5f2dd922..fca8eaa8 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp @@ -19,26 +19,27 @@ * */ class FleetInterfaceBase : public FIMSRcppInterfaceBase { - public: - static uint32_t id_g; /**< static id of the FleetInterfaceBase object */ - uint32_t id; /**< local id of the FleetInterfaceBase object */ - // live objects in C++ are objects that have been created and live in memory - static std::map live_objects; /**< +public: + static uint32_t id_g; /**< static id of the FleetInterfaceBase object */ + uint32_t id; /**< local id of the FleetInterfaceBase object */ + // live objects in C++ are objects that have been created and live in memory + static std::map live_objects; /**< map relating the ID of the FleetInterfaceBase to the FleetInterfaceBase objects */ - FleetInterfaceBase() { - this->id = FleetInterfaceBase::id_g++; - /* Create instance of map: key is id and value is pointer to - FleetInterfaceBase */ - FleetInterfaceBase::live_objects[this->id] = this; - FIMSRcppInterfaceBase::fims_interface_objects.push_back(this); - } + FleetInterfaceBase() { + this->id = FleetInterfaceBase::id_g++; + /* Create instance of map: key is id and value is pointer to + FleetInterfaceBase */ + FleetInterfaceBase::live_objects[this->id] = this; + FIMSRcppInterfaceBase::fims_interface_objects.push_back(this); + } - virtual ~FleetInterfaceBase() {} + virtual ~FleetInterfaceBase() { + } - /** @brief get_id method for child fleet interface objects to inherit */ - virtual uint32_t get_id() = 0; + /** @brief get_id method for child fleet interface objects to inherit */ + virtual uint32_t get_id() = 0; }; uint32_t FleetInterfaceBase::id_g = 1; @@ -51,107 +52,267 @@ std::map FleetInterfaceBase::live_objects; * */ class FleetInterface : public FleetInterfaceBase { - int interface_selectivity_id_m = -999; /**< id of selectivity component*/ + int interface_selectivity_id_m = -999; /**< id of selectivity component*/ - public: - bool is_survey = false; /**< whether this is a survey fleet */ - int nages; /**< number of ages in the fleet data*/ - int nyears; /**< number of years in the fleet data */ - double log_q; /**< log of catchability for the fleet*/ - ParameterVector - log_Fmort; /**< log of fishing mortality rate for the fleet*/ - ParameterVector log_expected_index; /**< expected index of abundance for the survey */ - ParameterVector proportion_catch_numbers_at_age; /**< expected catch numbers at age for the fleet */ - bool estimate_q = false; /**< whether the parameter q should be estimated*/ - bool random_q = false; /**< whether q should be a random effect*/ +public: + std::string name = "NA"; + bool is_survey = false; /**< whether this is a survey fleet */ + int nages; /**< number of ages in the fleet data*/ + int nyears; /**< number of years in the fleet data */ + ParameterVector log_q; /**< log of catchability for the fleet*/ + ParameterVector + log_Fmort; /**< log of fishing mortality rate for the fleet*/ + ParameterVector log_expected_index; /**< expected index of abundance for the survey */ + ParameterVector proportion_catch_numbers_at_age; /**< expected catch numbers at age for the fleet */ + bool estimate_q = false; /**< whether the parameter q should be estimated*/ + bool random_q = false; /**< whether q should be a random effect*/ -Rcpp::NumericVector derived_cnaa; /**< derived quantity: catch numbers at age */ -Rcpp::NumericVector derived_cwaa; /**< derived quantity: catch weight at age */ -Rcpp::NumericVector derived_index; /**< derived quantity: expected index */ -Rcpp::NumericVector derived_age_composition; /**< derived quantity: expected catch */ + Rcpp::NumericVector derived_cnaa; /**< derived quantity: catch numbers at age */ + Rcpp::NumericVector derived_cwaa; /**< derived quantity: catch weight at age */ + Rcpp::NumericVector derived_index; /**< derived quantity: expected index */ + Rcpp::NumericVector derived_age_composition; /**< derived quantity: expected catch */ + FleetInterface() : FleetInterfaceBase() { + } - FleetInterface() : FleetInterfaceBase() {} + virtual ~FleetInterface() { + } - virtual ~FleetInterface() {} + /** @brief returns the id for the fleet interface */ + virtual uint32_t get_id() { + return this->id; + } - /** @brief returns the id for the fleet interface */ - virtual uint32_t get_id() { return this->id; } + /** + * @brief Set the unique id for the Selectivity object + * + * @param selectivity_id Unique id for the Selectivity object + */ + void SetSelectivity(int selectivity_id) { + interface_selectivity_id_m = selectivity_id; + } - /** - * @brief Set the unique id for the Selectivity object - * - * @param selectivity_id Unique id for the Selectivity object - */ - void SetSelectivity(int selectivity_id) { - interface_selectivity_id_m = selectivity_id; - } + virtual void finalize() { + + if (this->finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("Fleet " + fims::to_string(this->id) + " has been finalized already."); + } + + this->finalized = true; //indicate this has been called already + + std::shared_ptr > info = + fims_info::Information::GetInstance(); + + fims_info::Information::fleet_iterator it; + + + it = info->fleets.find(this->id); + + if (it == info->fleets.end()) { + FIMS_WARNING_LOG("Fleet " + fims::to_string(this->id) + " not found in Information."); + return; + } else { + + std::shared_ptr > fleet = + std::dynamic_pointer_cast >(it->second); + + + for (size_t i = 0; i < this->log_Fmort.size(); i++) { + if (this->log_Fmort[i].estimated_m) { + this->log_Fmort[i].estimated_value_m = fleet->log_Fmort[i]; + } else { + this->log_Fmort[i].estimated_value_m = this->log_Fmort[i].value_m; + } + } + + for (size_t i = 0; i < this->log_q.size(); i++) { + if (this->log_q[i].estimated_m) { + this->log_q[i].estimated_value_m = fleet->log_q[i]; + } else { + this->log_q[i].estimated_value_m = this->log_q[i].value_m; + } + } + + this->derived_cnaa = Rcpp::NumericVector(fleet->catch_numbers_at_age.size()); + for (size_t i = 0; i < this->derived_cnaa.size(); i++) { + this->derived_cnaa[i] = fleet->catch_numbers_at_age[i]; + } + + this->derived_cwaa = Rcpp::NumericVector(fleet->catch_weight_at_age.size()); + for (size_t i = 0; i < this->derived_cwaa.size(); i++) { + this->derived_cwaa[i] = fleet->catch_weight_at_age[i]; + } + + this->derived_age_composition = Rcpp::NumericVector(fleet->proportion_catch_numbers_at_age.size()); + for (size_t i = 0; i < this->derived_age_composition.size(); i++) { + this->derived_age_composition[i] = fleet->proportion_catch_numbers_at_age[i]; + } + + this->derived_index = Rcpp::NumericVector(fleet->expected_index.size()); + for (size_t i = 0; i < this->derived_index.size(); i++) { + this->derived_index[i] = fleet->expected_index[i]; + } + + } -#ifdef TMB_MODEL - template - bool add_to_fims_tmb_internal() { - std::shared_ptr > info = - fims_info::Information::GetInstance(); - - std::shared_ptr > fleet = - std::make_shared >(); - - // set relative info - fleet->id = this->id; - fleet->is_survey = this->is_survey; - fleet->nages = this->nages; - fleet->nyears = this->nyears; - fleet->fleet_selectivity_id_m = interface_selectivity_id_m; - - fleet->log_q = this->log_q; - if (this->estimate_q) { - info->RegisterParameterName("log_q"); - if (this->random_q) { - info->RegisterRandomEffect(fleet->log_q); - } else { - info->RegisterParameter(fleet->log_q); - } } - fleet->log_Fmort.resize(this->log_Fmort.size()); - for (size_t i = 0; i < log_Fmort.size(); i++) { - fleet->log_Fmort[i] = this->log_Fmort[i].value_m; + virtual std::string to_json() { + std::stringstream ss; + + ss << "\"module\" : {\n"; + ss << " \"name\" : \"Fleet\",\n"; + + ss << " \"type\" : \"fleet\",\n"; + ss << " \"tag\" : \"" << this->name << "\",\n"; + ss << " \"id\": " << this->id << ",\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"log_Fmort\",\n"; + ss << " \"id\":" << this->log_Fmort.id_m << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\": " << this->log_Fmort << "\n},\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"log_Fmort\",\n"; + ss << " \"id\":" << this->log_q.id_m << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\": " << this->log_q << "\n},\n"; - if (this->log_Fmort[i].estimated_m) { - info->RegisterParameterName("log_Fmort"); - if (this->log_Fmort[i].is_random_effect_m) { - info->RegisterRandomEffect(fleet->log_Fmort[i]); + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"cnaa\",\n"; + ss << " \"values\":["; + if (this->derived_cnaa.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->derived_cnaa.size() - 1; i++) { + ss << this->derived_cnaa[i] << ", "; + } + ss << this->derived_cnaa[this->derived_cnaa.size() - 1] << "]\n"; + } + ss << " },\n"; + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"cwaa\",\n"; + ss << " \"values\":["; + if (this->derived_cwaa.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->derived_cwaa.size() - 1; i++) { + ss << this->derived_cwaa[i] << ", "; + } + ss << this->derived_cwaa[this->derived_cwaa.size() - 1] << "]\n"; + } + ss << " },\n"; + + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"age_composition \",\n"; + ss << " \"values\":["; + if (this->derived_age_composition.size() == 0) { + ss << "]\n"; } else { - info->RegisterParameter(fleet->log_Fmort[i]); + for (size_t i = 0; i < this->derived_age_composition.size() - 1; i++) { + ss << this->derived_age_composition[i] << ", "; + } + ss << this->derived_age_composition[this->derived_age_composition.size() - 1] << "]\n"; } - } + ss << " },\n"; + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"index \",\n"; + ss << " \"values\":["; + if (this->derived_index.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->derived_index.size() - 1; i++) { + ss << this->derived_index[i] << ", "; + } + ss << this->derived_index[this->derived_index.size() - 1] << "]\n"; + } + ss << " },\n"; + + return ss.str(); + } - //add to variable_map - info->variable_map[this->log_Fmort.id_m] = &(fleet)->log_Fmort; - //exp_catch - fleet->log_expected_index.resize(nyears); // assume index is for all ages. - info->variable_map[this->log_expected_index.id_m] = &(fleet)->log_expected_index; - fleet->proportion_catch_numbers_at_age.resize(nyears * nages); - info->variable_map[this->proportion_catch_numbers_at_age.id_m] = &(fleet)->proportion_catch_numbers_at_age; - // add to Information - info->fleets[fleet->id] = fleet; +#ifdef TMB_MODEL + + template + bool add_to_fims_tmb_internal() { + std::shared_ptr > info = + fims_info::Information::GetInstance(); + + std::shared_ptr > fleet = + std::make_shared >(); + + // set relative info + fleet->id = this->id; + fleet->is_survey = this->is_survey; + fleet->nages = this->nages; + fleet->nyears = this->nyears; + fleet->fleet_selectivity_id_m = interface_selectivity_id_m; + + fleet->log_q.resize(this->log_q.size()); + for (size_t i = 0; i < this->log_q.size(); i++) { + fleet->log_q[i] = this->log_q[i].value_m; + + if (this->log_q[i].estimated_m) { + info->RegisterParameterName("log_q"); + if (this->log_q[i].is_random_effect_m) { + info->RegisterRandomEffect(fleet->log_q[i]); + } else { + info->RegisterParameter(fleet->log_q[i]); + } + } + } + - return true; - } + fleet->log_Fmort.resize(this->log_Fmort.size()); + for (size_t i = 0; i < log_Fmort.size(); i++) { + fleet->log_Fmort[i] = this->log_Fmort[i].value_m; + + if (this->log_Fmort[i].estimated_m) { + info->RegisterParameterName("log_Fmort"); + if (this->log_Fmort[i].is_random_effect_m) { + info->RegisterRandomEffect(fleet->log_Fmort[i]); + } else { + info->RegisterParameter(fleet->log_Fmort[i]); + } + } + } + //add to variable_map + info->variable_map[this->log_Fmort.id_m] = &(fleet)->log_Fmort; - /** @brief this adds the values to the TMB model object */ - virtual bool add_to_fims_tmb() { - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); + //exp_catch + fleet->log_expected_index.resize(nyears); // assume index is for all ages. + info->variable_map[this->log_expected_index.id_m] = &(fleet)->log_expected_index; + fleet->proportion_catch_numbers_at_age.resize(nyears * nages); + info->variable_map[this->proportion_catch_numbers_at_age.id_m] = &(fleet)->proportion_catch_numbers_at_age; - return true; - } + + // add to Information + info->fleets[fleet->id] = fleet; + + return true; + } + + /** @brief this adds the values to the TMB model object */ + virtual bool add_to_fims_tmb() { + FIMS_INFO_LOG("adding Fleet object to TMB"); + + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + + return true; + } #endif }; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp index 192c11a9..1647a71a 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_interface_base.hpp @@ -9,6 +9,7 @@ #ifndef FIMS_INTERFACE_RCPP_RCPP_OBJECTS_RCPP_INTERFACE_BASE_HPP #define FIMS_INTERFACE_RCPP_RCPP_OBJECTS_RCPP_INTERFACE_BASE_HPP +#include #include #include @@ -16,6 +17,8 @@ #include "../../../common/information.hpp" #include "../../interface.hpp" + + #define RCPP_NO_SUGAR #include @@ -27,7 +30,8 @@ class Parameter { public: static uint32_t id_g; /**< global id of the parameter */ uint32_t id_m; /**< id of the parameter */ - double value_m; /**< initial value of the parameter */ + double value_m = 0.0; /**< initial value of the parameter */ + double estimated_value_m = 0.0; /**< estimated value of the parameter */ double min_m = -std::numeric_limits::infinity(); /**< min value of the parameter; default is negative infinity*/ double max_m = @@ -48,6 +52,31 @@ class Parameter { : id_m(Parameter::id_g++), value_m(value), min_m(min), max_m(max), estimated_m(estimated) { } + Parameter(const Parameter& other) : + id_m(other.id_m), value_m(other.value_m), + estimated_value_m(other.estimated_value_m), + min_m(other.min_m), max_m(other.max_m), + is_random_effect_m(other.is_random_effect_m), + estimated_m(other.estimated_m), + random_m(other.random_m) { + } + + Parameter& operator=(const Parameter& right) { + // Check for self-assignment! + if (this == &right) // Same object? + return *this; // Yes, so skip assignment, and just return *this. + this->id_m = right.id_m; + this->value_m = right.value_m; + this->estimated_m = right.estimated_m; + this->random_m = right.random_m; + this->min_m = right.min_m; + this->max_m = right.max_m; + this->is_random_effect_m = right.is_random_effect_m; + return *this; + } + + + /** * @brief Constructor for initializing Parameter. * @details Inputs include value. @@ -70,7 +99,9 @@ class Parameter { uint32_t Parameter::id_g = 0; std::ostream& operator<<(std::ostream& out, const Parameter& p) { - out << "Parameter:{" << "id:" << p.id_m << ",value:" << p.value_m << ",min:" << p.min_m << ",max:" << p.max_m << ",estimated:" << p.estimated_m << "}"; + out << "Parameter:{" << "id:" << p.id_m << ",\nvalue:" << p.value_m + << ",\nestimated_value:" << p.estimated_value_m << ",\nmin:" + << p.min_m << ",\nmax:" << p.max_m << ",\nestimated:" << p.estimated_m << "\n}"; return out; } @@ -81,7 +112,7 @@ std::ostream& operator<<(std::ostream& out, const Parameter& p) { class ParameterVector { public: static uint32_t id_g; /**< global identifier*/ - Rcpp::List storage_m; /**< list of parameter objects*/ + std::shared_ptr > storage_m; /**< list of parameter objects*/ uint32_t id_m; /**< unique identifier*/ /** @@ -89,18 +120,25 @@ class ParameterVector { */ ParameterVector() { this->id_m = ParameterVector::id_g++; - Parameter p; - this->storage_m.push_back(Rcpp::wrap(p)); + this->storage_m = std::make_shared >(); + // Parameter p; + this->storage_m->resize(1); //push_back(Rcpp::wrap(p)); + } + + ParameterVector(const ParameterVector& other) : + storage_m(other.storage_m), id_m(other.id_m) { } /** * @brief constructor */ ParameterVector(size_t size) { + this->id_m = ParameterVector::id_g++; + this->storage_m = std::make_shared >(); + this->storage_m->resize(size); for (size_t i = 0; i < size; i++) { - Parameter p; - this->storage_m.push_back(Rcpp::wrap(p)); + storage_m->at(i) = Parameter(); } } @@ -111,18 +149,21 @@ class ParameterVector { */ ParameterVector(Rcpp::NumericVector x, size_t size) { this->id_m = ParameterVector::id_g++; + this->storage_m = std::make_shared >(); + this->resize(size); for (size_t i = 0; i < size; i++) { - Parameter p = x[i]; - this->storage_m.push_back(Rcpp::wrap(p)); + storage_m->at(i).value_m = x[i]; } } - - ParameterVector(const fims::Vector& v){ - this->id_m = ParameterVector::id_g++; + + ParameterVector(const fims::Vector& v) { + this->id_m = ParameterVector::id_g++; + this->storage_m = std::make_shared >(); + this->storage_m->resize(v.size()); for (size_t i = 0; i < v.size(); i++) { - Parameter p = v[i]; - this->storage_m.push_back(Rcpp::wrap(p)); + storage_m->at(i).value_m = v[i]; } + } /** @@ -136,8 +177,9 @@ class ParameterVector { * @brief Accessor. First index starts is zero. * @param pos return a Parameter at position "pos". */ - inline Parameter operator[](R_xlen_t pos) { - return this->storage_m[pos]; + inline Parameter& operator[](size_t pos) { + + return this->storage_m->at(pos); } /** @@ -145,18 +187,39 @@ class ParameterVector { * @param pos return a Parameter at position "pos". */ SEXP at(R_xlen_t pos) { - if (pos == 0 || pos > this->storage_m.size()) { - Rcpp::Rcout << "Index out of range.\n"; + if (pos == 0 || pos > this->storage_m->size()) { + Rcpp::Rcout << "ParameterVector: Index out of range.\n"; + FIMS_ERROR_LOG(fims::to_string(pos) + "!<" + fims::to_string(this->size())); return NULL; } - return this->storage_m[pos - 1]; + std::cout << "id:" << this->storage_m->at(pos - 1).id_m << "\n"; + std::cout << "value:" << this->storage_m->at(pos - 1).value_m << "\n"; + return Rcpp::wrap(this->storage_m->at(pos - 1)); + } + + /** + * @brief Internal accessor. First index is one. For calling from R. + * @param pos return a Parameter at position "pos". + */ + Parameter& get(size_t pos) { + if (pos >= this->storage_m->size()) { + std::cout << "ParameterVector: Index out of range.\n"; + throw std::invalid_argument("ParameterVector: Index out of range"); + } + + std::cout << this->storage_m->at(pos) << "\n"; + return (this->storage_m->at(pos)); + } + + void set(size_t pos, const Parameter& p) { + this->storage_m->at(pos) = p; } /** * @brief returns vector length */ size_t size() { - return this->storage_m.size(); + return this->storage_m->size(); } /** @@ -164,24 +227,7 @@ class ParameterVector { * @param size new length of vector to be resized */ void resize(size_t size) { - size_t n = this->storage_m.size(); - - if (size > n) { - size_t m = size - n; - - for (size_t i = 0; i < m; i++) { - Parameter p; - this->storage_m.push_back(Rcpp::wrap(p)); - } - } else if (n > size) { - size_t m = size; - Rcpp::List l(m); - for (size_t i = 0; i < m; i++) { - l[i] = this->storage_m[i]; - } - this->storage_m = l; - } - + this->storage_m->resize(size); } /** @@ -190,10 +236,9 @@ class ParameterVector { * @param estimable Boolean; if true, all parameters are set to be estimated in the model */ void set_all_estimable(bool estimable) { - for (R_xlen_t i = 0; i < this->storage_m.size(); i++) { - Parameter p = Rcpp::as(this->storage_m[i]); - p.estimated_m = estimable; - this->storage_m[i] = Rcpp::wrap(p); + Rcpp::Rcout << this->storage_m->data() << "\n"; + for (size_t i = 0; i < this->storage_m->size(); i++) { + storage_m->at(i).estimated_m = estimable; } } @@ -203,10 +248,8 @@ class ParameterVector { * @param random Boolean; if true, all parameters are set to be random effects in the model */ void set_all_random(bool random) { - for (R_xlen_t i = 0; i < this->storage_m.size(); i++) { - Parameter p = Rcpp::as(this->storage_m[i]); - p.random_m = random; - this->storage_m[i] = Rcpp::wrap(p); + for (size_t i = 0; i < this->storage_m->size(); i++) { + storage_m->at(i).random_m = random; } } @@ -216,10 +259,8 @@ class ParameterVector { * @param value The value to be assigned */ void fill(double value) { - for (R_xlen_t i = 0; i < this->storage_m.size(); i++) { - Parameter p = Rcpp::as(this->storage_m[i]); - p.value_m = value; - this->storage_m[i] = Rcpp::wrap(p); + for (size_t i = 0; i < this->storage_m->size(); i++) { + storage_m->at(i).value_m = value; } } @@ -229,10 +270,8 @@ class ParameterVector { * @param value The value to be assigned */ void fill_min(double value) { - for (int i = 0; i < this->storage_m.size(); i++) { - Parameter p = Rcpp::as(this->storage_m[i]); - p.min_m = value; - this->storage_m[i] = Rcpp::wrap(p); + for (size_t i = 0; i < this->storage_m->size(); i++) { + storage_m->at(i).min_m = value; } } @@ -242,23 +281,29 @@ class ParameterVector { * @param value The value to be assigned */ void fill_max(double value) { - for (int i = 0; i < this->storage_m.size(); i++) { - Parameter p = Rcpp::as(this->storage_m[i]); - p.max_m = value; - this->storage_m[i] = Rcpp::wrap(p); + for (size_t i = 0; i < this->storage_m->size(); i++) { + storage_m->at(i).max_m = value; + } + } + + void show() { + std::cout << this->storage_m->data() << "\n"; + + for (size_t i = 0; i < this->storage_m->size(); i++) { + std::cout << storage_m->at(i) << " "; } } }; uint32_t ParameterVector::id_g = 0; -std::ostream& operator<<(std::ostream& out, ParameterVector& v) { - out<<"ParameterVector:{"; +std::ostream& operator<<(std::ostream& out, ParameterVector& v) { + out << "["; size_t size = v.size(); - for(size_t i = 0; i < size-1; i++){ - out< MaturityInterfaceBase::live_objects; */ class LogisticMaturityInterface : public MaturityInterfaceBase { public: - Parameter + ParameterVector inflection_point; /**< the index value at which the response reaches .5 */ - Parameter slope; /**< the width of the curve at the inflection_point */ - - double estimated_inflection_point; /**< estimmated result of the index value at which the response reaches .5 */ - double estimated_slope; /**< estimmated result of the width of the curve at the inflection_point */ + ParameterVector slope; /**< the width of the curve at the inflection_point */ LogisticMaturityInterface() : MaturityInterfaceBase() { } @@ -86,36 +83,78 @@ class LogisticMaturityInterface : public MaturityInterfaceBase { virtual double evaluate(double x) { fims_popdy::LogisticMaturity LogisticMat; LogisticMat.inflection_point.resize(1); - LogisticMat.inflection_point[0] = this->inflection_point.value_m; + LogisticMat.inflection_point[0] = this->inflection_point[0].value_m; LogisticMat.slope.resize(1); - LogisticMat.slope[0] = this->slope.value_m; + LogisticMat.slope[0] = this->slope[0].value_m; return LogisticMat.evaluate(x); } + virtual void finalize() { + + if (this->finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("Logistic Maturity " + fims::to_string(this->id) + " has been finalized already."); + } + + this->finalized = true; //indicate this has been called already + + + std::shared_ptr > info = + fims_info::Information::GetInstance(); + + fims_info::Information::maturity_models_iterator it; + + + //search for maturity in Information + it = info->maturity_models.find(this->id); + //if not found, just return + if (it == info->maturity_models.end()) { + FIMS_WARNING_LOG("Logistic Maturity " + fims::to_string(this->id) + " not found in Information."); + return; + } else { + std::shared_ptr > mat = + std::dynamic_pointer_cast >(it->second); + + for (size_t i = 0; i < inflection_point.size(); i++) { + if (this->inflection_point[i].estimated_m) { + this->inflection_point[i].estimated_value_m = mat->inflection_point[i]; + } else { + this->inflection_point[i].estimated_value_m = this->inflection_point[i].value_m; + } + + } + + for (size_t i = 0; i < slope.size(); i++) { + if (this->slope[i].estimated_m) { + this->slope[i].estimated_value_m = mat->slope[i]; + } else { + this->slope[i].estimated_value_m = this->slope[i].value_m; + } + + } + + + } + } + virtual std::string to_json() { std::stringstream ss; ss << "\"module\" : {\n"; ss << " \"name\": \"maturity\",\n"; - ss << " \"type\": \"Logistic\",\n"; + ss << " \"type\": \"logistic\",\n"; ss << " \"id\": " << this->id << ",\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"inflection_point\",\n"; ss << " \"id\":" << this->inflection_point.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->inflection_point.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_inflection_point << ",\n"; - ss << " \"is_estimated\":" << this->inflection_point.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->inflection_point.is_random_effect_m << "\n },\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->inflection_point << ",\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"slope\",\n"; ss << " \"id\":" << this->slope.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->slope.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_slope << ",\n"; - ss << " \"is_estimated\":" << this->slope.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->slope.is_random_effect_m << "\n }\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->slope << ",\n"; ss << "}"; @@ -135,24 +174,34 @@ class LogisticMaturityInterface : public MaturityInterfaceBase { // set relative info maturity->id = this->id; - maturity->inflection_point.resize(1); - maturity->inflection_point[0] = this->inflection_point.value_m; - if (this->inflection_point.estimated_m) { - info->RegisterParameterName("maturity inflection_point"); - if (this->inflection_point.is_random_effect_m) { - info->RegisterRandomEffect(maturity->inflection_point[0]); - } else { - info->RegisterParameter(maturity->inflection_point[0]); + std::stringstream ss; + maturity->inflection_point.resize(this->inflection_point.size()); + for (size_t i = 0; i < this->inflection_point.size(); i++) { + maturity->inflection_point[i] = this->inflection_point[i].value_m; + if (this->inflection_point[i].estimated_m) { + ss.str(""); + ss << "maturity.inflection_point." << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->inflection_point[i].is_random_effect_m) { + info->RegisterRandomEffect(maturity->inflection_point[i]); + } else { + info->RegisterParameter(maturity->inflection_point[i]); + } } } - maturity->slope.resize(1); - maturity->slope[0] = this->slope.value_m; - if (this->slope.estimated_m) { - info->RegisterParameterName("maturity slope"); - if (this->slope.is_random_effect_m) { - info->RegisterRandomEffect(maturity->slope[0]); - } else { - info->RegisterParameter(maturity->slope[0]); + + maturity->slope.resize(this->slope.size()); + for (size_t i = 0; i < this->slope.size(); i++) { + maturity->slope[i] = this->slope[i].value_m; + if (this->slope[i].estimated_m) { + ss.str(""); + ss << "maturity.slope_" << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->slope[i].is_random_effect_m) { + info->RegisterRandomEffect(maturity->slope[i]); + } else { + info->RegisterParameter(maturity->slope[i]); + } } } @@ -165,6 +214,7 @@ class LogisticMaturityInterface : public MaturityInterfaceBase { /** @brief this adds the parameter values and derivatives to the TMB model * object */ virtual bool add_to_fims_tmb() { + FIMS_INFO_LOG("adding Maturity object to TMB"); this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp index 144e72ac..51196d41 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp @@ -123,8 +123,7 @@ class PopulationInterface : public PopulationInterfaceBase { fims_popdy::Population population; return population.Evaluate(); } - - + /** * @brief finalize function. Extracts derived quantities back to * the Rcpp interface object from the Information object. @@ -236,85 +235,19 @@ class PopulationInterface : public PopulationInterfaceBase { ss << " \"name\": \"log_M\",\n"; ss << " \"id\":" << -999 << ",\n"; ss << " \"type\": \"vector\",\n"; - ss << " \"values\":["; - - if (this->log_M.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->log_M.size() - 1; i++) { - ss << this->log_M[i] << ", "; - } - ss << this->log_M[this->log_M.size() - 1] << "],\n"; - } - - ss << " \"estimated_values\":["; - if (this->estimated_log_M.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->estimated_log_M.size() - 1; i++) { - ss << this->estimated_log_M[i] << ", "; - } - ss << this->estimated_log_M[this->estimated_log_M.size() - 1] << "],\n"; - } - - ss << " \"is_estimated\":" << this->estimated_log_M << ",\n"; - ss << " \"is_random_effect\":" << 0 << "\n },\n"; + ss << " \"values\": " << this->log_M << "\n},\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"log_init_naa\",\n"; ss << " \"id\":" << -999 << ",\n"; ss << " \"type\": \"vector\",\n"; - ss << " \"values\":["; - - if (this->log_init_naa.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->log_init_naa.size() - 1; i++) { - ss << this->log_init_naa[i] << ", "; - } - ss << this->log_init_naa[this->log_init_naa.size() - 1] << "],\n"; - } - - ss << " \"estimated_values\":["; - if (this->estimated_log_init_naa.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->estimated_log_init_naa.size() - 1; i++) { - ss << this->estimated_log_init_naa[i] << ", "; - } - ss << this->estimated_log_init_naa[this->estimated_log_init_naa.size() - 1] << "],\n"; - } - - ss << " \"is_estimated\":" << this->estimated_log_init_naa << ",\n"; - ss << " \"is_random_effect\":" << 0 << "\n },\n"; + ss << " \"values\":" << this->log_init_naa << " \n},\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"proportion_female\",\n"; ss << " \"id\":" << -999 << ",\n"; ss << " \"type\": \"vector\",\n"; - ss << " \"values\":["; - - if (this->proportion_female.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->proportion_female.size() - 1; i++) { - ss << this->proportion_female[i] << ", "; - } - ss << this->proportion_female[this->proportion_female.size() - 1] << "],\n"; - } - - ss << " \"estimated_values\":["; - if (this->estimated_proportion_female.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->estimated_proportion_female.size() - 1; i++) { - ss << this->estimated_proportion_female[i] << ", "; - } - ss << this->estimated_proportion_female[this->estimated_proportion_female.size() - 1] << "],\n"; - } - - ss << " \"is_estimated\":" << this->estimate_prop_female << ",\n"; - ss << " \"is_random_effect\":" << 0 << "\n},\n"; + ss << " \"values\":" << this->proportion_female << "\n},\n"; ss << " \"derived_quantity\": {\n"; @@ -441,6 +374,7 @@ class PopulationInterface : public PopulationInterfaceBase { /** @brief this adds the parameter values and derivatives to the TMB model * object */ virtual bool add_to_fims_tmb() { + FIMS_INFO_LOG("adding Population object to TMB"); this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp index 78a07599..a41ea03a 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp @@ -65,9 +65,9 @@ RecruitmentInterfaceBase::live_objects; */ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { public: - Parameter logit_steep; /**< steepness or the productivity of the stock*/ - Parameter log_rzero; /**< recruitment at unfished biomass */ - Parameter + ParameterVector logit_steep; /**< steepness or the productivity of the stock*/ + ParameterVector log_rzero; /**< recruitment at unfished biomass */ + ParameterVector log_sigma_recruit; /**< the log of the stock recruit standard deviation */ ParameterVector log_devs; /**< log recruitment deviations*/ bool estimate_log_devs = false; /**< boolean describing whether to estimate */ @@ -90,14 +90,14 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { virtual double evaluate(double spawners, double ssbzero) { fims_popdy::SRBevertonHolt BevHolt; BevHolt.logit_steep.resize(1); - BevHolt.logit_steep[0] = this->logit_steep.value_m; - if (this->logit_steep.value_m == 1.0) { + BevHolt.logit_steep[0] = this->logit_steep[0].value_m; + if (this->logit_steep[0].value_m == 1.0) { warning( "Steepness is subject to a logit transformation, so its value is " "0.7848469. Fixing it at 1.0 is not currently possible."); } BevHolt.log_rzero.resize(1); - BevHolt.log_rzero[0] = this->log_rzero.value_m; + BevHolt.log_rzero[0] = this->log_rzero[0].value_m; return BevHolt.evaluate(spawners, ssbzero); } @@ -115,16 +115,6 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { fims_info::Information::GetInstance(); - this->estimated_logit_steep = this->logit_steep.value_m; - this->estimated_log_rzero = this->log_rzero.value_m; - this->estimated_log_sigma_recruit = this->log_sigma_recruit.value_m; - - this->estimated_log_devs = Rcpp::NumericVector(this->log_devs.size()); - - for (size_t i = 0; i < this->estimated_log_devs.size(); i++) { - this->estimated_log_devs[i] = log_devs[i].value_m; - } - fims_info::Information::recruitment_models_iterator it; it = info->recruitment_models.find(this->id); @@ -137,22 +127,36 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { std::shared_ptr > recr = std::dynamic_pointer_cast >(it->second); - if (this->logit_steep.estimated_m) { - this->estimated_logit_steep = recr->logit_steep[0]; + for (size_t i = 0; i < this->logit_steep.size(); i++) { + if (this->logit_steep[i].estimated_m) { + this->logit_steep[i].estimated_value_m = recr->logit_steep[i]; + } else { + this->logit_steep[i].estimated_value_m = this->logit_steep[i].value_m; + } } - if (log_rzero.estimated_m) { - this->estimated_log_rzero = recr->log_rzero[0]; + for (size_t i = 0; i < log_rzero.size(); i++) { + if (log_rzero[i].estimated_m) { + this->log_rzero[i].estimated_value_m = recr->log_rzero[i]; + } else { + this->log_rzero[i].estimated_value_m = this->log_rzero[i].value_m; + } } - if (log_sigma_recruit.estimated_m) { - this->estimated_log_sigma_recruit = recr->log_sigma_recruit[0]; - } + for (size_t i = 0; i < this->log_sigma_recruit.size(); i++) { + if (log_sigma_recruit[i].estimated_m) { + this->log_sigma_recruit[i].estimated_value_m = recr->log_sigma_recruit[i]; + } else { + this->log_sigma_recruit[i].estimated_value_m = this->log_sigma_recruit[i].value_m; + } + } - if (this->estimate_log_devs) { - for (size_t i = 0; i < this->estimated_log_devs.size(); i++) { - this->estimated_log_devs[i] = recr->log_recruit_devs[i]; + for (size_t i = 0; i < this->estimated_log_devs.size(); i++) { + if (this->log_devs[i].estimated_m) { + this->log_devs[i].estimated_value_m = recr->log_recruit_devs[i]; + } else { + this->log_devs[i].estimated_value_m = this->log_devs[i].value_m; } } } @@ -171,60 +175,27 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { ss << " \"parameter\": {\n"; ss << " \"name\": \"logit_steep\",\n"; ss << " \"id\":" << this->logit_steep.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->logit_steep.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_logit_steep << ",\n"; - ss << " \"is_estimated\":" << this->logit_steep.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->logit_steep.is_random_effect_m << "\n },\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->logit_steep << ",\n},\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"log_rzero\",\n"; ss << " \"id\":" << this->log_rzero.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->log_rzero.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_log_rzero << ",\n"; - ss << " \"is_estimated\":" << this->log_rzero.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->log_rzero.is_random_effect_m << "\n },\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->log_rzero << ",\n },\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"log_sigma_recruit\",\n"; ss << " \"id\":" << this->log_sigma_recruit.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->log_sigma_recruit.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_log_sigma_recruit << ",\n"; - ss << " \"is_estimated\":" << this->log_sigma_recruit.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->log_sigma_recruit.is_random_effect_m << "\n },\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->log_sigma_recruit << ",\n },\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"log_devs\",\n"; - ss << " \"id\":" << -999 << ",\n"; + ss << " \"id\":" << this->log_devs.id_m << ",\n"; ss << " \"type\": \"vector\",\n"; - ss << " \"values\":["; - - if (this->log_devs.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->log_devs.size() - 1; i++) { - ss << this->log_devs[i] << ", "; - } - ss << this->log_devs[this->log_devs.size() - 1] << "],\n"; - } - - ss << " \"estimated_values\":["; - if (this->estimated_log_devs.size() == 0) { - ss << "],\n"; - } else { - for (size_t i = 0; i < this->estimated_log_devs.size() - 1; i++) { - ss << this->estimated_log_devs[i] << ", "; - } - ss << this->estimated_log_devs[this->estimated_log_devs.size() - 1] << "],\n"; - } - - ss << " \"is_estimated\":" << this->estimate_log_devs << ",\n"; - ss << " \"is_random_effect\":" << 0 << "\n }\n"; - + ss << " \"values\":" << this->log_devs << ",\n },\n"; - ss << "}"; return ss.str(); } @@ -241,55 +212,76 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { // set relative info recruitment->id = this->id; + //set logit_steep - recruitment->logit_steep.resize(1); - recruitment->logit_steep[0] = this->logit_steep.value_m; - if (this->logit_steep.estimated_m) { - info->RegisterParameterName("logit_steep"); - if (this->logit_steep.is_random_effect_m) { - info->RegisterRandomEffect(recruitment->logit_steep[0]); - } else { - info->RegisterParameter(recruitment->logit_steep[0]); + recruitment->logit_steep.resize(this->logit_steep.size()); + for (size_t i = 0; i < this->logit_steep.size(); i++) { + + recruitment->logit_steep[i] = this->logit_steep[i].value_m; + + if (this->logit_steep[i].estimated_m) { + info->RegisterParameterName("logit_steep"); + if (this->logit_steep[i].is_random_effect_m) { + info->RegisterRandomEffect(recruitment->logit_steep[i]); + } else { + info->RegisterParameter(recruitment->logit_steep[i]); + } } + } + info->variable_map[this->logit_steep.id_m] = &(recruitment)->logit_steep; + //set log_rzero - recruitment->log_rzero.resize(1); - recruitment->log_rzero[0] = this->log_rzero.value_m; - if (this->log_rzero.estimated_m) { - info->RegisterParameterName("log_rzero"); - if (this->log_rzero.is_random_effect_m) { - info->RegisterRandomEffect(recruitment->log_rzero[0]); - } else { - info->RegisterParameter(recruitment->log_rzero[0]); + recruitment->log_rzero.resize(this->log_rzero.size()); + for (size_t i = 0; i < this->log_rzero.size(); i++) { + + recruitment->log_rzero[i] = this->log_rzero[i].value_m; + + if (this->log_rzero[i].estimated_m) { + info->RegisterParameterName("log_rzero"); + if (this->log_rzero[i].is_random_effect_m) { + info->RegisterRandomEffect(recruitment->log_rzero[i]); + } else { + info->RegisterParameter(recruitment->log_rzero[i]); + } } } + info->variable_map[this->log_rzero.id_m] = &(recruitment)->log_rzero; //set log_sigma_recruit - recruitment->log_sigma_recruit.resize(1); - recruitment->log_sigma_recruit[0] = this->log_sigma_recruit.value_m; - if (this->log_sigma_recruit.estimated_m) { - info->RegisterParameterName("log_sigma_recruit"); - if (this->log_sigma_recruit.is_random_effect_m) { - info->RegisterRandomEffect(recruitment->log_sigma_recruit[0]); - } else { - info->RegisterParameter(recruitment->log_sigma_recruit[0]); + recruitment->log_sigma_recruit.resize(this->log_sigma_recruit.size()); + for (size_t i = 0; i < this->log_sigma_recruit.size(); i++) { + + recruitment->log_sigma_recruit[i] = this->log_sigma_recruit[i].value_m; + if (this->log_sigma_recruit[i].estimated_m) { + info->RegisterParameterName("log_sigma_recruit"); + if (this->log_sigma_recruit[i].is_random_effect_m) { + info->RegisterRandomEffect(recruitment->log_sigma_recruit[i]); + } else { + info->RegisterParameter(recruitment->log_sigma_recruit[i]); + } } } + info->variable_map[this->log_sigma_recruit.id_m] = &(recruitment)->log_sigma_recruit; + + //set log_recruit_devs recruitment->log_recruit_devs.resize(this->log_devs.size()); - for (size_t i = 0; i < recruitment->log_recruit_devs.size(); i++) { + for (size_t i = 0; i < this->log_devs.size(); i++) { recruitment->log_recruit_devs[i] = this->log_devs[i].value_m; - if (this->estimate_log_devs) { + if (this->log_devs[i].estimated_m) { info->RegisterParameter(recruitment->log_recruit_devs[i]); } else { recruitment->estimate_log_recruit_devs = false; } + } + info->variable_map[this->log_devs.id_m] = &(recruitment)->log_recruit_devs; // add to Information @@ -301,6 +293,8 @@ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { /** @brief this adds the parameter values and derivatives to the TMB model * object */ virtual bool add_to_fims_tmb() { + FIMS_INFO_LOG("adding Recruitment object to TMB"); + this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp index c314c966..12e23c58 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp @@ -61,12 +61,9 @@ SelectivityInterfaceBase::live_objects; */ class LogisticSelectivityInterface : public SelectivityInterfaceBase { public: - Parameter + ParameterVector inflection_point; /**< the index value at which the response reaches .5 */ - Parameter slope; /**< the width of the curve at the inflection_point */ - - double estimated_inflection_point; /**< estimmated result of the index value at which the response reaches .5 */ - double estimated_slope; /**< estimmated result of the width of the curve at the inflection_point */ + ParameterVector slope; /**< the width of the curve at the inflection_point */ LogisticSelectivityInterface() : SelectivityInterfaceBase() { } @@ -86,9 +83,9 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { virtual double evaluate(double x) { fims_popdy::LogisticSelectivity LogisticSel; LogisticSel.inflection_point.resize(1); - LogisticSel.inflection_point[0] = this->inflection_point.value_m; + LogisticSel.inflection_point[0] = this->inflection_point[0].value_m; LogisticSel.slope.resize(1); - LogisticSel.slope[0] = this->slope.value_m; + LogisticSel.slope[0] = this->slope[0].value_m; return LogisticSel.evaluate(x); } @@ -105,12 +102,9 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { fims_info::Information::GetInstance(); - //set default values as initial values - this->estimated_inflection_point = this->inflection_point.value_m; - this->estimated_slope = this->slope.value_m; - fims_info::Information::selectivity_models_iterator it; + //search for maturity in Information it = info->selectivity_models.find(this->id); //if not found, just return @@ -121,14 +115,23 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { std::shared_ptr > sel = std::dynamic_pointer_cast >(it->second); - // if the parameter was estimated, set set the estimated value. - if (this->inflection_point.estimated_m) { - this->estimated_inflection_point = sel->inflection_point[0]; + for (size_t i = 0; i < inflection_point.size(); i++) { + if (this->inflection_point[i].estimated_m) { + this->inflection_point[i].estimated_value_m = sel->inflection_point[i]; + } else { + this->inflection_point[i].estimated_value_m = this->inflection_point[i].value_m; + } } - if (this->slope.estimated_m) { - this->estimated_slope = sel->slope[0]; + for (size_t i = 0; i < slope.size(); i++) { + if (this->slope[i].estimated_m) { + this->slope[i].estimated_value_m = sel->slope[i]; + } else { + this->slope[i].estimated_value_m = this->slope[i].value_m; + } + } + } } @@ -143,20 +146,14 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { ss << " \"parameter\": {\n"; ss << " \"name\": \"inflection_point\",\n"; ss << " \"id\":" << this->inflection_point.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->inflection_point.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_inflection_point << ",\n"; - ss << " \"is_estimated\":" << this->inflection_point.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->inflection_point.is_random_effect_m << "\n },\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->inflection_point << ",\n },\n"; ss << " \"parameter\": {\n"; ss << " \"name\": \"slope\",\n"; ss << " \"id\":" << this->slope.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->slope.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_slope << ",\n"; - ss << " \"is_estimated\":" << this->slope.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->slope.is_random_effect_m << "\n }\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->slope << ",\n}\n"; ss << "}"; @@ -174,30 +171,43 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { std::shared_ptr > selectivity = std::make_shared >(); - + std::stringstream ss; // set relative info selectivity->id = this->id; - selectivity->inflection_point.resize(1); - selectivity->inflection_point[0] = this->inflection_point.value_m; - if (this->inflection_point.estimated_m) { - info->RegisterParameterName("logistic selectivity inflection_point"); - if (this->inflection_point.is_random_effect_m) { - info->RegisterRandomEffect(selectivity->inflection_point[0]); - } else { - info->RegisterParameter(selectivity->inflection_point[0]); + selectivity->inflection_point.resize(this->inflection_point.size()); + for (size_t i = 0; i < this->inflection_point.size(); i++) { + selectivity->inflection_point[i] = this->inflection_point[i].value_m; + if (this->inflection_point[i].estimated_m) { + ss.str(""); + ss << "selectivity.inflection_point ." << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->inflection_point[i].is_random_effect_m) { + info->RegisterRandomEffect(selectivity->inflection_point[i]); + } else { + info->RegisterParameter(selectivity->inflection_point[i]); + } } } + info->variable_map[this->inflection_point.id_m] = &(selectivity)->inflection_point; - selectivity->slope.resize(1); - selectivity->slope[0] = this->slope.value_m; - if (this->slope.estimated_m) { - info->RegisterParameterName("logistic selectivity slope"); - if (this->slope.is_random_effect_m) { - info->RegisterRandomEffect(selectivity->slope[0]); - } else { - info->RegisterParameter(selectivity->slope[0]); + + + selectivity->slope.resize(this->slope.size()); + for (size_t i = 0; i < this->slope.size(); i++) { + selectivity->slope[i] = this->slope[i].value_m; + if (this->slope[i].estimated_m) { + ss.str(""); + ss << "selectivity.slope." << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->slope[i].is_random_effect_m) { + info->RegisterRandomEffect(selectivity->slope[i]); + } else { + info->RegisterParameter(selectivity->slope[i]); + } } } + + info->variable_map[this->slope.id_m] = &(selectivity)->slope; // add to Information @@ -209,6 +219,8 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { /** @brief this adds the parameter values and derivatives to the TMB model * object */ virtual bool add_to_fims_tmb() { + FIMS_INFO_LOG("adding Selectivity object to TMB"); + this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); this->add_to_fims_tmb_internal(); @@ -226,20 +238,14 @@ class LogisticSelectivityInterface : public SelectivityInterfaceBase { */ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { public: - Parameter inflection_point_asc; /**< the index value at which the response + ParameterVector inflection_point_asc; /**< the index value at which the response reaches .5 */ - Parameter slope_asc; /**< the width of the curve at the inflection_point */ - Parameter inflection_point_desc; /**< the index value at which the response + ParameterVector slope_asc; /**< the width of the curve at the inflection_point */ + ParameterVector inflection_point_desc; /**< the index value at which the response reaches .5 */ - Parameter slope_desc; /**< the width of the curve at the inflection_point */ + ParameterVector slope_desc; /**< the width of the curve at the inflection_point */ - - double estimated_inflection_point_asc; - double estimated_slope_asc; - double estimated_inflection_point_desc; - double estimated_slope_desc; - DoubleLogisticSelectivityInterface() : SelectivityInterfaceBase() { } @@ -258,18 +264,17 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { virtual double evaluate(double x) { fims_popdy::DoubleLogisticSelectivity DoubleLogisticSel; DoubleLogisticSel.inflection_point_asc.resize(1); - DoubleLogisticSel.inflection_point_asc[0] = this->inflection_point_asc.value_m; + DoubleLogisticSel.inflection_point_asc[0] = this->inflection_point_asc[0].value_m; DoubleLogisticSel.slope_asc.resize(1); - DoubleLogisticSel.slope_asc[0] = this->slope_asc.value_m; + DoubleLogisticSel.slope_asc[0] = this->slope_asc[0].value_m; DoubleLogisticSel.inflection_point_desc.resize(1); DoubleLogisticSel.inflection_point_desc[0] = - this->inflection_point_desc.value_m; + this->inflection_point_desc[0].value_m; DoubleLogisticSel.slope_desc.resize(1); - DoubleLogisticSel.slope_desc[0] = this->slope_desc.value_m; + DoubleLogisticSel.slope_desc[0] = this->slope_desc[0].value_m; return DoubleLogisticSel.evaluate(x); } - - + virtual void finalize() { if (this->finalized) { @@ -283,12 +288,7 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { fims_info::Information::GetInstance(); - //set default values as initial values - this->estimated_inflection_point_asc = this->inflection_point_asc.value_m; - this->estimated_slope_asc = this->slope_asc.value_m; - this->estimated_inflection_point_desc = this->inflection_point_desc.value_m; - this->estimated_slope_desc = this->slope_desc.value_m; - + fims_info::Information::selectivity_models_iterator it; //search for maturity in Information @@ -301,72 +301,82 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { std::shared_ptr > sel = std::dynamic_pointer_cast >(it->second); - // if the parameter was estimated, set set the estimated value. - if (this->inflection_point_asc.estimated_m) { - this->estimated_inflection_point_asc = sel->inflection_point_asc[0]; - } - if (this->slope_asc.estimated_m) { - this->estimated_slope_asc = sel->slope_asc[0]; - } + for (size_t i = 0; i < inflection_point_asc.size(); i++) { + if (this->inflection_point_asc[i].estimated_m) { + this->inflection_point_asc[i].estimated_value_m = sel->inflection_point_asc[i]; + } else { + this->inflection_point_asc[i].estimated_value_m = this->inflection_point_asc[i].value_m; + } - if (this->inflection_point_desc.estimated_m) { - this->estimated_inflection_point_desc = sel->inflection_point_desc[0]; } - if (this->slope_desc.estimated_m) { - this->estimated_slope_desc = sel->slope_desc[0]; - } + for (size_t i = 0; i < slope_asc.size(); i++) { + if (this->slope_asc[i].estimated_m) { + this->slope_asc[i].estimated_value_m = sel->slope_asc[i]; + } else { + this->slope_asc[i].estimated_value_m = this->slope_asc[i].value_m; + } - } - } + } - virtual std::string to_json() { - std::stringstream ss; + for (size_t i = 0; i < inflection_point_desc.size(); i++) { + if (this->inflection_point_desc[i].estimated_m) { + this->inflection_point_desc[i].estimated_value_m = sel->inflection_point_desc[i]; + } else { + this->inflection_point_desc[i].estimated_value_m = this->inflection_point_desc[i].value_m; + } - ss << "\"module\" : {\n"; - ss << " \"name\": \"selectivity\",\n"; - ss << " \"type\": \"DoubleLogistic\",\n"; - ss << " \"id\": " << this->id << ",\n"; + } - ss << " \"parameter\": {\n"; - ss << " \"name\": \"inflection_point_asc\",\n"; - ss << " \"id\":" << this->inflection_point_asc.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->inflection_point_asc.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_inflection_point_asc << ",\n"; - ss << " \"is_estimated\":" << this->inflection_point_asc.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->inflection_point_asc.is_random_effect_m << "\n },\n"; + for (size_t i = 0; i < slope_desc.size(); i++) { + if (this->slope_desc[i].estimated_m) { + this->slope_desc[i].estimated_value_m = sel->slope_desc[i]; + } else { + this->slope_desc[i].estimated_value_m = this->slope_desc[i].value_m; + } - ss << " \"parameter\": {\n"; - ss << " \"name\": \"slope_asc\",\n"; - ss << " \"id\":" << this->slope_asc.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->slope_asc.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_slope_asc << ",\n"; - ss << " \"is_estimated\":" << this->slope_asc.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->slope_asc.is_random_effect_m << "\n },\n"; + } - ss << " \"parameter\": {\n"; - ss << " \"name\": \"inflection_point_desc\",\n"; - ss << " \"id\":" << this->inflection_point_desc.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->inflection_point_desc.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_inflection_point_desc << ",\n"; - ss << " \"is_estimated\":" << this->inflection_point_desc.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->inflection_point_desc.is_random_effect_m << "\n },\n"; - ss << " \"parameter\": {\n"; - ss << " \"name\": \"slope_desc\",\n"; - ss << " \"id\":" << this->slope_desc.id_m << ",\n"; - ss << " \"type\": \"scalar\",\n"; - ss << " \"value\":" << this->slope_desc.value_m << ",\n"; - ss << " \"estimated_value\":" << this->estimated_slope_desc << ",\n"; - ss << " \"is_estimated\":" << this->slope_desc.estimated_m << ",\n"; - ss << " \"is_random_effect\":" << this->slope_desc.is_random_effect_m << "\n }\n"; + } + } - ss << "}"; + virtual std::string to_json() { + std::stringstream ss; + + ss << "\"module\" : {\n"; + ss << " \"name\": \"selectivity\",\n"; + ss << " \"type\": \"DoubleLogistic\",\n"; + ss << " \"id\": " << this->id << ",\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"inflection_point_asc\",\n"; + ss << " \"id\":" << this->inflection_point_asc.id_m << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->inflection_point_asc << ",\n},\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"slope_asc\",\n"; + ss << " \"id\":" << this->slope_asc.id_m << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->slope_asc << ",\n},\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"inflection_point_desc\",\n"; + ss << " \"id\":" << this->inflection_point_desc.id_m << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->inflection_point_desc << ",\n},\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"slope_desc\",\n"; + ss << " \"id\":" << this->slope_desc.id_m << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->slope_desc << ",\n}\n"; + + + ss << "}"; return ss.str(); } @@ -382,51 +392,79 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { std::shared_ptr > selectivity = std::make_shared >(); + std::stringstream ss; // set relative info selectivity->id = this->id; - selectivity->inflection_point_asc.resize(1); - selectivity->inflection_point_asc[0] = this->inflection_point_asc.value_m; - if (this->inflection_point_asc.estimated_m) { - info->RegisterParameterName("double logistic selectivity inflection_point_asc"); - if (this->inflection_point_asc.is_random_effect_m) { - info->RegisterRandomEffect(selectivity->inflection_point_asc[0]); - } else { - info->RegisterParameter(selectivity->inflection_point_asc[0]); + selectivity->inflection_point_asc.resize(this->inflection_point_asc.size()); + for (size_t i = 0; i < this->inflection_point_asc.size(); i++) { + selectivity->inflection_point_asc[i] = this->inflection_point_asc[i].value_m; + if (this->inflection_point_asc[i].estimated_m) { + ss.str(""); + ss << "selectivity.inflection_point_asc." << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->inflection_point_asc[i].is_random_effect_m) { + info->RegisterRandomEffect(selectivity->inflection_point_asc[i]); + } else { + info->RegisterParameter(selectivity->inflection_point_asc[i]); + } } } + info->variable_map[this->inflection_point_asc.id_m] = &(selectivity)->inflection_point_asc; - selectivity->slope_asc.resize(1); - selectivity->slope_asc[0] = this->slope_asc.value_m; - if (this->slope_asc.estimated_m) { - info->RegisterParameterName("double logistic selectivity slope_asc"); - if (this->slope_asc.is_random_effect_m) { - info->RegisterRandomEffect(selectivity->slope_asc[0]); - } else { - info->RegisterParameter(selectivity->slope_asc[0]); + + + selectivity->slope_asc.resize(this->slope_asc.size()); + for (size_t i = 0; i < this->slope_asc.size(); i++) { + selectivity->slope_asc[i] = this->slope_asc[i].value_m; + if (this->slope_asc[i].estimated_m) { + ss.str(""); + ss << "selectivity.slope_asc." << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->slope_asc[i].is_random_effect_m) { + info->RegisterRandomEffect(selectivity->slope_asc[i]); + } else { + info->RegisterParameter(selectivity->slope_asc[i]); + } } } + + info->variable_map[this->slope_asc.id_m] = &(selectivity)->slope_asc; - selectivity->inflection_point_desc.resize(1); - selectivity->inflection_point_desc[0] = this->inflection_point_desc.value_m; - if (this->inflection_point_desc.estimated_m) { - info->RegisterParameterName("double logistic selectivity inflection_point_desc"); - if (this->inflection_point_desc.is_random_effect_m) { - info->RegisterRandomEffect(selectivity->inflection_point_desc[0]); - } else { - info->RegisterParameter(selectivity->inflection_point_desc[0]); + + selectivity->inflection_point_desc.resize(this->inflection_point_desc.size()); + for (size_t i = 0; i < this->inflection_point_desc.size(); i++) { + selectivity->inflection_point_desc[i] = this->inflection_point_desc[i].value_m; + if (this->inflection_point_desc[i].estimated_m) { + ss.str(""); + ss << "selectivity.inflection_point_desc." << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->inflection_point_desc[i].is_random_effect_m) { + info->RegisterRandomEffect(selectivity->inflection_point_desc[i]); + } else { + info->RegisterParameter(selectivity->inflection_point_desc[i]); + } } } + info->variable_map[this->inflection_point_desc.id_m] = &(selectivity)->inflection_point_desc; - selectivity->slope_desc.resize(1); - selectivity->slope_desc[0] = this->slope_desc.value_m; - if (this->slope_desc.estimated_m) { - info->RegisterParameterName("double logistic selectivity slope_desc"); - if (this->slope_desc.is_random_effect_m) { - info->RegisterRandomEffect(selectivity->slope_desc[0]); - } else { - info->RegisterParameter(selectivity->slope_desc[0]); + + + selectivity->slope_desc.resize(this->slope_desc.size()); + for (size_t i = 0; i < this->slope_desc.size(); i++) { + selectivity->slope_desc[i] = this->slope_desc[i].value_m; + if (this->slope_asc[i].estimated_m) { + ss.str(""); + ss << "selectivity.slope_desc." << this->id << "." << i; + info->RegisterParameterName(ss.str()); + if (this->slope_desc[i].is_random_effect_m) { + info->RegisterRandomEffect(selectivity->slope_desc[i]); + } else { + info->RegisterParameter(selectivity->slope_desc[i]); + } } } + + info->variable_map[this->slope_desc.id_m] = &(selectivity)->slope_desc; // add to Information diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp index d09fdd73..a9815794 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp @@ -18,65 +18,65 @@ * */ class DistributionsInterfaceBase : public FIMSRcppInterfaceBase { - public: - static uint32_t - id_g; /**< static id of the DistributionsInterfaceBase object */ - uint32_t id_m; /**< local id of the DistributionsInterfaceBase object */ - std::vector key_m; /**< unique id for variable map that points to a fims::Vector */ - std::string input_type_m; /**< type of density input, options are: prior, re, or data */ - // live objects in C++ are objects that have been created and live in memory - static std::map live_objects; /**< +public: + static uint32_t + id_g; /**< static id of the DistributionsInterfaceBase object */ + uint32_t id_m; /**< local id of the DistributionsInterfaceBase object */ + std::vector key_m; /**< unique id for variable map that points to a fims::Vector */ + std::string input_type_m; /**< type of density input, options are: prior, re, or data */ + // live objects in C++ are objects that have been created and live in memory + static std::map live_objects; /**< map relating the ID of the DistributionsInterfaceBase to the DistributionsInterfaceBase objects */ -uint32_t interface_observed_data_id_m = - -999; /**< id of observed data object*/ + uint32_t interface_observed_data_id_m = + -999; /**< id of observed data object*/ + + DistributionsInterfaceBase() { + this->id_m = DistributionsInterfaceBase::id_g++; + /* Create instance of map: key is id and value is pointer to + DistributionsInterfaceBase */ + DistributionsInterfaceBase::live_objects[this->id_m] = this; + FIMSRcppInterfaceBase::fims_interface_objects.push_back(this); + } - DistributionsInterfaceBase() { - this->id_m = DistributionsInterfaceBase::id_g++; - /* Create instance of map: key is id and value is pointer to - DistributionsInterfaceBase */ - DistributionsInterfaceBase::live_objects[this->id_m] = this; - FIMSRcppInterfaceBase::fims_interface_objects.push_back(this); - } + virtual ~DistributionsInterfaceBase() { + } - virtual ~DistributionsInterfaceBase() {} + /** @brief get_id method for child distribution interface objects to inherit + */ + virtual uint32_t get_id() = 0; + + /** + * @brief set_distribution_links sets pointers for data observations, random effects, or priors + * + * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. + * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector + */ + virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids) { + return false; + } - /** @brief get_id method for child distribution interface objects to inherit - */ - virtual uint32_t get_id() = 0; + /** + * @brief Set the unique id for the Observed Data object + * + * @param observed_data_id Unique id for the Observed Age Comp Data + * object + */ + virtual bool set_observed_data(int observed_data_id) { + return false; + } -/** - * @brief set_distribution_links sets pointers for data observations, random effects, or priors - * - * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. - * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector - */ - virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids){ - return false; - } - - - /** - * @brief Set the unique id for the Observed Data object - * - * @param observed_data_id Unique id for the Observed Age Comp Data - * object - */ - virtual bool set_observed_data(int observed_data_id){ - return false; - } - - /** @brief evaluate method for child distribution interface objects to inherit - */ - virtual double evaluate() = 0; + /** @brief evaluate method for child distribution interface objects to inherit + */ + virtual double evaluate() = 0; }; uint32_t DistributionsInterfaceBase::id_g = - 1; /**< static id of the DistributionsInterfaceBase object */ + 1; /**< static id of the DistributionsInterfaceBase object */ std::map /**< local id of the DistributionsInterfaceBase object */ - DistributionsInterfaceBase::live_objects; /**< +DistributionsInterfaceBase +*> /**< local id of the DistributionsInterfaceBase object */ +DistributionsInterfaceBase::live_objects; /**< map relating the ID of the DistributionsInterfaceBase to the DistributionsInterfaceBase objects */ @@ -87,128 +87,191 @@ std::mapid_m; } - - virtual ~DnormDistributionsInterface() {} - - /** - * @brief Set the unique id for the Observed Data object - * - * @param observed_data_id Unique id for the Observed Age Comp Data - * object - */ - virtual bool set_observed_data(int observed_data_id) { - this->interface_observed_data_id_m = observed_data_id; - return true; - } - - /** - * @brief set_distribution_links sets pointers for data observations, random effects, or priors - * - * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. - * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector - */ - virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids){ - this->input_type_m = input_type; - this->key_m.resize(ids.size()); - for(int i; ikey_m[i] = ids[i]; - } - - return true; - } - - /** - * @brief Evaluate normal probability density function, default returns the - * log of the pdf - * - * @tparam T - * @return log pdf - */ - virtual double evaluate() { - fims_distributions::NormalLPDF dnorm; - dnorm.x.resize(this->x.size()); - dnorm.expected_values.resize(this->expected_values.size()); - dnorm.log_sd.resize(this->log_sd.size()); - for(size_t i=0; ix[i].value_m; - } - for(size_t i=0; iexpected_values[i].value_m; - } - for(size_t i=0; ilog_sd[i].value_m; - } - return dnorm.evaluate(); - } - - +public: + ParameterVector x; /**< observed data */ + ParameterVector expected_values; /**< mean of x for the normal distribution */ + ParameterVector log_sd; /**< sd of x for the normal distribution */ + Rcpp::NumericVector lpdf_vec; + + DnormDistributionsInterface() : DistributionsInterfaceBase() { + } + + virtual uint32_t get_id() { + return this->id_m; + } + + virtual ~DnormDistributionsInterface() { + } + + /** + * @brief Set the unique id for the Observed Data object + * + * @param observed_data_id Unique id for the Observed Age Comp Data + * object + */ + virtual bool set_observed_data(int observed_data_id) { + this->interface_observed_data_id_m = observed_data_id; + return true; + } + + /** + * @brief set_distribution_links sets pointers for data observations, random effects, or priors + * + * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. + * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector + */ + virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids) { + this->input_type_m = input_type; + this->key_m.resize(ids.size()); + for (int i; i < ids.size(); i++) { + this->key_m[i] = ids[i]; + } + + return true; + } + + /** + * @brief Evaluate normal probability density function, default returns the + * log of the pdf + * + * @tparam T + * @return log pdf + */ + virtual double evaluate() { + fims_distributions::NormalLPDF dnorm; + dnorm.x.resize(this->x.size()); + dnorm.expected_values.resize(this->expected_values.size()); + dnorm.log_sd.resize(this->log_sd.size()); + for (size_t i = 0; i < x.size(); i++) { + dnorm.x[i] = this->x[i].value_m; + } + for (size_t i = 0; i < expected_values.size(); i++) { + dnorm.expected_values[i] = this->expected_values[i].value_m; + } + for (size_t i = 0; i < log_sd.size(); i++) { + dnorm.log_sd[i] = this->log_sd[i].value_m; + } + return dnorm.evaluate(); + } + + virtual void finalize() { + if (this->finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("DnormDistribution " + fims::to_string(this->id_m) + " has been finalized already."); + } + + this->finalized = true; //indicate this has been called already + + + std::shared_ptr > info = + fims_info::Information::GetInstance(); + + fims_info::Information::density_components_iterator it; + + //search for density component in Information + it = info->density_components.find(this->id_m); + //if not found, just return + if (it == info->density_components.end()) { + FIMS_WARNING_LOG("DnormDistribution " + fims::to_string(this->id_m) + " not found in Information."); + return; + } else { + std::shared_ptr > dnorm = + std::dynamic_pointer_cast >(it->second); + + this->lpdf_vec = Rcpp::NumericVector(dnorm->lpdf_vec.size()); + for (size_t i = 0; i < this->lpdf_vec.size(); i++) { + this->lpdf_vec[i] = dnorm->lpdf_vec[i]; + } + + } + + + } + + virtual std::string to_json() { + std::stringstream ss; + + ss << "\"module\" : {\n"; + ss << " \"name\": \"DnormDistribution\",\n"; + ss << " \"type\": \"normal\",\n"; + ss << " \"id\": " << this->id_m << ",\n"; + + ss << " \"density_component\": {\n"; + ss << " \"name\": \"lpdf_vec\",\n"; + ss << " \"values\":["; + if (this->lpdf_vec.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->lpdf_vec.size() - 1; i++) { + ss << this->lpdf_vec[i] << ", "; + } + ss << this->lpdf_vec[this->lpdf_vec.size() - 1] << "]\n"; + } + ss << " }\n]"; + + return ss.str(); + } + + #ifdef TMB_MODEL - template - bool add_to_fims_tmb_internal() { - std::shared_ptr> info = - fims_info::Information::GetInstance(); - - std::shared_ptr> distribution = - std::make_shared>(); - - // interface to data/parameter value - - distribution->observed_data_id_m = - interface_observed_data_id_m; - distribution->input_type = this->input_type_m; - distribution->key.resize(this->key_m.size()); - for(size_t i=0; ikey_m.size(); i++){ - distribution->key[i] = this->key_m[i]; - } - distribution->id = this->id_m; - distribution->x.resize(this->x.size()); - for(size_t i=0; ix.size(); i++){ - distribution->x[i] = this->x[i].value_m; - } - // set relative info - distribution->expected_values.resize(this->expected_values.size()); - for(size_t i=0; iexpected_values.size(); i++){ - distribution->expected_values[i] = this->expected_values[i].value_m; - } - distribution->log_sd.resize(this->log_sd.size()); - for(size_t i=0; ilog_sd.size(); i++){ - distribution->log_sd[i] = this->log_sd[i].value_m; - if(this->log_sd[i].estimated_m){ - info->RegisterParameterName("normal log_sd"); - info->RegisterParameter(distribution->log_sd[i]); - } - if (this->log_sd[i].is_random_effect_m) { - error("standard deviations cannot be set to random effects"); - } - } - info->variable_map[this->log_sd.id_m] = &(distribution)->log_sd; - - info->density_components[distribution->id] = distribution; - - return true; - } - - /** - * @brief adds the dnorm distribution and its parameters to the TMB model - */ - virtual bool add_to_fims_tmb() { - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - - return true; - } + template + bool add_to_fims_tmb_internal() { + std::shared_ptr> info = + fims_info::Information::GetInstance(); + + std::shared_ptr> distribution = + std::make_shared> (); + + // interface to data/parameter value + + distribution->observed_data_id_m = + interface_observed_data_id_m; + distribution->input_type = this->input_type_m; + distribution->key.resize(this->key_m.size()); + for (size_t i = 0; ikey_m.size(); i++) { + distribution->key[i] = this->key_m[i]; + } + distribution->id = this->id_m; + distribution->x.resize(this->x.size()); + for (size_t i = 0; ix.size(); i++) { + distribution->x[i] = this->x[i].value_m; + } + // set relative info + distribution->expected_values.resize(this->expected_values.size()); + for (size_t i = 0; iexpected_values.size(); i++) { + distribution->expected_values[i] = this->expected_values[i].value_m; + } + distribution->log_sd.resize(this->log_sd.size()); + for (size_t i = 0; ilog_sd.size(); i++) { + distribution->log_sd[i] = this->log_sd[i].value_m; + if (this->log_sd[i].estimated_m) { + info->RegisterParameterName("normal log_sd"); + info->RegisterParameter(distribution->log_sd[i]); + } + if (this->log_sd[i].is_random_effect_m) { + error("standard deviations cannot be set to random effects"); + } + } + info->variable_map[this->log_sd.id_m] = &(distribution)->log_sd; + + info->density_components[distribution->id] = distribution; + + return true; + } + + /** + * @brief adds the dnorm distribution and its parameters to the TMB model + */ + virtual bool add_to_fims_tmb() { + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + + return true; + } #endif }; @@ -220,131 +283,196 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { * */ class DlnormDistributionsInterface : public DistributionsInterfaceBase { - public: - ParameterVector x; /**< observation */ - ParameterVector expected_values; /**< mean of the distribution of log(x) */ - ParameterVector log_logsd; /**< log standard deviation of the distribution of log(x) */ - Rcpp::String input_type; /**< character string indicating type of input: data, re, prior */ - - DlnormDistributionsInterface() : DistributionsInterfaceBase() {} - - virtual ~DlnormDistributionsInterface() {} - - /** - * @brief get the id of the Dlnorm distributions interface class object - */ - virtual uint32_t get_id() { return this->id_m; } - - /** - * @brief Set the unique id for the Observed Data object - * - * @param observed_data_id Unique id for the Observed Age Comp Data - * object - */ - virtual bool set_observed_data(int observed_data_id) { - this->interface_observed_data_id_m = observed_data_id; - - return true; - } - - /** - * @brief set_distribution_links sets pointers for data observations, random effects, or priors - * - * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. - * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector - */ - virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids){ - this->input_type_m = input_type; - this->key_m.resize(ids.size()); - for(int i; ikey_m[i] = ids[i]; - } - - return true; - } - - /** - * @brief Evaluate lognormal probability density function, default returns the - * log of the pdf - * - * @tparam T - * @return log pdf - */ - virtual double evaluate() { - fims_distributions::LogNormalLPDF dlnorm; - dlnorm.input_type = this->input_type; - dlnorm.x.resize(this->x.size()); - dlnorm.expected_values.resize(this->expected_values.size()); - dlnorm.log_logsd.resize(this->log_logsd.size()); - for(size_t i=0; ix[i].value_m; - } - for(size_t i=0; iexpected_values[i].value_m; - } - for(size_t i=0; ilog_logsd[i].value_m; - } - return dlnorm.evaluate(); - } +public: + ParameterVector x; /**< observation */ + ParameterVector expected_values; /**< mean of the distribution of log(x) */ + ParameterVector log_logsd; /**< log standard deviation of the distribution of log(x) */ + Rcpp::String input_type; /**< character string indicating type of input: data, re, prior */ + Rcpp::NumericVector lpdf_vec; + + DlnormDistributionsInterface() : DistributionsInterfaceBase() { + } + + virtual ~DlnormDistributionsInterface() { + } + + /** + * @brief get the id of the Dlnorm distributions interface class object + */ + virtual uint32_t get_id() { + return this->id_m; + } + + /** + * @brief Set the unique id for the Observed Data object + * + * @param observed_data_id Unique id for the Observed Age Comp Data + * object + */ + virtual bool set_observed_data(int observed_data_id) { + this->interface_observed_data_id_m = observed_data_id; + + return true; + } + + /** + * @brief set_distribution_links sets pointers for data observations, random effects, or priors + * + * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. + * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector + */ + virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids) { + this->input_type_m = input_type; + this->key_m.resize(ids.size()); + for (int i; i < ids.size(); i++) { + this->key_m[i] = ids[i]; + } + + return true; + } + + /** + * @brief Evaluate lognormal probability density function, default returns the + * log of the pdf + * + * @tparam T + * @return log pdf + */ + virtual double evaluate() { + fims_distributions::LogNormalLPDF dlnorm; + dlnorm.input_type = this->input_type; + dlnorm.x.resize(this->x.size()); + dlnorm.expected_values.resize(this->expected_values.size()); + dlnorm.log_logsd.resize(this->log_logsd.size()); + for (size_t i = 0; i < x.size(); i++) { + dlnorm.x[i] = this->x[i].value_m; + } + for (size_t i = 0; i < expected_values.size(); i++) { + dlnorm.expected_values[i] = this->expected_values[i].value_m; + } + for (size_t i = 0; i < log_logsd.size(); i++) { + dlnorm.log_logsd[i] = this->log_logsd[i].value_m; + } + return dlnorm.evaluate(); + } + + virtual void finalize() { + if (this->finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("LogNormalLPDF " + fims::to_string(this->id_m) + " has been finalized already."); + } + + this->finalized = true; //indicate this has been called already + + + std::shared_ptr > info = + fims_info::Information::GetInstance(); + + fims_info::Information::density_components_iterator it; + + //search for density component in Information + it = info->density_components.find(this->id_m); + //if not found, just return + if (it == info->density_components.end()) { + FIMS_WARNING_LOG("LogNormalLPDF " + fims::to_string(this->id_m) + " not found in Information."); + return; + } else { + std::shared_ptr > dnorm = + std::dynamic_pointer_cast >(it->second); + + this->lpdf_vec = Rcpp::NumericVector(dnorm->lpdf_vec.size()); + for (size_t i = 0; i < this->lpdf_vec.size(); i++) { + this->lpdf_vec[i] = dnorm->lpdf_vec[i]; + } + + } + + + } + + virtual std::string to_json() { + + std::stringstream ss; + ss << "\"module\" : {\n"; + ss << " \"name\": \"LogNormalLPDF\",\n"; + ss << " \"type\": \"log_normal\",\n"; + ss << " \"id\": " << this->id_m << ",\n"; + + ss << " \"density_component\": {\n"; + ss << " \"name\": \"lpdf_vec\",\n"; + ss << " \"values\":["; + if (this->lpdf_vec.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->lpdf_vec.size() - 1; i++) { + ss << this->lpdf_vec[i] << ", "; + } + ss << this->lpdf_vec[this->lpdf_vec.size() - 1] << "]\n"; + } + ss << " }\n]"; + + return ss.str(); + } + + #ifdef TMB_MODEL - template - bool add_to_fims_tmb_internal() { - std::shared_ptr> info = - fims_info::Information::GetInstance(); - - std::shared_ptr> distribution = - std::make_shared>(); - - // set relative info - distribution->id = this->id_m; - distribution->observed_data_id_m = - interface_observed_data_id_m; - distribution->input_type = this->input_type_m; - distribution->key.resize(this->key_m.size()); - for(size_t i=0; ikey_m.size(); i++){ - distribution->key[i] = this->key_m[i]; - } - distribution->x.resize(this->x.size()); - for(size_t i=0; ix.size(); i++){ - distribution->x[i] = this->x[i].value_m; - } - // set relative info - distribution->expected_values.resize(this->expected_values.size()); - for(size_t i=0; iexpected_values.size(); i++){ - distribution->expected_values[i] = this->expected_values[i].value_m; - } - distribution->log_logsd.resize(this->log_logsd.size()); - for(size_t i=0; ilog_logsd.size(); i++){ - distribution->log_logsd[i] = this->log_logsd[i].value_m; - if(this->log_logsd[i].estimated_m){ - info->RegisterParameterName("lognormal log_logsd"); - info->RegisterParameter(distribution->log_logsd[i]); - } - if (this->log_logsd[i].is_random_effect_m) { - error("standard deviations cannot be set to random effects"); - } - } - info->variable_map[this->log_logsd.id_m] = &(distribution)->log_logsd; - - info->density_components[distribution->id] = distribution; - - return true; - } - - /** - * @brief adds the dlnorm distribution and its parameters to the TMB model - */ - virtual bool add_to_fims_tmb() { - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - - return true; - } + template + bool add_to_fims_tmb_internal() { + std::shared_ptr> info = + fims_info::Information::GetInstance(); + + std::shared_ptr> distribution = + std::make_shared> (); + + // set relative info + distribution->id = this->id_m; + distribution->observed_data_id_m = + interface_observed_data_id_m; + distribution->input_type = this->input_type_m; + distribution->key.resize(this->key_m.size()); + for (size_t i = 0; ikey_m.size(); i++) { + distribution->key[i] = this->key_m[i]; + } + distribution->x.resize(this->x.size()); + for (size_t i = 0; ix.size(); i++) { + distribution->x[i] = this->x[i].value_m; + } + // set relative info + distribution->expected_values.resize(this->expected_values.size()); + for (size_t i = 0; iexpected_values.size(); i++) { + distribution->expected_values[i] = this->expected_values[i].value_m; + } + distribution->log_logsd.resize(this->log_logsd.size()); + for (size_t i = 0; ilog_logsd.size(); i++) { + distribution->log_logsd[i] = this->log_logsd[i].value_m; + if (this->log_logsd[i].estimated_m) { + info->RegisterParameterName("lognormal log_logsd"); + info->RegisterParameter(distribution->log_logsd[i]); + } + if (this->log_logsd[i].is_random_effect_m) { + error("standard deviations cannot be set to random effects"); + } + } + info->variable_map[this->log_logsd.id_m] = &(distribution)->log_logsd; + + info->density_components[distribution->id] = distribution; + + return true; + } + + /** + * @brief adds the dlnorm distribution and its parameters to the TMB model + */ + virtual bool add_to_fims_tmb() { + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + + return true; + } #endif }; @@ -358,116 +486,180 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { // template class DmultinomDistributionsInterface : public DistributionsInterfaceBase { - public: - ParameterVector x; /**< Vector of length K of integers */ - ParameterVector expected_values; /**< Vector of length K, specifying the probability +public: + ParameterVector x; /**< Vector of length K of integers */ + ParameterVector expected_values; /**< Vector of length K, specifying the probability for the K classes (note, unlike in R these must sum to 1). */ - Rcpp::NumericVector dims; /**< Dimensions of the number of rows and columns of the multivariate dataset */ - - DmultinomDistributionsInterface() : DistributionsInterfaceBase() {} - - virtual ~DmultinomDistributionsInterface() {} - - virtual uint32_t get_id() { return this->id_m; } - - /** - * @brief Set the unique id for the Observed Data object - * - * @param observed_data_id Unique id for the Observed Age Comp Data - * object - */ - virtual bool set_observed_data(int observed_data_id) { - this->interface_observed_data_id_m = observed_data_id; - - return true; - } - - /** - * @brief set_distribution_links sets pointers for data observations, random effects, or priors - * - * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. - * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector - */ - virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids){ - this->input_type_m = input_type; - this->key_m.resize(ids.size()); - for(int i; ikey_m[i] = ids[i]; - } - - return true; - } - - /** - * @brief Evaluate multinom probability density function, default returns the - * log of the pdf - * - * @tparam T - * @return log pdf - */ - virtual double evaluate() { - fims_distributions::MultinomialLPMF dmultinom; - // Declare TMBVector in this scope - dmultinom.x.resize(this->x.size()); - dmultinom.expected_values.resize(this->expected_values.size()); - for(size_t i=0; ix[i].value_m; - } - for(size_t i=0; iexpected_values[i].value_m; - } - dmultinom.dims.resize(2); - dmultinom.dims[0] = this->dims[0]; - dmultinom.dims[1] = this->dims[1]; - return dmultinom.evaluate(); - } + Rcpp::NumericVector dims; /**< Dimensions of the number of rows and columns of the multivariate dataset */ + Rcpp::NumericVector lpdf_vec; -#ifdef TMB_MODEL + DmultinomDistributionsInterface() : DistributionsInterfaceBase() { + } + + virtual ~DmultinomDistributionsInterface() { + } - template - bool add_to_fims_tmb_internal() { - std::shared_ptr> info = - fims_info::Information::GetInstance(); + virtual uint32_t get_id() { + return this->id_m; + } - std::shared_ptr> distribution = - std::make_shared>(); + /** + * @brief Set the unique id for the Observed Data object + * + * @param observed_data_id Unique id for the Observed Age Comp Data + * object + */ + virtual bool set_observed_data(int observed_data_id) { + this->interface_observed_data_id_m = observed_data_id; - distribution->id = this->id_m; - distribution->observed_data_id_m = - interface_observed_data_id_m; - distribution->input_type = this->input_type_m; - distribution->key.resize(this->key_m.size()); - for(size_t i=0; ikey_m.size(); i++){ - distribution->key[i] = this->key_m[i]; + return true; } - distribution->x.resize(this->x.size()); - for(size_t i=0; ix.size(); i++){ - distribution->x[i] = this->x[i].value_m; + + /** + * @brief set_distribution_links sets pointers for data observations, random effects, or priors + * + * @param input_type String that sets whether the distribution type is: priors, random_effects, or data. + * @param ids Vector of unique ids for each linked parameter/s, derived value/s, or observed data vector + */ + virtual bool set_distribution_links(std::string input_type, Rcpp::IntegerVector ids) { + this->input_type_m = input_type; + this->key_m.resize(ids.size()); + for (int i; i < ids.size(); i++) { + this->key_m[i] = ids[i]; + } + + return true; } - // set relative info - distribution->expected_values.resize(this->expected_values.size()); - for(size_t i=0; iexpected_values.size(); i++){ - distribution->expected_values[i] = this->expected_values[i].value_m; + + /** + * @brief Evaluate multinom probability density function, default returns the + * log of the pdf + * + * @tparam T + * @return log pdf + */ + virtual double evaluate() { + fims_distributions::MultinomialLPMF dmultinom; + // Declare TMBVector in this scope + dmultinom.x.resize(this->x.size()); + dmultinom.expected_values.resize(this->expected_values.size()); + for (size_t i = 0; i < x.size(); i++) { + dmultinom.x[i] = this->x[i].value_m; + } + for (size_t i = 0; i < expected_values.size(); i++) { + dmultinom.expected_values[i] = this->expected_values[i].value_m; + } + dmultinom.dims.resize(2); + dmultinom.dims[0] = this->dims[0]; + dmultinom.dims[1] = this->dims[1]; + return dmultinom.evaluate(); + } + + virtual void finalize() { + if (this->finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("MultinomialLPMF " + fims::to_string(this->id_m) + " has been finalized already."); + } + + this->finalized = true; //indicate this has been called already + + + std::shared_ptr > info = + fims_info::Information::GetInstance(); + + fims_info::Information::density_components_iterator it; + + //search for density component in Information + it = info->density_components.find(this->id_m); + //if not found, just return + if (it == info->density_components.end()) { + FIMS_WARNING_LOG("MultinomialLPMF " + fims::to_string(this->id_m) + " not found in Information."); + return; + } else { + std::shared_ptr > dnorm = + std::dynamic_pointer_cast >(it->second); + + this->lpdf_vec = Rcpp::NumericVector(dnorm->lpdf_vec.size()); + for (size_t i = 0; i < this->lpdf_vec.size(); i++) { + this->lpdf_vec[i] = dnorm->lpdf_vec[i]; + } + + } + + } - if(this->dims.size()>0){ - distribution->dims.resize(2); - distribution->dims[0] = this->dims[0]; - distribution->dims[1] = this->dims[1]; + + virtual std::string to_json() { + + std::stringstream ss; + ss << "\"module\" : {\n"; + ss << " \"name\": \"MultinomialLPMF\",\n"; + ss << " \"type\": \"multinomial\",\n"; + ss << " \"id\": " << this->id_m << ",\n"; + + ss << " \"density_component\": {\n"; + ss << " \"name\": \"lpdf_vec\",\n"; + ss << " \"values\":["; + if (this->lpdf_vec.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->lpdf_vec.size() - 1; i++) { + ss << this->lpdf_vec[i] << ", "; + } + ss << this->lpdf_vec[this->lpdf_vec.size() - 1] << "]\n"; + } + ss << " }\n]"; + + return ss.str(); } - info->density_components[distribution->id] = distribution; - return true; - } +#ifdef TMB_MODEL - virtual bool add_to_fims_tmb() { - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); - this->add_to_fims_tmb_internal(); + template + bool add_to_fims_tmb_internal() { + std::shared_ptr> info = + fims_info::Information::GetInstance(); + + std::shared_ptr> distribution = + std::make_shared> (); + + distribution->id = this->id_m; + distribution->observed_data_id_m = + interface_observed_data_id_m; + distribution->input_type = this->input_type_m; + distribution->key.resize(this->key_m.size()); + for (size_t i = 0; ikey_m.size(); i++) { + distribution->key[i] = this->key_m[i]; + } + distribution->x.resize(this->x.size()); + for (size_t i = 0; ix.size(); i++) { + distribution->x[i] = this->x[i].value_m; + } + // set relative info + distribution->expected_values.resize(this->expected_values.size()); + for (size_t i = 0; iexpected_values.size(); i++) { + distribution->expected_values[i] = this->expected_values[i].value_m; + } + if (this->dims.size() > 0) { + distribution->dims.resize(2); + distribution->dims[0] = this->dims[0]; + distribution->dims[1] = this->dims[1]; + } + + info->density_components[distribution->id] = distribution; + + return true; + } - return true; - } + virtual bool add_to_fims_tmb() { + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + this->add_to_fims_tmb_internal(); + + return true; + } #endif }; diff --git a/inst/include/population_dynamics/fleet/fleet.hpp b/inst/include/population_dynamics/fleet/fleet.hpp index 9ca0db58..bd120fde 100644 --- a/inst/include/population_dynamics/fleet/fleet.hpp +++ b/inst/include/population_dynamics/fleet/fleet.hpp @@ -18,153 +18,165 @@ namespace fims_popdy { -/** @brief Base class for all fleets. - * - * @tparam Type The type of the fleet object. - */ -template -struct Fleet : public fims_model_object::FIMSObject { - static uint32_t id_g; /*!< reference id for fleet object*/ - size_t nyears; /*!< the number of years in the model*/ - size_t nages; /*!< the number of ages in the model*/ - - // selectivity - int fleet_selectivity_id_m = -999; /*!< id of selectivity component*/ - std::shared_ptr> - selectivity; /*!< selectivity component*/ - - // Mortality and catchability - fims::Vector - log_Fmort; /*!< estimated parameter: log Fishing mortality*/ - Type log_q; /*!< estimated parameter: catchability of the fleet */ - - fims::Vector Fmort; /*!< transformed parameter: Fishing mortality*/ - Type q; /*!< transofrmed parameter: the catchability of the fleet */ - - // derived quantities - fims::Vector catch_at_age; /*! catch_index; /*! age_composition; /*! observed_catch_lpdf; /*! + struct Fleet : public fims_model_object::FIMSObject { + static uint32_t id_g; /*!< reference id for fleet object*/ + size_t nyears; /*!< the number of years in the model*/ + size_t nages; /*!< the number of ages in the model*/ + + // selectivity + int fleet_selectivity_id_m = -999; /*!< id of selectivity component*/ + std::shared_ptr> + selectivity; /*!< selectivity component*/ + + // Mortality and catchability + fims::Vector + log_Fmort; /*!< estimated parameter: log Fishing mortality*/ + fims::Vector log_q; /*!< estimated parameter: catchability of the fleet */ + + fims::Vector Fmort; /*!< transformed parameter: Fishing mortality*/ + fims::Vector q; /*!< transofrmed parameter: the catchability of the fleet */ + + // derived quantities + fims::Vector catch_at_age; /*! catch_index; /*! age_composition; /*! observed_catch_lpdf; /*! observed_index_lpdf; /*! observed_index_lpdf; /*! expected_catch; /*! expected_index; /*! log_expected_index; /*! expected_catch_lpdf; /*! expected_catch; /*! expected_index; /*! log_expected_index; /*! expected_catch_lpdf; /*! expected_index_lpdf; /*! expected_index_lpdf; /*! catch_numbers_at_age; /*! proportion_catch_numbers_at_age; /*! catch_weight_at_age; /*! catch_numbers_at_age; /*! proportion_catch_numbers_at_age; /*! catch_weight_at_age; /*! *of; + ::objective_function *of; #endif - /** - * @brief Constructor. - */ - Fleet() { this->id = Fleet::id_g++; } - - /** - * @brief Destructor. - */ - virtual ~Fleet() {} - - /** - * @brief Intialize Fleet Class - * @param nyears The number of years in the model. - * @param nages The number of ages in the model. - */ - void Initialize(int nyears, int nages) { - this->nyears = nyears; - this->nages = nages; - - catch_at_age.resize(nyears * nages); - catch_numbers_at_age.resize(nyears * nages); - catch_weight_at_age.resize(nyears * nages); - catch_index.resize(nyears); // assume index is for all ages. - expected_catch.resize(nyears); - expected_index.resize(nyears); - age_composition.resize(nyears * nages); - - log_Fmort.resize(nyears); - Fmort.resize(nyears); - } - - /** - * @brief Prepare to run the fleet module. Called at each model itartion, and - * used to exponentiate the log q and Fmort parameters prior to evaluation. - * - */ - void Prepare() { - // for(size_t fleet_ = 0; fleet_ <= this->nfleets; fleet_++) { - // this -> Fmort[fleet_] = fims_math::exp(this -> log_Fmort[fleet_]); - - // derived quantities - std::fill(catch_at_age.begin(), catch_at_age.end(), - 0); /**q = fims_math::exp(this->log_q); - for (size_t year = 0; year < this->nyears; year++) { - FLEET_LOG << "input F mort " << this->log_Fmort[year] << std::endl; - FLEET_LOG << "input q " << this->log_q << std::endl; - this->Fmort[year] = fims_math::exp(this->log_Fmort[year]); - } - } - - /** - * Evaluate the proportion of catch numbers at age. - */ - void evaluate_age_comp() { - for (size_t y = 0; y < this->nyears; y++) { - Type sum = 0.0; - for (size_t a = 0; a < this->nages; a++) { - size_t i_age_year = y * this->nages + a; - sum += this->catch_numbers_at_age[i_age_year]; - } - for (size_t a = 0; a < this->nages; a++) { - size_t i_age_year = y * this->nages + a; - this->proportion_catch_numbers_at_age[i_age_year] = this->catch_numbers_at_age[i_age_year] / sum; - - } - } - } - - /** - * Evaluate the log of the expected index. - */ - void evaluate_index() { - for(size_t i=0; iexpected_index.size(); i++){ - log_expected_index[i] = log(this->expected_index[i]); - } - } -}; - -// default id of the singleton fleet class -template -uint32_t Fleet::id_g = 0; - -} // end namespace fims_popdy + /** + * @brief Constructor. + */ + Fleet() { + this->id = Fleet::id_g++; + } + + /** + * @brief Destructor. + */ + virtual ~Fleet() { + } + + /** + * @brief Intialize Fleet Class + * @param nyears The number of years in the model. + * @param nages The number of ages in the model. + */ + void Initialize(int nyears, int nages) { + if(this->log_q.size()==0){ + this->log_q.resize(1); + this->log_q[0] = 0.0; + } + this->nyears = nyears; + this->nages = nages; + + catch_at_age.resize(nyears * nages); + catch_numbers_at_age.resize(nyears * nages); + catch_weight_at_age.resize(nyears * nages); + catch_index.resize(nyears); // assume index is for all ages. + expected_catch.resize(nyears); + expected_index.resize(nyears); + log_expected_index.resize(nyears); + age_composition.resize(nyears * nages); + q.resize(this->log_q.size()); + log_Fmort.resize(nyears); + Fmort.resize(nyears); + } + + /** + * @brief Prepare to run the fleet module. Called at each model itartion, and + * used to exponentiate the log q and Fmort parameters prior to evaluation. + * + */ + void Prepare() { + // for(size_t fleet_ = 0; fleet_ <= this->nfleets; fleet_++) { + // this -> Fmort[fleet_] = fims_math::exp(this -> log_Fmort[fleet_]); + + // derived quantities + std::fill(catch_at_age.begin(), catch_at_age.end(), + 0); /**log_q.size(); i++) { + this->q[i] = fims_math::exp(this->log_q[i]); + } + + for (size_t year = 0; year < this->nyears; year++) { + FLEET_LOG << "input F mort " << this->log_Fmort[year] << std::endl; + FLEET_LOG << "input q " << this->log_q << std::endl; + this->Fmort[year] = fims_math::exp(this->log_Fmort[year]); + } + } + + /** + * Evaluate the proportion of catch numbers at age. + */ + void evaluate_age_comp() { + for (size_t y = 0; y < this->nyears; y++) { + Type sum = 0.0; + for (size_t a = 0; a < this->nages; a++) { + size_t i_age_year = y * this->nages + a; + sum += this->catch_numbers_at_age[i_age_year]; + } + for (size_t a = 0; a < this->nages; a++) { + size_t i_age_year = y * this->nages + a; + this->proportion_catch_numbers_at_age[i_age_year] = this->catch_numbers_at_age[i_age_year] / sum; + + } + } + } + + /** + * Evaluate the log of the expected index. + */ + void evaluate_index() { + for (size_t i = 0; iexpected_index.size(); i++) { + log_expected_index[i] = log(this->expected_index[i]); + } + } + }; + + // default id of the singleton fleet class + template + uint32_t Fleet::id_g = 0; + +} // end namespace fims_popdy #endif /* FIMS_POPULATION_DYNAMICS_FLEET_HPP */ diff --git a/inst/include/population_dynamics/growth/functors/ewaa.hpp b/inst/include/population_dynamics/growth/functors/ewaa.hpp index 62ff5c5d..3ee420af 100644 --- a/inst/include/population_dynamics/growth/functors/ewaa.hpp +++ b/inst/include/population_dynamics/growth/functors/ewaa.hpp @@ -16,32 +16,39 @@ namespace fims_popdy { -/** - * @brief EWAAgrowth class that returns the EWAA function value. - */ -template -struct EWAAgrowth : public GrowthBase { - // add submodule class members here - // these include parameters of the submodule - // a map looks up values based on a reference key - // in this case, our key is age (first double), and - // the value is the weight at that age (second double) - std::map ewaa; /** + struct EWAAgrowth : public GrowthBase { + // add submodule class members here + // these include parameters of the submodule + // a map looks up values based on a reference key + // in this case, our key is age (first double), and + // the value is the weight at that age (second double) + std::map ewaa; /** */ + typedef typename std::map::iterator weight_iterator; - EWAAgrowth() : GrowthBase() {} + EWAAgrowth() : GrowthBase() { + } - virtual ~EWAAgrowth() {} + virtual ~EWAAgrowth() { + } - /** - * @brief Returns the weight at age a (in kg) from the input vector. - * - * @param a age of the fish, the age vector must start at zero - */ - virtual const Type evaluate(const double& a) { - Type ret = ewaa[a]; - return ret; - } -}; -} // namespace fims_popdy + /** + * @brief Returns the weight at age a (in kg) from the input vector. + * + * @param a age of the fish, the age vector must start at zero + */ + virtual const Type evaluate(const double& a) { + weight_iterator it = this->ewaa.find(a); + if(it == this->ewaa.end() ){ + return 0.0; + } + Type ret = (*it).second;//itewaa[a]; + return ret; + } + }; +} // namespace fims_popdy #endif /* POPULATION_DYNAMICS_GROWTH_EWAA_HPP */ \ No newline at end of file diff --git a/inst/include/population_dynamics/population/population.hpp b/inst/include/population_dynamics/population/population.hpp index 9d969415..9388becd 100644 --- a/inst/include/population_dynamics/population/population.hpp +++ b/inst/include/population_dynamics/population/population.hpp @@ -156,7 +156,7 @@ struct Population : public fims_model_object::FIMSObject { POPULATION_LOG << "nseasons: " << this->nseasons << std::endl; POPULATION_LOG << "nyears: " << this->nyears << std::endl; - for (size_t fleet = 0; fleet < this->nfleets; fleet++) { + for (size_t fleet = 0; fleet < this->fleets.size(); fleet++) { this->fleets[fleet]->Prepare(); } @@ -457,7 +457,7 @@ struct Population : public fims_model_object::FIMSObject { this->weight_at_age[age]; } else { POPULATION_LOG << "fleet " << fleet_ << " is a survey" << std::endl; - index_ = this->fleets[fleet_]->q * + index_ = this->fleets[fleet_]->q[0] * this->fleets[fleet_]->selectivity->evaluate(ages[age]) * this->numbers_at_age[i_age_year] * this->weight_at_age[age]; // this->weight_at_age[age]; diff --git a/tests/gtest/integration_test_population.cpp b/tests/gtest/integration_test_population.cpp index 606f22eb..a4d02345 100644 --- a/tests/gtest/integration_test_population.cpp +++ b/tests/gtest/integration_test_population.cpp @@ -218,10 +218,10 @@ namespace typename JsonObject::iterator fleet2_index; fleet2_index = it->second.GetObject().find("survey1"); JsonArray &fleet_index = (*fleet2_index).second.GetArray(); - EXPECT_EQ(pop.fleets[0]->q, 1.0); + EXPECT_EQ(pop.fleets[0]->q[0], 1.0); // Do not use EXPECT_EQ to compare floats or doubles // Use EXPECT_NEAR here - EXPECT_NEAR(pop.fleets[1]->q, fleet_q[0].GetDouble(), 1.0e-07); + EXPECT_NEAR(pop.fleets[1]->q[0], fleet_q[0].GetDouble(), 1.0e-07); if(pop.fleets[1]->is_survey){ for (int year = 0; year < pop.nyears; year++) diff --git a/tests/gtest/test_population_Index.cpp b/tests/gtest/test_population_Index.cpp index f1f0efbd..d70ecbe8 100644 --- a/tests/gtest/test_population_Index.cpp +++ b/tests/gtest/test_population_Index.cpp @@ -26,7 +26,7 @@ namespace // If not possible to test entire vector, test middle or second to last // than earlier years (collapses to mean in early years) expected_index[index_yf] += population.numbers_at_age[i_age_year]* - population.fleets[fleet_]->q* + population.fleets[fleet_]->q[0]* population.fleets[fleet_]->selectivity->evaluate(population.ages[age])* population.growth->evaluate(population.ages[age]); diff --git a/tests/gtest/test_population_dynamics_fleet_initialize_prepare.cpp b/tests/gtest/test_population_dynamics_fleet_initialize_prepare.cpp index f5c4c969..9d4b68d2 100644 --- a/tests/gtest/test_population_dynamics_fleet_initialize_prepare.cpp +++ b/tests/gtest/test_population_dynamics_fleet_initialize_prepare.cpp @@ -31,6 +31,7 @@ namespace fleet.expected_catch.resize(nyears); fleet.expected_index.resize(nyears); fleet.catch_numbers_at_age.resize(nyears * nages); + fleet.log_q.resize(1);//needs to be initialized here, size used by q in Initialize fleet.Initialize(nyears, nages); int seed = 1234; @@ -44,17 +45,19 @@ namespace double log_q_min = fims_math::log(0.1); double log_q_max = fims_math::log(1); std::uniform_real_distribution log_q_distribution(log_q_min, log_q_max); - fleet.log_q = log_q_distribution(generator); + + fleet.log_q[0] = log_q_distribution(generator); for(int i = 0; i < nyears; i++) { fleet.log_Fmort[i] = log_Fmort_distribution(generator); } + fleet.Prepare(); // Test fleet.Fmort and fleet.q std::vector Fmort(nyears, 0); - double q = fims_math::exp(fleet.log_q); - EXPECT_EQ(fleet.q, q); + double q = fims_math::exp(fleet.log_q[0]); + EXPECT_EQ(fleet.q[0], q); for (int i = 0; i < nyears; i++) { Fmort[i] = fims_math::exp(fleet.log_Fmort[i]); diff --git a/tests/gtest/test_population_test_fixture.hpp b/tests/gtest/test_population_test_fixture.hpp index 89a34c30..3f6272eb 100644 --- a/tests/gtest/test_population_test_fixture.hpp +++ b/tests/gtest/test_population_test_fixture.hpp @@ -2,333 +2,351 @@ #include "population_dynamics/population/population.hpp" + + + namespace { -// Use test fixture to reuse the same configuration of objects for -// several different tests. To use a test fixture, derive a class -// from testing::Test. -class PopulationInitializeTestFixture : public testing::Test { - // Make members protected and they can be accessed from - // sub-classes. - protected: - // Use SetUp function to prepare the objects for each test. - // Use override in C++11 to make sure SetUp (e.g., not Setup with - // a lowercase u) is spelled - // correctly. - void SetUp() override { - population.id_g = id_g; - population.nyears = nyears; - population.nseasons = nseasons; - population.nages = nages; - for (int i = 0; i < nfleets; i++) { - auto fleet = std::make_shared>(); - population.fleets.push_back(fleet); - } - } - - // Virtual void TearDown() will be called after each test is - // run. It needs to be defined if there is clearup work to - // do. Otherwise, it does not need to be provided. - virtual void TearDown() {} - - fims_popdy::Population population; - - // Use default values from the Li et al., 2021 - // https://github.com/Bai-Li-NOAA/Age_Structured_Stock_Assessment_Model_Comparison/blob/master/R/save_initial_input.R - int id_g = 0; - int nyears = 30; - int nseasons = 1; - int nages = 12; - int nfleets = 2; -}; - -class PopulationEvaluateTestFixture : public testing::Test { - protected: - void SetUp() override { - population.id_g = id_g; - population.nyears = nyears; - population.nseasons = nseasons; - population.nages = nages; - population.nfleets = nfleets; - - // C++ code to set up true values for log_naa, log_M, - // log_Fmort, and log_q: - int seed = 1234; - std::default_random_engine generator(seed); - - // log_Fmort - double log_Fmort_min = fims_math::log(0.1); - double log_Fmort_max = fims_math::log(2.3); - std::uniform_real_distribution log_Fmort_distribution( - log_Fmort_min, log_Fmort_max); - - // log_q - double log_q_min = fims_math::log(0.1); - double log_q_max = fims_math::log(1); - std::uniform_real_distribution log_q_distribution(log_q_min, - log_q_max); - - // Make a shared pointer to selectivity and fleet because - // fleet object needs a shared pointer in fleet.hpp - // (std::shared_ptr > selectivity;) - // and population object needs a shared pointer in population.hpp - // (std::vector > > fleets;) - - // Does Fmort need to be in side of the year loop like log_q? - for (int i = 0; i < nfleets; i++) { - auto fleet = std::make_shared>(); - auto selectivity = - std::make_shared>(); - selectivity->inflection_point.resize(1); - selectivity->inflection_point[0] = 7; - selectivity->slope.resize(1); - selectivity->slope[0] = 0.5; - - - fleet->expected_catch.resize(nyears); - fleet->expected_index.resize(nyears); - fleet->catch_numbers_at_age.resize(nyears * nages); - fleet->Initialize(nyears, nages); - fleet->selectivity = selectivity; - fleet->log_q = log_q_distribution(generator); - for (int year = 0; year < nyears; year++) { - fleet->log_Fmort[year] = log_Fmort_distribution(generator); - } - if (i == 0) { - fleet->is_survey = true; - } - fleet->Prepare(); - population.fleets.push_back(fleet); - } - population.numbers_at_age.resize((nyears + 1) * nages); - population.Initialize(nyears, nseasons, nages); - - for (int i = 0; i < nages; i++) { - population.ages[i] = i + 1; - } - - // log_naa - double log_init_naa_min = 10.0; - double log_init_naa_max = 12.0; - std::uniform_real_distribution log_naa_distribution( - log_init_naa_min, log_init_naa_max); - for (int i = 0; i < nages; i++) { - population.log_init_naa[i] = log_naa_distribution(generator); - } - - // prop_female - double prop_female_min = 0.1; - double prop_female_max = 0.9; - std::uniform_real_distribution prop_female_distribution( - prop_female_min, prop_female_max); - for (int i = 0; i < nages; i++) { - population.proportion_female[i] = prop_female_distribution(generator); - } - - // log_M - double log_M_min = fims_math::log(0.1); - double log_M_max = fims_math::log(0.3); - std::uniform_real_distribution log_M_distribution(log_M_min, - log_M_max); - for (int i = 0; i < nyears * nages; i++) { - population.log_M[i] = log_M_distribution(generator); - } - - // numbers_at_age - double numbers_at_age_min = fims_math::exp(10.0); - double numbers_at_age_max = fims_math::exp(12.0); - std::uniform_real_distribution numbers_at_age_distribution( - numbers_at_age_min, numbers_at_age_max); - for (int i = 0; i < (nyears + 1) * nages; i++) { - population.numbers_at_age[i] = numbers_at_age_distribution(generator); - } - - // weight_at_age - double weight_at_age_min = 0.5; - double weight_at_age_max = 12.0; - - std::shared_ptr> growth = - std::make_shared>(); - std::uniform_real_distribution weight_at_age_distribution( - weight_at_age_min, weight_at_age_max); - for (int i = 0; i < nages; i++) { - growth->ewaa[static_cast(population.ages[i])] = - weight_at_age_distribution(generator); - } - - population.growth = growth; - population.Prepare(); - - auto maturity = std::make_shared>(); - maturity->inflection_point.resize(1); - maturity->inflection_point[0] = 6; - maturity->slope.resize(1); - maturity->slope[0] = 0.15; - population.maturity = maturity; - - auto recruitment = std::make_shared>(); - recruitment->logit_steep.resize(1); - recruitment->log_rzero.resize(1); - recruitment->logit_steep[0] = fims_math::logit(0.2, 1.0, 0.75); - recruitment->log_rzero[0] = fims_math::log(1000000.0); - /*the log_recruit_dev vector does not include a value for year == 0 - and is of length nyears - 1 where the first position of the vector - corresponds to the second year of the time series.*/ - recruitment->log_recruit_devs.resize(nyears - 1); - for (int i = 0; i < recruitment->log_recruit_devs.size(); i++) { - recruitment->log_recruit_devs[i] = 0.0; - } - population.recruitment = recruitment; - - int year = 4; - int age = 6; - int i_age_year = year * population.nages + age; - int i_agem1_yearm1 = (year - 1) * population.nages + age - 1; - - population.CalculateMortality(i_age_year, year, age); - population.CalculateNumbersAA(i_age_year, i_agem1_yearm1, age); - } - - virtual void TearDown() {} - - fims_popdy::Population population; - int id_g = 0; - int nyears = 30; - int nseasons = 1; - int nages = 12; - int nfleets = 2; - - int year = 4; - int age = 6; - int i_age_year = year * nages + age; - int i_agem1_yearm1 = (year - 1) * nages + age - 1; -}; - -class PopulationPrepareTestFixture : public testing::Test { - protected: - void SetUp() override { - population.id_g = id_g; - population.nyears = nyears; - population.nseasons = nseasons; - population.nages = nages; - population.nfleets = nfleets; - - // C++ code to set up true values for log_Fmort, and log_q: - int seed = 1234; - std::default_random_engine generator(seed); - - // log_Fmort - double log_Fmort_min = fims_math::log(0.1); - double log_Fmort_max = fims_math::log(2.3); - std::uniform_real_distribution log_Fmort_distribution( - log_Fmort_min, log_Fmort_max); - - // log_q - double log_q_min = fims_math::log(0.1); - double log_q_max = fims_math::log(1); - std::uniform_real_distribution log_q_distribution(log_q_min, - log_q_max); - - // Make a shared pointer to selectivity and fleet because - // fleet object needs a shared pointer in fleet.hpp - // (std::shared_ptr > selectivity;) - // and population object needs a shared pointer in population.hpp - // (std::vector > > fleets;) - - for (int i = 0; i < nfleets; i++) { - auto fleet = std::make_shared>(); - auto selectivity = - std::make_shared>(); - selectivity->inflection_point.resize(1); - selectivity->slope.resize(1); - selectivity->inflection_point[0] = 7; - selectivity->slope[0] = 0.5; - - - fleet->expected_catch.resize(nyears); - fleet->expected_index.resize(nyears); - fleet->catch_numbers_at_age.resize(nyears * nages); - fleet->Initialize(nyears, nages); - fleet->selectivity = selectivity; - fleet->log_q = log_q_distribution(generator); - for (int year = 0; year < nyears; year++) { - fleet->log_Fmort[year] = log_Fmort_distribution(generator); - } - if (i == 0) { - fleet->is_survey = true; - } - fleet->Prepare(); - population.fleets.push_back(fleet); - } - - population.numbers_at_age.resize((nyears + 1) * nages); - population.Initialize(nyears, nseasons, nages); - - for (int i = 0; i < nages; i++) { - population.ages[i] = i + 1; - } - - // log_naa - double log_init_naa_min = 10.0; - double log_init_naa_max = 12.0; - std::uniform_real_distribution log_naa_distribution( - log_init_naa_min, log_init_naa_max); - for (int i = 0; i < nages; i++) { - population.log_init_naa[i] = log_naa_distribution(generator); - } - - // prop_female - double prop_female_min = 0.1; - double prop_female_max = 0.9; - std::uniform_real_distribution prop_female_distribution( - prop_female_min, prop_female_max); - for (int i = 0; i < nages; i++) { - population.proportion_female[i] = prop_female_distribution(generator); - } - - // log_M - double log_M_min = fims_math::log(0.1); - double log_M_max = fims_math::log(0.3); - std::uniform_real_distribution log_M_distribution(log_M_min, - log_M_max); - for (int i = 0; i < nyears * nages; i++) { - population.log_M[i] = log_M_distribution(generator); - } - - // numbers_at_age - double numbers_at_age_min = fims_math::exp(10.0); - double numbers_at_age_max = fims_math::exp(12.0); - std::uniform_real_distribution numbers_at_age_distribution( - numbers_at_age_min, numbers_at_age_max); - for (int i = 0; i < (nyears + 1) * nages; i++) { - population.numbers_at_age[i] = numbers_at_age_distribution(generator); - } - - // weight_at_age - double weight_at_age_min = 0.5; - double weight_at_age_max = 12.0; - - std::shared_ptr> growth = - std::make_shared>(); - std::uniform_real_distribution weight_at_age_distribution( - weight_at_age_min, weight_at_age_max); - for (int i = 0; i < nages; i++) { - growth->ewaa[static_cast(population.ages[i])] = - weight_at_age_distribution(generator); - } - - population.growth = growth; - - population.Prepare(); - } - - virtual void TearDown() {} - - fims_popdy::Population population; - int id_g = 0; - int nyears = 30; - int nseasons = 1; - int nages = 12; - int nfleets = 2; -}; -} // namespace \ No newline at end of file + // Use test fixture to reuse the same configuration of objects for + // several different tests. To use a test fixture, derive a class + // from testing::Test. + + class PopulationInitializeTestFixture : public testing::Test { + // Make members protected and they can be accessed from + // sub-classes. + protected: + // Use SetUp function to prepare the objects for each test. + // Use override in C++11 to make sure SetUp (e.g., not Setup with + // a lowercase u) is spelled + // correctly. + + void SetUp() override { + population.id_g = id_g; + population.nyears = nyears; + population.nseasons = nseasons; + population.nages = nages; + for (int i = 0; i < nfleets; i++) { + auto fleet = std::make_shared>(); + fleet->log_q.resize(1); + population.fleets.push_back(fleet); + } + } + + // Virtual void TearDown() will be called after each test is + // run. It needs to be defined if there is clearup work to + // do. Otherwise, it does not need to be provided. + + virtual void TearDown() { + } + + fims_popdy::Population population; + + // Use default values from the Li et al., 2021 + // https://github.com/Bai-Li-NOAA/Age_Structured_Stock_Assessment_Model_Comparison/blob/master/R/save_initial_input.R + int id_g = 0; + int nyears = 30; + int nseasons = 1; + int nages = 12; + int nfleets = 2; + }; + + class PopulationEvaluateTestFixture : public testing::Test { + protected: + + void SetUp() override { + population.id_g = id_g; + population.nyears = nyears; + population.nseasons = nseasons; + population.nages = nages; + population.nfleets = nfleets; + + // C++ code to set up true values for log_naa, log_M, + // log_Fmort, and log_q: + int seed = 1234; + std::default_random_engine generator(seed); + + // log_Fmort + double log_Fmort_min = fims_math::log(0.1); + double log_Fmort_max = fims_math::log(2.3); + std::uniform_real_distribution log_Fmort_distribution( + log_Fmort_min, log_Fmort_max); + + // log_q + double log_q_min = fims_math::log(0.1); + double log_q_max = fims_math::log(1); + std::uniform_real_distribution log_q_distribution(log_q_min, + log_q_max); + + // Make a shared pointer to selectivity and fleet because + // fleet object needs a shared pointer in fleet.hpp + // (std::shared_ptr > selectivity;) + // and population object needs a shared pointer in population.hpp + // (std::vector > > fleets;) + + // Does Fmort need to be in side of the year loop like log_q? + for (int i = 0; i < nfleets; i++) { + auto fleet = std::make_shared>(); + auto selectivity = + std::make_shared>(); + selectivity->inflection_point.resize(1); + selectivity->inflection_point[0] = 7; + selectivity->slope.resize(1); + selectivity->slope[0] = 0.5; + + + fleet->expected_catch.resize(nyears); + fleet->expected_index.resize(nyears); + fleet->catch_numbers_at_age.resize(nyears * nages); + fleet->log_q.resize(1); + fleet->Initialize(nyears, nages); + fleet->selectivity = selectivity; + fleet->log_q[0] = log_q_distribution(generator); + for (int year = 0; year < nyears; year++) { + fleet->log_Fmort[year] = log_Fmort_distribution(generator); + } + if (i == 0) { + fleet->is_survey = true; + } + fleet->Prepare(); + population.fleets.push_back(fleet); + } + population.numbers_at_age.resize((nyears + 1) * nages); + try { + population.Initialize(nyears, nseasons, nages); + } catch (std::exception& e) { + std::cout << e.what() << "\n"; + } + + for (int i = 0; i < nages; i++) { + population.ages[i] = i + 1; + } + + // log_naa + double log_init_naa_min = 10.0; + double log_init_naa_max = 12.0; + std::uniform_real_distribution log_naa_distribution( + log_init_naa_min, log_init_naa_max); + for (int i = 0; i < nages; i++) { + population.log_init_naa[i] = log_naa_distribution(generator); + } + + // prop_female + double prop_female_min = 0.1; + double prop_female_max = 0.9; + std::uniform_real_distribution prop_female_distribution( + prop_female_min, prop_female_max); + for (int i = 0; i < nages; i++) { + population.proportion_female[i] = prop_female_distribution(generator); + } + + // log_M + double log_M_min = fims_math::log(0.1); + double log_M_max = fims_math::log(0.3); + std::uniform_real_distribution log_M_distribution(log_M_min, + log_M_max); + for (int i = 0; i < nyears * nages; i++) { + population.log_M[i] = log_M_distribution(generator); + } + + // numbers_at_age + double numbers_at_age_min = fims_math::exp(10.0); + double numbers_at_age_max = fims_math::exp(12.0); + std::uniform_real_distribution numbers_at_age_distribution( + numbers_at_age_min, numbers_at_age_max); + for (int i = 0; i < (nyears + 1) * nages; i++) { + population.numbers_at_age[i] = numbers_at_age_distribution(generator); + } + + // weight_at_age + double weight_at_age_min = 0.5; + double weight_at_age_max = 12.0; + + std::shared_ptr> growth = + std::make_shared>(); + std::uniform_real_distribution weight_at_age_distribution( + weight_at_age_min, weight_at_age_max); + for (int i = 0; i < nages; i++) { + growth->ewaa[static_cast (population.ages[i])] = + weight_at_age_distribution(generator); + } + + population.growth = growth; + population.Prepare(); + + auto maturity = std::make_shared>(); + maturity->inflection_point.resize(1); + maturity->inflection_point[0] = 6; + maturity->slope.resize(1); + maturity->slope[0] = 0.15; + population.maturity = maturity; + + auto recruitment = std::make_shared>(); + recruitment->logit_steep.resize(1); + recruitment->log_rzero.resize(1); + recruitment->logit_steep[0] = fims_math::logit(0.2, 1.0, 0.75); + recruitment->log_rzero[0] = fims_math::log(1000000.0); + /*the log_recruit_dev vector does not include a value for year == 0 + and is of length nyears - 1 where the first position of the vector + corresponds to the second year of the time series.*/ + recruitment->log_recruit_devs.resize(nyears - 1); + for (int i = 0; i < recruitment->log_recruit_devs.size(); i++) { + recruitment->log_recruit_devs[i] = 0.0; + } + population.recruitment = recruitment; + + int year = 4; + int age = 6; + int i_age_year = year * population.nages + age; + int i_agem1_yearm1 = (year - 1) * population.nages + age - 1; + + population.CalculateMortality(i_age_year, year, age); + population.CalculateNumbersAA(i_age_year, i_agem1_yearm1, age); + } + + virtual void TearDown() { + } + + fims_popdy::Population population; + int id_g = 0; + int nyears = 30; + int nseasons = 1; + int nages = 12; + int nfleets = 2; + + int year = 4; + int age = 6; + int i_age_year = year * nages + age; + int i_agem1_yearm1 = (year - 1) * nages + age - 1; + }; + + class PopulationPrepareTestFixture : public testing::Test { + protected: + + void SetUp() override { + population.id_g = id_g; + population.nyears = nyears; + population.nseasons = nseasons; + population.nages = nages; + population.nfleets = nfleets; + + // C++ code to set up true values for log_Fmort, and log_q: + int seed = 1234; + std::default_random_engine generator(seed); + + // log_Fmort + double log_Fmort_min = fims_math::log(0.1); + double log_Fmort_max = fims_math::log(2.3); + std::uniform_real_distribution log_Fmort_distribution( + log_Fmort_min, log_Fmort_max); + + // log_q + double log_q_min = fims_math::log(0.1); + double log_q_max = fims_math::log(1); + std::uniform_real_distribution log_q_distribution(log_q_min, + log_q_max); + + // Make a shared pointer to selectivity and fleet because + // fleet object needs a shared pointer in fleet.hpp + // (std::shared_ptr > selectivity;) + // and population object needs a shared pointer in population.hpp + // (std::vector > > fleets;) + + for (int i = 0; i < nfleets; i++) { + auto fleet = std::make_shared>(); + auto selectivity = + std::make_shared>(); + selectivity->inflection_point.resize(1); + selectivity->slope.resize(1); + selectivity->inflection_point[0] = 7; + selectivity->slope[0] = 0.5; + + + fleet->expected_catch.resize(nyears); + fleet->expected_index.resize(nyears); + fleet->catch_numbers_at_age.resize(nyears * nages); + fleet->log_q.resize(1); + fleet->Initialize(nyears, nages); + fleet->selectivity = selectivity; + fleet->log_q[0] = log_q_distribution(generator); + for (int year = 0; year < nyears; year++) { + fleet->log_Fmort[year] = log_Fmort_distribution(generator); + } + if (i == 0) { + fleet->is_survey = true; + } + fleet->Prepare(); + population.fleets.push_back(fleet); + } + + population.numbers_at_age.resize((nyears + 1) * nages); + population.Initialize(nyears, nseasons, nages); + + for (int i = 0; i < nages; i++) { + population.ages[i] = i + 1; + } + + // log_naa + double log_init_naa_min = 10.0; + double log_init_naa_max = 12.0; + std::uniform_real_distribution log_naa_distribution( + log_init_naa_min, log_init_naa_max); + for (int i = 0; i < nages; i++) { + population.log_init_naa[i] = log_naa_distribution(generator); + } + + // prop_female + double prop_female_min = 0.1; + double prop_female_max = 0.9; + std::uniform_real_distribution prop_female_distribution( + prop_female_min, prop_female_max); + for (int i = 0; i < nages; i++) { + population.proportion_female[i] = prop_female_distribution(generator); + } + + // log_M + double log_M_min = fims_math::log(0.1); + double log_M_max = fims_math::log(0.3); + std::uniform_real_distribution log_M_distribution(log_M_min, + log_M_max); + for (int i = 0; i < nyears * nages; i++) { + population.log_M[i] = log_M_distribution(generator); + } + + // numbers_at_age + double numbers_at_age_min = fims_math::exp(10.0); + double numbers_at_age_max = fims_math::exp(12.0); + std::uniform_real_distribution numbers_at_age_distribution( + numbers_at_age_min, numbers_at_age_max); + for (int i = 0; i < (nyears + 1) * nages; i++) { + population.numbers_at_age[i] = numbers_at_age_distribution(generator); + } + + // weight_at_age + double weight_at_age_min = 0.5; + double weight_at_age_max = 12.0; + + std::shared_ptr> growth = + std::make_shared>(); + std::uniform_real_distribution weight_at_age_distribution( + weight_at_age_min, weight_at_age_max); + for (int i = 0; i < nages; i++) { + growth->ewaa[static_cast (population.ages[i])] = + weight_at_age_distribution(generator); + } + + population.growth = growth; + + population.Prepare(); + } + + virtual void TearDown() { + } + + fims_popdy::Population population; + int id_g = 0; + int nyears = 30; + int nseasons = 1; + int nages = 12; + int nfleets = 2; + }; +} // namespace \ No newline at end of file diff --git a/tests/integration/integration_class.hpp b/tests/integration/integration_class.hpp index fa79f8a3..50d27c9d 100644 --- a/tests/integration/integration_class.hpp +++ b/tests/integration/integration_class.hpp @@ -46,7 +46,11 @@ class IntegrationTest { good = false; } + if(good){ this->RunModelLoop(pop, input); + }else{ + throw std::invalid_argument("model not good!"); + } // if (!this->CheckModelOutput(pop, output)) { // good = false; @@ -152,6 +156,7 @@ class IntegrationTest { for (size_t i = 0; i < nfleets; i++) { std::shared_ptr > f = std::make_shared >(); + f->log_q.resize(1); f->Initialize(nyears, nages); // f->observed_index_data = std::make_shared >(nyears); // f->observed_agecomp_data = std::make_shared >(nyears, nages); @@ -261,7 +266,7 @@ class IntegrationTest { } - f->log_q = 0.0; + f->log_q[0] = 0.0; it = obj.find("f"); if ((*it).second.GetType() == JsonValueType::Array) { JsonArray f_values = (*it).second.GetArray(); @@ -302,6 +307,7 @@ class IntegrationTest { for (size_t i = 0; i < nsurveys; i++) { std::shared_ptr > s = std::make_shared >(); s->is_survey = true; + s->log_q.resize(1); s->Initialize(nyears, nages); // s->observed_index_data = std::make_shared >(nyears); // s->observed_agecomp_data = std::make_shared >(nyears, nages); @@ -389,7 +395,7 @@ class IntegrationTest { } - s->log_q = 0.0; + s->log_q[0] = 0.0; it = obj2.find("survey_q"); @@ -401,7 +407,7 @@ class IntegrationTest { if ((*qit).second.GetType() == JsonValueType::Array) { JsonArray a = (*qit).second.GetArray(); - s->log_q = fims_math::log(a[0].GetDouble()); + s->log_q[0] = fims_math::log(a[0].GetDouble()); if (this->print_statements) { std::cout << "q = " << a[0].GetDouble() << "\nlog(q) = " << s->log_q << "\n"; } diff --git a/tests/testthat/helper-integration-tests-setup.R b/tests/testthat/helper-integration-tests-setup.R index 10176b52..c1671f82 100644 --- a/tests/testthat/helper-integration-tests-setup.R +++ b/tests/testthat/helper-integration-tests-setup.R @@ -69,13 +69,13 @@ setup_and_run_FIMS <- function(iter_id, # place as appropriate. # set up log_rzero (equilibrium recruitment) - recruitment$log_rzero$value <- log(om_input$R0) - recruitment$log_rzero$is_random_effect <- FALSE - recruitment$log_rzero$estimated <- TRUE + recruitment$log_rzero[1]$value <- log(om_input$R0) + recruitment$log_rzero[1]$is_random_effect <- FALSE + recruitment$log_rzero[1]$estimated <- TRUE # set up logit_steep - recruitment$logit_steep$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) - recruitment$logit_steep$is_random_effect <- FALSE - recruitment$logit_steep$estimated <- FALSE + recruitment$logit_steep[1]$value <- -log(1.0 - om_input$h) + log(om_input$h - 0.2) + recruitment$logit_steep[1]$is_random_effect <- FALSE + recruitment$logit_steep[1]$estimated <- FALSE # turn on estimation of deviations # recruit deviations should enter the model in normal space. # The log is taken in the likelihood calculations @@ -122,31 +122,31 @@ setup_and_run_FIMS <- function(iter_id, # Maturity maturity <- new(LogisticMaturity) - maturity$inflection_point$value <- om_input$A50.mat - maturity$inflection_point$is_random_effect <- FALSE - maturity$inflection_point$estimated <- FALSE - maturity$slope$value <- om_input$slope - maturity$slope$is_random_effect <- FALSE - maturity$slope$estimated <- FALSE + maturity$inflection_point[1]$value <- om_input$A50.mat + maturity$inflection_point[1]$is_random_effect <- FALSE + maturity$inflection_point[1]$estimated <- FALSE + maturity$slope[1]$value <- om_input$slope + maturity$slope[1]$is_random_effect <- FALSE + maturity$slope[1]$estimated <- FALSE # Fleet # Create the fishing fleet fishing_fleet_selectivity <- new(LogisticSelectivity) - fishing_fleet_selectivity$inflection_point$value <- om_input$sel_fleet$fleet1$A50.sel1 - fishing_fleet_selectivity$inflection_point$is_random_effect <- FALSE + fishing_fleet_selectivity$inflection_point[1]$value <- om_input$sel_fleet$fleet1$A50.sel1 + fishing_fleet_selectivity$inflection_point[1]$is_random_effect <- FALSE # turn on estimation of inflection_point - fishing_fleet_selectivity$inflection_point$estimated <- TRUE - fishing_fleet_selectivity$slope$value <- om_input$sel_fleet$fleet1$slope.sel1 + fishing_fleet_selectivity$inflection_point[1]$estimated <- TRUE + fishing_fleet_selectivity$slope[1]$value <- om_input$sel_fleet$fleet1$slope.sel1 # turn on estimation of slope - fishing_fleet_selectivity$slope$is_random_effect <- FALSE - fishing_fleet_selectivity$slope$estimated <- TRUE + fishing_fleet_selectivity$slope[1]$is_random_effect <- FALSE + fishing_fleet_selectivity$slope[1]$estimated <- TRUE fishing_fleet <- new(Fleet) fishing_fleet$nages <- om_input$nages fishing_fleet$nyears <- om_input$nyr fishing_fleet$log_Fmort <- methods::new(ParameterVector, log(om_output$f), om_input$nyr) fishing_fleet$log_Fmort$set_all_estimable(TRUE) - fishing_fleet$log_q <- log(1.0) + fishing_fleet$log_q[1]$value <- log(1.0) fishing_fleet$estimate_q <- FALSE fishing_fleet$random_q <- FALSE fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) @@ -170,20 +170,20 @@ setup_and_run_FIMS <- function(iter_id, # Create the survey fleet survey_fleet_selectivity <- new(LogisticSelectivity) - survey_fleet_selectivity$inflection_point$value <- om_input$sel_survey$survey1$A50.sel1 - survey_fleet_selectivity$inflection_point$is_random_effect <- FALSE + survey_fleet_selectivity$inflection_point[1]$value <- om_input$sel_survey$survey1$A50.sel1 + survey_fleet_selectivity$inflection_point[1]$is_random_effect <- FALSE # turn on estimation of inflection_point - survey_fleet_selectivity$inflection_point$estimated <- TRUE - survey_fleet_selectivity$slope$value <- om_input$sel_survey$survey1$slope.sel1 - survey_fleet_selectivity$slope$is_random_effect <- FALSE + survey_fleet_selectivity$inflection_point[1]$estimated <- TRUE + survey_fleet_selectivity$slope[1]$value <- om_input$sel_survey$survey1$slope.sel1 + survey_fleet_selectivity$slope[1]$is_random_effect <- FALSE # turn on estimation of slope - survey_fleet_selectivity$slope$estimated <- TRUE + survey_fleet_selectivity$slope[1]$estimated <- TRUE survey_fleet <- new(Fleet) survey_fleet$is_survey <- TRUE survey_fleet$nages <- om_input$nages survey_fleet$nyears <- om_input$nyr - survey_fleet$log_q <- log(om_output$survey_q$survey1) + survey_fleet$log_q[1]$value <- log(om_output$survey_q$survey1) survey_fleet$estimate_q <- TRUE survey_fleet$random_q <- FALSE survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) diff --git a/tests/testthat/test-rcpp-fims.R b/tests/testthat/test-rcpp-fims.R index 58cd924a..f21e32b2 100644 --- a/tests/testthat/test-rcpp-fims.R +++ b/tests/testthat/test-rcpp-fims.R @@ -3,10 +3,10 @@ test_that("Rcpp interface works for modules", { expect_no_error(beverton_holt <- new(BevertonHoltRecruitment)) expect_no_error(logistic_selectivity <- new(LogisticSelectivity)) expect_no_error(ewaa_growth <- new(EWAAgrowth)) - logistic_selectivity$slope$value <- .7 - logistic_selectivity$inflection_point$value <- 5.0 + logistic_selectivity$slope[1]$value <- .7 + logistic_selectivity$inflection_point[1]$value <- 5.0 - expect_equal(logistic_selectivity$slope$value, 0.7) + expect_equal(logistic_selectivity$slope[1]$value, 0.7) expect_equal(logistic_selectivity$get_id(), 1) ewaa_growth$ages <- 1.0 ewaa_growth$weights <- 2.5 diff --git a/tests/testthat/test-rcpp-get_fixed.R b/tests/testthat/test-rcpp-get_fixed.R index fc76ab2a..3d026d5c 100644 --- a/tests/testthat/test-rcpp-get_fixed.R +++ b/tests/testthat/test-rcpp-get_fixed.R @@ -1,18 +1,18 @@ test_that("test get parameter vector", { # Create selectivity selectivity <- new(LogisticSelectivity) - selectivity$inflection_point$value <- 10.0 - selectivity$inflection_point$min <- 8.0 - selectivity$inflection_point$max <- 12.0 - selectivity$inflection_point$is_random_effect <- FALSE - selectivity$inflection_point$estimated <- TRUE - selectivity$slope$value <- 0.2 - selectivity$slope$is_random_effect <- FALSE - selectivity$slope$estimated <- TRUE + selectivity$inflection_point[1]$value <- 10.0 + selectivity$inflection_point[1]$min <- 8.0 + selectivity$inflection_point[1]$max <- 12.0 + selectivity$inflection_point[1]$is_random_effect <- FALSE + selectivity$inflection_point[1]$estimated <- TRUE + selectivity$slope[1]$value <- 0.2 + selectivity$slope[1]$is_random_effect <- FALSE + selectivity$slope[1]$estimated <- TRUE CreateTMBModel() p <- get_fixed() - sel_parm <- c(selectivity$inflection_point$value, selectivity$slope$value) + sel_parm <- c(selectivity$inflection_point[1]$value, selectivity$slope[1]$value) expect_equal(sel_parm, p) # test fims clear @@ -29,28 +29,28 @@ test_that("test get parameter vector", { p <- get_fixed() expect_equal(numeric(0), p) selectivity <- new(LogisticSelectivity) - selectivity$inflection_point$value <- 11.0 - selectivity$inflection_point$min <- 8.0 - selectivity$inflection_point$max <- 12.0 - selectivity$inflection_point$is_random_effect <- FALSE - selectivity$inflection_point$estimated <- TRUE - selectivity$slope$value <- 0.5 - selectivity$slope$is_random_effect <- FALSE - selectivity$slope$estimated <- TRUE - sel_parm <- c(selectivity$inflection_point$value, selectivity$slope$value) + selectivity$inflection_point[1]$value <- 11.0 + selectivity$inflection_point[1]$min <- 8.0 + selectivity$inflection_point[1]$max <- 12.0 + selectivity$inflection_point[1]$is_random_effect <- FALSE + selectivity$inflection_point[1]$estimated <- TRUE + selectivity$slope[1]$value <- 0.5 + selectivity$slope[1]$is_random_effect <- FALSE + selectivity$slope[1]$estimated <- TRUE + sel_parm <- c(selectivity$inflection_point[1]$value, selectivity$slope[1]$value) recruitment <- new(BevertonHoltRecruitment) h <- 0.75 r0 <- 1000000.0 spawns <- 9.55784 * 10^6 ssb0 <- 0.0102562 - recruitment$logit_steep$value <- -log(1.0 - h) + log(h - 0.2) - recruitment$logit_steep$min <- 0.21 - recruitment$logit_steep$max <- 1.0 - recruitment$logit_steep$is_random_effect <- FALSE - recruitment$logit_steep$estimated <- TRUE - recruitment$log_rzero$value <- log(r0) - recruitment$log_rzero$is_random_effect <- FALSE - recruitment$log_rzero$estimated <- TRUE + recruitment$logit_steep[1]$value <- -log(1.0 - h) + log(h - 0.2) + recruitment$logit_steep[1]$min <- 0.21 + recruitment$logit_steep[1]$max <- 1.0 + recruitment$logit_steep[1]$is_random_effect <- FALSE + recruitment$logit_steep[1]$estimated <- TRUE + recruitment$log_rzero[1]$value <- log(r0) + recruitment$log_rzero[1]$is_random_effect <- FALSE + recruitment$log_rzero[1]$estimated <- TRUE rec_parm <- c(-log(1.0 - h) + log(h - 0.2), log(r0)) CreateTMBModel() diff --git a/tests/testthat/test-rcpp-maturity-interface.R b/tests/testthat/test-rcpp-maturity-interface.R index 4ded4336..7d5ddbcd 100644 --- a/tests/testthat/test-rcpp-maturity-interface.R +++ b/tests/testthat/test-rcpp-maturity-interface.R @@ -2,20 +2,20 @@ test_that("Maturity input settings work as expected", { # Create maturity1 maturity1 <- new(LogisticMaturity) - maturity1$inflection_point$value <- 10.0 - maturity1$inflection_point$min <- 8.0 - maturity1$inflection_point$max <- 12.0 - maturity1$inflection_point$is_random_effect <- TRUE - maturity1$inflection_point$estimated <- TRUE - maturity1$slope$value <- 0.2 + maturity1$inflection_point[1]$value <- 10.0 + maturity1$inflection_point[1]$min <- 8.0 + maturity1$inflection_point[1]$max <- 12.0 + maturity1$inflection_point[1]$is_random_effect <- TRUE + maturity1$inflection_point[1]$estimated <- TRUE + maturity1$slope[1]$value <- 0.2 expect_equal(maturity1$get_id(), 1) - expect_equal(maturity1$inflection_point$value, 10.0) - expect_equal(maturity1$inflection_point$min, 8.0) - expect_equal(maturity1$inflection_point$max, 12.0) - expect_true(maturity1$inflection_point$is_random_effect) - expect_true(maturity1$inflection_point$estimated) - expect_equal(maturity1$slope$value, 0.2) + expect_equal(maturity1$inflection_point[1]$value, 10.0) + expect_equal(maturity1$inflection_point[1]$min, 8.0) + expect_equal(maturity1$inflection_point[1]$max, 12.0) + expect_true(maturity1$inflection_point[1]$is_random_effect) + expect_true(maturity1$inflection_point[1]$estimated) + expect_equal(maturity1$slope[1]$value, 0.2) expect_equal(maturity1$evaluate(10.0), 0.5) diff --git a/tests/testthat/test-rcpp-recruitment-interface.R b/tests/testthat/test-rcpp-recruitment-interface.R index 3d8242db..dadc7d27 100644 --- a/tests/testthat/test-rcpp-recruitment-interface.R +++ b/tests/testthat/test-rcpp-recruitment-interface.R @@ -7,21 +7,21 @@ test_that("Recruitment input settings work as expected", { spawns <- 9.55784 * 10^6 ssb0 <- 0.0102562 - recruitment$logit_steep$value <- -log(1.0 - h) + log(h - 0.2) - recruitment$logit_steep$min <- 0.21 - recruitment$logit_steep$max <- 1.0 - recruitment$logit_steep$is_random_effect <- TRUE - recruitment$logit_steep$estimated <- TRUE - recruitment$log_rzero$value <- log(r0) - recruitment$log_sigma_recruit$value <- log(0.7) + recruitment$logit_steep[1]$value <- -log(1.0 - h) + log(h - 0.2) + recruitment$logit_steep[1]$min <- 0.21 + recruitment$logit_steep[1]$max <- 1.0 + recruitment$logit_steep[1]$is_random_effect <- TRUE + recruitment$logit_steep[1]$estimated <- TRUE + recruitment$log_rzero[1]$value <- log(r0) + recruitment$log_sigma_recruit[1]$value <- log(0.7) expect_equal(recruitment$get_id(), 1) - expect_equal(recruitment$logit_steep$value, 0.78845736) - expect_equal(recruitment$logit_steep$min, 0.21) - expect_equal(recruitment$logit_steep$max, 1.0) - expect_true(recruitment$logit_steep$is_random_effect) - expect_true(recruitment$logit_steep$estimated) - expect_equal(recruitment$log_rzero$value, log(1000000.0)) + expect_equal(recruitment$logit_steep[1]$value, 0.78845736) + expect_equal(recruitment$logit_steep[1]$min, 0.21) + expect_equal(recruitment$logit_steep[1]$max, 1.0) + expect_true(recruitment$logit_steep[1]$is_random_effect) + expect_true(recruitment$logit_steep[1]$estimated) + expect_equal(recruitment$log_rzero[1]$value, log(1000000.0)) expect_equal(object = recruitment$evaluate(spawns, ssb0), expected = 1090802.68) diff --git a/tests/testthat/test-rcpp-selectivity-interface.R b/tests/testthat/test-rcpp-selectivity-interface.R index 398fe0c9..a58da95b 100644 --- a/tests/testthat/test-rcpp-selectivity-interface.R +++ b/tests/testthat/test-rcpp-selectivity-interface.R @@ -2,20 +2,20 @@ test_that("Selectivity input settings work as expected", { # Create selectivity1 selectivity1 <- new(LogisticSelectivity) - selectivity1$inflection_point$value <- 10.0 - selectivity1$inflection_point$min <- 8.0 - selectivity1$inflection_point$max <- 12.0 - selectivity1$inflection_point$is_random_effect <- TRUE - selectivity1$inflection_point$estimated <- TRUE - selectivity1$slope$value <- 0.2 + selectivity1$inflection_point[1]$value <- 10.0 + selectivity1$inflection_point[1]$min <- 8.0 + selectivity1$inflection_point[1]$max <- 12.0 + selectivity1$inflection_point[1]$is_random_effect <- TRUE + selectivity1$inflection_point[1]$estimated <- TRUE + selectivity1$slope[1]$value <- 0.2 expect_equal(selectivity1$get_id(), 1) - expect_equal(selectivity1$inflection_point$value, 10.0) - expect_equal(selectivity1$inflection_point$min, 8.0) - expect_equal(selectivity1$inflection_point$max, 12.0) - expect_true(selectivity1$inflection_point$is_random_effect) - expect_true(selectivity1$inflection_point$estimated) - expect_equal(selectivity1$slope$value, 0.2) + expect_equal(selectivity1$inflection_point[1]$value, 10.0) + expect_equal(selectivity1$inflection_point[1]$min, 8.0) + expect_equal(selectivity1$inflection_point[1]$max, 12.0) + expect_true(selectivity1$inflection_point[1]$is_random_effect) + expect_true(selectivity1$inflection_point[1]$estimated) + expect_equal(selectivity1$slope[1]$value, 0.2) expect_equal(selectivity1$evaluate(10.0), 0.5) @@ -26,14 +26,14 @@ test_that("Selectivity input settings work as expected", { # Test double logistic selectivity3 <- new(DoubleLogisticSelectivity) - selectivity3$inflection_point_asc$value <- 10.5 - selectivity3$slope_asc$value <- 0.2 - selectivity3$inflection_point_desc$value <- 15.0 - selectivity3$slope_desc$value <- 0.05 + selectivity3$inflection_point_asc[1]$value <- 10.5 + selectivity3$slope_asc[1]$value <- 0.2 + selectivity3$inflection_point_desc[1]$value <- 15.0 + selectivity3$slope_desc[1]$value <- 0.05 expect_equal(selectivity3$get_id(), 3) - expect_equal(selectivity3$inflection_point_asc$value, 10.5) - expect_equal(selectivity3$slope_asc$value, 0.2) + expect_equal(selectivity3$inflection_point_asc[1]$value, 10.5) + expect_equal(selectivity3$slope_asc[1]$value, 0.2) # R code that generates true value for the test # 1.0/(1.0+exp(-(34.5-10.5)*0.2)) * (1.0 - 1.0/(1.0+exp(-(34.5-15)*0.05))) = 0.2716494 expect_equal(selectivity3$evaluate(34.5), 0.2716494, tolerance = 0.0000001) diff --git a/tests/testthat/test-unit-rcpp-interface-variable-vector.R b/tests/testthat/test-unit-rcpp-interface-variable-vector.R index ab86c3a1..a107b2bc 100644 --- a/tests/testthat/test-unit-rcpp-interface-variable-vector.R +++ b/tests/testthat/test-unit-rcpp-interface-variable-vector.R @@ -13,188 +13,188 @@ expect_equal(v0$at(1)$value, 0) v1 <- new(ParameterVector, v_size) v1$fill(v1_value) for(i in 1:v_size){ - expect_equal(v1$at(i)$value, v1_value) + expect_equal(v1$get(i-1)$value, v1_value) } # Test that constructor that takes vector and size works. v2 <- new(ParameterVector, rep(v2_value, v_size), v_size) for(i in 1:v_size){ - expect_equal(v2$at(i)$value, v2_value) + expect_equal(v2$get(i-1)$value, v2_value) } #Test plus operator works. v_plus_test <- v1 + v2 for(i in 1:v_size){ - expect_equal(v_plus_test$at(i)$value, (v1[i]$value + v2[i]$value)) + expect_equal(v_plus_test$get(i-1)$value, (v1[i]$value + v2[i]$value)) } #Test minus operator works. v_minus_test <- v1 - v2 for(i in 1:v_size){ - expect_equal(v_minus_test$at(i)$value, (v1[i]$value - v2[i]$value)) + expect_equal(v_minus_test$get(i-1)$value, (v1[i]$value - v2[i]$value)) } #Test mult operator works. v_mult_test <- v1 * v2 for(i in 1:v_size){ - expect_equal(v_mult_test$at(i)$value, (v1[i]$value * v2[i]$value)) + expect_equal(v_mult_test$get(i-1)$value, (v1[i]$value * v2[i]$value)) } #Test div operator works. v_div_test <- v1 / v2 for(i in 1:v_size){ - expect_equal(v_div_test$at(i)$value, (v1[i]$value / v2[i]$value)) + expect_equal(v_div_test$get(i-1)$value, (v1[i]$value / v2[i]$value)) } #Test pre scalar plus operator works. v_plus_test_scalar <- v2_value + v1 for(i in 1:v_size){ - expect_equal(v_plus_test_scalar$at(i)$value, (v2_value+ v1[i]$value )) + expect_equal(v_plus_test_scalar$get(i-1)$value, (v2_value+ v1[i]$value )) } #Test pre scalar minus operator works. v_minus_test_scalar <- v2_value - v1 for(i in 1:v_size){ - expect_equal(v_minus_test_scalar$at(i)$value, (v2_value- v1[i]$value )) + expect_equal(v_minus_test_scalar$get(i-1)$value, (v2_value- v1[i]$value )) } #Test pre scalar mult operator works. v_mult_test_scalar <- v2_value * v1 for(i in 1:v_size){ - expect_equal(v_mult_test_scalar$at(i)$value, (v2_value* v1[i]$value )) + expect_equal(v_mult_test_scalar$get(i-1)$value, (v2_value* v1[i]$value )) } #Test pre scalar div operator works. v_div_test_scalar <- v2_value / v1 for(i in 1:v_size){ - expect_equal(v_div_test_scalar$at(i)$value, (v2_value/ v1[i]$value )) + expect_equal(v_div_test_scalar$get(i-1)$value, (v2_value/ v1[i]$value )) } #Test post scalar plus operator works. v_plus_test_scalar <- v1 + v2_value for(i in 1:v_size){ - expect_equal(v_plus_test_scalar$at(i)$value, (v1[i]$value + v2_value)) + expect_equal(v_plus_test_scalar$get(i-1)$value, (v1[i]$value + v2_value)) } #Test post scalar minus operator works. v_minus_test_scalar <- v1 - v2_value for(i in 1:v_size){ - expect_equal(v_minus_test_scalar$at(i)$value, (v1[i]$value - v2_value)) + expect_equal(v_minus_test_scalar$get(i-1)$value, (v1[i]$value - v2_value)) } #Test post scalar mult operator works. v_mult_test_scalar <- v1 * v2_value for(i in 1:v_size){ - expect_equal(v_mult_test_scalar$at(i)$value, (v1[i]$value * v2_value)) + expect_equal(v_mult_test_scalar$get(i-1)$value, (v1[i]$value * v2_value)) } #Test post scalar div operator works. v_div_test_scalar <- v1 / v2_value for(i in 1:v_size){ - expect_equal(v_div_test_scalar$at(i)$value, (v1[i]$value / v2_value)) + expect_equal(v_div_test_scalar$get(i-1)$value, (v1[i]$value / v2_value)) } #Test acos function works. v_acos_test <- acos(v1) for(i in 1:v_size){ - expect_equal(v_acos_test$at(i)$value, acos(v1_value)) + expect_equal(v_acos_test$get(i-1)$value, acos(v1_value)) } #Test asin function works. v_asin_test <- asin(v1) for(i in 1:v_size){ - expect_equal(v_asin_test$at(i)$value, asin(v1_value)) + expect_equal(v_asin_test$get(i-1)$value, asin(v1_value)) } #Test atan function works. v_atan_test <- atan(v1) for(i in 1:v_size){ - expect_equal(v_atan_test$at(i)$value, atan(v1_value)) + expect_equal(v_atan_test$get(i-1)$value, atan(v1_value)) } #Test cos function works. v_cos_test <- cos(v1) for(i in 1:v_size){ - expect_equal(v_cos_test$at(i)$value, cos(v1_value)) + expect_equal(v_cos_test$get(i-1)$value, cos(v1_value)) } #Test cosh function works. v_cosh_test <- cosh(v1) for(i in 1:v_size){ - expect_equal(v_cosh_test$at(i)$value, cosh(v1_value)) + expect_equal(v_cosh_test$get(i-1)$value, cosh(v1_value)) } #Test sin function works. v_sin_test <- sin(v1) for(i in 1:v_size){ - expect_equal(v_sin_test$at(i)$value, sin(v1_value)) + expect_equal(v_sin_test$get(i-1)$value, sin(v1_value)) } #Test sinh function works. v_sinh_test <- sinh(v1) for(i in 1:v_size){ - expect_equal(v_sinh_test$at(i)$value, sinh(v1_value)) + expect_equal(v_sinh_test$get(i-1)$value, sinh(v1_value)) } #Test tan function works. v_tan_test <- tan(v1) for(i in 1:v_size){ - expect_equal(v_tan_test$at(i)$value, tan(v1_value)) + expect_equal(v_tan_test$get(i-1)$value, tan(v1_value)) } #Test tanh function works. v_tanh_test <- tanh(v1) for(i in 1:v_size){ - expect_equal(v_tanh_test$at(i)$value, tanh(v1_value)) + expect_equal(v_tanh_test$get(i-1)$value, tanh(v1_value)) } #Test exp function works. v_exp_test <- exp(v1) for(i in 1:v_size){ - expect_equal(v_exp_test$at(i)$value, exp(v1_value)) + expect_equal(v_exp_test$get(i-1)$value, exp(v1_value)) } #Test log10 function works. v_log10_test <- log10(v1) for(i in 1:v_size){ - expect_equal(v_log10_test$at(i)$value, log10(v1_value)) + expect_equal(v_log10_test$get(i-1)$value, log10(v1_value)) } #Test sqrt function works. v_sqrt_test <- sqrt(v1) for(i in 1:v_size){ - expect_equal(v_sqrt_test$at(i)$value, sqrt(v1_value)) + expect_equal(v_sqrt_test$get(i-1)$value, sqrt(v1_value)) } #Test log function works. v_log_test <- log(v1) for(i in 1:v_size){ - expect_equal(v_log_test$at(i)$value, log(v1_value)) + expect_equal(v_log_test$get(i-1)$value, log(v1_value)) } clear() diff --git a/vignettes/fims-demo.Rmd b/vignettes/fims-demo.Rmd index 4919a6a8..07ef6fa1 100644 --- a/vignettes/fims-demo.Rmd +++ b/vignettes/fims-demo.Rmd @@ -110,12 +110,12 @@ Each variable of [Parameter class](https://github.com/NOAA-FIMS/FIMS/blob/main/i ```{r fleet_selectivity} methods::show(LogisticSelectivity) fishing_fleet_selectivity <- methods::new(LogisticSelectivity) -fishing_fleet_selectivity$inflection_point$value <- 2.0 -fishing_fleet_selectivity$inflection_point$is_random_effect <- FALSE -fishing_fleet_selectivity$inflection_point$estimated <- TRUE -fishing_fleet_selectivity$slope$value <- 1.0 -fishing_fleet_selectivity$slope$is_random_effect <- FALSE -fishing_fleet_selectivity$slope$estimated <- TRUE +fishing_fleet_selectivity$inflection_point[1]$value <- 2.0 +fishing_fleet_selectivity$inflection_point[1]$is_random_effect <- FALSE +fishing_fleet_selectivity$inflection_point[1]$estimated <- TRUE +fishing_fleet_selectivity$slope[1]$value <- 1.0 +fishing_fleet_selectivity$slope[1]$is_random_effect <- FALSE +fishing_fleet_selectivity$slope[1]$estimated <- TRUE ``` #### Creating the Fleet Object @@ -147,7 +147,7 @@ fishing_fleet$log_Fmort <- new(ParameterVector, log(c( # Turn on estimation for F fishing_fleet$log_Fmort$set_all_estimable(TRUE) # Set value for log_q -fishing_fleet$log_q <- log(1.0) +fishing_fleet$log_q[1]$value <- log(1.0) fishing_fleet$estimate_q <- FALSE fishing_fleet$random_q <- FALSE fishing_fleet$SetSelectivity(fishing_fleet_selectivity$get_id()) @@ -199,12 +199,12 @@ survey_fleet_age_comp$age_comp_data <- survey_agecomp * 200 # unit: number at ag ```{r survey-selectivity} survey_fleet_selectivity <- new(LogisticSelectivity) -survey_fleet_selectivity$inflection_point$value <- 1.5 -survey_fleet_selectivity$inflection_point$is_random_effect <- FALSE -survey_fleet_selectivity$inflection_point$estimated <- TRUE -survey_fleet_selectivity$slope$value <- 2.0 -survey_fleet_selectivity$slope$is_random_effect <- FALSE -survey_fleet_selectivity$slope$estimated <- TRUE +survey_fleet_selectivity$inflection_point[1]$value <- 1.5 +survey_fleet_selectivity$inflection_point[1]$is_random_effect <- FALSE +survey_fleet_selectivity$inflection_point[1]$estimated <- TRUE +survey_fleet_selectivity$slope[1]$value <- 2.0 +survey_fleet_selectivity$slope[1]$is_random_effect <- FALSE +survey_fleet_selectivity$slope[1]$estimated <- TRUE ``` @@ -217,7 +217,7 @@ survey_fleet$nages <- nages survey_fleet$nyears <- nyears # survey_fleet$estimate_F <- FALSE # survey_fleet$random_F <- FALSE -survey_fleet$log_q <- log(3.315143e-07) +survey_fleet$log_q[1]$value <- log(3.315143e-07) survey_fleet$estimate_q <- TRUE survey_fleet$random_q <- FALSE survey_fleet$SetSelectivity(survey_fleet_selectivity$get_id()) @@ -257,12 +257,12 @@ methods::show(BevertonHoltRecruitment) There are three parameters we need to set-up: *log_sigma_recruit*, *log_rzero*, and *logit_steep*. ```{r set-up-recruitment} -recruitment$log_rzero$value <- log(1e+06) # unit: log(number) -recruitment$log_rzero$is_random_effect <- FALSE -recruitment$log_rzero$estimated <- TRUE -recruitment$logit_steep$value <- -log(1.0 - 0.75) + log(0.75 - 0.2) -recruitment$logit_steep$is_random_effect <- FALSE -recruitment$logit_steep$estimated <- FALSE +recruitment$log_rzero[1]$value <- log(1e+06) # unit: log(number) +recruitment$log_rzero[1]$is_random_effect <- FALSE +recruitment$log_rzero[1]$estimated <- TRUE +recruitment$logit_steep[1]$value <- -log(1.0 - 0.75) + log(0.75 - 0.2) +recruitment$logit_steep[1]$is_random_effect <- FALSE +recruitment$logit_steep[1]$estimated <- FALSE recruitment$log_devs <- new(ParameterVector, c( 0.08904850, 0.43787763, -0.13299042, -0.43251973, 0.64861200, 0.50640852, -0.06958319, 0.30246260, @@ -311,12 +311,12 @@ Each population will also need a maturity model. Here we define a logistic matur ```{r maturity} # Maturity maturity <- new(LogisticMaturity) -maturity$inflection_point$value <- 2.25 -maturity$inflection_point$is_random_effect <- FALSE -maturity$inflection_point$estimated <- FALSE -maturity$slope$value <- 3 -maturity$slope$is_random_effect <- FALSE -maturity$slope$estimated <- FALSE +maturity$inflection_point[1]$value <- 2.25 +maturity$inflection_point[1]$is_random_effect <- FALSE +maturity$inflection_point[1]$estimated <- FALSE +maturity$slope[1]$value <- 3 +maturity$slope[1]$is_random_effect <- FALSE +maturity$slope[1]$estimated <- FALSE ``` Now that our life history sub-models are defined, lets define the actual population. diff --git a/vignettes/fims-path-maturity.Rmd b/vignettes/fims-path-maturity.Rmd index 744a6513..61672a56 100644 --- a/vignettes/fims-path-maturity.Rmd +++ b/vignettes/fims-path-maturity.Rmd @@ -57,12 +57,12 @@ library(FIMS) # Create a new maturity model maturity <- new(LogisticMaturity) # Populate the maturity module with parameter values. -maturity$inflection_point$value <- 10 -maturity$inflection_point$is_random_effect <- FALSE -maturity$inflection_point$estimated <- FALSE -maturity$slope$value <- 0.2 -maturity$slope$is_random_effect <- FALSE -maturity$slope$estimated <- FALSE +maturity$inflection_point[1]$value <- 10 +maturity$inflection_point[1]$is_random_effect <- FALSE +maturity$inflection_point[1]$estimated <- FALSE +maturity$slope[1]$value <- 0.2 +maturity$slope[1]$is_random_effect <- FALSE +maturity$slope[1]$estimated <- FALSE ``` After that, we need to create and set up a population and link the maturity module to the population: