// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! // There are three debug levels (o to 2). 0 will show now degugging output. #define DEBUG 2 #if DEBUG == 1 #define DEBUG_OUTER #endif #if DEBUG == 2 #define DEBUG_INNER #define DEBUG_OUTER #endif // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! // Macro for printing debug output #define PRINTDBG(x) Rcout << __FUNCTION__ << "> " << #x << std::endl << x << std::endl #define PRINTDBG_MSG(x) Rcout << __FUNCTION__ << "> " << x << std::endl #include using namespace Rcpp; // [[Rcpp::export()]] NumericVector cpp_rvonmises(int n, double mu, double kappa) { /* Harry Southworth, January 2005. GENERATE RANDOM NUMBERS FROM A vonMises DISTRIBUTION. FOLLOWING THE ALGORITHM ON PAGE 43 OF MARDIA AND JUPP, DUE TO BEST AND FISHER. NOTE TAKEN OF ULRIC LUND'S IMPLEMENTATION IN S. THE ALGORITHM IS A REJECTION ALGORITHM AND RUNS SLOWLY IN S. Ported the C to R by Claudio Agostinelli, August 2006. Modified by James Campbell on 2023-06-22 for use in Rcpp code. Function written by James Campbell */ if(kappa == 0){ return(runif(n, -1*M_PI, M_PI)); } // init output vector NumericVector x(n); double U1, U2, U3; double a, b, r, z, f, c; /* SET a, b AND r */ a = 1 + sqrt(1+4* kappa * kappa); b = (a - sqrt(2*a))/(2* kappa); r = (1 + b*b)/(2*b); int i = 0; GetRNGstate(); do{ /* GENERATE U(0,1) RANDOM NUMBERS WITH THE R ROUTINE */ U1 = unif_rand(); z = cos(M_PI * U1); f = (1. + r * z)/(r + z); c = kappa * (r - f); U2 = unif_rand(); if(c * (2 - c) - U2 > 0) { U3 = unif_rand(); if (U3 > 0.50) x[i] = acos(f) + mu; else x[i] = -acos(f) + mu; i++; } else { if(log(c/U2) + 1 - c >= 0.) { U3 = unif_rand(); if (U3 > 0.50) x[i] = acos(f) + mu; else x[i] = -acos(f) + mu; i++; } } } while(i < n); PutRNGstate(); return(x); } enum TRANSFORM{NONE = 1, LOG = 2, COS = 3}; enum STEP_COL{X = 0, Y = 1, LEN = 2, BEARING = 3, TA = 4}; // [[Rcpp::export()]] double angle(double x){ if(x > 2*M_PI){ return(fmod(x, 2*M_PI)); } if(x < -2*M_PI){ return(fmod(x, 2*M_PI) + 2*M_PI); } if(x < 0){ return(x + 2*M_PI); } return(x); } NumericVector angle(NumericVector &x){ for(int i = 0; i < x.size(); i++){ x(i) = angle(x(i)); } return(x); } // Calculate random step vector // // Return a matrix where each row is a pseudo step. Columns 1:2 are the x/y // coords, and columns 3:4 are the step length and bearing. Bearing is calculate // from the vector difference between two step positions. // // [[Rcpp::export()]] void r_steps( NumericMatrix &new_steps, NumericVector ¤t_step, double shape, double scale, double kappa){ int n = new_steps.nrow(); NumericMatrix::Column col_x = new_steps( _ , STEP_COL::X); NumericMatrix::Column col_y = new_steps( _ , STEP_COL::Y); NumericMatrix::Column col_step = new_steps( _ , STEP_COL::LEN); NumericMatrix::Column col_bear = new_steps( _ , STEP_COL::BEARING); NumericMatrix::Column col_ta = new_steps( _ , STEP_COL::TA); // step lengths col_step = rgamma(n, shape, scale); // bearing NumericVector tmp; tmp = cpp_rvonmises(n, current_step(STEP_COL::BEARING), kappa); col_bear = angle(tmp); // turning angle tmp = tmp - current_step(STEP_COL::BEARING); col_ta = angle(tmp); // x/y position // Here, we're treating angles like geographic system, where north is 0 // degrees, pointing upwards (+y) and east is clockwise (+x). col_x = -sin(-col_bear) * col_step + current_step[STEP_COL::X]; col_y = cos(-col_bear) * col_step + current_step[STEP_COL::Y]; #ifdef DEBUG_INNER PRINTDBG(new_steps); #endif } // Extract covariate data from RasterLayer // [[Rcpp::export()]] double extract_covariate(List &iSSF_raster, NumericVector &step){ // Raster attributes NumericVector dims = iSSF_raster["dims"]; NumericVector extent = iSSF_raster["extent"]; bool directional = iSSF_raster["directional"]; NumericMatrix values = iSSF_raster["values"]; // Check if out of bounds if( step[0] < extent[0] | step[0] > extent[1] | step[1] < extent[2] | step[1] > extent[3] ){ return R_NaN; } // Get raster index to extract int icol = (int) floor(((step[0] - extent[0]) / (extent[1] - extent[0])) * (dims[1])); int irow = (int) floor(((extent[3] - step[1]) / (extent[3] - extent[2])) * (dims[0])); double out = values(irow, icol); // Check if covariate data is directional if(directional == true){ // Convert covariate into bearing difference out = angle(angle(out) - angle(step[STEP_COL::BEARING])); if(out > M_PI){ out -= 2*M_PI; }else if(out < -M_PI){ out += 2*M_PI; } out = cos(out); } return out; } // Get varaible value for formula variable double design_var(List &iSSF_raster_list, NumericVector &step, NumericVector ¤t_step, int ras_idx, int ras_trns, int ras_prev){ List iSSF_raster; double var; // Extract variable value switch(ras_idx) { case -1: // Extract step length var = step[2]; break; case -2: // Extract turning angle var = step[4]; break; default: // Extract raster value iSSF_raster = iSSF_raster_list[ras_idx]; if(ras_prev){ var = extract_covariate(iSSF_raster, current_step); }else{ var = extract_covariate(iSSF_raster, step); } } // Apply transformations switch (ras_trns) { case TRANSFORM::LOG: var = log(var); break; case TRANSFORM::COS: var = cos(var); break; } // #ifdef DEBUG_INNER // PRINTDBG(step); // PRINTDBG(var); // PRINTDBG(ras_idx); // PRINTDBG(ras_trns); // #endif return(var); } // Get value for foruma linear term double design_term(List &iSSF_raster_list, NumericVector &step, NumericVector ¤t_step, IntegerVector &ras_idx, IntegerVector &ras_trns, IntegerVector &ras_prev){ int n_var = ras_idx.size(); // Extract first variable from term double ret = design_var(iSSF_raster_list, step, current_step, ras_idx[0], ras_trns[0], ras_prev[0]); // Apply remaining interaction variables for(int i = 1; i < n_var; i++){ ret *= design_var(iSSF_raster_list, step, current_step, ras_idx[i], ras_trns[i], ras_prev[i]); } return ret; } // Combine extracted raster covariates to make design matrix. // This step is required for interaction variables. NumericMatrix design_matrix( List &iSSF_raster_list, NumericMatrix &steps, NumericVector current_step, List &ras_idxs, List &ras_trnsfs, List &ras_prev){ int n_term = ras_idxs.size(); int n_steps = steps.nrow(); NumericVector step(4); List term; NumericMatrix X(n_steps, n_term); IntegerVector idx, trns, prev; // loop steps for(int i_step = 0; i_step < n_steps; i_step++){ step = steps(i_step, _); // Loop and extract linear terms for(int i_term = 0; i_term < n_term; i_term++){ idx = ras_idxs[i_term]; trns = ras_trnsfs[i_term]; prev = ras_prev[i_term]; X(i_step, i_term) = design_term(iSSF_raster_list, step, current_step, idx, trns, prev); } } #ifdef DEBUG_INNER PRINTDBG(X); #endif return(X); } // Calculate step probability form covairate data // [[Rcpp::export()]] NumericVector step_probability_rast(NumericMatrix &X, NumericVector &beta_ras){ int n_step = X.nrow(); int n_param = beta_ras.size(); NumericVector numerator(n_step); NumericVector p(n_step); double tmp; // Loop steps for(int i_step = 0; i_step < n_step; i_step++){ // Set init probability tmp = 0; // Loop beta parameters for(int i_par = 0; i_par < n_param; i_par ++){ // Set probability to 0 if no covariate values if( R_isnancpp(X(i_step, i_par)) ){ tmp = R_NegInf; break; } tmp += X(i_step, i_par) * beta_ras[i_par]; } numerator(i_step) = exp(tmp); } double denominator = sum(numerator); if(denominator == 0){ p.fill(0); }else{ p = numerator / denominator; } #ifdef DEBUG_INNER PRINTDBG(p); #endif return(p); } //' Simulate an iSSF agent //' //' This function will simulated the movement of an agent over a provided //' landscape of covariate rasters. //' //' @param position_start Coordinate the agent will start at. //' @param bearing_start Initial bearing of the agent. //' @param n_steps The number of steps to simulate per track. //' @param len_track The number of tracks to simulate. //' @param ras_lst A list of raster covariate data, returned from `iSSF_raster_list()`. //' @param beta_par A vector holding the beta values for the model. //' @param beta_ras_idx A list holding indices of the raster layers for each term in `beta_par`. //' @param scale Scale parameter for gamma step length distribution. //' @param shape Shape parameter for gamma step length distribution. //' @param kappa kappa parameter for von Mises turning angle distribution. //' // [[Rcpp::export()]] List cpp_iSSF_simulation( NumericVector position_start, double bearing_start, int n_steps, int len_track, List iSSF_raster_list, NumericVector coef, List ras_idx, List ras_trnsfs, List ras_prev, double scale, double shape, double kappa ){ // Get number of raster layers int n_rasters = iSSF_raster_list.size(); // Get number of model parameters int n_coef = coef.size(); // Init agent data objects NumericMatrix positions(len_track + 1, 2); positions(0,_) = position_start; NumericVector bearing(len_track + 1); bearing(0) = bearing_start; NumericVector step_length(len_track + 1); step_length(0) = 0; NumericVector turning_angle(len_track + 1); turning_angle(0) = 0; NumericVector current_step(5); current_step[seq(0,1)] = position_start; current_step[2] = 0; current_step[3] = bearing_start; current_step[4] = 0; // Init step selection probability vector double p_vals[n_steps]; #ifdef DEBUG_OUTER PRINTDBG(position_start); PRINTDBG(bearing_start); PRINTDBG(n_steps); PRINTDBG(len_track); PRINTDBG(scale); PRINTDBG(shape); PRINTDBG(kappa); #endif // temp variables for holding agent position and possible new steps NumericMatrix new_steps(n_steps, 5); // Matrix holding raster covariates (columns) for each possible new step (rows) NumericMatrix ras_vals(n_steps, n_rasters); // Design matrix NumericMatrix X(n_steps, n_rasters); NumericVector p(n_steps); // rmultinom requires a pointer to store the return value. int drawn[n_steps]; int select; // Make objects for holding list elements List ras_lst_element; IntegerVector beta_ras_idx_element; // ============= loop time steps and simulate agent movement ================= for(int i = 1; i <= len_track; i++){ #ifdef DEBUG_OUTER PRINTDBG(i); #endif #ifdef DEBUG_OUTER PRINTDBG_MSG("Make new steps"); #endif // ------ Populate new_step matrix // Make new random steps r_steps(new_steps, current_step, shape, scale, kappa); #ifdef DEBUG_OUTER PRINTDBG_MSG("Make design matrix"); #endif // ------ Get linear design matrix // i.e. combine values for interacting covariate layers. X = design_matrix(iSSF_raster_list, new_steps, current_step, ras_idx ,ras_trnsfs, ras_prev); #ifdef DEBUG_OUTER PRINTDBG_MSG("Get new step probabilities"); #endif // ------ Calculate step probabilities: Raster covariates p = step_probability_rast(X, coef); // If all steps are out of bounds, stay still, skip to next step if(sum(p) == 0){ positions(i, 0) = positions(i - 1, 0); bearing(i) = bearing(i-1); turning_angle(i) = turning_angle(i-1); step_length(i) = R_NaN; continue; } // Draw new step // rmunlinom docs: // https://dirk.eddelbuettel.com/code/rcpp/html/Rmath_8h_source.html#l00102 // https://stackoverflow.com/questions/63368639/how-to-use-multinomial-function-in-c-with-rmath-h rmultinom(1, p.begin(), n_steps, drawn); for(int j = 0; j < n_steps; j++){ if(drawn[j] != 0){ select = j; break; } } // update current step current_step = new_steps(select, _); // write output positions(i, 0) = current_step[0]; positions(i, 1) = current_step[1]; step_length[i] = current_step[2]; #ifdef DEBUG_OUTER PRINTDBG(select); #endif } // ------ Return an agent List agent = List::create( _["track"] = positions, _["step_length"] = step_length, _["bearing"] = bearing, _["turning_angle"] = turning_angle); return(agent); }