Adding migration process equivalent to Stock Synthesis 3

BuildSystem/Version.iss

@@ -1,2 +1,2 @@
 // WARNING: THIS FILE IS AUTOMATICALLY GENERATED BY doBuild version. DO NOT EDIT THIS FILE
-AppVersion='v24.06 (2024-06-26)'
+AppVersion='v24.10 (2024-10-07)'

CASAL2/source/Selectivities/Common/MigrationSS3.Test.cpp

@@ -0,0 +1,124 @@
+/**
+ * @file MigrationSS3.Test.cpp
+ * @author  S. Datta
+ * @version 1.0
+ * @date 2024
+ * @section LICENSE
+ *
+ * Copyright Casal2 Project 2024 - https://github.com/Casal2/
+ *
+ * $Date: 2024-10-02 12:19:00 +1300 (Wed, 02 Oct 2024) $
+ */
+#ifdef TESTMODE
+
+// Headers
+#include <gmock/gmock.h>
+#include <gtest/gtest.h>
+
+#include <boost/lexical_cast.hpp>
+#include <string>
+
+#include "../../TestResources/MockClasses/Model.h"
+#include "../../TestResources/MockClasses/ModelLength.h"
+#include "MigrationSS3.h"
+
+// Namespaces
+namespace niwa {
+
+using ::testing::Return;
+using ::testing::ReturnRef;
+
+/**
+ * Test the results of our selectivity are correct
+
+# R function
+
+migration_stock_synthesis <- function(x,L_, H_, y1_, y2_,
+		alt1_ = F, alt1_rate_ = 0,
+		alt2_ = F, alt2_rate_ = 0,
+		r_ = 1e-4){
+		#
+		# migration rate for ages between L & H interpolated
+		# rates at L is y1, at H is y2
+		# rate gounger ages than L = y1 (unless alt1 is false, then alt1_rate)
+		# similarly for age older than H
+		#
+	y1dash = y1_ + r_
+	y2dash = y2_ + r_
+	k      = log(y2dash/y1dash)/(H_ - L_)
+	if(alt1_) yyoung = alt1_rate_ else yyoung = y1dash 
+	if(alt2_) yold   = alt2_rate_ else yold   = y2dash 
+
+	rate   = rep(0,length(x))
+	younger= x < L_
+	if(length(younger[younger]) > 0) rate[younger] <- yyoung
+	old    = x > H_
+	if(length(old[old]) > 0) rate[old] <- yold
+
+	IN     = x >= L_ & x <= H_
+	if(length(IN[IN]) > 0) rate[IN] <- y1dash * exp(k*(x[IN] - L_))
+
+	return(data.frame(x=x,y=rate))
+}
+
+y <- migration_stock_synthesis(x,5,12,0.3,0.05,T,0.05,T,0.02)  #usual case
+plot(y$x,y$y,type="l",lwd=3,ylim=c(0,0.35),las=1,ylab="Migration rate",xlab="Age")
+text(c(5,12),rep(0.01,2),labels=c("L","H"))
+abline(v=c(5,12),lty=2)
+
+y$y
+
+ */
+
+TEST(Selectivities, Migration_Age_ss3) {
+  shared_ptr<MockModel> model = shared_ptr<MockModel>(new MockModel());
+  EXPECT_CALL(*model, min_age()).WillRepeatedly(Return(1));
+  EXPECT_CALL(*model, max_age()).WillRepeatedly(Return(20));
+  EXPECT_CALL(*model, age_spread()).WillRepeatedly(Return(20));
+  EXPECT_CALL(*model, partition_type()).WillRepeatedly(Return(PartitionType::kAge));
+  niwa::selectivities::MigrationSS3 migration_stock_synthesis(model);
+
+  migration_stock_synthesis.parameters().Add(PARAM_LABEL, "unit_test_migration_ss3", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_TYPE, "not needed in test", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_L, "5", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_H, "12", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_Y1, "0.3", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_Y2, "0.05", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_ALT_RATE_LOWER_SWITCH, "1", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_ALT_RATE_UPPER_SWITCH, "1", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_ALT_RATE_LOWER, "0.05", __FILE__, __LINE__);
+  migration_stock_synthesis.parameters().Add(PARAM_ALT_RATE_UPPER, "0.02", __FILE__, __LINE__);
+
+
+  migration_stock_synthesis.Validate();
+  migration_stock_synthesis.Build();
+  migration_stock_synthesis.Reset();
+
+  ASSERT_THROW(migration_stock_synthesis.GetAgeResult(0, nullptr), std::string);      
+  // values from R
+  EXPECT_NEAR(0.05000000, migration_stock_synthesis.GetAgeResult(1, nullptr), 1e-4); 
+  EXPECT_NEAR(0.05000000, migration_stock_synthesis.GetAgeResult(2, nullptr), 1e-4); 
+  EXPECT_NEAR(0.05000000, migration_stock_synthesis.GetAgeResult(3, nullptr), 1e-4); 
+  EXPECT_NEAR(0.05000000, migration_stock_synthesis.GetAgeResult(4, nullptr), 1e-4); 
+  EXPECT_NEAR(0.30010000, migration_stock_synthesis.GetAgeResult(5, nullptr), 1e-4); 
+  EXPECT_NEAR(0.23238325, migration_stock_synthesis.GetAgeResult(6, nullptr), 1e-4); 
+  EXPECT_NEAR(0.17994659, migration_stock_synthesis.GetAgeResult(7, nullptr), 1e-4); 
+  EXPECT_NEAR(0.13934213, migration_stock_synthesis.GetAgeResult(8, nullptr), 1e-4); 
+  EXPECT_NEAR(0.10789996, migration_stock_synthesis.GetAgeResult(9, nullptr), 1e-4); 
+  EXPECT_NEAR(0.08355262, migration_stock_synthesis.GetAgeResult(10, nullptr), 1e-4); 
+  EXPECT_NEAR(0.06469920, migration_stock_synthesis.GetAgeResult(11, nullptr), 1e-4);   
+  EXPECT_NEAR(0.05010000, migration_stock_synthesis.GetAgeResult(12, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(13, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(14, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(15, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(16, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(17, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(18, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(19, nullptr), 1e-4); 
+  EXPECT_NEAR(0.02000000, migration_stock_synthesis.GetAgeResult(20, nullptr), 1e-4); 
+
+}
+
+}  // namespace niwa
+
+#endif /* ifdef TESTMODE */

CASAL2/source/Selectivities/Common/MigrationSS3.cpp

@@ -0,0 +1,120 @@
+/**
+ * @file MigrationSS3.cpp
+ * @author  S Datta
+ * @version 1.0
+ * @date 2024
+ * @section LICENSE
+ *
+ * Copyright Casal2 Project 2024 - https://github.com/Casal2/
+ *
+ * $Date: 2024-10-02 12:19:00 +1300 (Wed, 02 Oct 2024) $
+ */
+
+// Headers
+#include "MigrationSS3.h"
+
+#include "../../Model/Model.h"
+
+// Namespaces
+namespace niwa {
+namespace selectivities {
+
+/**
+ * Default constructor
+ */
+MigrationSS3::MigrationSS3(shared_ptr<Model> model) : Selectivity(model) {
+  parameters_.Bind<Double>(PARAM_L, &min_set_age_, "Lower age to specifiy migration rate for.", "")->set_lower_bound(0.0, false);
+  parameters_.Bind<Double>(PARAM_H, &max_set_age_, "Upper age to specifiy migration rate for.", "");
+  parameters_.Bind<Double>(PARAM_Y1, &migration_rate_at_L_, "Migration rate at lower age, between 0 and 1", "")->set_range(0.0, 1.0, true, true);
+  parameters_.Bind<Double>(PARAM_Y2, &migration_rate_at_H_, "Migration rate at upper age, between 0 and 1", "")->set_range(0.0, 1.0, true, true);
+  parameters_.Bind<Double>(PARAM_ALT_RATE_LOWER_SWITCH, &alt_rate_lower_switch_, "Switch for setting different value below lower age.", "")->set_range(0.0, 1.0, true, true);
+  parameters_.Bind<Double>(PARAM_ALT_RATE_UPPER_SWITCH, &alt_rate_upper_switch_, "Switch for setting different value above upper age.", "")->set_range(0.0, 1.0, true, true);
+  parameters_.Bind<Double>(PARAM_ALT_RATE_LOWER, &alt_rate_lower_, "Value of migration below lower age.", "")->set_range(0.0, 1.0, true, true);
+  parameters_.Bind<Double>(PARAM_ALT_RATE_UPPER, &alt_rate_upper_, "Value of migration below lower age.", "")->set_range(0.0, 1.0, true, true);
+
+
+  RegisterAsAddressable(PARAM_L, &min_set_age_);
+  RegisterAsAddressable(PARAM_H, &max_set_age_);
+  RegisterAsAddressable(PARAM_Y1, &migration_rate_at_L_);
+  RegisterAsAddressable(PARAM_Y2, &migration_rate_at_H_);
+  RegisterAsAddressable(PARAM_ALT_RATE_LOWER_SWITCH, &alt_rate_lower_switch_);
+  RegisterAsAddressable(PARAM_ALT_RATE_UPPER_SWITCH, &alt_rate_upper_switch_);
+  RegisterAsAddressable(PARAM_ALT_RATE_LOWER, &alt_rate_lower_);
+  RegisterAsAddressable(PARAM_ALT_RATE_UPPER, &alt_rate_upper_);
+  allowed_length_based_in_age_based_model_ = false;
+
+}
+
+/**
+ * Validate this selectivity. This will load the
+ * values that were passed in from the configuration
+ * file and assign them to the local variables.
+ *
+ * Then do some basic checks on the local
+ * variables to ensure they are within the business
+ * rules for the model.
+ */
+void MigrationSS3::DoValidate() {
+
+}
+
+void MigrationSS3::DoReset() {
+  r_ = 1e-4; // small amount to add to values
+  y1dash_ = migration_rate_at_L_ + r_;
+  y2dash_ = migration_rate_at_H_ + r_;
+  k_ = log(y2dash_/y1dash_)/(max_set_age_ - min_set_age_);
+  if (alt_rate_lower_switch_ > 0) {
+	yyoung_ = alt_rate_lower_;
+  } else {
+	yyoung_ = migration_rate_at_L_;  
+  }
+  if (alt_rate_upper_switch_ > 0) {
+	yold_ = alt_rate_upper_;
+  } else {
+	yold_ = migration_rate_at_H_;  
+  }
+}
+/**
+ * The core function
+ */
+Double MigrationSS3::get_value(Double value) {
+  Double sel_val = 0.0;
+  // NOTE: NEED to get the model min and max age or min and max size, and the bin width (assumed to be 1 below)
+  // min age is popn.state.col_min or get_int("min_age")
+  // max age is popn.state.col_max or get_int("max_age")
+  // length bins is popn.state.class_mins; min size is the value in the first bin and max size is the value in the last bin
+
+  // from Ian Doonan, sent on 2024-10-02
+  if (value < min_set_age_) {
+	sel_val = yyoung_;
+  } else if (value > max_set_age_) {
+	sel_val = yold_;
+  } else {
+	sel_val = y1dash_ * exp(k_*(value - min_set_age_));
+  }
+  return sel_val;
+}
+
+/**
+ * The core function
+ */
+Double MigrationSS3::get_value(unsigned value) {
+ Double sel_val = 0.0;
+  // NOTE: NEED to get the model min and max age or min and max size, and the bin width (assumed to be 1 below)
+  // min age is popn.state.col_min or get_int("min_age")
+  // max age is popn.state.col_max or get_int("max_age")
+  // length bins is popn.state.class_mins; min size is the value in the first bin and max size is the value in the last bin
+
+  // from Ian Doonan, sent on 2024-10-02
+  if (value < min_set_age_) {
+	sel_val = yyoung_;
+  } else if (value > max_set_age_) {
+	sel_val = yold_;
+  } else {
+	sel_val = y1dash_ * exp(k_*(value - min_set_age_));
+  }
+  return sel_val;
+}
+
+} /* namespace selectivities */
+} /* namespace niwa */

CASAL2/source/Selectivities/Common/MigrationSS3.h

@@ -0,0 +1,61 @@
+/**
+ * @file MigrationSS3.h
+ * @author  S. Datta
+ * @version 1.0
+ * @date 2024
+ * @section LICENSE
+ *
+ * Copyright Casal2 Project 2024 - https://github.com/Casal2/
+ *
+ * @section DESCRIPTION
+ *
+ * The time class represents a moment of time.
+ *
+ * $Date: 2024-10-02 12:19:00 +1300 (Wed, 02 Oct 2024) $
+ */
+#ifndef MIGRATIONSS3_H_
+#define MIGRATIONSS3_H_
+
+// Headers
+#include "../../Selectivities/Selectivity.h"
+
+// Namespaces
+namespace niwa {
+namespace selectivities {
+
+/**
+ * Class definition
+ */
+class MigrationSS3 : public niwa::Selectivity {
+public:
+  // Methods
+  explicit MigrationSS3(shared_ptr<Model> model);
+  virtual ~MigrationSS3() = default;
+  void DoValidate() override final;
+  void DoReset() override final;
+
+  Double get_value(Double value)  override final; 
+  Double get_value(unsigned value)  override final; 
+
+private:
+  // Members
+  Double min_set_age_;
+  Double max_set_age_;
+  Double migration_rate_at_L_;
+  Double migration_rate_at_H_;
+  Double alt_rate_lower_switch_;
+  Double alt_rate_upper_switch_;
+  Double alt_rate_lower_;
+  Double alt_rate_upper_;
+  // transformed parameters
+  Double r_;
+  Double y1dash_;
+  Double y2dash_;
+  Double k_;
+  Double yyoung_;
+  Double yold_;
+};
+
+} /* namespace selectivities */
+} /* namespace niwa */
+#endif /* MIGRATIONSS3_H_ */

CASAL2/source/Selectivities/Factory.cpp

@@ -30,6 +30,7 @@
 #include "../Selectivities/Common/CompoundAll.h"
 #include "../Selectivities/Common/DoubleNormalPlateau.h"
 #include "../Selectivities/Common/DoubleNormalSS3.h"
+#include "../Selectivities/Common/MigrationSS3.h"
 #include "../Selectivities/Common/Increasing.h"
 #include "../Selectivities/Common/InverseLogistic.h"
 #include "../Selectivities/Common/MultiSelectivity.h"
@@ -76,6 +77,8 @@ Selectivity* Factory::Create(shared_ptr<Model> model, const string& object_type,
       result = new DoubleNormalPlateau(model);
     else if (sub_type == PARAM_DOUBLE_NORMAL_SS3)
       result = new DoubleNormalSS3(model);
+    else if (sub_type == PARAM_MIGRATION_SS3)
+      result = new MigrationSS3(model);
     else if (sub_type == PARAM_INCREASING)
       result = new Increasing(model);
     else if (sub_type == PARAM_INVERSE_LOGISTIC)

CASAL2/source/Translations/English_UK.h

@@ -146,6 +146,10 @@
 #define PARAM_ALGORITHM							  "algorithm"
 #define PARAM_ALL_VALUES_BOUNDED                  "all_values_bounded"
 #define PARAM_ALL_VALUES                          "all_values"
+#define PARAM_ALT_RATE_LOWER_SWITCH               "switch_for_alternate_lower_rate"
+#define PARAM_ALT_RATE_UPPER_SWITCH				  "switch_for_alternate_upper_rate"
+#define PARAM_ALT_RATE_LOWER                      "alternate_lower_rate"
+#define PARAM_ALT_RATE_UPPER                      "alternate_upper_rate"
 #define PARAM_ALPHA                               "alpha"
 #define PARAM_ANNUAL_DURATION                     "annual_duration"
 #define PARAM_ANNUAL_MORTALITY_RATE               "annual_mortality_rate"
@@ -472,6 +476,7 @@
 #define PARAM_METROPOLIS_HASTINGS                 "metropolis_hastings"
 #define PARAM_METHOD_OF_REMOVAL                   "method_of_removal"
 #define PARAM_MIGRATION_MOVEMENT                  "migration"
+#define PARAM_MIGRATION_SS3                       "migration_stock_synthesis"
 #define PARAM_MIN_AGE                             "min_age"
 #define PARAM_MIN_BOBYQA                          "min_bobyqa"
 #define PARAM_MIN_BOX_CONSTRAINED                 "min_box_constrained"

CASAL2/source/Version.h

@@ -1,13 +1,13 @@
 // WARNING: THIS FILE IS AUTOMATICALLY GENERATED BY doBuild version. DO NOT EDIT THIS FILE
 #ifndef VERSION_H_
 #define VERSION_H_
-#define VERSION "24.06 (2024-06-26)"
-#define SOURCE_REPOS "git@github.com:casal2/casal2.git:master"
-#define VERSION_NUMBER "24.06"
-#define SOURCE_CONTROL_REVISION efe547d0eac16aefb2c97b7f7fa3cfd8590d0d4a
-#define SOURCE_CONTROL_DATE "2024-06-26"
+#define VERSION "24.10 (2024-10-07)"
+#define SOURCE_REPOS "git@github.com:SamikDatta/CASAL2.git:add_mig_ss3"
+#define VERSION_NUMBER "24.10"
+#define SOURCE_CONTROL_REVISION 9ac0f6ba2dfbbe83cc656d33bb6b6c86ecd57078
+#define SOURCE_CONTROL_DATE "2024-10-07"
 #define SOURCE_CONTROL_YEAR "2024"
-#define SOURCE_CONTROL_MONTH "June"
-#define SOURCE_CONTROL_TIME "09:08:16"
-#define SOURCE_CONTROL_VERSION "2024-06-26 09:08:16 UTC (rev. efe547d0e)"
+#define SOURCE_CONTROL_MONTH "October"
+#define SOURCE_CONTROL_TIME "23:11:56"
+#define SOURCE_CONTROL_VERSION "2024-10-07 23:11:56 UTC (rev. 9ac0f6ba2)"
 #endif

Documentation/UserManual/Version.tex

@@ -1,10 +1,10 @@
 % WARNING: THIS FILE IS AUTOMATICALLY GENERATED BY doBuild version. DO NOT EDIT THIS FILE
-\newcommand{\Version}{24.06 (2024-06-26)}
-\newcommand{\SourceRepos}{git@github.com:casal2/casal2.git:master}
-\newcommand{\VersionNumber}{24.06}
-\newcommand{\SourceControlRevision}{efe547d0eac16aefb2c97b7f7fa3cfd8590d0d4a}
-\newcommand{\SourceControlDateDoc}{2024-06-26}
+\newcommand{\Version}{24.10 (2024-10-07)}
+\newcommand{\SourceRepos}{git@github.com:SamikDatta/CASAL2.git:add_mig_ss3}
+\newcommand{\VersionNumber}{24.10}
+\newcommand{\SourceControlRevision}{9ac0f6ba2dfbbe83cc656d33bb6b6c86ecd57078}
+\newcommand{\SourceControlDateDoc}{2024-10-07}
 \newcommand{\SourceControlYearDoc}{2024}
-\newcommand{\SourceControlMonthDoc}{June}
-\newcommand{\SourceControlTimeDoc}{09:08:16}
-\newcommand{\SourceControlVersion}{2024-06-26 09:08:16 UTC (rev. efe547d0e)}
+\newcommand{\SourceControlMonthDoc}{October}
+\newcommand{\SourceControlTimeDoc}{23:11:56}
+\newcommand{\SourceControlVersion}{2024-10-07 23:11:56 UTC (rev. 9ac0f6ba2)}