Skip to content

Commit 4dbc264

Browse files
alistairdunn1alistairdunn1
committed
FIX: MH MCMC ; Tidy AgeLength & TagByLength code; Added error warning on HMC as we need to check; UPD clang-format for Craig
1 parent abe13d2 commit 4dbc264

File tree

6 files changed

+206
-218
lines changed

6 files changed

+206
-218
lines changed

‎.clang-format

+1-1
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ AllowShortCaseLabelsOnASingleLine: false
1919
AllowShortFunctionsOnASingleLine: Inline
2020
AllowShortLambdasOnASingleLine: All
2121
AllowShortIfStatementsOnASingleLine: Never
22-
AllowShortLoopsOnASingleLine: true
22+
AllowShortLoopsOnASingleLine: false
2323
AlwaysBreakAfterDefinitionReturnType: None
2424
AlwaysBreakAfterReturnType: None
2525
AlwaysBreakBeforeMultilineStrings: true

‎CASAL2/source/AgeLengths/AgeLength.cpp

+155-166
Large diffs are not rendered by default.

‎CASAL2/source/AgeLengths/AgeLength.h

+6-3
Original file line numberDiff line numberDiff line change
@@ -60,8 +60,10 @@ class AgeLength : public niwa::base::Object {
6060
vector<int>& map_length_bin_ndx); // overloaded for the case with no selectivity and class has bespoke length bins
6161
void populate_numbers_at_age_with_length_based_exploitation(vector<Double>& numbers_at_age, vector<Double>& numbers_at_age_with_exploitation, Double& exploitation_by_length,
6262
unsigned model_length_bin_ndx, Selectivity* selectivity); //
63-
void populate_age_length_matrix(vector<Double> numbers_at_age, vector<vector<Double>>& numbers_by_age_length); // overloaded for the case with no selectivity and class has bespoke length bins
64-
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
63+
void populate_age_length_matrix(vector<Double> numbers_at_age,
64+
vector<vector<Double>>& numbers_by_age_length); // overloaded for the case with no selectivity and class has bespoke length bins
65+
void populate_age_length_matrix(vector<Double> numbers_at_age, vector<vector<Double>>& numbers_by_age_length,
66+
Selectivity* selectivity); // overloaded for the case with no selectivity and class has bespoke length bins
6567

6668
// For reporting in the AgeLength
6769
void FillReportCache(ostringstream& cache);
@@ -100,7 +102,8 @@ class AgeLength : public niwa::base::Object {
100102
vector<vector<Double>> mean_weight_by_timestep_age_; // mean_weight_by_timestep_age_[time_step][age]
101103
vector<unsigned> age_length_matrix_years_;
102104
map<unsigned, unsigned> age_length_matrix_year_key_; // [year, dimension]
103-
// 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_
105+
// because this may not be in sequential order (see method BuildAgeLengthMatrixForTheseYears) we have a key which maps the right year with first dimension of
106+
// age_length_transition_matrix_
104107
vector<vector<Double>> numbers_by_age_length_transition_; // age x length used as a temporarey container
105108
vector<vector<vector<vector<Double>>>> age_length_transition_matrix_; // dims years x timesteps x age x length
106109
// these members are used in the functions

‎CASAL2/source/MCMCs/Common/HamiltonianMonteCarlo.cpp

+9-7
Original file line numberDiff line numberDiff line change
@@ -35,15 +35,16 @@ using niwa::utilities::Vector_Debug;
3535
using niwa::utilities::Vector_Scale;
3636

3737
/**
38-
* @brief norm2. Return cummulative squard values of
38+
* @brief norm2. Return cumulative squard values of
3939
* the vector
4040
*
4141
* @param target
4242
* @return double
4343
*/
4444
double norm2(const std::vector<double>& target) {
4545
double result = 0.0;
46-
for (auto d : target) result += d * d;
46+
for (auto d : target)
47+
result += d * d;
4748
return result;
4849
}
4950

@@ -104,12 +105,12 @@ void HamiltonianMonteCarlo::DoExecute(shared_ptr<ThreadPool> thread_pool) {
104105

105106
estimate_lower_bounds_.resize(estimates_.size(), 0.0);
106107
estimate_upper_bounds_.resize(estimates_.size(), 0.0);
107-
for(unsigned i = 0; i < estimates_.size(); ++i) {
108+
for (unsigned i = 0; i < estimates_.size(); ++i) {
108109
estimate_lower_bounds_[i] = estimates_[i]->lower_bound();
109110
estimate_upper_bounds_[i] = estimates_[i]->upper_bound();
110111
}
111-
//estimate_lower_bounds_ = model_->managers()->estimate()->lower_bounds();
112-
//estimate_upper_bounds_ = model_->managers()->estimate()->upper_bounds();
112+
// estimate_lower_bounds_ = model_->managers()->estimate()->lower_bounds();
113+
// estimate_upper_bounds_ = model_->managers()->estimate()->upper_bounds();
113114

114115
// Lambda to increment vector e by scaled g
115116
auto increment_vector = [](vector<double> e, vector<double> g, double delta) {
@@ -123,7 +124,8 @@ void HamiltonianMonteCarlo::DoExecute(shared_ptr<ThreadPool> thread_pool) {
123124

124125
// Lambda to easily run a model run with some iterations
125126
auto quick_run = [this](const vector<double>& candidates) {
126-
for (unsigned i = 0; i < candidates.size(); ++i) estimates_[i]->set_value(candidates[i]);
127+
for (unsigned i = 0; i < candidates.size(); ++i)
128+
estimates_[i]->set_value(candidates[i]);
127129

128130
model_->FullIteration();
129131
model_->objective_function().CalculateScore();
@@ -178,7 +180,7 @@ void HamiltonianMonteCarlo::DoExecute(shared_ptr<ThreadPool> thread_pool) {
178180
if (q_score >= previous_score) {
179181
ratio = exp(previous_score - q_score);
180182
}
181-
183+
LOG_WARNING() << " HMC: LOOKS WRONG: Don't this rejection should be saving the last chain value. NEED TO CHECK";
182184
// Check if we accept this jump
183185
if (math::IsEqual(ratio, 1.0) || rng_uniform < ratio) {
184186
LOG_MEDIUM() << "Accepting Jump (ratio: " << ratio << ", rng: " << rng_uniform << ")" << endl;

‎CASAL2/source/MCMCs/Common/RandomWalkMetropolisHastings.cpp

+17-32
Original file line numberDiff line numberDiff line change
@@ -52,7 +52,6 @@ void RandomWalkMetropolisHastings::DoExecute(shared_ptr<ThreadPool> thread_pool)
5252

5353
mcmc::ChainLink last_link;
5454

55-
bool accept_jump = false;
5655
double score = 0.0;
5756
int invalid_counter = 0;
5857

@@ -64,14 +63,14 @@ void RandomWalkMetropolisHastings::DoExecute(shared_ptr<ThreadPool> thread_pool)
6463
UpdateStepSize();
6564
UpdateCovarianceMatrix();
6665
GenerateNewCandidates();
67-
accept_jump = false;
68-
score = 0.0;
66+
score = 0.0;
6967

7068
if (WithinBounds()) {
7169
++jumps_;
7270
++jumps_since_adapt_;
7371

74-
for (unsigned i = 0; i < candidates_.size(); ++i) estimates_[i]->set_value(candidates_[i]);
72+
for (unsigned i = 0; i < candidates_.size(); ++i)
73+
estimates_[i]->set_value(candidates_[i]);
7574

7675
model_->FullIteration();
7776
obj_function.CalculateScore();
@@ -86,8 +85,6 @@ void RandomWalkMetropolisHastings::DoExecute(shared_ptr<ThreadPool> thread_pool)
8685
// double rng_uniform = rng.uniform();
8786
if (math::IsEqual(ratio, 1.0) || rng.uniform() < ratio) {
8887
LOG_MEDIUM() << "Accept: Possible. Iteration = " << jumps_ << ", score = " << score << " Previous score " << previous_score;
89-
// Accept this jump
90-
accept_jump = true;
9188
successful_jumps_++;
9289
successful_jumps_since_adapt_++;
9390
previous_score = score;
@@ -102,31 +99,20 @@ void RandomWalkMetropolisHastings::DoExecute(shared_ptr<ThreadPool> thread_pool)
10299
mcmc_state_ = PARAM_MCMC;
103100
else
104101
mcmc_state_ = PARAM_BURN_IN;
105-
if (accept_jump) {
106-
mcmc::ChainLink new_link{.iteration_ = jumps_,
107-
.mcmc_state_ = mcmc_state_,
108-
.score_ = obj_function.score(),
109-
.likelihood_ = obj_function.likelihoods(),
110-
.prior_ = obj_function.priors(),
111-
.penalty_ = obj_function.penalties(),
112-
.additional_priors_ = obj_function.additional_priors(),
113-
.jacobians_ = obj_function.jacobians(),
114-
.acceptance_rate_ = double(successful_jumps_) / double(jumps_),
115-
.acceptance_rate_since_adapt_ = double(successful_jumps_since_adapt_) / double(jumps_since_adapt_),
116-
.step_size_ = step_size_,
117-
.values_ = candidates_};
118-
119-
chain_.push_back(new_link);
120-
} else {
121-
// Copy the last chain accepted link to the end of the vector
122-
auto temp = *chain_.rbegin();
123-
temp.iteration_ = jumps_;
124-
temp.mcmc_state_ = mcmc_state_;
125-
temp.acceptance_rate_ = double(successful_jumps_) / double(jumps_);
126-
temp.acceptance_rate_since_adapt_ = double(successful_jumps_since_adapt_) / double(jumps_since_adapt_);
127-
chain_.push_back(temp);
128-
}
129-
102+
mcmc::ChainLink new_link{.iteration_ = jumps_,
103+
.mcmc_state_ = mcmc_state_,
104+
.score_ = obj_function.score(),
105+
.likelihood_ = obj_function.likelihoods(),
106+
.prior_ = obj_function.priors(),
107+
.penalty_ = obj_function.penalties(),
108+
.additional_priors_ = obj_function.additional_priors(),
109+
.jacobians_ = obj_function.jacobians(),
110+
.acceptance_rate_ = double(successful_jumps_) / double(jumps_),
111+
.acceptance_rate_since_adapt_ = double(successful_jumps_since_adapt_) / double(jumps_since_adapt_),
112+
.step_size_ = step_size_,
113+
.values_ = candidates_};
114+
115+
chain_.push_back(new_link);
130116
// LOG_MEDIUM() << "Storing: Successful Jumps " << successful_jumps_ << " Jumps : " << jumps_;
131117
model_->managers()->report()->Execute(model_, State::kIterationComplete);
132118
}
@@ -136,7 +122,6 @@ void RandomWalkMetropolisHastings::DoExecute(shared_ptr<ThreadPool> thread_pool)
136122
LOG_MEDIUM() << "Reject: Bounds. Iteration = " << jumps_ << ", score = " << score << " Previous score " << previous_score << ". There were " << invalid_counter
137123
<< " invalid jumps";
138124
}
139-
140125
} while (jumps_ < length_);
141126
}
142127

‎CASAL2/source/Processes/Age/TagByLength.cpp

+18-9
Original file line numberDiff line numberDiff line change
@@ -110,7 +110,8 @@ void TagByLength::DoValidate() {
110110
}
111111
}
112112
}
113-
for (auto& split_category : split_category_labels) split_from_category_labels_.push_back(split_category);
113+
for (auto& split_category : split_category_labels)
114+
split_from_category_labels_.push_back(split_category);
114115
} else
115116
split_from_category_labels_.push_back(category);
116117
}
@@ -273,7 +274,8 @@ void TagByLength::DoBuild() {
273274
numbers_at_age_by_category_[i].resize(model_->age_spread(), 0.0);
274275
}
275276

276-
for (unsigned i = 0; i < model_->age_spread(); ++i) numbers_at_age_and_length_[i].resize(n_length_bins_, 0.0);
277+
for (unsigned i = 0; i < model_->age_spread(); ++i)
278+
numbers_at_age_and_length_[i].resize(n_length_bins_, 0.0);
277279

278280
if (penalty_label_ != "")
279281
penalty_ = model_->managers()->penalty()->GetPenalty(penalty_label_);
@@ -298,7 +300,8 @@ void TagByLength::DoBuild() {
298300

299301
for (unsigned year_ndx = 0; year_ndx < years_.size(); ++year_ndx) {
300302
proportion_by_length_[year_ndx].resize(n_length_bins_, 0.0);
301-
for (unsigned length_ndx = 0; length_ndx < n_length_bins_; ++length_ndx) tagged_fish_by_year_[year_ndx] += numbers_[years_[year_ndx]][length_ndx];
303+
for (unsigned length_ndx = 0; length_ndx < n_length_bins_; ++length_ndx)
304+
tagged_fish_by_year_[year_ndx] += numbers_[years_[year_ndx]][length_ndx];
302305
for (unsigned length_ndx = 0; length_ndx < n_length_bins_; ++length_ndx)
303306
proportion_by_length_[year_ndx][length_ndx] = numbers_[years_[year_ndx]][length_ndx] / tagged_fish_by_year_[year_ndx];
304307
}
@@ -363,7 +366,8 @@ void TagByLength::DoExecute() {
363366

364367
for (unsigned length_ndx = 0; length_ndx < n_length_bins_; ++length_ndx) {
365368
sum_age = 0.0;
366-
for (unsigned age_ndx = 0; age_ndx < model_->age_spread(); ++age_ndx) sum_age += numbers_at_age_and_length_[age_ndx][length_ndx];
369+
for (unsigned age_ndx = 0; age_ndx < model_->age_spread(); ++age_ndx)
370+
sum_age += numbers_at_age_and_length_[age_ndx][length_ndx];
367371
for (unsigned age_ndx = 0; age_ndx < model_->age_spread(); ++age_ndx)
368372
exploitation_by_age_[age_ndx] += proportion_by_length_[year_ndx][length_ndx] * (numbers_at_age_and_length_[age_ndx][length_ndx] / sum_age);
369373
}
@@ -381,7 +385,8 @@ void TagByLength::DoExecute() {
381385
from_category_iter = 0;
382386

383387
for (; from_iter != from_partition_.end(); from_iter++, from_category_iter++) {
384-
for (unsigned age_ndx = 0; age_ndx < (*from_iter)->data_.size(); ++age_ndx) vulnerable_fish_by_age_[age_ndx] += (*from_iter)->data_[age_ndx];
388+
for (unsigned age_ndx = 0; age_ndx < (*from_iter)->data_.size(); ++age_ndx)
389+
vulnerable_fish_by_age_[age_ndx] += (*from_iter)->data_[age_ndx];
385390
}
386391
LOG_FINE() << "vulnerable to tag";
387392
// for (unsigned age_ndx = 0; age_ndx < model_->age_spread(); ++age_ndx)
@@ -551,23 +556,27 @@ void TagByLength::FillReportCache(ostringstream& cache) {
551556
for (unsigned category_ndx = 0; category_ndx < to_category_labels_.size(); ++category_ndx) {
552557
cache << "tags-after-initial_mortality-" << to_category_labels_[category_ndx] << " " << REPORT_R_DATAFRAME_ROW_LABELS << "\n";
553558
cache << "year ";
554-
for (unsigned age = min_age_; age <= max_age_; ++age) cache << age << " ";
559+
for (unsigned age = min_age_; age <= max_age_; ++age)
560+
cache << age << " ";
555561
cache << REPORT_EOL;
556562
for (unsigned year_ndx = 0; year_ndx < years_.size(); ++year_ndx) {
557563
cache << years_[year_ndx] << " ";
558-
for (unsigned age_ndx = 0; age_ndx < age_spread_; ++age_ndx) cache << AS_DOUBLE(tagged_fish_after_init_mort_[year_ndx][category_ndx][age_ndx]) << " ";
564+
for (unsigned age_ndx = 0; age_ndx < age_spread_; ++age_ndx)
565+
cache << AS_DOUBLE(tagged_fish_after_init_mort_[year_ndx][category_ndx][age_ndx]) << " ";
559566
cache << REPORT_EOL;
560567
}
561568
}
562569

563570
for (unsigned category_ndx = 0; category_ndx < to_category_labels_.size(); ++category_ndx) {
564571
cache << "tag-releases-" << to_category_labels_[category_ndx] << " " << REPORT_R_DATAFRAME_ROW_LABELS << "\n";
565572
cache << "year ";
566-
for (unsigned age = min_age_; age <= max_age_; ++age) cache << age << " ";
573+
for (unsigned age = min_age_; age <= max_age_; ++age)
574+
cache << age << " ";
567575
cache << REPORT_EOL;
568576
for (unsigned year_ndx = 0; year_ndx < years_.size(); ++year_ndx) {
569577
cache << years_[year_ndx] << " ";
570-
for (unsigned age_ndx = 0; age_ndx < age_spread_; ++age_ndx) cache << AS_DOUBLE(actual_tagged_fish_to_[year_ndx][category_ndx][age_ndx]) << " ";
578+
for (unsigned age_ndx = 0; age_ndx < age_spread_; ++age_ndx)
579+
cache << AS_DOUBLE(actual_tagged_fish_to_[year_ndx][category_ndx][age_ndx]) << " ";
571580
cache << REPORT_EOL;
572581
}
573582
}

0 commit comments

Comments
 (0)