Merge pull request #91 from James-Thorson/development

Includes V4.1.0 which adds feature to compute the GMRF normalization constant outside of the inner optimizer.

R/Build_TMB_Fn.R

@@ -103,7 +103,13 @@ function( TmbData, Version, Q_Config=TRUE, CovConfig=TRUE,
   # Build object
   dyn.load( paste0(RunDir,"/",TMB::dynlib(Version)) ) # random=Random,
   Obj <- MakeADFun(data=TmbData, parameters=Parameters, hessian=FALSE, map=Map, random=Random, inner.method="newton", DLL=Version)  #
-  Obj$control <- list(trace=1, parscale=1, REPORT=1, reltol=1e-12, maxit=100)
+  Obj$control <- list(parscale=1, REPORT=1, reltol=1e-12, maxit=100)
+
+  # Add normalization in
+  if(Version %in% c("VAST_v4_1_0") & TmbData$Options['normalize_GMRF_in_CPP']==FALSE ){
+    message("Normalizing GMRF in R using `TMB::normalize` feature")
+    Obj = TMB::normalize(Obj, flag="include_data", value=FALSE)
+  }
 
   # Diagnostic functions (optional)
   if( !is.null(DiagnosticDir) ){

R/Data_Fn.R

@@ -54,8 +54,8 @@ function( Version, FieldConfig, OverdispersionConfig=c("eta1"=0,"eta2"=0), ObsMo
   Options=c('SD_site_logdensity'=0,'Calculate_Range'=0,'Calculate_effective_area'=0,'Calculate_Cov_SE'=0,'Calculate_Synchrony'=0,'Calculate_proportion'=0) ){
 
   # Specify default values for `Options`
-  Options2use = c('SD_site_density'=0,'SD_site_logdensity'=0,'Calculate_Range'=0,'Calculate_evenness'=0,'Calculate_effective_area'=0,
-    'Calculate_Cov_SE'=0,'Calculate_Synchrony'=0,'Calculate_Coherence'=0,'Calculate_proportion'=0)
+  Options2use = c('SD_site_density'=0, 'SD_site_logdensity'=0, 'Calculate_Range'=0, 'Calculate_evenness'=0, 'Calculate_effective_area'=0,
+    'Calculate_Cov_SE'=0, 'Calculate_Synchrony'=0, 'Calculate_Coherence'=0, 'Calculate_proportion'=0, 'normalize_GMRF_in_CPP'=TRUE)
 
   # Replace defaults for `Options` with provided values (if any)
   for( i in 1:length(Options)){
@@ -192,6 +192,11 @@ function( Version, FieldConfig, OverdispersionConfig=c("eta1"=0,"eta2"=0), ObsMo
     }
   }
 
+  # Check for incompatible settings
+  if( ncol(t_iz)>=2 & ( RhoConfig[["Beta1"]]!=0 | RhoConfig[["Beta2"]]!=0 ) ){
+    stop("Temporal structure on intercepts is not implemented for seasonal models")
+  }
+
   # switch defaults if necessary
   if( Method=="Grid" ){
     Aniso = 0
@@ -232,6 +237,9 @@ function( Version, FieldConfig, OverdispersionConfig=c("eta1"=0,"eta2"=0), ObsMo
   if(Version%in%c("VAST_v4_0_0")){
     Return = list( "n_i"=n_i, "n_s"=c(MeshList$anisotropic_spde$n.spde,n_x)[Options_vec['Method']+1], "n_x"=n_x, "n_t"=n_t, "n_c"=n_c, "n_e"=n_e, "n_j"=n_j, "n_p"=n_p, "n_k"=n_k, "n_v"=n_v, "n_l"=n_l, "n_m"=ncol(Z_xm), "Options_vec"=Options_vec, "FieldConfig"=FieldConfig_input, "OverdispersionConfig"=OverdispersionConfig_input, "ObsModel_ez"=ObsModel_ez, "Options"=Options2use, "yearbounds_zz"=yearbounds_zz, "b_i"=b_i, "a_i"=a_i, "c_iz"=c_iz, "e_i"=e_i, "s_i"=s_i, "t_iz"=t_iz-min(t_iz,na.rm=TRUE), "v_i"=match(v_i,sort(unique(v_i)))-1, "PredTF_i"=PredTF_i, "a_xl"=a_xl, "X_xj"=X_xj, "X_xtp"=X_xtp, "Q_ik"=Q_ik, "t_yz"=t_yz, "Z_xm"=Z_xm, "spde"=list(), "spde_aniso"=list(), "M0"=GridList$M0, "M1"=GridList$M1, "M2"=GridList$M2 )
   }
+  if(Version%in%c("VAST_v4_1_0")){
+    Return = list( "n_i"=n_i, "n_s"=c(MeshList$anisotropic_spde$n.spde,n_x)[Options_vec['Method']+1], "n_x"=n_x, "n_t"=n_t, "n_c"=n_c, "n_e"=n_e, "n_p"=n_p, "n_v"=n_v, "n_l"=n_l, "n_m"=ncol(Z_xm), "Options_vec"=Options_vec, "FieldConfig"=FieldConfig_input, "OverdispersionConfig"=OverdispersionConfig_input, "ObsModel_ez"=ObsModel_ez, "include_data"=TRUE, "Options"=Options2use, "yearbounds_zz"=yearbounds_zz, "b_i"=b_i, "a_i"=a_i, "c_iz"=c_iz, "e_i"=e_i, "s_i"=s_i, "t_iz"=t_iz-min(t_iz,na.rm=TRUE), "v_i"=match(v_i,sort(unique(v_i)))-1, "PredTF_i"=PredTF_i, "a_xl"=a_xl, "X_xj"=X_xj, "X_xtp"=X_xtp, "Q_ik"=Q_ik, "t_yz"=t_yz, "Z_xm"=Z_xm, "spde"=list(), "spde_aniso"=list(), "M0"=GridList$M0, "M1"=GridList$M1, "M2"=GridList$M2 )
+  }
   if( "spde" %in% names(Return) ) Return[['spde']] = MeshList$isotropic_spde$param.inla[c("M0","M1","M2")]
   if( "spde_aniso" %in% names(Return) ) Return[['spde_aniso']] = list("n_s"=MeshList$anisotropic_spde$n.spde, "n_tri"=nrow(MeshList$anisotropic_mesh$graph$tv), "Tri_Area"=MeshList$Tri_Area, "E0"=MeshList$E0, "E1"=MeshList$E1, "E2"=MeshList$E2, "TV"=MeshList$TV-1, "G0"=MeshList$anisotropic_spde$param.inla$M0, "G0_inv"=INLA::inla.as.dgTMatrix(solve(MeshList$anisotropic_spde$param.inla$M0)) )
 

R/Make_Map.R

@@ -85,8 +85,8 @@ function( DataList, TmbParams, CovConfig=TRUE, DynCovConfig=TRUE, Q_Config=TRUE,
     Map[["delta_i"]] = factor(Map[["delta_i"]])
   }
 
-  # Change beta1_ct if
-  if( any(DataList$ObsModel_ez[,2]%in%c(3)) ){
+  # Change beta1_ct if 100% encounters (not designed to work with seasonal models
+  if( any(DataList$ObsModel_ez[,2]%in%c(3)) & ncol(DataList$t_iz)==1 ){
     Tmp_ct = tapply(ifelse(DataList$b_i>0,1,0), INDEX=list(factor(DataList$c_iz[,1],levels=sort(unique(DataList$c_iz[,1]))),DataList$t_iz[,1]), FUN=mean)
     Map[["beta1_ct"]] = array( 1:prod(dim(Tmp_ct)), dim=dim(Tmp_ct) )
     Map[["beta1_ct"]][which(is.na(Tmp_ct) | Tmp_ct==1)] = NA
@@ -150,8 +150,12 @@ function( DataList, TmbParams, CovConfig=TRUE, DynCovConfig=TRUE, Q_Config=TRUE,
   }
   # fix betas and/or epsilons for missing years if betas are fixed-effects
   #YearNotInData = !( (1:DataList$n_t) %in% (unique(DataList$t_i)+1) )
-  Num_ct = tapply( DataList$b_i, INDEX=list(factor(DataList$c_iz[,1],levels=1:DataList$n_c-1),factor(DataList$t_iz[,1],levels=1:DataList$n_t-1)), FUN=function(vec){sum(!is.na(vec))} )
-  Num_ct = ifelse( is.na(Num_ct), 0, Num_ct )
+  Num_ct = array(0, dim=c(DataList$n_c,DataList$n_t))
+  for( tz in 1:ncol(DataList$t_iz) ){
+    Temp = tapply( DataList$b_i, INDEX=list(factor(DataList$c_iz[,1],levels=1:DataList$n_c-1),factor(DataList$t_iz[,tz],levels=1:DataList$n_t-1)), FUN=function(vec){sum(!is.na(vec))} )
+    Temp = ifelse( is.na(Temp), 0, Temp )
+    Num_ct = Num_ct + Temp
+  }
   if( sum(Num_ct==0)>0 ){
     # Beta1 -- Fixed
     if( RhoConfig["Beta1"]==0){
@@ -193,15 +197,20 @@ function( DataList, TmbParams, CovConfig=TRUE, DynCovConfig=TRUE, Q_Config=TRUE,
     # (Re)start map for intercepts
     if( !("beta1_ct" %in% names(Map)) ){
       Map[["beta1_ct"]] = 1:prod(dim(TmbParams[["beta1_ct"]]))
-    }else{ Map[["beta1_ct"]] = as.numeric(Map[["beta1_ct"]]) }
+    }else{
+      Map[["beta1_ct"]] = as.numeric(Map[["beta1_ct"]])
+    }
     if( !("beta2_ct" %in% names(Map)) ){
       Map[["beta2_ct"]] = 1:prod(dim(TmbParams[["beta2_ct"]]))
-    }else{ Map[["beta2_ct"]] = as.numeric(Map[["beta2_ct"]]) }
+    }else{
+      Map[["beta2_ct"]] = as.numeric(Map[["beta2_ct"]])
+    }
     # Add fixed values for lowest value of 2nd and higher columns
     for( zI in 2:ncol(DataList$t_iz) ){
       Which2Fix = min( DataList$t_iz[,zI] )
-      Map[["beta1_ct"]][Which2Fix+1] = NA
-      Map[["beta2_ct"]][Which2Fix+1] = NA
+      Which2Fix = matrix( 1:(TmbData$n_c*TmbData$n_t), ncol=TmbData$n_t, nrow=TmbData$n_c )[,Which2Fix+1]
+      Map[["beta1_ct"]][Which2Fix] = NA
+      Map[["beta2_ct"]][Which2Fix] = NA
     }
     # Remake as factor
     Map[["beta1_ct"]] = factor(Map[["beta1_ct"]])
@@ -228,7 +237,7 @@ function( DataList, TmbParams, CovConfig=TRUE, DynCovConfig=TRUE, Q_Config=TRUE,
 
   # Static covariates
   Var_j = apply( DataList[["X_xj"]], MARGIN=2, FUN=var )
-  Map[["gamma1_j"]] = Map[["gamma2_j"]] = 1:DataList$n_j
+  Map[["gamma1_j"]] = Map[["gamma2_j"]] = 1:ncol(DataList$X_xj)
   for(j in 1:length(Var_j)){
     if( Var_j[j]==0 || sum(CovConfig)==0 ){
       Map[["gamma1_j"]][j] = NA
@@ -240,7 +249,7 @@ function( DataList, TmbParams, CovConfig=TRUE, DynCovConfig=TRUE, Q_Config=TRUE,
 
   # Catchability variables
   Var_k = apply( DataList[["Q_ik"]], MARGIN=2, FUN=var )
-  Map[["lambda1_k"]] = Map[["lambda2_k"]] = 1:DataList$n_k
+  Map[["lambda1_k"]] = Map[["lambda2_k"]] = 1:ncol(DataList$Q_ik)
   for(k in 1:length(Var_k)){
     if( Var_k[k]==0 || sum(Q_Config)==0 ){
       Map[["lambda1_k"]][k] = NA