Use protection

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: poppr
 Type: Package
 Title: Genetic Analysis of Populations with Mixed Reproduction
-Version: 2.3.0.99-61
-Date: 2017-04-05
+Version: 2.3.0.99-64
+Date: 2017-04-06
 Authors@R: c(person(c("Zhian", "N."), "Kamvar", role = c("cre", "aut"),
     email = "zkamvar@gmail.com"),
     person(c("Javier", "F."), "Tabima", role = "aut",

NEWS

@@ -54,12 +54,20 @@ BUG FIX
 * A bug where round-robin allele frequencies calculated with by_pop = TRUE were
   inaccurate for all but the first population was fixed. For details, see
   https://github.com/grunwaldlab/poppr/issues/132.
+* A potential integer overflow was fixed in `SEXP association_index_haploid`.
+  This was a ghost from https://github.com/grunwaldlab/poppr/issues/100.
+* PROTECT statements were placed around allocation statements. For details, see
+  https://github.com/grunwaldlab/poppr/issues/133.
 
 MISC
 ----
 
 * The documentation for `bitwise.dist()` clarifies the role of the
   `differences_only` flag (see https://github.com/grunwaldlab/poppr/issues/119).
+* Interruptions in C code is now handled gracefully via `R_CheckUserInterrupt()`.
+  The benefit is that long-running calculations are interrupted near instantly, 
+  but at the cost of a few more milliseconds of computation time.
+  (see https://github.com/grunwaldlab/poppr/issues/86)
 
 poppr 2.3.0
 ===========

README.md

@@ -1,7 +1,7 @@
 # Poppr version 2 <img src="vignettes/small_logo.png" align="right"/>
 
-[![Build Status](https://travis-ci.org/grunwaldlab/poppr.svg?branch=master)](https://travis-ci.org/grunwaldlab/poppr)
-[![Coverage Status](https://coveralls.io/repos/grunwaldlab/poppr/badge.svg?branch=master)](https://coveralls.io/r/grunwaldlab/poppr?branch=master)
+[![Build Status](https://travis-ci.org/grunwaldlab/poppr.svg?branch=use-protection)](https://travis-ci.org/grunwaldlab/poppr)
+[![Coverage Status](https://coveralls.io/repos/grunwaldlab/poppr/badge.svg?branch=use-protection)](https://coveralls.io/r/grunwaldlab/poppr?branch=use-protection)
 [![CRAN version](http://www.r-pkg.org/badges/version/poppr)](https://cran.r-project.org/package=poppr)
 <!-- 
 [![Downloads from Rstudio mirror per month](https://cranlogs.r-pkg.org/badges/poppr)](http://www.r-pkg.org/pkg/poppr)

src/bitwise_distance.c

@@ -36,6 +36,7 @@
 
 #include <stdio.h>
 #include <Rinternals.h>
+#include <R_ext/Utils.h>
 #include <Rdefines.h>
 #include <R.h>
 #include <math.h>
@@ -227,6 +228,7 @@ SEXP bitwise_distance_haploid(SEXP genlight, SEXP missing, SEXP requested_thread
   // Loop through every genotype 
   for(i = 0; i < num_gens; i++)
   {
+    R_CheckUserInterrupt();
     // Set R_chr1_1 to be genlight@gen[[i]]@snp[[1]], aka a raw list
     // representing the entire first set of chromosomes in this genotype
     R_chr1_1 = VECTOR_ELT(getAttrib(VECTOR_ELT(R_gen,i),R_chr_symbol),0);
@@ -549,6 +551,7 @@ SEXP bitwise_distance_diploid(SEXP genlight, SEXP missing, SEXP differences_only
   // Loop through every genotype 
   for(i = 0; i < num_gens; i++)
   {
+    R_CheckUserInterrupt();
     // Set R_chr1_1 to be genlight@gen[[i]]@snp[[1]], 
     // aka a raw list representing the entire first set of chromosomes in this genotype
     R_chr1_1 = VECTOR_ELT(getAttrib(VECTOR_ELT(R_gen,i),R_chr_symbol),0); 
@@ -811,8 +814,8 @@ SEXP association_index_haploid(SEXP genlight, SEXP missing, SEXP requested_threa
   double* vars; // Variance at each locus
   double* M;  // Sum of distances at each locus
   double* M2; // Sum of squared distances at each locus
-  int D;   // Sum of distances between each sample
-  int D2;  // Sum of squared distances between each sample
+  long int D;   // Sum of distances between each sample
+  long int D2;  // Sum of squared distances between each sample
   double Vo; // Observed variance
   double Ve; // Expected variance
   double Nc2;  // num_gens choose 2
@@ -916,6 +919,7 @@ SEXP association_index_haploid(SEXP genlight, SEXP missing, SEXP requested_threa
   #endif
   for(i = 0; i < num_chunks; i++)
   {
+    R_CheckUserInterrupt();
     // Loop through all samples
     for(j = 0; j < num_gens; j++)
     {
@@ -1041,7 +1045,10 @@ SEXP association_index_haploid(SEXP genlight, SEXP missing, SEXP requested_threa
       D2 += INTEGER(R_dists)[i + j*num_gens]*INTEGER(R_dists)[i + j*num_gens];
     }
   }
-
+  if (D2 < 0)
+  {
+    warning("\nAn integer overflow has occured and the resulting index will not be accurate.\nPlease consider using a smaller sample.\n");
+  }  
   // Calculate (num_gens choose 2), which will always be (n*n-n)/2 
   Nc2 = (num_gens*num_gens - num_gens)/2.0;
   // Calculate the observed variance using D and D2
@@ -1070,6 +1077,7 @@ SEXP association_index_haploid(SEXP genlight, SEXP missing, SEXP requested_threa
   #endif
   for(i = 0; i < num_loci; i++)
   {
+    R_CheckUserInterrupt();
     for(j = i+1; j < num_loci; j++)  
     {
       if(i != j)
@@ -1292,6 +1300,7 @@ SEXP association_index_diploid(SEXP genlight, SEXP missing, SEXP differences_onl
   #endif
   for(i = 0; i < num_chunks; i++)
   {
+    R_CheckUserInterrupt();
     // Loop through all samples
     for(j = 0; j < num_gens; j++)
     {
@@ -1469,6 +1478,7 @@ SEXP association_index_diploid(SEXP genlight, SEXP missing, SEXP differences_onl
   #endif
   for(i = 0; i < num_loci; i++)
   {
+    R_CheckUserInterrupt();
     for(j = i+1; j < num_loci; j++)  
     {
       if(i != j)
@@ -1564,6 +1574,7 @@ SEXP get_pgen_matrix_genind(SEXP genind, SEXP freqs, SEXP pops, SEXP npop)
   // Loop through all samples
   for (i = 0; i < num_gens; i++)
   {
+    R_CheckUserInterrupt();
     // Get the index for the allele frequencies
     pop = INTEGER(pops)[i] - 1;
     // h depends on ploidy. 

src/mlg_clustering.c

@@ -36,6 +36,7 @@
 
 #include <stdio.h>
 #include <Rinternals.h>
+#include <R_ext/Utils.h>
 #include <Rdefines.h>
 #include <R.h>
 #include <math.h>
@@ -234,6 +235,7 @@ SEXP neighbor_clustering(SEXP dist, SEXP mlg, SEXP threshold, SEXP algorithm, SE
   min_cluster_distance = -1;
   while(min_cluster_distance < thresh && num_clusters > 1)
   {
+    R_CheckUserInterrupt();
     min_cluster_distance = -1;
     closest_pair[0] = -1;
     closest_pair[1] = -1;

src/mlg_counter.c

@@ -37,6 +37,7 @@
 #include <string.h>
 #include <stdlib.h>
 #include <Rinternals.h>
+#include <R_ext/Utils.h>
 #include <R.h>
 
 int mlg_round_robin_cmpr (const void *a, const void *b);
@@ -185,6 +186,7 @@ SEXP mlg_round_robin(SEXP mat)
   // locus in the genotype_matrix, effectively masking the j + 1 locus. 
   for (j = 0; j < cols; j++)
   {
+    R_CheckUserInterrupt();
     qsort(mask_matrix, rows, sizeof(struct mask), mlg_round_robin_cmpr);
 
     for (i = 0; i < rows; i++)
@@ -282,6 +284,7 @@ SEXP genotype_curve_internal(SEXP mat, SEXP iter, SEXP maxloci, SEXP report)
   // Step 1: Loop over the number of loci
   while (nloci < nmax + 1)
   {
+    R_CheckUserInterrupt();
     // Initialize the global variable to be the number of loci we want to 
     // compare. 
     NLOCI = nloci*sizeof(int);

src/msn.c

@@ -36,6 +36,7 @@
 
 #include <stdio.h>
 #include <Rinternals.h>
+#include <R_ext/Utils.h>
 #include <Rdefines.h>
 #include <R.h>
 #include <math.h>
@@ -92,6 +93,7 @@ SEXP msn_tied_edges(SEXP mst, SEXP bclone, SEXP epsi)
   num_vertices = INTEGER(getAttrib(bclone,R_DimSymbol))[1];
   for(int i = 0; i < num_vertices; i++)
   {
+    R_CheckUserInterrupt();
     // Find the shortest path out of this vertex
     mn = -1;
     for(int j = 0; j < num_vertices; j++)

src/permut_shuffler.c

@@ -34,6 +34,7 @@
 #
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
 #include <Rinternals.h>
+#include <R_ext/Utils.h>
 #include <R.h>
 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 A slightly faster method of permuting alleles at a locus. 
@@ -63,12 +64,14 @@ SEXP permute_shuff(SEXP locus, SEXP alleles, SEXP ploidy)
 	rows = INTEGER(Rdim)[0]; 
 	cols = INTEGER(Rdim)[1]; 
 	PROTECT(Rout = allocMatrix(INTSXP, rows, cols));
-	alleles = coerceVector(alleles, INTSXP);
-	ploidy = coerceVector(ploidy, INTSXP);
+	PROTECT(alleles = coerceVector(alleles, INTSXP));
+	PROTECT(ploidy = coerceVector(ploidy, INTSXP));
 	ploid = INTEGER(ploidy);
 	int* inmat = INTEGER(locus);
 	int* outmat = INTEGER(Rout);
 	int* alle = INTEGER(alleles);
+	R_CheckUserInterrupt(); // This function doesn't run very long so, checking
+	                        // at this point is fine.
 	for(i = 0; i < rows; i++)
 	{
 		// loop through all columns first and initialize
@@ -97,7 +100,7 @@ SEXP permute_shuff(SEXP locus, SEXP alleles, SEXP ploidy)
 			}
 		}
 	}
-	UNPROTECT(1);
+	UNPROTECT(3); // Rout; alleles; ploidy
 	return Rout;
 }
 
@@ -124,6 +127,7 @@ SEXP expand_indices(SEXP indices, SEXP length) {
 	int max;
 	int min = 1;
 	rows = INTEGER(length)[0];
+	R_CheckUserInterrupt();
 	PROTECT(res = allocVector(VECSXP, rows));
 	for (i = 0; i < rows; i++)
 	{
@@ -135,9 +139,9 @@ SEXP expand_indices(SEXP indices, SEXP length) {
 			INTEGER(tempvec)[j] = min + j;
 		}
 		SET_VECTOR_ELT(res, i, tempvec);
-		UNPROTECT(1);
+		UNPROTECT(1); // for tempvec
 		min = ind[i] + 1;
 	}
-	UNPROTECT(1);
+	UNPROTECT(1); // for res
 	return res;
 }

src/poppr_distance.c

@@ -37,6 +37,7 @@
 #include <Rinternals.h>
 #include <Rdefines.h>
 #include <R.h>
+#include <R_ext/Utils.h>
 #include <math.h>
 #include <time.h>
 #include <string.h>
@@ -82,17 +83,18 @@ SEXP pairwise_covar(SEXP pair_vec)
 	int count;
 	SEXP Rout;
 	I = length(pair_vec);
-	pair_vec = coerceVector(pair_vec, REALSXP);
+	PROTECT(pair_vec = coerceVector(pair_vec, REALSXP));
 	PROTECT(Rout = allocVector(REALSXP, (I*(I-1)/2) ));
 	count = 0;
 	for(i = 0; i < I-1; i++)
 	{
+		R_CheckUserInterrupt();
 		for(j = i+1; j < I; j++)
 		{
 			REAL(Rout)[count++] = sqrt(REAL(pair_vec)[i] * REAL(pair_vec)[j]);
 		}
 	}
-	UNPROTECT(1);
+	UNPROTECT(2); // pair_vec; Rout
 	return Rout;
 }
 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -123,6 +125,7 @@ SEXP pairdiffs(SEXP freq_mat)
 
 	for(i = 0; i < rows-1; i++)
 	{
+		R_CheckUserInterrupt();
 		for(j = i+1; j < rows; j++)
 		{
 			val = 0;
@@ -225,7 +228,7 @@ SEXP bruvo_distance(SEXP bruvo_mat, SEXP permutations, SEXP alleles, SEXP m_add,
 	ploidy = INTEGER(coerceVector(alleles, INTSXP))[0];
 	loss = asLogical(m_loss);
 	add = asLogical(m_add);
-	bruvo_mat = coerceVector(bruvo_mat, INTSXP);
+	PROTECT(bruvo_mat = coerceVector(bruvo_mat, INTSXP));
 	perm = INTEGER(coerceVector(permutations, INTSXP));
 	PROTECT(Rval = allocMatrix(REALSXP, rows*(rows-1)/2, cols/ploidy));
 	PROTECT(pair_matrix = allocVector(INTSXP, 2*ploidy));
@@ -235,6 +238,7 @@ SEXP bruvo_distance(SEXP bruvo_mat, SEXP permutations, SEXP alleles, SEXP m_add,
 	{
 		for(i = 0; i < rows - 1; i++)
 		{
+			R_CheckUserInterrupt(); // in case the user wants to quit
 			for(allele = 0; allele < ploidy; allele++) 
 			{
 				clm = (allele + locus)*rows;
@@ -251,7 +255,7 @@ SEXP bruvo_distance(SEXP bruvo_mat, SEXP permutations, SEXP alleles, SEXP m_add,
 			}
 		}
 	}
-	UNPROTECT(2);
+	UNPROTECT(3); // bruvo_mat; Rval; pair_matrix
 	return Rval;
 }
 
@@ -365,6 +369,7 @@ polysat_bruvo() == poppr_bruvo()
 ==============================================================================*/
 double bruvo_dist(int *in, int *nall, int *perm, int *woo, int *loss, int *add)
 {
+	// R_CheckUserInterrupt();
 	int i; 
 	int j; 
 	int counter = 0;    // counter used for building arrays
@@ -662,6 +667,7 @@ void genome_add_calc(int perms, int alleles, int *perm, double **dist,
 	int zeroes, int *zero_ind, int curr_zero, int miss_ind, int *replacement, 
 	int inds, int curr_ind, double *genome_add_sum, int *tracker)
 {
+	// R_CheckUserInterrupt();
 	int i;
 	int j;
 	//==========================================================================
@@ -749,6 +755,7 @@ void genome_loss_calc(int *genos, int nalleles, int *perm_array, int woo,
 		int miss_ind, int curr_allele, double *genome_loss_sum, 
 		int *loss_tracker)
 {
+	// R_CheckUserInterrupt();
 	int i; 
 	int full_ind;
 	full_ind = 1 + (0 - miss_ind);
@@ -800,6 +807,7 @@ void fill_short_geno(int *genos, int nalleles, int *perm_array, int *woo,
 		int miss_ind, int *replacement, int inds, int curr_ind, double *res, 
 		int *tracker)
 {
+	// R_CheckUserInterrupt();
 	int i; //full_ind;
 	genos[miss_ind*nalleles + zero_ind[curr_zero]] = 
 		genos[miss_ind*nalleles + replacement[curr_ind]];