gadget-framework / gadget3

TMB-based gadget implemtation
GNU General Public License v2.0
8 stars 5 forks source link

Weightloss action #174

Closed lentinj closed 1 month ago

lentinj commented 2 months ago

Besides pupping / spawning weightloss, GoB ringed seals has a ~constant seasonal weightloss from moulting, etc.

Mouting weightloss will likely be proportional

vbartolino commented 1 month ago

Do I understand correctly that now we can have weight loss both as part of the spawning and as an independent process? I've started to work on the implementation in the ringed seal model and in brief I'm proceeding with the following: timestep1 -> spawning incl an absolute wgt loss timestep2 -> additional physiological wgt loss timestep3-4 -> feeding with regain of wgt

Blubber thickness data, which inform about the seasonal weight changes, show a minimum in the step2. I'm not sure about the order of calculation, would the weight loss in step2 happen before the data fitting in step2?

At the moment g3a_weightloss and g3l_sparsesample to include blubber thickness data are in separate branches, right? Would it make sense to merge them?

lentinj commented 1 month ago

Blubber thickness data, which inform about the seasonal weight changes, show a minimum in the step2. I'm not sure about the order of calculation, would the weight loss in step2 happen before the data fitting in step2?

We regard g3a_weightloss as naturalmortality for the sake of gadget order of actions (you can see that in the function's defaults https://github.com/gadget-framework/gadget3/blob/master/man/action_weightloss.Rd#L19).

Spawning weightloss we append to the spawn stage in order of actions

So, referring to https://github.com/gadget-framework/gadget3/blob/master/R/action_order.R the order within a step should be:

Is this workable, or do we need to re-arrange some of this?

At the moment g3a_weightloss and g3l_sparsesample to include blubber thickness data are in separate branches, right?

I've merged g3a_weightloss into master, but not g3l_sparsesample. I'll tidy this up now and let you know.

vbartolino commented 1 month ago

Is this workable, or do we need to re-arrange some of this?

it seems good to me, also thanks to point to the links

At the moment, the argument min_weight in g3a_weightlossis useful so we can have some control and bound the weight loss. To make full sense of this, I'd like to have a min_weight by length. Does it make sense to use a formula like below?

    # ADDITIONAL PHYSIOLOGICAL WEIGHT LOSS 
    g3a_weightloss(mat_stock,
             # 10% of body weight physiological loss to the breeding and lactation costs (i.e., moulting)
             rel_loss = g3_parameterized("moult.weightloss", value = 0.1),
             min_weight = g3_formula(
                 wmin.a * stock__midlen^wmin.b,
                 wmin.a = g3_parameterized("wmin.a", by_stock = TRUE),
                 wmin.b = g3_parameterized("wmin.b", by_stock = TRUE),
                 end = NULL ),
             run_step = 2),
lentinj commented 1 month ago

I've rebased likelihood-individual on top of master, so now if you re-install that branch you should have all you need.

It needs some more work before I can merge it into master, but that shouldn't stop you if you're happy with the interface.

lentinj commented 1 month ago

Does it make sense to use a formula like below?

Looks good to me. I'll add it to the examples & make sure it works.

vbartolino commented 1 month ago

I'm having problem with

    g3a_weightloss(mat_stock,
             abs_loss = g3_parameterized("spawn.weightloss", value = 15),
             ...)

but not with

    g3a_weightloss(mat_stock,
             rel_loss = g3_parameterized("moult.weightloss", value = 0.1),
             ...)

this is the error message

> > obj.fun <- gadget3::g3_tmb_adfun(tmb_model, tmb_param)
g++ -std=gnu++14 -I"/usr/local/lib/R412/lib/R/include" -DNDEBUG -I"/usr/local/lib/R412/lib/R/library/TMB/include" -I"/usr/local/lib/R412/lib/R/library/RcppEigen/include"  -DTMB_SAFEBOUNDS -DTMB_EIGEN_DISABLE_WARNINGS -DLIB_UNLOAD=R_unload_g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b  -DTMB_LIB_INIT=R_init_g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b  -DTMBAD_FRAMEWORK  -I/usr/local/include   -fpic  -std=gnu++11 -Wno-ignored-attributes -DEIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS -O3 -flto=auto -march=native -c /tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp -o /tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.o
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp: In instantiation of ‘Type objective_function<Type>::operator()() [with Type = TMBad::global::ad_aug]’:
/usr/local/lib/R412/lib/R/library/TMB/include/tmb_core.hpp:1277:7:   required from here
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:1231:97: error: no match for call to ‘(objective_function<Type>::operator()() [with Type = TMBad::global::ad_aug]::<lambda(tmbutils::vector<TMBad::global::ad_aug>, TMBad::global::ad_aug)>) (tmbutils::array<TMBad::global::ad_aug>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<TMBad::global::ad_aug, TMBad::global::ad_aug>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<TMBad::global::ad_aug>, const Eigen::Array<TMBad::global::ad_aug, -1, 1, 0, -1, 1> >, const Eigen::CwiseBinaryOp<Eigen::internal::scalar_pow_op<TMBad::global::ad_aug, TMBad::global::ad_aug>, const Eigen::Array<TMBad::global::ad_aug, -1, 1, 0, -1, 1>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<TMBad::global::ad_aug>, const Eigen::Array<TMBad::global::ad_aug, -1, 1, 0, -1, 1> > > >&)’
                     rin_mat__wgt.col(rin_mat__age_idx).col(rin_mat__area_idx) = logspace_add_vec((rin_mat__wgt.col(rin_mat__age_idx).col(rin_mat__area_idx) - abs_loss)*(double)(1e+07), min_weight) / (double)(1e+07);
                                                                                 ~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:1231:97: note: candidate: tmbutils::vector<TMBad::global::ad_aug> (*)(tmbutils::vector<TMBad::global::ad_aug>, TMBad::global::ad_aug) <conversion>
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:1231:97: note:   candidate expects 3 arguments, 3 provided
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:484:59: note: candidate: objective_function<Type>::operator()()::<lambda(tmbutils::vector<Type>, Type)> [with Type = TMBad::global::ad_aug]
     auto logspace_add_vec = [](vector<Type> a, Type b) -> vector<Type> {
                                                           ^~~~~~~~~~~~
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:484:59: note:   no known conversion for argument 2 from ‘Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<TMBad::global::ad_aug, TMBad::global::ad_aug>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<TMBad::global::ad_aug>, const Eigen::Array<TMBad::global::ad_aug, -1, 1, 0, -1, 1> >, const Eigen::CwiseBinaryOp<Eigen::internal::scalar_pow_op<TMBad::global::ad_aug, TMBad::global::ad_aug>, const Eigen::Array<TMBad::global::ad_aug, -1, 1, 0, -1, 1>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<TMBad::global::ad_aug>, const Eigen::Array<TMBad::global::ad_aug, -1, 1, 0, -1, 1> > > >’ to ‘TMBad::global::ad_aug’
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp: In instantiation of ‘Type objective_function<Type>::operator()() [with Type = double]’:
/usr/local/lib/R412/lib/R/library/TMB/include/tmb_core.hpp:2031:7:   required from here
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:1231:97: error: no match for call to ‘(objective_function<Type>::operator()() [with Type = double]::<lambda(tmbutils::vector<double>, double)>) (tmbutils::array<double>, Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, const Eigen::Array<double, -1, 1> >, const Eigen::CwiseBinaryOp<Eigen::internal::scalar_pow_op<double, double>, const Eigen::Array<double, -1, 1>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, const Eigen::Array<double, -1, 1> > > >&)’
                     rin_mat__wgt.col(rin_mat__age_idx).col(rin_mat__area_idx) = logspace_add_vec((rin_mat__wgt.col(rin_mat__age_idx).col(rin_mat__area_idx) - abs_loss)*(double)(1e+07), min_weight) / (double)(1e+07);
                                                                                 ~~~~~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:1231:97: note: candidate: tmbutils::vector<double> (*)(tmbutils::vector<double>, double) <conversion>
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:1231:97: note:   candidate expects 3 arguments, 3 provided
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:484:59: note: candidate: objective_function<Type>::operator()()::<lambda(tmbutils::vector<Type>, Type)> [with Type = double]
     auto logspace_add_vec = [](vector<Type> a, Type b) -> vector<Type> {
                                                           ^~~~~~~~~~~~
/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.cpp:484:59: note:   no known conversion for argument 2 from ‘Eigen::CwiseBinaryOp<Eigen::internal::scalar_product_op<double, double>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, const Eigen::Array<double, -1, 1> >, const Eigen::CwiseBinaryOp<Eigen::internal::scalar_pow_op<double, double>, const Eigen::Array<double, -1, 1>, const Eigen::CwiseNullaryOp<Eigen::internal::scalar_constant_op<double>, const Eigen::Array<double, -1, 1> > > >’ to ‘double’
/usr/local/lib/R412/lib/R/etc/Makeconf:177: recipe for target '/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.o' failed
make: *** [/tmp/RtmptO7Vnj/g30_12_1_999_tmb1_9_14_42412448b477a81f67c7878f37afa59f583c807b.o] Error 1
Error in (function (file, flags = "", safebounds = TRUE, safeunload = TRUE,  : 
  Compilation failed
lentinj commented 1 month ago

Bah, I only tried that in the examples, which don't compile to TMB.

Do you need a quick fix, or would you be okay with having a single min_weight value for a day? Ideally I need to re-jig logspace_add_vec so we don't keep needing infinite variations of it.

vbartolino commented 1 month ago

Do you need a quick fix, or would you be okay with having a single min_weight value for a day? Ideally I need to re-jig logspace_add_vec so we don't keep needing infinite variations of it.

no need of quick fix, take your time

lentinj commented 1 month ago

Right, I think I've (finally) solved this. Short version: Your code with length-dependent min_length should work as above.

Long version: The problem is that we had different versions for logspace_add and similar functions, depending if the provided variables were scalar, vector & array. The combinatorial explosion of functions required was getting out of hand, and meant you had to choose which to use (thus the problem here, I'd expected min_weight to be scalar, and used the wrong function for your case). I've added dif_pmax, dif_pmin, dif_pminmax. These use C++ templates so the _vec and _arr variations aren't necessary. I've also re-written avoid_zero to use them, and removed avoid_zero_vec.

Others could be re-written too, but I've not got there yet.

vbartolino commented 1 month ago

Right, I think I've (finally) solved this. Short version: Your code with length-dependent min_length should work as above.

I haven't checked the actual outcomes of an optimisation, but it runs both abs_loss in Q1 and rel_lossin Q2 with no complain, thanks!