From 4a8a78da334b06d6f8cb1647d006e0790f4da68a Mon Sep 17 00:00:00 2001 From: Jim Thorson Date: Mon, 12 Nov 2018 08:19:22 -0800 Subject: [PATCH] Fixed bug in V5.3.0... ... which caused R to crash when FieldConfig[1] != FieldConfig[2] due to wrong initialization of length of Epsilon_rho1_f --- DESCRIPTION | 4 +- R/Param_Fn.R | 4 +- inst/extdata/EBS_3species/VAST_v5_2_0.cpp | 1443 --------------------- 3 files changed, 4 insertions(+), 1447 deletions(-) delete mode 100644 inst/extdata/EBS_3species/VAST_v5_2_0.cpp diff --git a/DESCRIPTION b/DESCRIPTION index eb15f51..422b7c8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: VAST Type: Package Title: Vector-autoregressive spatio-temporal (VAST) model -Version: 2.0.0 -Date: 2018-10-31 +Version: 2.0.1 +Date: 2018-11-12 Authors@R: person("James","Thorson", email="James.Thorson@noaa.gov", role=c("aut","cre")) Maintainer: James Thorson Description: VAST is an R package for conducting spatio-temporal analysis diff --git a/R/Param_Fn.R b/R/Param_Fn.R index fdac644..06f1c88 100644 --- a/R/Param_Fn.R +++ b/R/Param_Fn.R @@ -119,8 +119,8 @@ function( Version, DataList, RhoConfig=c("Beta1"=0,"Beta2"=0,"Epsilon1"=0,"Epsil if( "L_omega2_z" %in% names(Return)) Return = Add_factor( List=Return, n_c=DataList$n_c, n_f=DataList$FieldConfig[3], n_i=DataList$n_s, list_names=c("L_omega2_z","Omegainput2_sf"), sd=0 ) if( "L_epsilon2_z" %in% names(Return)) Return = Add_factor( List=Return, n_c=DataList$n_c, n_f=DataList$FieldConfig[4], n_i=DataList$n_s, n_t=DataList$n_t, list_names=c("L_epsilon2_z","Epsiloninput2_sft"), sd=0 ) # Autocorrelation - if( "Epsilon_rho1_f" %in% names(Return)) Return[["Epsilon_rho1_f"]] = rep(0, ncol(Return[["Omegainput1_sf"]])) - if( "Epsilon_rho2_f" %in% names(Return)) Return[["Epsilon_rho2_f"]] = rep(0, ncol(Return[["Omegainput2_sf"]])) + if( "Epsilon_rho1_f" %in% names(Return)) Return[["Epsilon_rho1_f"]] = rep(0, dim(Return[["Epsiloninput1_sft"]])[2]) + if( "Epsilon_rho2_f" %in% names(Return)) Return[["Epsilon_rho2_f"]] = rep(0, dim(Return[["Epsiloninput2_sft"]])[2]) # Initial values if( all(DataList$ObsModel_ez[,2]%in%c(0,3)) ){ Return[["beta1_ct"]] = qlogis(0.01*0.99*tapply(ifelse(DataList$b_i>0,1,0),INDEX=factor(DataList$c_iz[,1],levels=sort(unique(DataList$c_iz[,1]))),FUN=mean)) %o% rep(1,DataList$n_t) diff --git a/inst/extdata/EBS_3species/VAST_v5_2_0.cpp b/inst/extdata/EBS_3species/VAST_v5_2_0.cpp deleted file mode 100644 index 112d570..0000000 --- a/inst/extdata/EBS_3species/VAST_v5_2_0.cpp +++ /dev/null @@ -1,1443 +0,0 @@ -#include -#include - -// Needed for returning SparseMatrix -template -Eigen::SparseMatrix Q_network( Type log_theta, int n_s, vector parent_s, vector child_s, vector dist_s ){ - Eigen::SparseMatrix Q( n_s, n_s ); - Type theta = exp( log_theta ); - for(int s=0; s -bool isNA(Type x){ - return R_IsNA(asDouble(x)); -} - -// Posfun -template -Type posfun(Type x, Type lowerlimit, Type &pen){ - pen += CppAD::CondExpLt(x,lowerlimit,Type(0.01)*pow(x-lowerlimit,2),Type(0)); - return CppAD::CondExpGe(x,lowerlimit,x,lowerlimit/(Type(2)-x/lowerlimit)); -} - -// Variance -template -Type var( array vec ){ - Type vec_mod = vec - (vec.sum()/vec.size()); - Type res = pow(vec_mod, 2).sum() / vec.size(); - return res; -} - -// dlnorm -template -Type dlnorm(Type x, Type meanlog, Type sdlog, int give_log=0){ - //return 1/(sqrt(2*M_PI)*sd)*exp(-.5*pow((x-mean)/sd,2)); - Type logres = dnorm( log(x), meanlog, sdlog, true) - log(x); - if(give_log) return logres; else return exp(logres); -} - -// Generate loadings matrix -template -matrix loadings_matrix( vector L_val, int n_rows, int n_cols ){ - matrix L_rc(n_rows, n_cols); - int Count = 0; - for(int r=0; r=c){ - L_rc(r,c) = L_val(Count); - Count++; - }else{ - L_rc(r,c) = 0.0; - } - }} - return L_rc; -} - -// IN: eta1_vf; L1_z -// OUT: jnll_comp; eta1_vc -template -matrix overdispersion_by_category_nll( int n_f, int n_v, int n_c, matrix eta_vf, vector L_z, Type &jnll_pointer, objective_function* of){ - using namespace density; - matrix eta_vc(n_v, n_c); - vector Tmp_c; - // Turn off - if(n_f<0){ - eta_vc.setZero(); - } - // AR1 structure - if( n_f==0 ){ - for(int v=0; v::value && of->do_simulate){ - SCALE( AR1(L_z(1)), exp(L_z(0)) ).simulate(Tmp_c); - eta_vf.row(v) = Tmp_c; - } - } - eta_vc = eta_vf; - } - // Factor analysis structure - if( n_f>0 ){ - // Assemble the loadings matrix - matrix L_cf = loadings_matrix( L_z, n_c, n_f ); - // Probability of overdispersion - for(int v=0; v::value && of->do_simulate){ - eta_vf(v,f) = rnorm( Type(0.0), Type(1.0) ); - } - }} - // Multiply out overdispersion - eta_vc = eta_vf * L_cf.transpose(); - } - return eta_vc; -} - -template // -matrix convert_upper_cov_to_cor( matrix cov ){ - int nrow = cov.row(0).size(); - for( int i=0; i // -matrix gmrf_by_category_nll( int n_f, int method, int timing, int n_s, int n_c, Type logkappa, array gmrf_input_sf, array gmrf_mean_sf, vector L_z, density::GMRF_t gmrf_Q, Type &jnll_pointer, objective_function* of){ - using namespace density; - matrix gmrf_sc(n_s, n_c); - vector gmrf_s(n_s); - matrix Cov_cc(n_c,n_c); - array diff_gmrf_sc(n_s, n_c); // Requires an array - Type logtau; - if(method==0) logtau = log( 1 / (exp(logkappa) * sqrt(4*M_PI)) ); - if(method==1) logtau = log( 1 / sqrt(1-exp(logkappa*2)) ); - if(method==2) logtau = Type(0.0); - // IID - if(n_f == -2){ - for( int c=0; c::value && of->do_simulate) { - gmrf_Q.simulate(gmrf_s); - gmrf_input_sf.col(c) = gmrf_s + gmrf_mean_sf.col(c); - } - // Rescale - gmrf_sc.col(c) = gmrf_input_sf.col(c) / exp(logtau) * L_z(c); // Rescaling from comp_index_v1d.cpp - } - } - // Turn off - if(n_f == -1){ - gmrf_sc.setZero(); - } - // AR1 structure - if(n_f==0){ - jnll_pointer += SEPARABLE( AR1(L_z(1)), gmrf_Q )(gmrf_input_sf - gmrf_mean_sf); - // Simulate new values when using obj.simulate() - if(isDouble::value && of->do_simulate) { - SEPARABLE( AR1(L_z(1)), gmrf_Q ).simulate(gmrf_input_sf); - gmrf_input_sf += gmrf_input_sf; - } - // Rescale - logtau = L_z(0) - logkappa; // - gmrf_sc = gmrf_input_sf / exp(logtau); // Rescaling from comp_index_v1d.cpp - } - // Factor analysis structure - if(n_f>0){ - // PDF if density-dependence/interactions occurs prior to correlated dynamics - if( timing==0 ){ - for( int f=0; f::value && of->do_simulate) { - gmrf_Q.simulate(gmrf_s); - gmrf_input_sf.col(f) = gmrf_s + gmrf_mean_sf.col(f); - } - } - // Rescale - matrix L_cf = loadings_matrix( L_z, n_c, n_f ); - gmrf_sc = (gmrf_input_sf.matrix() * L_cf.transpose()) / exp(logtau); - } - // PDF if density-dependence/interactions occurs after correlated dynamics (Only makes sense if n_f == n_c) - if( timing==1 ){ - // Calculate difference without rescaling - gmrf_sc = gmrf_input_sf.matrix(); - for( int s=0; s L_cf = loadings_matrix( L_z, n_c, n_f ); - Cov_cc = L_cf * L_cf.transpose(); - jnll_pointer += SCALE(SEPARABLE(MVNORM(Cov_cc), gmrf_Q), exp(-logtau))( diff_gmrf_sc ); - //gmrf_sc = gmrf_sc / exp(logtau); - // Simulate new values when using obj.simulate() - if(isDouble::value && of->do_simulate) { - SEPARABLE(MVNORM(Cov_cc), gmrf_Q).simulate( diff_gmrf_sc ); - gmrf_sc = gmrf_mean_sf + diff_gmrf_sc/exp(logtau); - } - } - } - return gmrf_sc; -} - -// CMP distribution -template -Type dCMP(Type x, Type mu, Type nu, int give_log=0, int iter_max=30, int break_point=10){ - // Explicit - Type ln_S_1 = nu*mu - ((nu-1)/2)*log(mu) - ((nu-1)/2)*log(2*M_PI) - 0.5*log(nu); - // Recursive - vector S_i(iter_max); - S_i(0) = 1; - for(int i=1; i -Type dPoisGam( Type x, Type shape, Type scale, Type intensity, vector &diag_z, int maxsum=50, int minsum=1, int give_log=0 ){ - // Maximum integration constant to prevent numerical overflow, but capped at value for maxsum to prevent numerical underflow when subtracting by a higher limit than is seen in the sequence - Type max_log_wJ, z1, maxJ_bounded; - if( x==0 ){ - diag_z(0) = 1; - max_log_wJ = 0; - diag_z(1) = 0; - }else{ - z1 = log(intensity) + shape*log(x/scale) - shape*log(shape) + 1; - diag_z(0) = exp( (z1 - 1) / (1 + shape) ); - maxJ_bounded = CppAD::CondExpGe(diag_z(0), Type(maxsum), Type(maxsum), diag_z(0)); - max_log_wJ = maxJ_bounded*log(intensity) + (maxJ_bounded*shape)*log(x/scale) - lgamma(maxJ_bounded+1) - lgamma(maxJ_bounded*shape); - diag_z(1) = diag_z(0)*log(intensity) + (diag_z(0)*shape)*log(x/scale) - lgamma(diag_z(0)+1) - lgamma(diag_z(0)*shape); - } - // Integration constant - Type W = 0; - Type log_w_j; - //Type pos_penalty; - for( int j=minsum; j<=maxsum; j++ ){ - Type j2 = j; - //W += pow(intensity,j) * pow(x/scale, j2*shape) / exp(lgamma(j2+1)) / exp(lgamma(j2*shape)) / exp(max_log_w_j); - log_w_j = j2*log(intensity) + (j2*shape)*log(x/scale) - lgamma(j2+1) - lgamma(j2*shape); - //W += exp( posfun(log_w_j, Type(-30), pos_penalty) ); - W += exp( log_w_j - max_log_wJ ); - if(j==minsum) diag_z(2) = log_w_j; - if(j==maxsum) diag_z(3) = log_w_j; - } - // Loglikelihood calculation - Type loglike = 0; - if( x==0 ){ - loglike = -intensity; - }else{ - loglike = -x/scale - intensity - log(x) + log(W) + max_log_wJ; - } - // Return - if(give_log) return loglike; else return exp(loglike); -} - -// Calculate B_cc -template -matrix calculate_B( int method, int n_f, int n_r, matrix Chi_fr, matrix Psi_fr, Type &jnll_pointer ){ - matrix B_ff( n_f, n_f ); - matrix BplusI_ff( n_f, n_f ); - matrix Chi_rf = Chi_fr.transpose(); - matrix Psi_rf = Psi_fr.transpose(); - matrix Identity_ff( n_f, n_f ); - Identity_ff.setIdentity(); - - // No interactions (default) - if( method==0 ){ - B_ff.setZero(); - } - // Simple co-integration -- complex unbounded eigenvalues - if( method==1 ){ - B_ff = Chi_fr * Psi_rf; - } - // Real eigenvalues - if( method==2 ){ - matrix Chi_ff( n_f, n_f ); - Chi_ff = Identity_ff; - // Make Chi_ff - vector colnorm_r( n_r ); - colnorm_r.setZero(); - for(int f=0; f Psi_ff( n_f, n_f ); - Psi_ff = Identity_ff; - for(int f=n_r; f L_ff(n_f, n_f); - L_ff.setZero(); - for(int r=0; r invChi_ff = atomic::matinv( Chi_ff ); - matrix trans_Psi_ff = Psi_ff.transpose(); - matrix trans_invPsi_ff = atomic::matinv( Psi_ff ).transpose(); - B_ff = Chi_ff * trans_Psi_ff; - B_ff = B_ff * L_ff; - B_ff = B_ff * trans_invPsi_ff; - B_ff = B_ff * invChi_ff; - // Penalize colnorm_r - jnll_pointer += ( log(colnorm_r)*log(colnorm_r) ).sum(); - } - // Complex bounded eigenvalues - if( method==3 ){ - BplusI_ff = Chi_fr * Psi_rf + Identity_ff; - // Extract eigenvalues - vector< std::complex > eigenvalues_B_ff = B_ff.eigenvalues(); - vector real_eigenvalues_B_ff = eigenvalues_B_ff.real(); - vector imag_eigenvalues_B_ff = eigenvalues_B_ff.imag(); - vector mod_eigenvalues_B_ff( n_f ); - // Calculate maximum eigenvalues - Type MaxEigen = 1; - for(int f=0; f -Type objective_function::operator() () -{ - using namespace R_inla; - using namespace Eigen; - using namespace density; - - // Dimensions - DATA_INTEGER(n_i); // Number of observations (stacked across all years) - DATA_INTEGER(n_s); // Number of "strata" (i.e., vectices in SPDE mesh) - DATA_INTEGER(n_x); // Number of real "strata" (i.e., k-means locations) - DATA_INTEGER(n_t); // Number of time-indices - DATA_INTEGER(n_c); // Number of categories (e.g., length bins) - DATA_INTEGER(n_e); // Number of error distributions - DATA_INTEGER(n_p); // Number of dynamic covariates - DATA_INTEGER(n_v); // Number of tows/vessels (i.e., levels for the factor explaining overdispersion) - DATA_INTEGER(n_l); // Number of indices to post-process - DATA_INTEGER(n_m); // Number of range metrics to use (probably 2 for Eastings-Northings) - - // Config - DATA_IVECTOR( Options_vec ); - // Slot 0 -- Aniso: 0=No, 1=Yes - // Slot 1 -- DEPRECATED - // Slot 2 -- AR1 on beta1 (year intercepts for 1st linear predictor) to deal with missing years: 0=No, 1=Yes - // Slot 3 -- AR1 on beta2 (year intercepts for 2nd linear predictor) to deal with missing years: 0=No, 1=Yes - // Slot 4 -- DEPRECATED - // Slot 5 -- Upper limit constant of integration calculation for infinite-series density functions (Conway-Maxwell-Poisson and Tweedie) - // Slot 6 -- Breakpoint in CMP density function - // Slot 7 -- Whether to use SPDE or 2D-AR1 hyper-distribution for spatial process: 0=SPDE; 1=2D-AR1; 2=Stream-network - // Slot 8 -- Whether to use F_ct or ignore it for speedup - DATA_IVECTOR(FieldConfig); // Input settings (vector, length 4) - DATA_IVECTOR(OverdispersionConfig); // Input settings (vector, length 2) - DATA_IMATRIX(ObsModel_ez); // Observation model - // Column 0: Probability distribution for data for each level of e_i - // Column 1: Link function for linear predictors for each level of c_i - // NOTE: nlevels(c_i) must be <= nlevels(e_i) - DATA_IVECTOR(VamConfig); - // Slot 0 -- method for calculating n_c-by-n_c interaction matrix, B_ff - // Slot 1 -- rank of interaction matrix B_ff - // Current implementation only makes sense when (1) intercepts are constant among years; (2) using a Poisson-link delta model; (3) n_f=n_c for spatio-temporal variation; (4) starts near equilibrium manifold - DATA_INTEGER(include_data); // Always use TRUE except for internal usage to extract GRMF normalization when turn off GMRF normalization in CPP - DATA_IVECTOR(Options); // Reporting options - // Slot 0: Calculate SE for Index_xctl - // Slot 1: Calculate SE for log(Index_xctl) - // Slot 2: Calculate mean_Z_ctm (i.e., center-of-gravity) - // Slot 3: Calculate SE for D_i (expected density for every observation) - // Slot 4: Calculate mean_D_tl and effective_area_tl - // Slot 5: Calculate standard errors for Covariance and Correlation among categories using factor-analysis parameterization - // Slot 6: Calculate synchrony for different periods specified via yearbounds_zz - // Slot 7: Calculate coherence and variance for Epsilon1_sct and Epsilon2_sct - // Slot 8: Calculate proportions and SE - // Slot 9: Include normalization in GMRF PDF - DATA_IMATRIX(yearbounds_zz); - // Two columns, and 1+ rows, specifying first and last t for each period used in calculating synchrony - - // Data vectors - DATA_VECTOR(b_i); // Response (biomass) for each observation - DATA_VECTOR(a_i); // Area swept for each observation (km^2) - DATA_IMATRIX(c_iz); // Category for each observation - DATA_IVECTOR(e_i); // Error distribution for each observation - DATA_IVECTOR(s_i); // Station for each observation - DATA_IMATRIX(t_iz); // Time-indices (year, season, etc.) for each observation - DATA_IVECTOR(v_i); // tows/vessels for each observation (level of factor representing overdispersion) - DATA_VECTOR(PredTF_i); // vector indicating whether an observatino is predictive (1=used for model evaluation) or fitted (0=used for parameter estimation) - DATA_MATRIX(a_xl); // Area for each "real" stratum(km^2) in each stratum - DATA_MATRIX(X_xj); // Covariate design matrix (strata x covariate) - DATA_ARRAY(X_xtp); // Covariate design matrix (strata x covariate) - DATA_MATRIX(Q_ik); // Catchability matrix (observations x variable) - DATA_IMATRIX(t_yz); // Matrix for time-indices of calculating outputs (abundance index and "derived-quantity") - DATA_MATRIX(Z_xm); // Derived quantity matrix - DATA_MATRIX(F_ct); // Matrix of annual fishing mortality for each category - - // Spatial network inputs - DATA_IVECTOR(parent_s); // Columns: 0=Parent index, 1=Child index, 2=Distance from parent to child - DATA_IVECTOR(child_s); // Columns: 0=Parent index, 1=Child index, 2=Distance from parent to child - DATA_VECTOR(dist_s); // Columns: 0=Parent index, 1=Child index, 2=Distance from parent to child - - // SPDE objects - DATA_STRUCT(spde,spde_t); - - // Aniso objects - DATA_STRUCT(spde_aniso,spde_aniso_t); - - // Sparse matrices for precision matrix of 2D AR1 process - // Q = M0*(1+rho^2)^2 + M1*(1+rho^2)*(-rho) + M2*rho^2 - DATA_SPARSE_MATRIX(M0); - DATA_SPARSE_MATRIX(M1); - DATA_SPARSE_MATRIX(M2); - - // Parameters - PARAMETER_VECTOR(ln_H_input); // Anisotropy parameters - PARAMETER_MATRIX(Chi_fr); // error correction responses - PARAMETER_MATRIX(Psi_fr); // error correction loadings, B_ff = Chi_fr %*% t(Psi_fr) - - // -- presence/absence fixed effects - PARAMETER_MATRIX(beta1_ct); // Year effect - PARAMETER_VECTOR(gamma1_j); // Static covariate effect - PARAMETER_ARRAY(gamma1_ctp); // Dynamic covariate effect - PARAMETER_VECTOR(lambda1_k); // Catchability coefficients - PARAMETER_VECTOR(L1_z); // Overdispersion parameters - PARAMETER_VECTOR(L_omega1_z); - PARAMETER_VECTOR(L_epsilon1_z); - PARAMETER(logkappa1); - PARAMETER(Beta_mean1); // mean-reversion for beta1_t - PARAMETER(logsigmaB1); // SD of beta1_t (default: not included in objective function) - PARAMETER(Beta_rho1); // AR1 for positive catch Epsilon component, Default=0 - PARAMETER(Epsilon_rho1); // AR1 for presence/absence Epsilon component, Default=0 - PARAMETER_VECTOR(log_sigmaratio1_z); // Ratio of variance for columns of t_iz - - // -- presence/absence random effects - PARAMETER_MATRIX(eta1_vf); - PARAMETER_ARRAY(Omegainput1_sf); // Expectation - PARAMETER_ARRAY(Epsiloninput1_sft); // Annual variation - - // -- positive catch rates fixed effects - PARAMETER_MATRIX(beta2_ct); // Year effect - PARAMETER_VECTOR(gamma2_j); // Covariate effect - PARAMETER_ARRAY(gamma2_ctp); // Dynamic covariate effect - PARAMETER_VECTOR(lambda2_k); // Catchability coefficients - PARAMETER_VECTOR(L2_z); // Overdispersion parameters - PARAMETER_VECTOR(L_omega2_z); - PARAMETER_VECTOR(L_epsilon2_z); - PARAMETER(logkappa2); - PARAMETER(Beta_mean2); // mean-reversion for beta2_t - PARAMETER(logsigmaB2); // SD of beta2_t (default: not included in objective function) - PARAMETER(Beta_rho2); // AR1 for positive catch Epsilon component, Default=0 - PARAMETER(Epsilon_rho2); // AR1 for positive catch Epsilon component, Default=0 - PARAMETER_VECTOR(log_sigmaratio2_z); // Ratio of variance for columns of t_iz - - // Error distribution parameters - PARAMETER_ARRAY(logSigmaM); - // Columns: 0=CV, 1=[usually not used], 2=[usually not used] - // Rows: Each level of e_i and/or c_i - // SigmaM[,0] indexed by e_i, e.g., SigmaM(e_i(i),0) - // SigmaM[,1] and SigmaM[,2] indexed by c_i, e.g., SigmaM(c_i(i),2) - - // -- positive catch rates random effects - PARAMETER_VECTOR(delta_i); - PARAMETER_MATRIX(eta2_vf); - PARAMETER_ARRAY(Omegainput2_sf); // Expectation - PARAMETER_ARRAY(Epsiloninput2_sft); // Annual variation - - // Indices -- i=Observation; j=Covariate; v=Vessel; t=Year; s=Stratum - int i,t,c; - - // Objective function - vector jnll_comp(14); - // Slot 0 -- spatial, encounter - // Slot 1 -- spatio-temporal, encounter - // Slot 2 -- spatial, positive catch - // Slot 3 -- spatio-temporal, positive catch - // Slot 4 -- tow/vessel overdispersion, encounter - // Slot 5 -- tow/vessel overdispersion, positive catch - // Slot 8 -- penalty on beta, encounter - // Slot 9 -- penalty on beta, positive catch - // Slot 10 -- likelihood of data, encounter - // Slot 11 -- likelihood of data, positive catch - // Slot 12 -- Likelihood of Lognormal-Poisson overdispersion delta_i - // Slot 13 -- penalty on estimate_B structure - jnll_comp.setZero(); - Type jnll = 0; - - // Derived parameters - Type Range_raw1, Range_raw2; - if( Options_vec(7)==0 ){ - Range_raw1 = sqrt(8) / exp( logkappa1 ); // Range = approx. distance @ 10% correlation - Range_raw2 = sqrt(8) / exp( logkappa2 ); // Range = approx. distance @ 10% correlation - } - if( (Options_vec(7)==1) | (Options_vec(7)==2) ){ - Range_raw1 = log(0.1) / logkappa1; // Range = approx. distance @ 10% correlation - Range_raw2 = log(0.1) / logkappa2; // Range = approx. distance @ 10% correlation - } - array SigmaM( n_e, 3 ); - SigmaM = exp( logSigmaM ); - - // Anisotropy elements - matrix H(2,2); - H(0,0) = exp(ln_H_input(0)); - H(1,0) = ln_H_input(1); - H(0,1) = ln_H_input(1); - H(1,1) = (1+ln_H_input(1)*ln_H_input(1)) / exp(ln_H_input(0)); - - //////////////////////// - // Calculate joint likelihood - //////////////////////// - - // Define interaction matrix for Epsilon1, and also the imapct of F_ct on intercepts - int n_f; - n_f = Epsiloninput1_sft.col(0).cols(); - matrix B_ff( n_f, n_f ); // Interactions among factors - B_ff = calculate_B( VamConfig(0), n_f, VamConfig(1), Chi_fr, Psi_fr, jnll_comp(13) ); - matrix iota_ct( n_c, n_t ); // Cumulative impact of fishing mortality F_ct in years <= current year t - // Calculate interaction matrix B_cc for categories if feasible - if( n_c==n_f ){ - matrix L_epsilon1_cf = loadings_matrix( L_epsilon1_z, n_c, n_f ); - // Assemble interaction matrix - matrix B_cc( n_c, n_c ); // Interactions among categories - B_cc = B_ff; - for( int c=0; c Btemp_cc( n_c, n_c ); - Btemp_cc = L_epsilon1_cf * B_cc; - B_cc = Btemp_cc * L_epsilon1_cf.inverse(); - } - REPORT( B_cc ); - REPORT( L_epsilon1_cf ); - ADREPORT( B_cc ); - // Define impact of F_ct on intercepts - if( Options_vec(8)==0 ){ - iota_ct.setZero(); - }else{ - // Use F_ct in first year as initial condition - if( Options_vec(8)==1 ){ - iota_ct.col(0) = -1 * F_ct.col(0); - } - // Use median of stationary distribution given F_ct in first year as initial condition - if( Options_vec(8)==2 ){ - vector iota0_c( n_c ); - matrix I_cc( n_c, n_c ); - matrix sumF_cc( n_c, n_c ); - I_cc.setIdentity(); - I_cc = I_cc - B_cc; - sumF_cc = I_cc.inverse(); - iota_ct.col(0) -= sumF_cc * F_ct.col(0); - REPORT( sumF_cc ); - } - for( int t=1; t Q1( n_s, n_s ); - Eigen::SparseMatrix Q2( n_s, n_s ); - GMRF_t gmrf_Q; - if( (Options_vec(7)==0) & (Options_vec(0)==0) ){ - Q1 = Q_spde(spde, exp(logkappa1)); - Q2 = Q_spde(spde, exp(logkappa2)); - } - if( (Options_vec(7)==0) & (Options_vec(0)==1) ){ - Q1 = Q_spde(spde_aniso, exp(logkappa1), H); - Q2 = Q_spde(spde_aniso, exp(logkappa2), H); - } - if( Options_vec(7)==1 ){ - Q1 = M0*pow(1+exp(logkappa1*2),2) + M1*(1+exp(logkappa1*2))*(-exp(logkappa1)) + M2*exp(logkappa1*2); - Q2 = M0*pow(1+exp(logkappa2*2),2) + M1*(1+exp(logkappa2*2))*(-exp(logkappa2)) + M2*exp(logkappa2*2); - } - if( Options_vec(7)==2 ){ - Q1 = Q_network( logkappa1, n_s, parent_s, child_s, dist_s ); - Q2 = Q_network( logkappa2, n_s, parent_s, child_s, dist_s ); - } - // Probability of encounter - gmrf_Q = GMRF( Q1, bool(Options(9)) ); - // Omega1 - n_f = Omegainput1_sf.cols(); - array Omegamean1_sf(n_s, n_f); - Omegamean1_sf.setZero(); - array Omega1_sc(n_s, n_c); - Omega1_sc = gmrf_by_category_nll(FieldConfig(0), Options_vec(7), VamConfig(2), n_s, n_c, logkappa1, Omegainput1_sf, Omegamean1_sf, L_omega1_z, gmrf_Q, jnll_comp(0), this); - // Epsilon1 - n_f = Epsiloninput1_sft.col(0).cols(); - array Epsilonmean1_sf(n_s, n_f); - // PDF for Epsilon1 - array Epsilon1_sct(n_s, n_c, n_t); - for(t=0; t=1){ - // Prediction for spatio-temporal component - // Default, and also necessary whenever VamConfig(2)==1 & n_f!=n_c - if( (VamConfig(0)==0) | ((n_f!=n_c) & (VamConfig(2)==1)) ){ - // If no interactions, then just autoregressive for factors - Epsilonmean1_sf = Epsilon_rho1 * Epsiloninput1_sft.col(t-1); - }else{ - // Impact of interactions, B_ff - Epsilonmean1_sf.setZero(); - for(int s=0; s Omegamean2_sf(n_s, n_f); - Omegamean2_sf.setZero(); - array Omega2_sc(n_s, n_c); - Omega2_sc = gmrf_by_category_nll(FieldConfig(2), Options_vec(7), VamConfig(2), n_s, n_c, logkappa2, Omegainput2_sf, Omegamean2_sf, L_omega2_z, gmrf_Q, jnll_comp(2), this); - // Epsilon2 - n_f = Epsiloninput2_sft.col(0).cols(); - array Epsilonmean2_sf(n_s, n_f); - // PDF for Epsilon1 - array Epsilon2_sct(n_s, n_c, n_t); - for(t=0; t=1){ - // Prediction for spatio-temporal component - Epsilonmean2_sf = Epsilon_rho2 * Epsiloninput2_sft.col(t-1); - // Hyperdistribution for spatio-temporal component - Epsilon2_sct.col(t) = gmrf_by_category_nll(FieldConfig(3), Options_vec(7), VamConfig(2), n_s, n_c, logkappa2, Epsiloninput2_sft.col(t), Epsilonmean2_sf, L_epsilon2_z, gmrf_Q, jnll_comp(3), this); - } - } - - // Normalization of GMRFs to normalize during outer-optimization step in R - Type jnll_GMRF = jnll_comp(0) + jnll_comp(1) + jnll_comp(2) + jnll_comp(3); - if( include_data == 0 ){ - return( jnll_GMRF ); - } - - ////// Probability of correlated overdispersion among bins - // IN: eta1_vf; n_f; L1_z - // OUT: jnll_comp; eta1_vc - matrix eta1_vc(n_v, n_c); - eta1_vc = overdispersion_by_category_nll( OverdispersionConfig(0), n_v, n_c, eta1_vf, L1_z, jnll_comp(4), this ); - matrix eta2_vc(n_v, n_c); - eta2_vc = overdispersion_by_category_nll( OverdispersionConfig(1), n_v, n_c, eta2_vf, L2_z, jnll_comp(5), this ); - - // Possible structure on betas - if( Options_vec(2)!=0 ){ - for(c=0; c eta1_x = X_xj * gamma1_j.matrix(); - vector zeta1_i = Q_ik * lambda1_k.matrix(); - vector eta2_x = X_xj * gamma2_j.matrix(); - vector zeta2_i = Q_ik * lambda2_k.matrix(); - array eta1_xct(n_x, n_c, n_t); - array eta2_xct(n_x, n_c, n_t); - eta1_xct.setZero(); - eta2_xct.setZero(); - for(int x=0; x var_i(n_i); - Type tmp_calc1; - Type tmp_calc2; - // Linear predictor (pre-link) for presence/absence component - matrix P1_iz(n_i,c_iz.row(0).size()); - // Response predictor (post-link) - // ObsModel_ez(e,0) = 0:4 or 11:12: probability ("phi") that data is greater than zero - // ObsModel_ez(e,0) = 5 (ZINB): phi = 1-ZeroInflation_prob -> Pr[D=0] = NB(0|mu,var)*phi + (1-phi) -> Pr[D>0] = phi - NB(0|mu,var)*phi - vector R1_i(n_i); - vector log_one_minus_R1_i(n_i); - vector log_R1_i(n_i); - vector LogProb1_i(n_i); - // Linear predictor (pre-link) for positive component - matrix P2_iz(n_i,c_iz.row(0).size()); - // Response predictor (post-link) - // ObsModel_ez(e,0) = 0:3, 11:12: expected value of data, given that data is greater than zero -> E[D] = mu*phi - // ObsModel_ez(e,0) = 4 (ZANB): expected value ("mu") of neg-bin PRIOR to truncating Pr[D=0] -> E[D] = mu/(1-NB(0|mu,var))*phi ALSO Pr[D] = NB(D|mu,var)/(1-NB(0|mu,var))*phi - // ObsModel_ez(e,0) = 5 (ZINB): expected value of data for non-zero-inflation component -> E[D] = mu*phi - vector R2_i(n_i); - vector LogProb2_i(n_i); - vector maxJ_i(n_i); - vector diag_z(4); - matrix diag_iz(n_i,4); - diag_iz.setZero(); // Used to track diagnostics for Tweedie distribution (columns: 0=maxJ; 1=maxW; 2=lowerW; 3=upperW) - P1_iz.setZero(); - P2_iz.setZero(); - - // Likelihood contribution from observations - LogProb1_i.setZero(); - LogProb2_i.setZero(); - for(int i=0; i=0) & (c_iz(i,zc)=0) & (t_iz(i,zt)=0) & (c_iz(i,zc) 0 ){ - LogProb1_i(i) = log_R1_i(i); - }else{ - LogProb1_i(i) = log_one_minus_R1_i(i); - } - }else{ - if( b_i(i) > 0 ){ - LogProb1_i(i) = log( R1_i(i) ); - }else{ - LogProb1_i(i) = log( 1-R1_i(i) ); - } - } - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rbinom( Type(1), R1_i(i) ); - } - // Positive density likelihood -- models with continuous positive support - if( b_i(i) > 0 ){ // 1e-500 causes overflow on laptop - if(ObsModel_ez(e_i(i),0)==0){ - LogProb2_i(i) = dnorm(b_i(i), R2_i(i), SigmaM(e_i(i),0), true); - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rnorm( R2_i(i), SigmaM(e_i(i),0) ); - } - } - if(ObsModel_ez(e_i(i),0)==1){ - LogProb2_i(i) = dlnorm(b_i(i), log(R2_i(i))-pow(SigmaM(e_i(i),0),2)/2, SigmaM(e_i(i),0), true); // log-space - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = exp(rnorm( log(R2_i(i))-pow(SigmaM(e_i(i),0),2)/2, SigmaM(e_i(i),0) )); - } - } - if(ObsModel_ez(e_i(i),0)==2){ - LogProb2_i(i) = dgamma(b_i(i), 1/pow(SigmaM(e_i(i),0),2), R2_i(i)*pow(SigmaM(e_i(i),0),2), true); // shape = 1/CV^2, scale = mean*CV^2 - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rgamma( 1/pow(SigmaM(e_i(i),0),2), R2_i(i)*pow(SigmaM(e_i(i),0),2) ); - } - } - }else{ - LogProb2_i(i) = 0; - } - } - // Likelihood for Tweedie model with continuous positive support - if(ObsModel_ez(e_i(i),0)==8){ - LogProb1_i(i) = 0; - //dPoisGam( Type x, Type shape, Type scale, Type intensity, Type &max_log_w_j, int maxsum=50, int minsum=1, int give_log=0 ) - LogProb2_i(i) = dPoisGam( b_i(i), SigmaM(e_i(i),0), R2_i(i), R1_i(i), diag_z, Options_vec(5), Options_vec(6), true ); - diag_iz.row(i) = diag_z; - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = 0; // Option not available - } - } - // Likelihood #2 for Tweedie model with continuous positive support - if(ObsModel_ez(e_i(i),0)==10){ - // Packaged code - LogProb1_i(i) = 0; - // dtweedie( Type y, Type mu, Type phi, Type p, int give_log=0 ) - // R1*R2 = mean - LogProb2_i(i) = dtweedie( b_i(i), R1_i(i)*R2_i(i), R1_i(i), invlogit(SigmaM(e_i(i),0))+1.0, true ); - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = 0; // Option not available - } - } - ///// Likelihood for models with discrete support - // Zero-inflated negative binomial (not numerically stable!) - if(ObsModel_ez(e_i(i),0)==5){ - var_i(i) = R2_i(i)*(1.0+SigmaM(e_i(i),0)) + pow(R2_i(i),2.0)*SigmaM(c_iz(i,0),1); - if( b_i(i)==0 ){ - //LogProb2_i(i) = log( (1-R1_i(i)) + dnbinom2(Type(0.0), R2_i(i), var_i(i), false)*R1_i(i) ); // Pr[X=0] = 1-phi + NB(X=0)*phi - LogProb2_i(i) = logspace_add( log(1-R1_i(i)), dnbinom2(Type(0.0),R2_i(i),var_i(i),true)+log(R1_i(i)) ); // Pr[X=0] = 1-phi + NB(X=0)*phi - }else{ - LogProb2_i(i) = dnbinom2(b_i(i), R2_i(i), var_i(i), true) + log(R1_i(i)); // Pr[X=x] = NB(X=x)*phi - } - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rbinom( Type(1), R1_i(i) ); - if( b_i(i)>0 ){ - b_i(i) = rnbinom2( R2_i(i), var_i(i) ); - } - } - } - // Conway-Maxwell-Poisson - if(ObsModel_ez(e_i(i),0)==6){ - LogProb2_i(i) = dCMP(b_i(i), R2_i(i), exp(P1_iz(i,0)), true, Options_vec(5)); - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = 0; // Option not available - } - } - // Zero-inflated Poisson - if(ObsModel_ez(e_i(i),0)==7){ - if( b_i(i)==0 ){ - //LogProb2_i(i) = log( (1-R1_i(i)) + dpois(Type(0.0), R2_i(i), false)*R1_i(i) ); // Pr[X=0] = 1-phi + Pois(X=0)*phi - LogProb2_i(i) = logspace_add( log(1-R1_i(i)), dpois(Type(0.0),R2_i(i),true)+log(R1_i(i)) ); // Pr[X=0] = 1-phi + Pois(X=0)*phi - }else{ - LogProb2_i(i) = dpois(b_i(i), R2_i(i), true) + log(R1_i(i)); // Pr[X=x] = Pois(X=x)*phi - } - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rbinom( Type(1), R1_i(i) ); - if( b_i(i)>0 ){ - b_i(i) = rpois( R2_i(i) ); - } - } - } - // Binned Poisson (for REEF data: 0=none; 1=1; 2=2-10; 3=>11) - /// Doesn't appear stable given spatial or spatio-temporal variation - if(ObsModel_ez(e_i(i),0)==9){ - vector logdBinPois(4); - logdBinPois(0) = logspace_add( log(1-R1_i(i)), dpois(Type(0.0), R2_i(i), true) + log(R1_i(i)) ); // Pr[X=0] = 1-phi + Pois(X=0)*phi - logdBinPois(1) = dpois(Type(1.0), R2_i(i), true) + log(R1_i(i)); // Pr[X | X>0] = Pois(X)*phi - logdBinPois(2) = dpois(Type(2.0), R2_i(i), true) + log(R1_i(i)); // SUM_J( Pr[X|X>0] ) = phi * SUM_J( Pois(J) ) - for(int j=3; j<=10; j++){ - logdBinPois(2) += logspace_add( logdBinPois(2), dpois(Type(j), R2_i(i), true) + log(R1_i(i)) ); - } - logdBinPois(3) = logspace_sub( log(Type(1.0)), logdBinPois(0) ); - logdBinPois(3) = logspace_sub( logdBinPois(3), logdBinPois(1) ); - logdBinPois(3) = logspace_sub( logdBinPois(3), logdBinPois(2) ); - if( b_i(i)==0 ) LogProb2_i(i) = logdBinPois(0); - if( b_i(i)==1 ) LogProb2_i(i) = logdBinPois(1); - if( b_i(i)==2 ) LogProb2_i(i) = logdBinPois(2); - if( b_i(i)==3 ) LogProb2_i(i) = logdBinPois(3); - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = 0; // Option not available - } - } - // Zero-inflated Lognormal Poisson - if(ObsModel_ez(e_i(i),0)==11){ - if( b_i(i)==0 ){ - //LogProb2_i(i) = log( (1-R1_i(i)) + dpois(Type(0.0), R2_i(i), false)*R1_i(i) ); // Pr[X=0] = 1-phi + Pois(X=0)*phi - LogProb2_i(i) = logspace_add( log(1-R1_i(i)), dpois(Type(0.0),R2_i(i)*exp(SigmaM(e_i(i),0)*delta_i(i)-0.5*pow(SigmaM(e_i(i),0),2)),true)+log(R1_i(i)) ); // Pr[X=0] = 1-phi + Pois(X=0)*phi - }else{ - LogProb2_i(i) = dpois(b_i(i), R2_i(i)*exp(SigmaM(e_i(i),0)*delta_i(i)-0.5*pow(SigmaM(e_i(i),0),2)), true) + log(R1_i(i)); // Pr[X=x] = Pois(X=x)*phi - } - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rbinom( Type(1), R1_i(i) ); - if( b_i(i)>0 ){ - b_i(i) = rpois( R2_i(i)*exp(SigmaM(e_i(i),0)*delta_i(i)-0.5*pow(SigmaM(e_i(i),0),2)) ); - } - } - } - // Non-zero-inflated Poisson using log link from 1st linear predictor - if(ObsModel_ez(e_i(i),0)==12){ - LogProb2_i(i) = dpois(b_i(i), R1_i(i), true); - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rpois( R1_i(i) ); - } - } - // Non-zero-inflated Bernoulli using cloglog link from 1st lilnear predict - if(ObsModel_ez(e_i(i),0)==13){ - if( b_i(i)==0 ){ - LogProb2_i(i) = dpois(Type(0), R1_i(i), true); - }else{ - LogProb2_i(i) = logspace_sub( log(Type(1.0)), dpois(Type(0), R1_i(i), true) ); - } - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rpois( R1_i(i) ); - if( b_i(i)>0 ){ - b_i(i) = 1; - } - } - } - // Non-zero-inflated Lognormal-Poisson using log link from 1st linear predictor - if(ObsModel_ez(e_i(i),0)==14){ - LogProb2_i(i) = dpois(b_i(i), R1_i(i)*exp(SigmaM(e_i(i),0)*delta_i(i)-0.5*pow(SigmaM(e_i(i),0),2)), true); - // Simulate new values when using obj.simulate() - SIMULATE{ - b_i(i) = rpois( R1_i(i)*exp(SigmaM(e_i(i),0)*delta_i(i)-0.5*pow(SigmaM(e_i(i),0),2)) ); - } - } - } - } - REPORT( diag_iz ); - - // Joint likelihood - jnll_comp(10) = -1 * (LogProb1_i * (Type(1.0)-PredTF_i)).sum(); - jnll_comp(11) = -1 * (LogProb2_i * (Type(1.0)-PredTF_i)).sum(); - jnll = jnll_comp.sum(); - Type pred_jnll = -1 * ( LogProb1_i*PredTF_i + LogProb2_i*PredTF_i ).sum(); - REPORT( pred_jnll ); - REPORT( tmp_calc1 ); - REPORT( tmp_calc2 ); - - //////////////////////// - // Calculate outputs - //////////////////////// - - // Number of output-years - int n_y = t_yz.col(0).size(); - - // Predictive distribution -- ObsModel_ez(e,0)==4 isn't implemented (it had a bug previously) - Type a_average = a_i.sum()/a_i.size(); - array P1_xcy(n_x, n_c, n_y); - array R1_xcy(n_x, n_c, n_y); - array P2_xcy(n_x, n_c, n_y); - array R2_xcy(n_x, n_c, n_y); - array D_xcy(n_x, n_c, n_y); - for(int c=0; c=0) & (t_yz(y,z) Index_xcyl(n_x, n_c, n_y, n_l); - array Index_cyl(n_c, n_y, n_l); - array ln_Index_cyl(n_c, n_y, n_l); - Index_cyl.setZero(); - for(int c=0; c mean_Z_cym(n_c, n_y, n_m); - if( Options(2)==1 ){ - mean_Z_cym.setZero(); - int report_summary_TF = false; - for(int c=0; c mean_D_cyl(n_c, n_y, n_l); - array log_mean_D_cyl(n_c, n_y, n_l); - mean_D_cyl.setZero(); - for(int c=0; c effective_area_cyl(n_c, n_y, n_l); - array log_effective_area_cyl(n_c, n_y, n_l); - effective_area_cyl = Index_cyl / (mean_D_cyl/1000); // Correct for different units of Index and density - log_effective_area_cyl = log( effective_area_cyl ); - REPORT( effective_area_cyl ); - ADREPORT( effective_area_cyl ); - ADREPORT( log_effective_area_cyl ); - } - - // Reporting and standard-errors for covariance and correlation matrices - if( Options(5)==1 ){ - if( FieldConfig(0)>0 ){ - matrix L1_omega_cf = loadings_matrix( L_omega1_z, n_c, FieldConfig(0) ); - matrix lowercov_uppercor_omega1 = L1_omega_cf * L1_omega_cf.transpose(); - lowercov_uppercor_omega1 = convert_upper_cov_to_cor( lowercov_uppercor_omega1 ); - REPORT( lowercov_uppercor_omega1 ); - ADREPORT( lowercov_uppercor_omega1 ); - } - if( FieldConfig(1)>0 ){ - matrix L1_epsilon_cf = loadings_matrix( L_epsilon1_z, n_c, FieldConfig(1) ); - matrix lowercov_uppercor_epsilon1 = L1_epsilon_cf * L1_epsilon_cf.transpose(); - lowercov_uppercor_epsilon1 = convert_upper_cov_to_cor( lowercov_uppercor_epsilon1 ); - REPORT( lowercov_uppercor_epsilon1 ); - ADREPORT( lowercov_uppercor_epsilon1 ); - } - if( FieldConfig(2)>0 ){ - matrix L2_omega_cf = loadings_matrix( L_omega2_z, n_c, FieldConfig(2) ); - matrix lowercov_uppercor_omega2 = L2_omega_cf * L2_omega_cf.transpose(); - lowercov_uppercor_omega2 = convert_upper_cov_to_cor( lowercov_uppercor_omega2 ); - REPORT( lowercov_uppercor_omega2 ); - ADREPORT( lowercov_uppercor_omega2 ); - } - if( FieldConfig(3)>0 ){ - matrix L2_epsilon_cf = loadings_matrix( L_epsilon2_z, n_c, FieldConfig(3) ); - matrix lowercov_uppercor_epsilon2 = L2_epsilon_cf * L2_epsilon_cf.transpose(); - lowercov_uppercor_epsilon2 = convert_upper_cov_to_cor( lowercov_uppercor_epsilon2 ); - REPORT( lowercov_uppercor_epsilon2 ); - ADREPORT( lowercov_uppercor_epsilon2 ); - } - } - - // Synchrony - if( Options(6)==1 ){ - int n_z = yearbounds_zz.col(0).size(); - // Density ("D") or area-expanded total biomass ("B") for each category (use B when summing across sites) - matrix D_xy( n_x, n_y ); - matrix B_cy( n_c, n_y ); - vector B_y( n_y ); - D_xy.setZero(); - B_cy.setZero(); - B_y.setZero(); - // Sample variance in category-specific density ("D") and biomass ("B") - array varD_xcz( n_x, n_c, n_z ); - array varD_xz( n_x, n_z ); - array varB_cz( n_c, n_z ); - vector varB_z( n_z ); - vector varB_xbar_z( n_z ); - vector varB_cbar_z( n_z ); - vector ln_varB_z( n_z ); - vector ln_varB_xbar_z( n_z ); - vector ln_varB_cbar_z( n_z ); - array maxsdD_xz( n_x, n_z ); - array maxsdB_cz( n_c, n_z ); - vector maxsdB_z( n_z ); - varD_xcz.setZero(); - varD_xz.setZero(); - varB_cz.setZero(); - varB_z.setZero(); - varB_xbar_z.setZero(); - varB_cbar_z.setZero(); - maxsdD_xz.setZero(); - maxsdB_cz.setZero(); - maxsdB_z.setZero(); - // Proportion of total biomass ("P") for each location or each category - matrix propB_xz( n_x, n_z ); - matrix propB_cz( n_c, n_z ); - propB_xz.setZero(); - propB_cz.setZero(); - // Synchrony indices - matrix phi_xz( n_x, n_z ); - matrix phi_cz( n_c, n_z ); - vector phi_xbar_z( n_z ); - vector phi_cbar_z( n_z ); - vector phi_z( n_z ); - phi_xbar_z.setZero(); - phi_cbar_z.setZero(); - phi_z.setZero(); - // Calculate total biomass for different categories - for( int y=0; y CovHat( n_c, n_c ); - matrix CovHat( n_c, n_c ); - CovHat.setIdentity(); - CovHat *= pow(0.0001, 2); - if( FieldConfig(1)>0 ) CovHat += loadings_matrix(L_epsilon1_z, n_c, FieldConfig(1)) * loadings_matrix(L_epsilon1_z, n_c, FieldConfig(1)).transpose(); - if( FieldConfig(3)>0 ) CovHat += loadings_matrix(L_epsilon2_z, n_c, FieldConfig(3)) * loadings_matrix(L_epsilon2_z, n_c, FieldConfig(3)).transpose(); - // Coherence ranges from 0 (all factors are equal) to 1 (first factor explains all variance) - SelfAdjointEigenSolver > es(CovHat); - vector eigenvalues_c = es.eigenvalues(); // Ranked from lowest to highest for some reason - Type psi = 0; - for(int c=0; c diag_CovHat( n_c ); - vector log_diag_CovHat( n_c ); - for(int c=0; c PropIndex_cyl(n_c, n_y, n_l); - array ln_PropIndex_cyl(n_c, n_y, n_l); - Type sumtemp; - for(int y=0; y D_i( n_i ); - D_i = R1_i * R2_i; - ADREPORT( D_i ); - } - - return jnll; - -}