diff --git a/CMakeLists.txt b/CMakeLists.txt index 64d98dd30..008dea9ed 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 3c7c0bd84..105b76b7f 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 a3f6e85c5..c03c68071 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 79ddf8ad3..80e4b1602 100644 --- a/inst/include/common/def.hpp +++ b/inst/include/common/def.hpp @@ -16,12 +16,71 @@ #include #include #include +//#include "fims_log.hpp" + + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +#if defined(linux) || defined(__linux) || defined(__linux__) +#define FIMS_LINUX +#elif defined(__FreeBSD__) || defined(__NetBSD__) || defined(__OpenBSD__) || defined(__DragonFly__) +#define FIMS_BSD +#elif defined(sun) || defined(__sun) +#define FIMS_SOLARIS +#elif defined(__sgi) +#define FIMS_IRIX +#elif defined(__hpux) +#define FIMS_HPUX +#elif defined(__CYGWIN__) +#define FIMS_CYGWIN +#elif defined(_WIN32) || defined(__WIN32__) || defined(WIN32) +#define FIMS_WIN32 +#elif defined(_WIN64) || defined(__WIN64__) || defined(WIN64) +#define FIMS_WIN64 +#elif defined(__BEOS__) +#define FIMS_BEOS +#elif defined(macintosh) || defined(__APPLE__) || defined(__APPLE_CC__) +#define FIMS_MACOS +#elif defined(__IBMCPP__) || defined(_AIX) +#define FIMS_AIX +#elif defined(__amigaos__) +#define FIMS_AMIGAOS +#elif defined(__QNXNTO__) +#define FIMS_QNXNTO +#endif + +#if defined(FIMS_WIN32) || defined(FIMS_WIN64) +#define FIMS_WINDOWS +#endif + +#ifdef FIMS_WINDOWS +#include +//#define __PRETTY_FUNCTION__ __FUNCSIG__ +#endif + +#if !defined(__PRETTY_FUNCTION__) && !defined(__GNUC__) +#ifdef FIMS_WINDOWS +#define __PRETTY_FUNCTION__ __FUNCTION__ +#endif +#endif // The following rows initialize default log files for outputing model progress // 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 +109,459 @@ std::ofstream DEBUG_LOG( namespace fims { -/** - * A static class for FIMS logging. - */ + /** + * Log entry. + */ + struct LogEntry { + std::string timestamp; + std::string message; + std::string level; + size_t rank; + std::string user; + std::string wd; + std::string file; + std::string func; + int line; -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()); - } + /** + * Convert this object to a string. + */ + std::string to_string() { + std::stringstream ss; + ss << "\"timestamp\" : " << "\"" << this->timestamp << "\"" << ",\n"; + ss << "\"level\" : " << "\"" << this->level << "\",\n"; + ss << "\"message\" : " << "\"" << this->message << "\",\n"; + ss << "\"id\" : " << "\"" << this->rank << "\",\n"; + ss << "\"user\" : " << "\"" << this->user << "\",\n"; + ss << "\"wd\" : " << "\"" << this->wd << "\",\n"; + ss << "\"file\" : " << "\"" << this->file << "\",\n"; + ss << "\"routine\" : " << "\"" << this->func << "\",\n"; + ss << "\"line\" : " << "\"" << this->line << "\"\n"; + return ss.str(); + } + + }; + + /** + * FIMS logging class. + */ + class FIMSLog { + std::vector entries; + std::vector log_entries; + size_t entry_number = 0; + std::string path = "fims.log"; + size_t warning_count = 0; + size_t error_count = 0; + + /** + * Get username. + * + * @return username. + */ + std::string get_user() { + char * user; + std::string user_ret = "UNKOWN_USER"; + +#ifdef FIMS_WINDOWS + user = getenv("username"); + user_ret = std::string(user); +#endif +#ifdef FIMS_LINUX + user = getenv("USER"); + user_ret = std::string(user); +#endif + +#ifdef FIMS_MACOS + user = getenv("USER"); + user_ret = std::string(user); +#endif + + return user_ret; + } + public: + bool write_on_exit = true; + bool throw_on_error = false; + static std::shared_ptr fims_log; + + /** + * Default constructor. + */ + FIMSLog() { + + } + + /** + * Destructor. If write_on_exit is set to true, + * the log will be written to the disk in JSON format. + */ + ~FIMSLog() { + if (this->write_on_exit) { + std::ofstream log(this->path); + log << this->get_log(); + log.close(); + } + } + + /** + * Set a path for the log file. + * + * @param path + */ + void set_path(std::string path) { + this->path = path; + } + + /** + * Get the path for the log file. + * + * @return + */ + std::string get_path() { + return this->path; + } + + /** + * Add a "info" level message to the log. + * + * @param str + * @param line + * @param file + * @param func + */ + void info_message(std::string str, int line, const char* file, const char* func) { + + std::filesystem::path cwd = std::filesystem::current_path(); + std::stringstream ss; + auto now = std::chrono::system_clock::now(); + std::time_t now_time = std::chrono::system_clock::to_time_t(now); + std::string ctime_no_newline = strtok(ctime(&now_time), "\n"); + + LogEntry l; + l.timestamp = ctime_no_newline; + l.message = str; + l.level = "info"; + l.rank = this->log_entries.size(); + l.user = this->get_user(); + l.wd = cwd.string(); + l.file = file; + l.line = line; + l.func = func; + this->log_entries.push_back(l); + + } + + /** + * Add a "debug" level message to the log. + * + * @param str + * @param line + * @param file + * @param func + */ + void debug_message(std::string str, int line, const char* file, const char* func) { + + std::filesystem::path cwd = std::filesystem::current_path(); + std::stringstream ss; + auto now = std::chrono::system_clock::now(); + std::time_t now_time = std::chrono::system_clock::to_time_t(now); + std::string ctime_no_newline = strtok(ctime(&now_time), "\n"); + + LogEntry l; + l.timestamp = ctime_no_newline; + l.message = str; + l.level = "debug"; + l.rank = this->log_entries.size(); + l.user = this->get_user(); + l.wd = cwd.string(); + l.file = file; + l.line = line; + l.func = func; + this->log_entries.push_back(l); + + } + + /** + * Add a "error" level message to the log. + * + * @param str + * @param line + * @param file + * @param func + */ + void error_message(std::string str, int line, const char* file, const char* func) { + this->error_count++; + std::filesystem::path cwd = std::filesystem::current_path(); + + std::stringstream ss; + auto now = std::chrono::system_clock::now(); + std::time_t now_time = std::chrono::system_clock::to_time_t(now); + std::string ctime_no_newline = strtok(ctime(&now_time), "\n"); + + LogEntry l; + l.timestamp = ctime_no_newline; + l.message = str; + l.level = "error"; + l.rank = this->log_entries.size(); + l.user = this->get_user(); + l.wd = cwd.string(); + l.file = file; + l.line = line; + l.func = func; + this->log_entries.push_back(l); + + if (this->throw_on_error) { + std::stringstream ss; + ss << "\n\n" << l.to_string() << "\n\n"; + throw std::runtime_error(ss.str().c_str()); + } + + } + + /** + * Add a "warning" level message to the log. + * + * @param str + * @param line + * @param file + * @param func + */ + void warning_message(std::string str, int line, const char* file, const char* func) { + this->warning_count++; + std::filesystem::path cwd = std::filesystem::current_path(); + + std::stringstream ss; + auto now = std::chrono::system_clock::now(); + std::time_t now_time = std::chrono::system_clock::to_time_t(now); + std::string ctime_no_newline = strtok(ctime(&now_time), "\n"); + + LogEntry l; + l.timestamp = ctime_no_newline; + l.message = str; + l.level = "warning"; + l.rank = this->log_entries.size(); + l.user = this->get_user(); + l.wd = cwd.string(); + l.file = file; + l.line = line; + l.func = func; + this->log_entries.push_back(l); + + } + + /** + * Get the log as a string. + * + * @return + */ + std::string get_log() { + std::stringstream ss; + if (log_entries.size() == 0) { + ss << "[\n]"; + } else { + ss << "[\n"; + for (size_t i = 0; i < log_entries.size() - 1; i++) { + ss << "{\n" << this->log_entries[i].to_string() << "},\n"; + + } + ss << "{\n" << this->log_entries[log_entries.size() - 1].to_string() << "}\n]"; + } + return ss.str(); + } + + /** + * Return only error entries from the log. + * + * @return + */ + std::string get_errors() { + std::stringstream ss; + std::vector errors; + for (size_t i = 0; i < log_entries.size(); i++) { + if (log_entries[i].level == "error") { + errors.push_back(this->log_entries[i]); + } + } + + if (errors.size() == 0) { + ss << "[\n]"; + } else { + ss << "[\n"; + for (size_t i = 0; i < errors.size() - 1; i++) { - return fims_log::FIMS_LOGS[l]; - } -}; + ss << "{\n" << errors[i].to_string() << "},\n"; -std::map fims_log::FIMS_LOGS; + } + + ss << "{\n" << errors[errors.size() - 1].to_string() << "}\n]"; + + } + return ss.str(); + } + + /** + * Return only warning entries from the log. + * + * @return + */ + std::string get_warnings() { + std::stringstream ss; + std::vector warnings; + for (size_t i = 0; i < log_entries.size(); i++) { + if (log_entries[i].level == "warning") { + warnings.push_back(this->log_entries[i]); + } + } + + if (warnings.size() == 0) { + ss << "[\n]"; + } else { + ss << "[\n"; + for (size_t i = 0; i < warnings.size() - 1; i++) { + + ss << "{\n" << warnings[i].to_string() << "},\n"; + + } + + ss << "{\n" << warnings[warnings.size() - 1].to_string() << "}\n]"; + + } + return ss.str(); + } + + /** + * Return only info entries from the log. + * + * @return + */ + std::string get_info() { + std::stringstream ss; + std::vector info; + for (size_t i = 0; i < log_entries.size(); i++) { + if (log_entries[i].level == "info") { + info.push_back(this->log_entries[i]); + } + } + + if (info.size() == 0) { + ss << "[\n]"; + } else { + ss << "[\n"; + for (size_t i = 0; i < info.size() - 1; i++) { + + ss << "{\n" << info[i].to_string() << "},\n"; + + } + + ss << "{\n" << info[info.size() - 1].to_string() << "}\n]"; + + } + return ss.str(); + } + + /** + * Query the log by module. + * + * @param module + * @return + */ + std::string get_module(const std::string& module) { + std::stringstream ss; + std::vector info; + for (size_t i = 0; i < log_entries.size(); i++) { + if (log_entries[i].file.find(module) != std::string::npos) { + info.push_back(this->log_entries[i]); + } + } + + if (info.size() == 0) { + ss << "[\n]"; + } else { + ss << "[\n"; + for (size_t i = 0; i < info.size() - 1; i++) { + + ss << "{\n" << info[i].to_string() << "},\n"; + + } + + ss << "{\n" << info[info.size() - 1].to_string() << "}\n]"; + + } + return ss.str(); + } + + size_t get_error_count() const { + return error_count; + } + + size_t get_warning_count() const { + return warning_count; + } + + + + }; + + + std::shared_ptr FIMSLog::fims_log = std::make_shared(); + +} // namespace fims + + + + + +#ifdef FIMS_DEBUG + +#define FIMS_DEBUG_LOG(MESSAGE) FIMSLog::fims_log->debug_message(MESSAGE, __LINE__, __FILE__, __PRETTY_FUNCTION__); + +#else + +#define FIMS_DEBUG_LOG(MESSAGE) + +#endif + +#define FIMS_INFO_LOG(MESSAGE) fims::FIMSLog::fims_log->info_message(MESSAGE, __LINE__, __FILE__, __PRETTY_FUNCTION__); + +#define FIMS_WARNING_LOG(MESSAGE) fims::FIMSLog::fims_log->warning_message(MESSAGE, __LINE__, __FILE__, __PRETTY_FUNCTION__); + +#define FIMS_ERROR_LOG(MESSAGE) fims::FIMSLog::fims_log->error_message(MESSAGE, __LINE__, __FILE__, __PRETTY_FUNCTION__); + +#define FIMS_STR(s) #s + + +namespace fims { + + /** + * Signal intercept function. Writes the log to the disk before + * a crash occurs. + * + * @param sig + */ + void WriteAtExit(int sig) { + + + if (FIMSLog::fims_log->write_on_exit) { + std::ofstream log(FIMSLog::fims_log->get_path()); + log << FIMSLog::fims_log->get_log(); + log.close(); + } + std::signal(sig, SIG_DFL); + raise(sig); + } + + /** + * Converts an object T to a string. + * + * @param v + * @return + */ + template + std::string to_string(T v) { + std::stringstream ss; + ss << v; + return ss.str(); + } -} // namespace fims +} #endif /* TRAITS_HPP */ diff --git a/inst/include/common/fims_math.hpp b/inst/include/common/fims_math.hpp index dcdf42d06..6be340cf6 100644 --- a/inst/include/common/fims_math.hpp +++ b/inst/include/common/fims_math.hpp @@ -18,8 +18,11 @@ // preprocessing macros //#include "def.hpp" #include +#include +#include #include "../interface/interface.hpp" +#include "fims_vector.hpp" namespace fims_math { #ifdef STD_LIB @@ -45,6 +48,21 @@ template inline const Type log(const Type &x) { return std::log(x); } + +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 @@ -87,6 +105,36 @@ 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 /** @@ -220,6 +268,420 @@ 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; + } + + + } // namespace fims_math #endif /* FIMS_MATH_HPP */ diff --git a/inst/include/common/fims_vector.hpp b/inst/include/common/fims_vector.hpp index 0892846d4..3d619ff23 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,19 @@ 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); } } else { @@ -293,39 +298,40 @@ 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); } // 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"); } } else { valid_model = false; ERROR_LOG << "Error: No data input for density " << d->id << std::endl; - exit(1); } } // end set data @@ -335,6 +341,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 +356,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 +367,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 +383,6 @@ class Information { "population " << p->id << ", recruitment function " << recruitment_uint << std::endl; - exit(1); } } else { @@ -384,7 +392,6 @@ class Information { << ". FIMS requires recruitment functions be defined for all " "populations." << std::endl; - exit(1); } INFO_LOG << "Checking for available growth function." << std::endl; @@ -393,6 +400,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 +412,10 @@ 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); } } else { @@ -417,12 +425,12 @@ class Information { << ". FIMS requires growth functions be defined for all " "populations." << std::endl; - 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 +444,8 @@ class Information { ERROR_LOG << "Error: Expected maturity function not defined for population " << p->id << ", maturity function " << maturity_uint << std::endl; - exit(1); } - + FIMS_INFO_LOG("segment"); } else { valid_model = false; ERROR_LOG << "Error: No maturity function defined for population " @@ -446,13 +453,13 @@ class Information { << ". FIMS requires maturity functions be defined for all " "populations." << std::endl; - 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/multinomial_lpmf.hpp b/inst/include/distributions/functors/multinomial_lpmf.hpp index 87b5035b8..3927d040f 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]); 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] = dmultinom((vector)x_vector, (vector) prob_vector, true); lpdf += this->lpdf_vec[i]; /* if (this->simulate_flag) diff --git a/inst/include/interface/rcpp/rcpp_interface.hpp b/inst/include/interface/rcpp/rcpp_interface.hpp index e84805245..9ae3addc4 100644 --- a/inst/include/interface/rcpp/rcpp_interface.hpp +++ b/inst/include/interface/rcpp/rcpp_interface.hpp @@ -21,10 +21,20 @@ #include "rcpp_objects/rcpp_recruitment.hpp" #include "rcpp_objects/rcpp_selectivity.hpp" #include "rcpp_objects/rcpp_tmb_distribution.hpp" + + +SEXP FIMS_objective_function; +SEXP FIMS_gradient_function; +double FIMS_function_value = 0; +Rcpp::NumericVector FIMS_function_gradient; +double FIMS_mgc_value = 0; +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(); @@ -53,6 +63,103 @@ bool CreateTMBModel() { return true; } +void SetFIMSFunctions(SEXP fn, SEXP gr) { + FIMS_objective_function = fn; + FIMS_gradient_function = gr; +} + +/** + * @brief Extracts derived quantities from model objects. + */ +void Finalize(Rcpp::NumericVector p) { + FIMS_finalized = true; + std::shared_ptr> information = + fims_info::Information::GetInstance(); + + std::shared_ptr> model = + fims_model::Model::GetInstance(); + + for (size_t i = 0; i < information->fixed_effects_parameters.size(); i++) { + *information->fixed_effects_parameters[i] = p[i]; + } + + model->Evaluate(); + + Rcpp::Function f = Rcpp::as(FIMS_objective_function); + Rcpp::Function g = Rcpp::as(FIMS_gradient_function); + double ret = Rcpp::as(f(p)); + Rcpp::NumericVector grad = Rcpp::as(g(p)); + + FIMS_function_value = ret; + FIMS_function_gradient = grad; + std::cout << "Final value = " << FIMS_function_value << "\nGradient: \n"; + double maxgc = -9999; + for (size_t i = 0; i < FIMS_function_gradient.size(); i++) { + if (std::fabs(FIMS_function_gradient[i]) > maxgc) { + maxgc = std::fabs(FIMS_function_gradient[i]); + } + } + FIMS_mgc_value = maxgc; + // std::cout<<"mgc = "<::population_iterator pit; + // for (pit = information->populations.begin(); pit != information->populations.end(); ++pit) { + // pit->second->Prepare(); + // pit->second->Evaluate(); + // } + // + // fims_info::Information < double>::fleet_iterator fit; + // for (fit = information->fleets.begin(); fit != information->fleets.end(); ++fit) { + // fit->second->Prepare(); + // fit->second->Evaluate(); + //// double ac = fit->second->evaluate_age_comp_nll(); + //// double i = fit->second->evaluate_index_nll(); + // } + + + for (size_t i = 0; i < FIMSRcppInterfaceBase::fims_interface_objects.size(); + i++) { + FIMSRcppInterfaceBase::fims_interface_objects[i]->finalize(); + } +} + +/** + * @brief Extracts derived quantities from model objects. + */ +std::string ToJSON() { + auto now = std::chrono::system_clock::now(); + std::time_t now_time = std::chrono::system_clock::to_time_t(now); + std::string ctime_no_newline = strtok(ctime(&now_time), "\n"); + std::shared_ptr> info = + fims_info::Information::GetInstance(); + std::stringstream ss; + ss << "{\n"; + ss << "\"timestamp\": \"" << ctime_no_newline << "\",\n"; + ss << "\"nyears\":" << info->nyears << ",\n"; + ss << "\"nseasons\":" << info->nseasons << ",\n"; + ss << "\"nages\":" << info->nages << ",\n"; + ss << "\"finalized\":" << FIMS_finalized << ",\n"; + ss << "\"objective_function_value\": " << FIMS_function_value << ",\n"; + ss << "\"max_gradient_component\": " << FIMS_mgc_value << ",\n"; + ss << "\"final_gradient\": ["; + 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 << "],"; + } + 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"; + } + + ss << FIMSRcppInterfaceBase::fims_interface_objects[length - 1]->to_json() << "\n}"; + return ss.str(); +} + Rcpp::NumericVector get_fixed_parameters_vector() { // base model std::shared_ptr> d0 = @@ -114,10 +221,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); } /** @@ -339,12 +446,139 @@ void clear() { clear_internal(); clear_internal(); clear_internal(); + + FIMS_finalized = false; +} + +/** + * Returns the entire log as a string in JSON format. + */ +std::string get_log() { + return fims::FIMSLog::fims_log->get_log(); +} + +/** + * Returns only error entries from log as a string in JSON format. + */ +std::string get_log_errors() { + return fims::FIMSLog::fims_log->get_errors(); +} + +/** + * Returns only warning entries from log as a string in JSON format. + */ +std::string get_log_warnings() { + return fims::FIMSLog::fims_log->get_warnings(); +} + +/** + * Returns only info entries from log as a string in JSON format. + */ +std::string get_log_info() { + return fims::FIMSLog::fims_log->get_info(); +} + +/** + * Returns log entries by module as a string in JSON format. + */ +std::string get_log_module(const std::string& module) { + return fims::FIMSLog::fims_log->get_module(module); +} + +/** + * If true, writes the log on exit . + */ +void write_log(bool write) { + FIMS_INFO_LOG("Setting FIMS write log: " + fims::to_string(write)); + fims::FIMSLog::fims_log->write_on_exit = write; +} + +/** + * Sets the path for the log file to written. + */ +void set_log_path(const std::string& path) { + FIMS_INFO_LOG("Setting FIMS log path: " + path); + fims::FIMSLog::fims_log->set_path(path); +} + +/** + * If true, throws a runtime exception when an error is logged . + */ +void set_log_throw_on_error(bool throw_on_error) { + fims::FIMSLog::fims_log->throw_on_error = throw_on_error; +} + +/** + * Initializes the logging system, sets all signal handling. + */ +void init_logging() { + + FIMS_INFO_LOG("Initializing FIMS logging system."); + std::signal(SIGSEGV, &fims::WriteAtExit); + std::signal(SIGINT, &fims::WriteAtExit); + std::signal(SIGABRT, &fims::WriteAtExit); + std::signal(SIGFPE, &fims::WriteAtExit); + std::signal(SIGILL, &fims::WriteAtExit); + std::signal(SIGTERM, &fims::WriteAtExit); +} + +/** + * Add log info entry from R. + */ +void log_info(std::string log_entry) { + fims::FIMSLog::fims_log->info_message(log_entry, -1, "R_env", "R_script_entry"); +} + +/** + * Add log warning entry from R. + */ +void log_warning(std::string log_entry) { + fims::FIMSLog::fims_log->warning_message(log_entry, -1, "R_env", "R_script_entry"); +} + +/** + * Add log error entry from R. + */ +void log_error(std::string log_entry) { + + std::stringstream ss; + ss << "capture.output(traceback(4))"; + SEXP expression, result; + ParseStatus status; + + PROTECT(expression = R_ParseVector(Rf_mkString(ss.str().c_str()), 1, &status, R_NilValue)); + if (status != PARSE_OK) { + std::cout << "Error parsing expression" << std::endl; + UNPROTECT(1); + } + Rcpp::Rcout << "before call."; + PROTECT(result = Rf_eval(VECTOR_ELT(expression, 0), R_GlobalEnv)); + Rcpp::Rcout << "after call."; + UNPROTECT(2); + std::stringstream ss_ret; + ss_ret << "traceback:\n"; + for (int j = 0; j < LENGTH(result); j++) { + std::string str(CHAR(STRING_ELT(result, j))); + ss_ret << str << "\n"; + } + + std::string ret = ss_ret.str(); //"find error";//Rcpp::as(result); + + + // Rcpp::Environment base = Rcpp::Environment::global_env(); + // Rcpp::Function f = base["traceback"]; + // std::string ret = Rcpp::as(f()); + fims::FIMSLog::fims_log->error_message(log_entry, -1, "R_env", ret.c_str()); } RCPP_EXPOSED_CLASS(Parameter) RCPP_EXPOSED_CLASS(ParameterVector) + RCPP_MODULE(fims) { Rcpp::function("CreateTMBModel", &CreateTMBModel); + Rcpp::function("SetFIMSFunctions", &SetFIMSFunctions); + Rcpp::function("Finalize", &Finalize); + Rcpp::function("ToJSON", &ToJSON); Rcpp::function("get_fixed", &get_fixed_parameters_vector); Rcpp::function("get_random", &get_random_parameters_vector); Rcpp::function("get_parameter_names", &get_parameter_names); @@ -362,12 +596,24 @@ RCPP_MODULE(fims) { Rcpp::function("clear_maturity_log", clear_maturity_log); Rcpp::function("clear_selectivity_log", clear_selectivity_log); Rcpp::function("clear_debug_log", clear_debug_log); - + Rcpp::function("get_log", get_log); + Rcpp::function("get_log_errors", get_log_errors); + Rcpp::function("get_log_warnings", get_log_warnings); + Rcpp::function("get_log_info", get_log_info); + Rcpp::function("get_log_module", get_log_module); + Rcpp::function("write_log", write_log); + Rcpp::function("set_log_path", set_log_path); + Rcpp::function("init_logging", init_logging); + Rcpp::function("set_log_throw_on_error", set_log_throw_on_error); + Rcpp::function("log_info", log_info); + Rcpp::function("log_warning", log_warning); + Rcpp::function("log_error", log_error); Rcpp::class_("Parameter", "FIMS Parameter Class") .constructor() .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") @@ -378,7 +624,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_data.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_data.hpp index 3f25e3f7e..fa870ff61 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_data.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_data.hpp @@ -80,6 +80,26 @@ class AgeCompDataInterface : public DataInterfaceBase { /** @brief get the ID of the interface base object */ virtual uint32_t get_id() { return this->id; } + + + virtual std::string to_json() { + std::stringstream ss; + + ss << "\"module\" : {\n"; + ss << " \"name\": \"data\",\n"; + ss << " \"type\" : \"AgeComp\",\n"; + ss << " \"id\":" << this->id << ",\n"; + ss << " \"rank\": " << 2 << ",\n"; + ss << " \"dimensions\": [" << this->ymax << "," << this->amax << "],\n"; + ss << " \"values\": ["; + for (size_t i = 0; i < age_comp_data.size() - 1; i++) { + ss << age_comp_data[i] << ", "; + } + ss << age_comp_data[age_comp_data.size() - 1] << "]\n"; + ss << "}"; + return ss.str(); + } + #ifdef TMB_MODEL @@ -143,6 +163,25 @@ class IndexDataInterface : public DataInterfaceBase { /** @brief get the ID of the interface base object */ virtual uint32_t get_id() { return this->id; } + + + virtual std::string to_json() { + std::stringstream ss; + + ss << "\"module\" : {\n"; + ss << " \"name\": \"data\",\n"; + ss << " \"type\": \"Index\",\n"; + ss << " \"id\": " << this->id << ",\n"; + ss << " \"rank\": " << 1 << ",\n"; + ss << " \"dimensions\": [" << this->ymax << "],\n"; + ss << " \"values\": ["; + for (size_t i = 0; i < index_data.size() - 1; i++) { + ss << index_data[i] << ", "; + } + ss << index_data[index_data.size() - 1] << "]\n"; + ss << "}"; + return ss.str(); + } #ifdef TMB_MODEL diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp index b5dc3600f..17b29c324 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_fleet.hpp @@ -53,11 +53,12 @@ std::map FleetInterfaceBase::live_objects; class FleetInterface : public FleetInterfaceBase { int interface_selectivity_id_m = -999; /**< id of selectivity component*/ - public: +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 */ - double log_q; /**< log of catchability for the fleet*/ + 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 */ @@ -65,6 +66,12 @@ class FleetInterface : public FleetInterfaceBase { 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 */ + FleetInterface() : FleetInterfaceBase() {} virtual ~FleetInterface() {} @@ -81,7 +88,156 @@ class FleetInterface : public FleetInterfaceBase { 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]; + } + + } + + } + + 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"; + + + 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 { + 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(); + + } + + + #ifdef TMB_MODEL + template bool add_to_fims_tmb_internal() { std::shared_ptr > info = @@ -97,16 +253,21 @@ class FleetInterface : public FleetInterfaceBase { fleet->nyears = this->nyears; fleet->fleet_selectivity_id_m = interface_selectivity_id_m; - fleet->log_q = this->log_q; - if (this->estimate_q) { + 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->random_q) { - info->RegisterRandomEffect(fleet->log_q); - } else { - info->RegisterParameter(fleet->log_q); + if (this->log_q[i].is_random_effect_m) { + info->RegisterRandomEffect(fleet->log_q[i]); + } else { + info->RegisterParameter(fleet->log_q[i]); + } } } + 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; @@ -138,6 +299,7 @@ class FleetInterface : public FleetInterfaceBase { /** @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(); diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_growth.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_growth.hpp index 8b60a2c4e..ea98e36af 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_growth.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_growth.hpp @@ -108,6 +108,32 @@ vectors have been set */ EWAAGrowth.ewaa = this->ewaa; return EWAAGrowth.evaluate(age); } + + + virtual std::string to_json() { + std::stringstream ss; + ss << "\"module\" : {\n"; + ss << " \"name\": \"growth\",\n"; + ss << " \"type\" : \"EWAA\",\n"; + ss << " \"id\":" << this->id << ",\n"; + ss << " \"rank\": " << 1 << ",\n"; + ss << " \"dimensions\": [" << this->weights.size() << "],\n"; + + ss << " \"ages\": ["; + for (size_t i = 0; i < ages.size() - 1; i++) { + ss << ages[i] << ", "; + } + ss << ages[ages.size() - 1] << "],\n"; + + ss << " \"values\": ["; + for (size_t i = 0; i < weights.size() - 1; i++) { + ss << weights[i] << ", "; + } + ss << weights[weights.size() - 1] << "]\n"; + ss << "}"; + return ss.str(); + } + #ifdef TMB_MODEL template 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 01b256246..90220cf9e 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 @@ -27,7 +28,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 = @@ -47,14 +49,39 @@ class Parameter { Parameter(double value, double min, double max, bool estimated) : id_m(Parameter::id_g++), value_m(value), min_m(min), max_m(max), estimated_m(estimated) {} - /** - * @brief Constructor for initializing Parameter. - * @details Inputs include value. - */ - Parameter(double value) { - value_m = value; - id_m = Parameter::id_g++; - } + 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. + */ + Parameter(double value) { + value_m = value; + id_m = Parameter::id_g++; + } /** * @brief Constructor for initializing Parameter. @@ -67,37 +94,50 @@ class Parameter { uint32_t Parameter::id_g = 0; +std::ostream& operator<<(std::ostream& out, const Parameter& p) { + 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; +} /** * @brief Rcpp representation of a Parameter vector * interface between R and cpp. */ 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*/ - /** * @brief default constructor */ 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++; - for(size_t i =0; i < size; i++){ - Parameter p; - this->storage_m.push_back(Rcpp::wrap(p)); + this->storage_m = std::make_shared >(); + this->storage_m->resize(size); + for (size_t i = 0; i < size; i++) { + storage_m->at(i) = Parameter(); } } + /** * @brief vector constructor * @param x numeric vector @@ -105,10 +145,21 @@ class ParameterVector{ */ ParameterVector(Rcpp::NumericVector x, size_t size){ this->id_m = ParameterVector::id_g++; - for(size_t i =0; i < size; i++){ - Parameter p = x[i]; - this->storage_m.push_back(Rcpp::wrap(p)); + this->storage_m = std::make_shared >(); + this->resize(size); + for (size_t i = 0; i < size; i++) { + storage_m->at(i).value_m = x[i]; + } + } + + 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++) { + storage_m->at(i).value_m = v[i]; } + } /** @@ -120,51 +171,55 @@ 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); } /** * @brief Accessor. First index is one. For calling from R. * @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(); + size_t size() { + return this->storage_m->size(); } /** * @brief resize to length "size" * @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; - } - + void resize(size_t size) { + this->storage_m->resize(size); } /** @@ -173,10 +228,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; } } @@ -186,10 +240,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; } } @@ -199,10 +251,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; } } @@ -212,10 +262,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; } } @@ -225,24 +273,42 @@ 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 << "["; + size_t size = v.size(); + for (size_t i = 0; i < size - 1; i++) { + out << v[i] << ", "; + } + out << v[size - 1] << "]"; + return out; +} /** *@brief Base class for all interface objects */ class FIMSRcppInterfaceBase { - public: - /**< FIMS interface object vectors */ - static std::vector fims_interface_objects; +public: + + bool finalized = false; + + /**< FIMS interface object vectors */ + static std::vector fims_interface_objects; /** @brief virtual method to inherit to add objects to the TMB model */ virtual bool add_to_fims_tmb() { @@ -250,6 +316,13 @@ class FIMSRcppInterfaceBase { "implemented.\n"; return false; } + virtual void finalize() { + + } + + virtual std::string to_json() { + return ""; + } }; std::vector FIMSRcppInterfaceBase::fims_interface_objects; diff --git a/inst/include/interface/rcpp/rcpp_objects/rcpp_maturity.hpp b/inst/include/interface/rcpp/rcpp_objects/rcpp_maturity.hpp index a845f33b2..fcf606eb1 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_maturity.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_maturity.hpp @@ -59,12 +59,13 @@ std::map MaturityInterfaceBase::live_objects; * instantiate from R: logistic_maturity <- new(logistic_maturity) */ class LogisticMaturityInterface : public MaturityInterfaceBase { - public: - Parameter +public: + ParameterVector inflection_point; /**< the index value at which the response reaches .5 */ - Parameter slope; /**< the width of the curve at the inflection_point */ + ParameterVector slope; /**< the width of the curve at the inflection_point */ - LogisticMaturityInterface() : MaturityInterfaceBase() {} + LogisticMaturityInterface() : MaturityInterfaceBase() { + } virtual ~LogisticMaturityInterface() {} @@ -78,12 +79,85 @@ 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 << " \"id\": " << this->id << ",\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"inflection_point\",\n"; + ss << " \"id\":" << this->inflection_point.id_m << ",\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\": \"vector\",\n"; + ss << " \"values\":" << this->slope << ",\n"; + + + ss << "}"; + + return ss.str(); + } + #ifdef TMB_MODEL template @@ -96,26 +170,36 @@ 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]); - } - } - 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]); - } - } + 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(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]); + } + } + } // add to Information info->maturity_models[maturity->id] = maturity; @@ -126,6 +210,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 07885098d..a6f3a0c26 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_population.hpp @@ -21,7 +21,7 @@ */ class PopulationInterfaceBase : public FIMSRcppInterfaceBase { public: - static uint32_t id_g; /**< static id of the population interface base*/ + static uint32_t id_g; /**< static id of the population interface base*/ uint32_t id; /**< id of the population interface base */ // live objects in C++ are objects that have been created and live in memory static std::map @@ -68,6 +68,18 @@ class PopulationInterface : public PopulationInterfaceBase { of female individuals */ bool estimate_prop_female; /**finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("Population " + 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(); + + + + + this->estimated_log_M = Rcpp::NumericVector(this->log_M.size()); + for (size_t i = 0; i < this->log_M.size(); i++) { + this->estimated_log_M[i] = this->log_M[i].value_m; + } + + this->estimated_log_init_naa = Rcpp::NumericVector(this->log_init_naa.size()); + for (size_t i = 0; i < this->log_init_naa.size(); i++) { + this->estimated_log_init_naa[i] = this->log_init_naa[i].value_m; + } + + this->estimated_proportion_female = Rcpp::NumericVector(this->proportion_female.size()); + for (size_t i = 0; i < this->proportion_female.size(); i++) { + this->estimated_proportion_female[i] = this->proportion_female[i]; + } + + fims_info::Information::population_iterator it; + + + it = info->populations.find(this->id); + + std::shared_ptr > pop = + info->populations[this->id]; + it = info->populations.find(this->id); + if (it == info->populations.end()) { + FIMS_WARNING_LOG("Population " + fims::to_string(this->id) + " not found in Information."); + return; + } else { + + if (this->estimated_log_M) { + for (size_t i = 0; i < this->log_M.size(); i++) { + this->estimated_log_M[i] = pop->log_M[i]; + } + } + + if (this->estimated_log_init_naa) { + for (size_t i = 0; i < this->log_init_naa.size(); i++) { + this->estimated_log_init_naa[i] = pop->log_init_naa[i]; + } + } + + if (this->estimate_prop_female) { + for (size_t i = 0; i < this->proportion_female.size(); i++) { + this->estimated_proportion_female[i] = pop->proportion_female[i]; + } + } + + this->derived_naa = Rcpp::NumericVector(pop->numbers_at_age.size()); + this->derived_ssb = Rcpp::NumericVector(pop->spawning_biomass.size()); + this->derived_biomass = Rcpp::NumericVector(pop->biomass.size()); + this->derived_recruitment = Rcpp::NumericVector(pop->expected_recruitment.size()); + + //set naa from Information/ + for (size_t i = 0; i < this->derived_naa.size(); i++) { + this->derived_naa[i] = pop->numbers_at_age[i]; + } + + //set ssb from Information/ + for (size_t i = 0; i < this->derived_ssb.size(); i++) { + this->derived_ssb[i] = pop->spawning_biomass[i]; + } + + //set biomass from Information + for (size_t i = 0; i < this->derived_biomass.size(); i++) { + this->derived_biomass[i] = pop->biomass[i]; + } + + //set recruitment from Information/ + for (size_t i = 0; i < this->derived_recruitment.size(); i++) { + this->derived_recruitment[i] = pop->expected_recruitment[i]; + } + + } + + } + + virtual std::string to_json() { + std::stringstream ss; + + ss << "\"module\" : {\n"; + ss << " \"name\" : \"Population\",\n"; + + ss << " \"type\" : \"population\",\n"; + ss << " \"tag\" : \"" << this->name << "\",\n"; + ss << " \"id\": " << this->id << ",\n"; + ss << " \"recruitment_id\": " << this->recruitment_id << ",\n"; + ss << " \"growth_id\": " << this->growth_id << ",\n"; + ss << " \"maturity_id\": " << this->maturity_id << ",\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"log_M\",\n"; + ss << " \"id\":" << -999 << ",\n"; + ss << " \"type\": \"vector\",\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\":" << this->log_init_naa << " \n},\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"proportion_female\",\n"; + ss << " \"id\":" << -999 << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->proportion_female << "\n},\n"; + + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"ssb\",\n"; + ss << " \"values\":["; + if (this->derived_ssb.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->derived_ssb.size() - 1; i++) { + ss << this->derived_ssb[i] << ", "; + } + ss << this->derived_ssb[this->derived_ssb.size() - 1] << "]\n"; + } + ss << " },\n"; + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"naa\",\n"; + ss << " \"values\":["; + if (this->derived_naa.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->derived_naa.size() - 1; i++) { + ss << this->derived_naa[i] << ", "; + } + ss << this->derived_naa[this->derived_naa.size() - 1] << "]\n"; + } + ss << " },\n"; + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"biomass\",\n"; + ss << " \"values\":["; + if (this->derived_biomass.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->derived_biomass.size() - 1; i++) { + ss << this->derived_biomass[i] << ", "; + } + ss << this->derived_biomass[this->derived_biomass.size() - 1] << "]\n"; + } + ss << " },\n"; + + ss << " \"derived_quantity\": {\n"; + ss << " \"name\": \"recruitment\",\n"; + ss << " \"values\":["; + if (this->derived_recruitment.size() == 0) { + ss << "]\n"; + } else { + for (size_t i = 0; i < this->derived_recruitment.size() - 1; i++) { + ss << this->derived_recruitment[i] << ", "; + } + ss << this->derived_recruitment[this->derived_recruitment.size() - 1] << "]\n"; + } + ss << " }\n"; + + ss << "}"; + + return ss.str(); + } + + #ifdef TMB_MODEL template @@ -133,8 +329,8 @@ class PopulationInterface : public PopulationInterfaceBase { for (size_t i = 0; i < log_M.size(); i++) { population->log_M[i] = this->log_M[i].value_m; if (this->log_M[i].estimated_m) { - info->RegisterParameterName("log_M"); - info->RegisterParameter(population->log_M[i]); + info->RegisterParameterName("log_M"); + info->RegisterParameter(population->log_M[i]); } } info->variable_map[this->log_M.id_m] = &(population)->log_M; @@ -169,6 +365,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 37506b4d7..84e48fbb1 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_recruitment.hpp @@ -64,11 +64,16 @@ std::map */ class BevertonHoltRecruitmentInterface : public RecruitmentInterfaceBase { public: - Parameter logit_steep; /**< steepness or the productivity of the stock*/ - Parameter log_rzero; /**< recruitment at unfished biomass */ + ParameterVector logit_steep; /**< steepness or the productivity of the stock*/ + ParameterVector log_rzero; /**< recruitment at unfished biomass */ ParameterVector log_devs; /**< log recruitment deviations*/ bool estimate_log_devs = false; /**< boolean describing whether to estimate */ + double estimated_logit_steep; /**< estimated steepness or the productivity of the stock*/ + double estimated_log_rzero; /**< estimated recruitment at unfished biomass */ + double estimated_log_sigma_recruit; /**< estimated log of the stock recruit standard deviation */ + Rcpp::NumericVector estimated_log_devs; /**< estimated log recruitment deviations*/ + BevertonHoltRecruitmentInterface() : RecruitmentInterfaceBase() {} virtual ~BevertonHoltRecruitmentInterface() {} @@ -78,18 +83,115 @@ 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); } + virtual void finalize() { + + if (this->finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("Beverton-Holt Recruitment " + 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::recruitment_models_iterator it; + + it = info->recruitment_models.find(this->id); + + if (it == info->recruitment_models.end()) { + FIMS_WARNING_LOG("Beverton-Holt Recruitment " + fims::to_string(this->id) + " not found in Information."); + return; + } else { + + std::shared_ptr > recr = + std::dynamic_pointer_cast >(it->second); + + 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; + } + } + + 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; + } + } + + 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; + } + + } + + 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; + } + } + } + + + } + + virtual std::string to_json() { + std::stringstream ss; + + ss << "\"module\" : {\n"; + ss << " \"name\": \"recruitment\",\n"; + ss << " \"type\": \"Beverton-Holt\",\n"; + ss << " \"id\": " << this->id << ",\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"logit_steep\",\n"; + ss << " \"id\":" << this->logit_steep.id_m << ",\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\": \"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\": \"vector\",\n"; + ss << " \"values\":" << this->log_sigma_recruit << ",\n },\n"; + + ss << " \"parameter\": {\n"; + ss << " \"name\": \"log_devs\",\n"; + ss << " \"id\":" << this->log_devs.id_m << ",\n"; + ss << " \"type\": \"vector\",\n"; + ss << " \"values\":" << this->log_devs << ",\n },\n"; + + + return ss.str(); + } #ifdef TMB_MODEL @@ -104,41 +206,54 @@ 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]); - } - } - 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]); - } - } - info->variable_map[this->log_rzero.id_m] = &(recruitment)->log_rzero; - - //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++) { - recruitment->log_recruit_devs[i] = this->log_devs[i].value_m; - if (this->estimate_log_devs) { - info->RegisterParameter(recruitment->log_recruit_devs[i]); - } else { - recruitment->estimate_log_recruit_devs = false; - } - } + 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(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_recruit_devs + recruitment->log_recruit_devs.resize(this->log_devs.size()); + for (size_t i = 0; i < this->log_devs.size(); i++) { + recruitment->log_recruit_devs[i] = this->log_devs[i].value_m; + 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 @@ -150,6 +265,7 @@ 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 9e5e08094..894a8ca61 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_selectivity.hpp @@ -60,9 +60,9 @@ std::map */ 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 */ + ParameterVector slope; /**< the width of the curve at the inflection_point */ LogisticSelectivityInterface() : SelectivityInterfaceBase() {} @@ -78,12 +78,85 @@ 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); } + virtual void finalize() { + + if (this->finalized) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("Logistic Selectivity " + 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::selectivity_models_iterator it; + + + //search for maturity in Information + it = info->selectivity_models.find(this->id); + //if not found, just return + if (it == info->selectivity_models.end()) { + FIMS_WARNING_LOG("Logistic Selectivity " + fims::to_string(this->id) + " not found in Information."); + return; + } else { + std::shared_ptr > sel = + 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 = sel->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 = sel->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\":\"selectivity\",\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\": \"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\": \"vector\",\n"; + ss << " \"values\":" << this->slope << ",\n}\n"; + + + ss << "}"; + + return ss.str(); + } + + #ifdef TMB_MODEL template @@ -93,30 +166,39 @@ 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 @@ -128,6 +210,7 @@ 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(); @@ -145,12 +228,13 @@ 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 */ + DoubleLogisticSelectivityInterface() : SelectivityInterfaceBase() {} @@ -166,17 +250,124 @@ 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) { + //log warning that finalize has been called more than once. + FIMS_WARNING_LOG("Double Logistic Selectivity " + 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::selectivity_models_iterator it; + + //search for maturity in Information + it = info->selectivity_models.find(this->id); + //if not found, just return + if (it == info->selectivity_models.end()) { + FIMS_WARNING_LOG("Double Logistic Selectivity " + fims::to_string(this->id) + " not found in Information."); + return; + } else { + std::shared_ptr > sel = + std::dynamic_pointer_cast >(it->second); + + + 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; + } + + } + + 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; + } + + } + + 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; + } + + } + + 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; + } + + } + + + + } + } + + 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(); + } + + #ifdef TMB_MODEL template @@ -187,51 +378,75 @@ class DoubleLogisticSelectivityInterface : public SelectivityInterfaceBase { std::shared_ptr > selectivity = std::make_shared >(); - // 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]); - } - } + std::stringstream ss; + // set relative info + selectivity->id = this->id; + 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 d09fdd738..df52f5d5f 100644 --- a/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp +++ b/inst/include/interface/rcpp/rcpp_objects/rcpp_tmb_distribution.hpp @@ -45,7 +45,7 @@ uint32_t interface_observed_data_id_m = */ 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. @@ -88,17 +88,18 @@ 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 @@ -147,9 +148,66 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { 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 @@ -177,8 +235,8 @@ class DnormDistributionsInterface : public DistributionsInterfaceBase { } // 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; + 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++){ @@ -225,7 +283,8 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { 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() {} @@ -288,6 +347,66 @@ class DlnormDistributionsInterface : public DistributionsInterfaceBase { 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 @@ -363,6 +482,7 @@ class DmultinomDistributionsInterface : public DistributionsInterfaceBase { 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 */ + Rcpp::NumericVector lpdf_vec; DmultinomDistributionsInterface() : DistributionsInterfaceBase() {} diff --git a/inst/include/population_dynamics/fleet/fleet.hpp b/inst/include/population_dynamics/fleet/fleet.hpp index 9ca0db58e..cba3a26a2 100644 --- a/inst/include/population_dynamics/fleet/fleet.hpp +++ b/inst/include/population_dynamics/fleet/fleet.hpp @@ -36,10 +36,10 @@ struct Fleet : public fims_model_object::FIMSObject { // Mortality and catchability fims::Vector log_Fmort; /*!< estimated parameter: log Fishing mortality*/ - Type log_q; /*!< estimated parameter: catchability of the fleet */ + fims::Vector log_q; /*!< estimated parameter: catchability of the fleet */ fims::Vector Fmort; /*!< transformed parameter: Fishing mortality*/ - Type q; /*!< transofrmed parameter: the catchability of the fleet */ + fims::Vector q; /*!< transofrmed parameter: the catchability of the fleet */ // derived quantities fims::Vector catch_at_age; /*! { * @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; @@ -92,8 +96,9 @@ struct Fleet : public fims_model_object::FIMSObject { 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); } @@ -125,7 +130,11 @@ struct Fleet : public fims_model_object::FIMSObject { 0); /**q = fims_math::exp(this->log_q); + + for (size_t i = 0; i < this->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; diff --git a/inst/include/population_dynamics/growth/functors/ewaa.hpp b/inst/include/population_dynamics/growth/functors/ewaa.hpp index 62ff5c5d8..fba443353 100644 --- a/inst/include/population_dynamics/growth/functors/ewaa.hpp +++ b/inst/include/population_dynamics/growth/functors/ewaa.hpp @@ -28,6 +28,7 @@ struct EWAAgrowth : public GrowthBase { // the value is the weight at that age (second double) std::map ewaa; /** */ + typedef typename std::map::iterator weight_iterator; EWAAgrowth() : GrowthBase() {} @@ -39,7 +40,11 @@ struct EWAAgrowth : public GrowthBase { * @param a age of the fish, the age vector must start at zero */ virtual const Type evaluate(const double& a) { - Type ret = ewaa[a]; + weight_iterator it = this->ewaa.find(a); + if(it == this->ewaa.end() ){ + return 0.0; + } + Type ret = (*it).second;//itewaa[a]; return ret; } }; diff --git a/inst/include/population_dynamics/population/population.hpp b/inst/include/population_dynamics/population/population.hpp index 9d9694154..9388becd9 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/src/Makevars b/src/Makevars index a3b36318d..506bceeb3 100644 --- a/src/Makevars +++ b/src/Makevars @@ -1,4 +1,4 @@ -CXX_STD = CXX11 -PKG_CXXFLAGS = -DTMB_MODEL -DTMB_EIGEN_DISABLE_WARNINGS -w -CXX17STD = -std=c++11 -w +CXX_STD = CXX17 +PKG_CXXFLAGS = -std=c++17 -DTMB_MODEL -DTMB_EIGEN_DISABLE_WARNINGS -w +CXX17STD = -std=c++17 -w diff --git a/src/Makevars.win b/src/Makevars.win index ab39ecce2..8598441dd 100644 --- a/src/Makevars.win +++ b/src/Makevars.win @@ -1,6 +1,6 @@ -CXX_STD = CXX11 +CXX_STD = CXX17 PKG_CXXFLAGS = -DTMB_MODEL -DTMB_EIGEN_DISABLE_WARNINGS -CXX17STD = -std=c++11 +CXX17STD = -std=c++17 CXX14FLAGS = Wa, -mbig-obj -O3 CXX17FLAGS = Wa, -mbig-obj -O3 diff --git a/tests/gtest/integration_test_population.cpp b/tests/gtest/integration_test_population.cpp index 606f22eb2..a4d023458 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 f1f0efbdc..d70ecbe8d 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 f5c4c9693..9d4b68d2a 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 89a34c30d..3dcfd4533 100644 --- a/tests/gtest/test_population_test_fixture.hpp +++ b/tests/gtest/test_population_test_fixture.hpp @@ -2,6 +2,9 @@ #include "population_dynamics/population/population.hpp" + + + namespace { // Use test fixture to reuse the same configuration of objects for @@ -22,6 +25,7 @@ class PopulationInitializeTestFixture : public testing::Test { 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); } } @@ -84,13 +88,13 @@ class PopulationEvaluateTestFixture : public testing::Test { 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 = log_q_distribution(generator); + fleet->log_q[0] = log_q_distribution(generator); for (int year = 0; year < nyears; year++) { fleet->log_Fmort[year] = log_Fmort_distribution(generator); } @@ -101,7 +105,11 @@ class PopulationEvaluateTestFixture : public testing::Test { population.fleets.push_back(fleet); } population.numbers_at_age.resize((nyears + 1) * nages); - population.Initialize(nyears, nseasons, 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; @@ -243,14 +251,14 @@ class PopulationPrepareTestFixture : public testing::Test { 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 = log_q_distribution(generator); + fleet->log_q[0] = log_q_distribution(generator); for (int year = 0; year < nyears; year++) { fleet->log_Fmort[year] = log_Fmort_distribution(generator); } diff --git a/tests/integration/integration_class.hpp b/tests/integration/integration_class.hpp index e17092bff..1957a6c22 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 10176b528..c1671f822 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 58cd924a8..f21e32b20 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 fc76ab2a9..3d026d5cb 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 4ded4336c..7d5ddbcdd 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 5c3903dc3..f8c7af539 100644 --- a/tests/testthat/test-rcpp-recruitment-interface.R +++ b/tests/testthat/test-rcpp-recruitment-interface.R @@ -7,20 +7,20 @@ 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$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) 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 398fe0c9d..a58da95bf 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 ab86c3a1d..a107b2bcf 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 334d11a1f..4f72db9fa 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) Set up the following parameters: *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 744a65137..61672a56c 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: