[Refactor]: make use of join functions more consistent

[Rick-Methot-NOAA]

Refactor request

A user noticed that a forecast result did not match expectation output (see issue #485 ). This is due to the control rule's JOIN not being steep enough, which allows values >1.0 to occur near the transition

The logistic join function is used to smoothly connect two functions without using an IF statement. It is used in double normal selectivity, control rule inflection, hockey stick SRR, and other places.
If the join is not steep enough, then the resultant value of the join can be > 1.0 close to the transition point.

The function shown below is not used universally to achieve the join. It has a very steep transition (1000). Other custom join implementations use 10, 20, 30, 40, 50, 100. 10 is the value in the control rule join:

FUNCTION dvariable Join_Fxn(const prevariable& MinPoss, const prevariable& MaxPoss, const prevariable& Inflec, const prevariable& Xvar, const prevariable& Y1, const prevariable& Y2)
  {
  RETURN_ARRAYS_INCREMENT();
  dvariable Yresult;
  dvariable join;
  join = 1.000 / (1.000 + mfexp(1000.0 * (Xvar - Inflec) / (MaxPoss - MinPoss))); //  steep joiner at the inflection
  Yresult = Y1 * (join) + Y2 * (1.000 - join);
  RETURN_ARRAYS_DECREMENT();
  return Yresult;
  }

This should be cleaned up after an evaluation of a good default steepness (suggest 100).

• create test case to guide selection of optimal steepness for the joiner
• rewrite all joiner uses to use the function: Join_Fxn

Expected behavior

more consistent code

[Rick-Methot-NOAA]

@iantaylor-NOAA - here is the issue regarding joiners.

[Rick-Methot-NOAA]

`C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_benchfore.tpl
2001,5: join1 = 1. / (1. + mfexp(30. * (Fcast_Fmult - max_harvest_rate)));
2632,13: join1 = 1. / (1. + mfexp(10. * (SSB_current - H4010_bot * temp)));
2831,23: join1 = 1. / (1. + mfexp(30. * (temp1 - max_harvest_rate)));
2857,23: join1 = 1. / (1. + mfexp(30. * (temp1 - max_harvest_rate)));
2992,23: join1 = 1. / (1. + mfexp(30. * (temp - 0.95 * max_harvest_rate)));
3031,23: join1 = 1. / (1. + mfexp(30. * (temp - 0.95 * max_harvest_rate)));
3561,17: join1 = 1. / (1. + mfexp(1000. * (temp - 1.0))); // steep logistic joiner at adjustment of 1.0
3590,17: join1 = 1. / (1. + mfexp(1000. * (temp - 1.0))); // steep logistic joiner at adjustment of 1.0

C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_biofxn.tpl
621,23: join1 = 1.0 / (1.0 + mfexp(-(50. * t2 / (1.0 + fabs(t2))))); // note the logit transform is not perfect, so growth near Linf will not be exactly same as with native growth function
741,23: join1 = 1.0 / (1.0 + mfexp(-(50. * t2 / (1.0 + fabs(t2))))); // note the logit transform is not perfect, so growth near Linf will not be exactly same as with native growth function
871,17: join1 = 1.0 / (1.0 + mfexp(-(50. * t2 / (1.0 + fabs(t2))))); // note the logit transform is not perfect, so growth near Linf will not be exactly same as with native growth function
953,17: join1 = 1.0 / (1.0 + mfexp(-(50. * t2 / (1.0 + fabs(t2))))); // note the logit transform is not perfect, so growth near Linf will not be exactly same as with native growth function

C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_popdyn.tpl
1141,19: join1 = 1. / (1. + mfexp(30. * (temp - 0.95))); // steep logistic joiner at harvest rate of 0.95
1222,25: join1 = 1. / (1. + mfexp(30. * (temp - 0.95 * max_harvest_rate)));
1896,21: join1 = 1. / (1. + mfexp(40.0 * (crashtemp - max_harvest_rate))); // steep joiner logistic curve at limit

C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_selex.tpl
177,19: join1 = 1.0 / (1.0 + mfexp(-(20. * t1 / (1.0 + fabs(t1))))); // note the logit transform on t1 causes range of mfexp to be over -20 to 20
378,21: join1 = 1. / (1. + mfexp(10. * (len_bins_m(j) - sp(1))));
533,19: join1 = 1.0 / (1.0 + mfexp(-(20. * t1 / (1.0 + fabs(t1))))); // note the logit transform on t1 causes range of mfexp to be over -20 to 20
1460,19: join1 = 1. / (1. + mfexp(-(20. / (1. + fabs(t1))) * t1));

JOIN2
C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_benchfore.tpl
2633,13: join2 = 1. / (1. + mfexp(10. * (SSB_current - H4010_top * temp)));

C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_selex.tpl
178,19: join2 = 1.0 / (1.0 + mfexp(-(20. * t2 / (1.0 + fabs(t2)))));
379,21: join2 = 1. / (1. + mfexp(10. * (len_bins_m(j) - (sp(1) + sp(8)))));
534,19: join2 = 1.0 / (1.0 + mfexp(-(20. * t2 / (1.0 + fabs(t2)))));
1461,19: join2 = 1. / (1. + mfexp(-(20. / (1. + fabs(t2))) * t2));

JOIN3
C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_selex.tpl
380,21: join3 = 1. / (1. + mfexp(10. * (len_bins_m(j) - sel_maxL)));

JOIN
C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_miscfxn.tpl
14,3: join = 1.000 / (1.000 + mfexp(1000.0 * (Xvar - Inflec) / (MaxPoss - MinPoss))); // steep joiner at the inflection

C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_recruit.tpl
138,7: join = 1. / (1. + mfexp(100 * (SRZ_surv - 1.)));

C:\Users\Richard.Methot\Documents\GitHub\StockSynthesis_git\stock-synthesis\SS_readcontrol_330.tpl
2544,9: Equ_F_joiner = 10; // defaults
2575,11: Equ_F_joiner = (log(1. / max_harvest_rate - 1.)) / (max_harvest_rate - 0.2); // used to spline the harvest rate

`

[iantaylor-NOAA]

@Rick-Methot-NOAA, thanks for the additional detail on this issue and sorry for failing to remember this issue and #485 earlier.

I think one of the challenging in making a consistent join function is that the units of the dependent variable being joined could differ, from small fractions representing harvest rates or fraction unfished to larger values associated with fish length used in selectivity functions.

A large number like 1000 would presumably be more than steep enough for all these cases, but it would be good to see if large numbers cause an impact in parameter estimation by making changes in the likelihood surface more abrupt.

[Rick-Methot-NOAA]

I think I will start by making all that are <50 be equal to 50. I noticed that 30 was enough to fix the petrale case.
Sorry for not getting to this a year ago when the petrale situation was first noticed.

[Rick-Methot-NOAA]

I created constant joinsteep so I could recompile with a range of options. No obvious sensitivity of overall petrale model results to values of 20, 50, or 200. Next need to compare variance.
effect_of_joinsteep.txt

[Rick-Methot-NOAA]

attached find ss_summary from base (3.30.22.1) and from joinsteep=20 and joinsteep=200 and joinsteep=50

I conclude that joinsteep 200 is too extreme and would like to settle on 50. Uncomfortable that all users will see some small numerical differences.
ss_summary_joinsteep.zip