Age size observation 202206

.github/workflows/Casal2_testsuite_modelrunner_archive.yml

@@ -128,7 +128,38 @@ jobs:
         with:
           name: Casal2-linux-unittest
           path: BuildSystem/Casal2/unittests.log
-
+  Base-R-Ubuntu:
+    needs: Compile-Casal2-linux
+    runs-on: ubuntu-20.04
+    steps:
+      - uses: actions/checkout@v2
+      - uses: r-lib/actions/setup-r@v2
+      - uses: r-lib/actions/setup-pandoc@v2
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          working-directory: R-libraries/casal2
+          extra-packages: |
+            any::devtools
+            any::bayesplot
+            any::purrr
+            any::DHARMa
+            any::mvtnorm
+          dependencies: '"hard"'
+      - name: Install Casal2 base
+        env:
+           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # https://github.com/settings/tokens
+        run: Rscript -e 'devtools::install_github("https://github.com/NIWAFisheriesModelling/CASAL2", subdir="R-libraries/casal2", ref = "HEAD")'
+      ## download model runs
+      - name: Download TestModels artifact
+        uses: actions/[email protected]
+        with:
+          # Artifact name
+          name: Casal2-linux-modelrunner # optional
+          # Destination path
+          path: CheckTestModels # optional
+      ## now run an R script that will check R-library compatible with exe output
+      - name: Run Casal2 R test suite
+        run: Rscript -e 'source(file.path("R-libraries","run_rlibrary_tests.R"))'          
   Compile-Casal2-Windows:
     # https://docs.github.com/en/actions/using-github-hosted-runners/about-github-hosted-runners#supported-runners-and-hardware-resources
     # https://github.com/actions/virtual-environments/blob/main/images/win/Windows2019-Readme.md

BuildSystem/Version.iss

@@ -1,2 +1,2 @@
 // WARNING: THIS FILE IS AUTOMATICALLY GENERATED BY doBuild version. DO NOT EDIT THIS FILE
-AppVersion='v22.06 (2022-06-28)'
+AppVersion='v22.07 (2022-07-20)'

CASAL2/source/AddressableInputLoader/AddressableInputLoader.cpp

@@ -125,6 +125,8 @@ void AddressableInputLoader::LoadValues(unsigned index) {
     auto estimate = model_->managers()->estimate()->GetEstimate(iter.first);
     if (estimate != nullptr) {
       estimate->set_value(*iter.second);
+      estimate->set_initial_value(*iter.second);
+      estimate->flag_value_has_been_initialised();
       ++estimate_count;
     }
   }

CASAL2/source/AgeLengths/AgeLength.cpp

@@ -169,6 +169,14 @@ void AgeLength::Build() {
   UpdateYearContainers();
   // Build Age-length matrix
   PopulateAgeLengthMatrix();
+  // used in get_pdf_with_sized_based_selectivity
+  quantiles_.resize(5,0);
+  inverse_values_.resize(5,0);
+  quantiles_[0]=0.1;
+  quantiles_[1]=0.3;
+  quantiles_[2]=0.5;
+  quantiles_[3]=0.7;
+  quantiles_[4]=0.9;
 }
 
 /**
@@ -848,6 +856,89 @@ void AgeLength::populate_age_length_matrix(vector<Double> numbers_at_age, vector
     }
   }
 }
+/**
+ * @details get_pdf
+ * @param age actual age
+ * @param length length
+ * @param year the year (needed for Data type)
+ * @param time_step needed for the correct inter annual interpolation.
+ * Given an age and length, figure the pdf of the size-at-age distribution for partition element [row,age]
+ */
+Double AgeLength::get_pdf(unsigned age, Double length, unsigned year, unsigned time_step) {
+  Double mu    = calculate_mean_length(year, time_step, age);
+  Double cv    = cvs_[year - year_offset_][time_step - time_step_offset_][age - age_offset_];
+  Double sigma = cv * mu;
+  //LOG_FINE() << "year: " << year << "; time_step " << time_step << "; age: " << age << "; mu: " << mu << "; cv: " << cv << "; sigma: " << sigma;
+  if (distribution_ == Distribution::kLogNormal) {
+    // Transform parameters in to log space
+    Double Lvar  = log(cv * cv + 1.0);
+    mu    = log(mu) - Lvar / 2.0;
+    sigma = sqrt(Lvar);
+    return utilities::math::dlognorm(length, mu, sigma);
+  } else {
+    return utilities::math::dnorm(length, mu, sigma);
+  }
+  return 0.0;
+}
+
+/**
+ * @details get_pdf_with_sized_based_selectivity
+ * @param age actual age
+ * @param length length
+ * @param year the year (needed for Data type)
+ * @param time_step needed for the correct inter annual interpolation.
+ * @param selectivity selectivity pointer assumed to be length-based selectivity
+ * Given an age and length, figure the pdf of the size-at-age distribution for partition element [age]
+ * The pdf is modified by a size-based selectivity.
+ * f_new(size) = Sel(size)f_old(size) / integral_over_size=s of Sel(s)f_old(s).
+ *  We approximate the integral by a discrete approximation over 5 appropriately spaced points.
+ */
+Double AgeLength::get_pdf_with_sized_based_selectivity(unsigned age, Double length, unsigned year, unsigned time_step, Selectivity* selectivity) {
+  Double numerator = get_pdf(age, length, year, time_step) * selectivity->GetResult(length);
+  if (numerator == 0) return 0;
+  Double denominator = 0.0;
+  // clear this inverse_values_ container
+  fill(inverse_values_.begin(), inverse_values_.end(), 0.0);
+  //
+  get_cdf_inverse(age, year, time_step, inverse_values_);
+  // now evaluate the selectivity over these integration points
+  for(unsigned i = 0; i < inverse_values_.size(); ++i) {
+    denominator += selectivity->GetResult(inverse_values_[i]);
+  }
+  // the integral is over
+  LOG_FINE() << "numerator = " << numerator << " denominator " << denominator;
+  if (denominator == 0) return 0;
+  return numerator / (denominator * 0.2);
+}
+
+/**
+ * @details get_cdf_inverse
+ * @param age actual age
+ * @param year the year (needed for Data type)
+ * @param time_step needed for the correct inter annual interpolation.
+ * at the specified time, figure the inverse cdf of the size-at-age
+ * distribution for partition element [age]
+ * used by get_pdf_with_sized_based_selectivity
+ */
+
+void AgeLength::get_cdf_inverse(unsigned age, unsigned year, unsigned time_step, vector<Double>& inverse_result) {
+  Double mu    = calculate_mean_length(year, time_step, age);
+  Double cv    = cvs_[year - year_offset_][time_step - time_step_offset_][age - age_offset_];
+  Double sigma = cv * mu;
+  if(inverse_result.size() != quantiles_.size())
+    LOG_CODE_ERROR() << "inverse_result.size() != quantiles_.size()";
+  for(unsigned i = 0; i < inverse_result.size(); ++i) {
+    if (distribution_ == Distribution::kLogNormal) {
+      // Transform parameters in to log space
+      Double Lvar  = log(cv * cv + 1.0);
+      mu    = log(mu) - Lvar / 2.0;
+      sigma = sqrt(Lvar);
+      inverse_result[i] =  utilities::math::qlognorm(quantiles_[i],mu,sigma);
+    } else {
+      inverse_result[i] =  utilities::math::qnorm(quantiles_[i],mu,sigma);
+    }
+  }
+}
 
 /**
  * @details FillReportCache
@@ -899,4 +990,7 @@ void AgeLength::FillReportCache(ostringstream& cache) {
     }
   }
 }
+
+
+
 } /* namespace niwa */

CASAL2/source/AgeLengths/AgeLength.h

@@ -72,6 +72,9 @@ class AgeLength : public niwa::base::Object {
                                   vector<vector<Double>>& numbers_by_age_length);  // overloaded for the case with no selectivity and class has bespoke length bins
   void populate_age_length_matrix(vector<Double> numbers_at_age, vector<vector<Double>>& numbers_by_age_length,
                                   Selectivity* selectivity);  // overloaded for the case with no selectivity and class has bespoke length bins
+  Double get_pdf(unsigned age, Double length, unsigned year, unsigned time_step);
+  Double get_pdf_with_sized_based_selectivity(unsigned age, Double length, unsigned year, unsigned time_step, Selectivity* selectivity);
+  void get_cdf_inverse(unsigned age, unsigned year, unsigned time_step, vector<Double>& inverse_result);
 
   // For reporting in the AgeLength
   void FillReportCache(ostringstream& cache);
@@ -110,6 +113,8 @@ class AgeLength : public niwa::base::Object {
   vector<vector<Double>>  mean_weight_by_timestep_age_;  // mean_weight_by_timestep_age_[time_step][age]
   vector<unsigned>        age_length_matrix_years_;
   map<unsigned, unsigned> age_length_matrix_year_key_;  // [year, dimension]
+  vector<Double>          quantiles_;
+  vector<Double>          inverse_values_;
   // because this may not be in sequential order (see method BuildAgeLengthMatrixForTheseYears) we have a key which maps the right year with first dimension of
   // age_length_transition_matrix_
   vector<vector<Double>>                 numbers_by_age_length_transition_;  // age x length used as a temporarey container

CASAL2/source/ConfigurationLoader/Loader.cpp

@@ -446,6 +446,7 @@ void Loader::ParseBlock(shared_ptr<Model> model, vector<FileLine>& block) {
  *
  */
 void Loader::HandleInlineDefinitions() {
+  LOG_FINEST() << "Enter: HandleInlineDefinitions";
   vector<string>           replacements = {"one", "two", "three", "four", "five", "six", "seven", "eight", "nine"};
   vector<string>           temp;  // temp to old outputs of splits we only want for a short period
   vector<vector<FileLine>> new_blocks;

CASAL2/source/Estimates/Estimate.cpp

@@ -49,7 +49,8 @@ void Estimate::Validate() {
  * Build the estimate
  */
 void Estimate::Build() {
-
+  if(!value_been_initialised_)
+    set_initial_value(value());
   Reset();
 }
 
@@ -63,6 +64,7 @@ void Estimate::Reset() {
    */
   if (utilities::math::IsEqual(lower_bound_, upper_bound_))
     set_value(value());
+
 }
 
 /**

CASAL2/source/Estimates/Estimate.h

@@ -65,7 +65,10 @@ class Estimate : public niwa::base::Object {
   void                   set_estimated_in_phasing(bool new_value) { estimated_in_phasing_ = new_value; }
   void                   set_mcmc_fixed(bool new_value) { mcmc_fixed_ = new_value; }
   Double                 value() { return *target_; }
+  Double                 get_initial_value() { return initial_value_; }
   void                   set_value(Double new_value);
+  void                   set_initial_value(Double new_value) {initial_value_ = new_value;}; // currently used by AddressableInputLoader if users use -i
+  void                   flag_value_has_been_initialised() {value_been_initialised_ = true;};
   bool                   mcmc_fixed() const { return mcmc_fixed_; }
   void                   set_in_objective_function(bool value) { in_objective_ = value; }
   bool                   in_objective_function() const { return in_objective_; }
@@ -89,6 +92,8 @@ class Estimate : public niwa::base::Object {
   bool            estimated_    = true;
   bool            estimated_in_phasing_               = true;
   bool            in_objective_ = true;
+  Double          initial_value_;
+  bool            value_been_initialised_ = false;
 };
 } /* namespace niwa */
 #endif /* ESTIMATE_H_ */

CASAL2/source/Likelihoods/Common/Bernoulli.cpp

@@ -0,0 +1,81 @@
+/**
+ * @file Bernoulli.cpp
+ * @author  C.Marsh
+ * @version 1.0
+ * @date 2022
+ * @section LICENSE
+ *
+ * Copyright NIWA Science 2022 - www.niwa.co.nz
+ *
+ */
+
+// Headers
+#include "Bernoulli.h"
+
+#include "../../Utilities/Math.h"
+#include "../../Utilities/RandomNumberGenerator.h"
+
+// Namespaces
+namespace niwa {
+namespace likelihoods {
+
+namespace math = niwa::utilities::math;
+
+/**
+ * Adjust the error value based on the process error
+ *
+ * @param process_error The observations process_error
+ * @param error_value The observations error_value
+ * @return An adjusted error value
+ */
+Double Bernoulli::AdjustErrorValue(const Double process_error, const Double error_value) {
+  return 1.0;
+}
+
+/**
+ * Calculate the scores
+ *
+ * @param comparisons A collection of comparisons passed by the observation
+ *
+ */
+void Bernoulli::GetScores(map<unsigned, vector<observations::Comparison> >& comparisons) {
+  for (auto year_iterator = comparisons.begin(); year_iterator != comparisons.end(); ++year_iterator) {
+    for (observations::Comparison& comparison : year_iterator->second) {
+      Double error_value = AdjustErrorValue(comparison.process_error_, comparison.error_value_);
+      if(comparison.observed_ == 1.0) {
+        Double score = log(math::ZeroFun(comparison.expected_, comparison.delta_));
+        comparison.adjusted_error_ = error_value;
+        comparison.score_          = -score;
+      } else if (comparison.observed_ == 0.0) {
+        Double score = log(math::ZeroFun(1.0 - comparison.expected_, comparison.delta_));
+        comparison.adjusted_error_ = error_value;
+        comparison.score_          = -score;
+      } else {
+        LOG_FATAL() << "found an observed value = " << comparison.observed_  << " that wasn't '1' or '0' which is expected for the Bernoulli likelihood. This is not allowed";
+      }
+    }
+  }
+}
+
+/**
+ * Simulate observed values
+ *
+ * @param comparisons A collection of comparisons passed by the observation
+ */
+void Bernoulli::SimulateObserved(map<unsigned, vector<observations::Comparison> >& comparisons) {
+  utilities::RandomNumberGenerator& rng = utilities::RandomNumberGenerator::Instance();
+  auto   iterator    = comparisons.begin();
+  for (; iterator != comparisons.end(); ++iterator) {
+    LOG_FINE() << "Simulating values for year: " << iterator->first;
+    for (observations::Comparison& comparison : iterator->second) {
+      if(rng.uniform() <= comparison.expected_) {
+        comparison.observed_ = 1.0;
+      } else {
+        comparison.observed_ = 0.0;
+      }
+    }
+  }
+}
+
+} /* namespace likelihoods */
+} /* namespace niwa */

CASAL2/source/Likelihoods/Common/Bernoulli.h

@@ -0,0 +1,40 @@
+/**
+ * @file Bernoulli.h
+ * @author  C.Marsh
+ * @version 1.0
+ * @date 2022
+ * @section LICENSE
+ *
+ * Copyright NIWA Science 2022 - www.niwa.co.nz
+ *
+ * @section DESCRIPTION
+ *
+ * The time class represents a moment of time.
+ */
+#ifndef LIKELIHOODS_BERNOULLI_H_
+#define LIKELIHOODS_BERNOULLI_H_
+
+// Headers
+#include "../../Likelihoods/Likelihood.h"
+
+// Namespaces
+namespace niwa {
+namespace likelihoods {
+
+/**
+ * Class definition
+ */
+class Bernoulli : public niwa::Likelihood {
+public:
+  // Methods
+  Bernoulli(shared_ptr<Model> model) : Likelihood(model){};
+  virtual ~Bernoulli() = default;
+  void   DoValidate() override final{};
+  Double AdjustErrorValue(const Double process_error, const Double error_value) override final;
+  void   SimulateObserved(map<unsigned, vector<observations::Comparison> >& comparisons) override final;
+  void   GetScores(map<unsigned, vector<observations::Comparison> >& comparisons) override final;
+};
+
+} /* namespace likelihoods */
+} /* namespace niwa */
+#endif /* LIKELIHOODS_BERNOULLI_H_ */

CASAL2/source/Likelihoods/Factory.cpp

@@ -5,14 +5,15 @@
  * @date 22/03/2013
  * @section LICENSE
  *
- * Copyright NIWA Science �2013 - www.niwa.co.nz
+ * Copyright NIWA Science �2013 - www.niwa.co.nz
  *
  * $Date: 2008-03-04 16:33:32 +1300 (Tue, 04 Mar 2008) $
  */
 
 // Headers
 #include "Factory.h"
 
+#include "../Likelihoods/Common/Bernoulli.h"
 #include "../Likelihoods/Common/Binomial.h"
 #include "../Likelihoods/Common/BinomialApprox.h"
 #include "../Likelihoods/Common/Dirichlet.h"
@@ -40,6 +41,8 @@ Likelihood* Factory::Create(shared_ptr<Model> model, const string& object_type,
 
   if (sub_type == PARAM_BINOMIAL)
     result = new Binomial(model);
+  else if (sub_type == PARAM_BERNOULLI)
+    result = new Bernoulli(model);
   else if (sub_type == PARAM_BINOMIAL_APPROX)
     result = new BinomialApprox(model);
   else if (sub_type == PARAM_DIRICHLET)