test with wham

can't sort out what is not working; disable calls to sptlvb and recompile

NAMESPACE

@@ -7,3 +7,4 @@ export(confint_Deriv)
 export(get_Breaks)
 export(refit_Growth)
 import(mgcv)
+useDynLib(wham)

src/sptlVB_Sel_Sigma.cpp

@@ -0,0 +1,116 @@
+#define TMB_LIB_INIT R_init_growthbreaks
+#include <TMB.hpp>
+#include <iostream>
+// https://kaskr.github.io/adcomp/register_atomic_8cpp-example.html#a3
+// Similar to example 'adaptive_integration' using CppAD Romberg integration. REGISTER_ATOMIC is used to reduce tape size.
+// #include <TMB.hpp> // per https://github.com/mlysy/TMBtools this shouldn't be here
+template <class Type> Type selex(Type x){return 1/(1+exp( 52.976 - x)) ;}
+
+
+template<class Type>
+struct univariate {
+  Type Linf, k, Age, t0, Sigma; // variables in integrand,
+  Type operator() (Type u) {
+    Type yfit = 0;
+    Type selPred = 0;
+    Type ans = 0;
+    yfit = Linf*(1-exp(-k*(Age - t0))); // generate predicted length
+    selPred = selex(yfit); //1.0 for testing
+    // ans = selPred * exp(-(yfit - u)/(2*pow(Sigma*Age,2.0))); // u is what is plugged a and b
+    ans = selPred * dnorm(yfit,u,Sigma,false); // just hardcode the individuals
+    return ans;
+  }
+};
+
+template<class Type>
+vector<Type> mydenom(vector<Type> input) {
+  Type Linf  = input[0], k  = input[1];         // Data
+  Type Age = input[2],  t0 = input[3]; Type Sigma = input[4];         // Parameters
+  univariate<Type> f = {Linf, k, Age, t0, Sigma};
+  Type a = 0, b = 200; // max obs was 160
+  vector<Type> res(1); // a place to store this output
+  res[0] = romberg::integrate(f, a, b);
+  return res;
+}
+REGISTER_ATOMIC(mydenom)
+  template<class Type>
+  Type objective_function<Type>::operator() ()
+  {
+
+    DATA_VECTOR(Length_cm); // observed lengths
+    DATA_VECTOR(Age); // observed ages
+    DATA_VECTOR(Sel); // selectivity at observed length
+    DATA_IVECTOR(selType); // turns selectivity on/off for DFO data
+    DATA_IVECTOR(DES); // to organize data
+    DATA_INTEGER(nStrata);
+    DATA_INTEGER(a2); // fixed, for reporing L1/L2
+
+    PARAMETER_VECTOR(t0); // t0 can be negative
+
+    PARAMETER_VECTOR(log_Sigma);
+    vector<Type> Sigma(nStrata); //nstrata
+    Sigma = exp(log_Sigma);
+    // Type Sigma = exp(log_Sigma); // use exponentiated here
+
+    PARAMETER_VECTOR(log_Linf);
+    vector<Type> Linf(nStrata); //nstrata
+    Linf = exp(log_Linf);
+
+    PARAMETER_VECTOR(log_k);
+    vector<Type> k(nStrata);
+    k = exp(log_k);
+
+    // things to calculate
+    vector<Type> ypreds(Age.rows());
+    // Type selPred = 0.0;
+    Type selObs = 0.0;
+    Type yfit = 0;
+    Type ans = 0;
+    // Type tiny = 0.0; // Set to 1e-12 for robustness
+
+    Type numerator = 0.0; // will be the same for everyone
+    Type denominator = 0.0;
+
+    vector<Type> L1(nStrata);
+    vector<Type> L2(nStrata);
+
+    for(int i=0; i < Age.rows(); i++) {
+      yfit = Linf(DES(i))*(1-exp(-k(DES(i))*(Age(i) - t0(DES(i)))));
+
+
+      if(selType(i) != 2){
+        selObs = 1.0; // coerce to 1.0 for non DFO regions
+        numerator = selObs*dnorm(Length_cm(i),yfit,Sigma(DES(i)),false) ; // typical dnorm, sel obs is 1
+        denominator = 1.0;
+      }
+      if(selType(i) == 2){
+        selObs = selex(Length_cm(i)); // selex of OBSERVED value
+        numerator = selObs*dnorm(Length_cm(i),yfit,Sigma(DES(i)),false); // typical dnorm, sel obs is 1
+
+        vector<Type> input(5);
+        input << Linf(DES(i)), k(DES(i)), Age(DES(i)), t0(DES(i)), Sigma(DES(i)); // bundle the current estimates
+        denominator = mydenom(input)[0]; // calc denominator based on curr estimates
+
+        // ans -= log( mydenom(input)[0] + tiny );
+      }
+      ans -= log(numerator + 1e-5)-log(denominator + 1e-5); // will just be numerator for denom = 1
+      ypreds(i) = yfit; // store estimated length
+
+
+      L1(DES(i)) =  Linf(DES(i))*(1-exp(-k(DES(i))*(4.0 - t0(DES(i))))); // change 0, 0.5
+      L2(DES(i))  =  Linf(DES(i))*(1-exp(-k(DES(i))*(a2 - t0(DES(i)))));
+
+
+    } // end rows
+
+    REPORT(ypreds);
+    ADREPORT(ypreds);
+    REPORT(denominator)
+      ADREPORT(t0);
+    ADREPORT(k);
+    ADREPORT(Linf);
+    ADREPORT(L1);
+    ADREPORT(L2);
+    ADREPORT(Sigma);
+    return(ans);
+  }