RF: ever_active is now 1-based indices to be more R like

jonathan-taylor · jonathan-taylor · commit 93a960275317 · 2017-08-20T16:05:30.000-07:00
diff --git a/selectiveInference/R/funs.fixed.R b/selectiveInference/R/funs.fixed.R
@@ -348,7 +348,7 @@ InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL, kkt
   if (is.null(soln_result)) {
      soln = rep(0, p)
      ever_active = rep(0, p)
-     ever_active[1] = i-1             # 0-based
+     ever_active[1] = i      # 1-based
      ever_active = as.integer(ever_active)
      nactive = as.integer(1)
      if (use_QP) {
@@ -373,7 +373,7 @@ InverseLinftyOneRow <- function (Sigma, i, mu, maxiter=50, soln_result=NULL, kkt
   if (use_QP) {
       result = solve_QP(Sigma, mu, maxiter, soln, linear_func, gradient, ever_active, nactive, kkt_tol, objective_tol) 
   } else {
-      result = find_one_row_debiasingM(Sigma, i-1, mu, maxiter, soln, gradient, ever_active, nactive, kkt_tol, objective_tol) # C function uses 0-based indexing
+      result = find_one_row_debiasingM(Sigma, i, mu, maxiter, soln, gradient, ever_active, nactive, kkt_tol, objective_tol) # C function uses 1-based indexing for active set
   }
 
   # Check feasibility
diff --git a/selectiveInference/src/debias.c b/selectiveInference/src/debias.c
@@ -11,12 +11,12 @@
 // Update one coordinate 
 
 double objective(double *Sigma_ptr,       /* A covariance matrix: X^TX/n */
-		 int *ever_active_ptr,    /* Ever active set: 0-based */ 
-		 int *nactive_ptr,    /* Size of ever active set */
-		 int nrow,            /* how many rows in Sigma */
-		 int row,             /* which row: 0-based */
-		 double bound,        /* Lagrange multipler for \ell_1 */
-		 double *theta)       /* current value */
+		 int *ever_active_ptr,    /* Ever active set: 1-based */ 
+		 int *nactive_ptr,        /* Size of ever active set */
+		 int nrow,                /* how many rows in Sigma */
+		 int row,                 /* which row: 1-based */
+		 double bound,            /* Lagrange multipler for \ell_1 */
+		 double *theta)           /* current value */
 {
   int irow, icol;
   double value = 0;
@@ -32,13 +32,13 @@ double objective(double *Sigma_ptr,       /* A covariance matrix: X^TX/n */
   for (irow=0; irow<nactive; irow++) {
 
     active_row_ptr = ((int *) ever_active_ptr + irow);
-    active_row = *active_row_ptr;
+    active_row = *active_row_ptr - 1;    // Ever-active set is 1-based
     theta_row_ptr = ((double *) theta + active_row);
 
     for (icol=0; icol<nactive; icol++) {
       
       active_col_ptr = ((int *) ever_active_ptr + icol);
-      active_col = *active_col_ptr;
+      active_col = *active_col_ptr - 1;    // Ever-active set is 1-based
       theta_col_ptr = ((double *) theta + active_col);
 
       Sigma_ptr_tmp = ((double *) Sigma_ptr + nrow * active_col + active_row); // Matrices are column-major order
@@ -48,14 +48,15 @@ double objective(double *Sigma_ptr,       /* A covariance matrix: X^TX/n */
     value = value + bound * fabs((*theta_row_ptr)); // the \ell_1 term
   }
 
-  theta_row_ptr = ((double *) theta + row);
+  theta_row_ptr = ((double *) theta + (row - 1));     // row is 1-based index
   value -= (*theta_row_ptr); // the elementary basis vector term
 
   return(value);
 }
 
 // Check if active and add it to active list if necessary
-
+// Ever active set is stored as 1-based indices
+// coord is 0-based
 int update_ever_active(int coord,
 		       int *ever_active_ptr,
 		       int *nactive_ptr) {
@@ -67,7 +68,7 @@ int update_ever_active(int coord,
   for (iactive=0; iactive<nactive; iactive++) {
     ever_active_ptr_tmp = ((int *) ever_active_ptr + iactive);
     active_var = (*ever_active_ptr_tmp);
-    if (active_var == coord) {
+    if (active_var - 1 == coord) {       // Indices in active set are 1-based
       return(1);
     }
   }
@@ -79,7 +80,7 @@ int update_ever_active(int coord,
   // number of active variables
 
   ever_active_ptr_tmp = ((int *) ever_active_ptr + *nactive_ptr);
-  *ever_active_ptr_tmp = coord;
+  *ever_active_ptr_tmp = coord + 1;     // Indices are 1-based
   *nactive_ptr += 1;
 
   return(0);
@@ -88,7 +89,7 @@ int update_ever_active(int coord,
 int check_KKT(double *theta,        /* current theta */
 	      double *gradient_ptr, /* Sigma times theta */
 	      int nrow,             /* how many rows in Sigma */
-	      int row,              /* which row: 0-based */
+	      int row,              /* which row: 1-based */
 	      double bound,         /* Lagrange multipler for \ell_1 */
 	      double tol)           /* precision for checking KKT conditions */
 {
@@ -108,7 +109,7 @@ int check_KKT(double *theta,        /* current theta */
 
     // For the basis vector
 
-    if (row == irow) {
+    if (row - 1 == irow) {                // Row is a 1-based index
        gradient -= 1;
     }
 
@@ -139,7 +140,7 @@ double update_one_coord(double *Sigma_ptr,           /* A covariance matrix: X^T
 			int nrow,                    /* How many rows in Sigma */
 			double bound,                /* feasibility parameter */
 			double *theta,               /* current value */
-			int row,                     /* which row: 0-based */
+			int row,                     /* which row: 1-based */
 			int coord,                   /* which coordinate to update: 0-based */
 			int is_active)               /* Is this part of ever_active */     
 {
@@ -156,8 +157,6 @@ double update_one_coord(double *Sigma_ptr,           /* A covariance matrix: X^T
   double *quadratic_ptr = ((double *) Sigma_diag_ptr + coord);
   double quadratic_term = *quadratic_ptr;
 
-  // int *ever_active_ptr_tmp;
-
   gradient_ptr_tmp = ((double *) gradient_ptr + coord);
   linear_term = *gradient_ptr_tmp;
 
@@ -167,9 +166,10 @@ double update_one_coord(double *Sigma_ptr,           /* A covariance matrix: X^T
   // The coord entry of gradient_ptr term has a diagonal term in it:
   // Sigma[coord, coord] * theta[coord]
   // This removes it. 
+
   linear_term -= quadratic_term * old_value;
 
-  if (row == coord) {
+  if (row - 1 == coord) {              // Row is 1-based
     linear_term -= 1;
   }
 
@@ -220,13 +220,13 @@ double update_one_coord(double *Sigma_ptr,           /* A covariance matrix: X^T
 int find_one_row_(double *Sigma_ptr,          /* A covariance matrix: X^TX/n */
 		  double *Sigma_diag_ptr,     /* Diagonal entry of covariance matrix */
 		  double *gradient_ptr,       /* Sigma times theta */
-		  int *ever_active_ptr,       /* Ever active set: 0-based */ 
+		  int *ever_active_ptr,       /* Ever active set: 1-based */ 
 		  int *nactive_ptr,           /* Size of ever active set */
 		  int nrow,                   /* How many rows in Sigma */
 		  double bound,               /* feasibility parameter */
 		  double *theta,              /* current value */
 		  int maxiter,                /* how many iterations */
-		  int row,                    /* which coordinate to update: 0-based */
+		  int row,                    /* which coordinate to solve: 1-based */
 		  double kkt_tol,             /* precision for checking KKT conditions */
 		  double objective_tol)       /* precision for checking relative decrease in objective value */
 {
@@ -267,7 +267,7 @@ int find_one_row_(double *Sigma_ptr,          /* A covariance matrix: X^TX/n */
 		       bound,
 		       theta,
 		       row,
-		       *active_ptr,
+		       *active_ptr-1,  // Ever active set is 1-based
 		       1);
       active_ptr++;
     }
diff --git a/selectiveInference/src/quadratic_program.c b/selectiveInference/src/quadratic_program.c
@@ -34,13 +34,13 @@ double objective_qp(double *Sigma_ptr,       /* A covariance matrix: X^TX/n */
   for (irow=0; irow<nactive; irow++) {
 
     active_row_ptr = ((int *) ever_active_ptr + irow);
-    active_row = *active_row_ptr;
+    active_row = *active_row_ptr - 1;          // Ever-active is 1-based
     theta_row_ptr = ((double *) theta + active_row);
 
     for (icol=0; icol<nactive; icol++) {
       
       active_col_ptr = ((int *) ever_active_ptr + icol);
-      active_col = *active_col_ptr;
+      active_col = *active_col_ptr - 1;          // Ever-active is 1-based
       theta_col_ptr = ((double *) theta + active_col);
 
       Sigma_ptr_tmp = ((double *) Sigma_ptr + nrow * active_col + active_row); // Matrices are column-major order
@@ -59,6 +59,8 @@ double objective_qp(double *Sigma_ptr,       /* A covariance matrix: X^TX/n */
   return(value);
 }
 
+// Ever-active is 1-based
+// coord is 0-based
 int update_ever_active_qp(int coord,
 			  int *ever_active_ptr,
 			  int *nactive_ptr) {
@@ -70,7 +72,7 @@ int update_ever_active_qp(int coord,
   for (iactive=0; iactive<nactive; iactive++) {
     ever_active_ptr_tmp = ((int *) ever_active_ptr + iactive);
     active_var = *ever_active_ptr_tmp;
-    if (active_var == coord) {
+    if (active_var - 1 == coord) {          // Ever-active is 1-based
       return(1);
     }
   }
@@ -82,7 +84,7 @@ int update_ever_active_qp(int coord,
   // number of active variables
 
   ever_active_ptr_tmp = ((int *) ever_active_ptr + *nactive_ptr);
-  *ever_active_ptr_tmp = coord;
+  *ever_active_ptr_tmp = coord + 1;         // Ever-active is 1-based
   *nactive_ptr += 1;
 
   return(0);
@@ -130,13 +132,13 @@ double update_one_coord_qp(double *Sigma_ptr,           /* A covariance matrix:
 			   double *linear_func_ptr,     /* Linear term in objective */
 			   double *Sigma_diag_ptr,      /* Diagonal entries of Sigma */
 			   double *gradient_ptr,        /* Sigma times theta */
-			   int *ever_active_ptr,        /* Ever active set: 0-based */ 
+			   int *ever_active_ptr,        /* Ever active set: 1-based */ 
 			   int *nactive_ptr,            /* Size of ever active set */
 			   int nrow,                    /* How many rows in Sigma */
 			   double bound,                /* feasibility parameter */
 			   double *theta,               /* current value */
 			   int coord,                   /* which coordinate to update: 0-based */
-			   int is_active)               /* Is this part of ever_active */     
+			   int is_active)               /* Is this coord in ever_active */     
 {
 
   double delta;
@@ -211,7 +213,7 @@ int solve_qp(double *Sigma_ptr,          /* A covariance matrix: X^TX/n */
 	     double *linear_func_ptr,    /* Linear term in objective */
 	     double *Sigma_diag_ptr,     /* Diagonal entry of covariance matrix */
 	     double *gradient_ptr,       /* Sigma times theta */
-	     int *ever_active_ptr,       /* Ever active set: 0-based */ 
+	     int *ever_active_ptr,       /* Ever active set: 1-based */ 
 	     int *nactive_ptr,           /* Size of ever active set */
 	     int nrow,                   /* How many rows in Sigma */
 	     double bound,               /* feasibility parameter */
@@ -259,7 +261,7 @@ int solve_qp(double *Sigma_ptr,          /* A covariance matrix: X^TX/n */
 			  nrow,
 			  bound,
 			  theta,
-			  *active_ptr,
+			  *active_ptr - 1,   // Ever-active is 1-based
 			  1);
       active_ptr++;
     }
