diff --git a/DESCRIPTION b/DESCRIPTION
index 8e928ec..8260019 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: designmatch
 Type: Package
 Title: Matched Samples that are Balanced and Representative by Design
-Version: 0.5.1
-Date: 2023-06-21
+Version: 0.5.3
+Date: 2023-07-10
 Author: Jose R. Zubizarreta <zubizarreta@hcp.med.harvard.edu>, Cinar Kilcioglu <ckilcioglu16@gsb.columbia.edu>, Juan P. Vielma <jvielma@mit.edu>, Eric R. Cohn <ericcohn@g.harvard.edu>
 Maintainer: Jose R. Zubizarreta <zubizarreta@hcp.med.harvard.edu>
 Depends: R (>= 3.2), lattice, MASS, slam, highs
@@ -10,7 +10,6 @@ Enhances: gurobi, Rcplex, Rglpk, Rmosek, Rsymphony
 License: GPL-2 | GPL-3
 Description: Includes functions for the construction of matched samples that are balanced and representative by design.  Among others, these functions can be used for matching in observational studies with treated and control units, with cases and controls, in related settings with instrumental variables, and in discontinuity designs.  Also, they can be used for the design of randomized experiments, for example, for matching before randomization.  By default, 'designmatch' uses the 'highs' optimization solver, but its performance is greatly enhanced by the 'Gurobi' optimization solver and its associated R interface.  For their installation, please follow the instructions at <https://www.gurobi.com/documentation/quickstart.html> and <https://www.gurobi.com/documentation/7.0/refman/r_api_overview.html>.  We have also included directions in the gurobi_installation file in the inst folder.
 NeedsCompilation: no
-Packaged: 2023-06-26 19:11:34 UTC; eric
+Packaged: 2023-07-10 18:04:58 UTC; eric
 Repository: CRAN
-Date/Publication: 2023-06-26 22:10:02 UTC
-RoxygenNote: 7.2.3
+Date/Publication: 2023-07-10 22:40:02 UTC
diff --git a/MD5 b/MD5
index ef86d6c..d1b7827 100644
--- a/MD5
+++ b/MD5
@@ -1,10 +1,10 @@
-33685adfb5060742d93ab1f1bfa579c2 *DESCRIPTION
+55a37634505f6b52c530e90598041930 *DESCRIPTION
 45cba158a6e3de9c2156a518f2f24628 *NAMESPACE
 4035589e4087aa8c0dfc97e90db340f8 *R/absstddif.r
 7d1a85944bcd2594e3ab13ffdd80bd1e *R/addcalip.r
 9017e1fba5bece7acf56712cef12c2d8 *R/bmatch.r
 9dc5166095b7987dbf4e742345004242 *R/cardmatch.r
-9e77b9b929f571ca8ac335edda5e67f1 *R/constraintmatrix.r
+ca1d392d51c97a4cfd362f0a8f8147f2 *R/constraintmatrix.r
 ddd36bec26100d3064a1bf106c36740c *R/distmat.r
 692b745a194a270376faf4842b958a02 *R/distmatch.r
 f562ca6fc19a2951c9fcaaa01cae81dd *R/ecdfplot.r
@@ -12,10 +12,11 @@ d41d8cd98f00b204e9800998ecf8427e *R/errorhandling.r
 a7c14294e78b350eca6b9b036f9f0817 *R/finetab.r
 c98b06b08fd9add6edc0464b9d8d3443 *R/loveplot.r
 d3df5386a7329cd2c93dedac69a83db4 *R/meantab.r
-09803a6063703ebc862e15f714aeda64 *R/nmatch.r
+69aaf15b43d4c72ff3606b2d2ece217f *R/nmatch.r
+61da7fc7a11b84f586c6a093459df04d *R/onattach.r
 c445bfdc992990fc1986f8d2bf524573 *R/oneprob_profmatch.r
 98f291ae461a0b1c18ea5d8858b3eb22 *R/pairsplot.r
-36d380e4bd2fb8bbd551af173cc12056 *R/problemparameters.r
+0ea3c1641bba57aa254dd2834f3c5554 *R/problemparameters.r
 90e08362114fb82097302f29ef213616 *R/problemparameters_cardmatch.r
 60d3a3e224db7c5adf9e7789fa981e5e *R/problemparameters_profmatch.R
 4ef6dd019eb18dbc4f060d44a730d440 *R/profmatch.r
diff --git a/R/constraintmatrix.r b/R/constraintmatrix.r
index 16cb8ab..b8efb1a 100644
--- a/R/constraintmatrix.r
+++ b/R/constraintmatrix.r
@@ -10,14 +10,14 @@
                             far_covs, far_pairs, far_groups,
                             use_controls,
                             approximate) {
-  
+
   #! Number of treated units, number of controls
   n_t = sum(t_ind)
-  n_c = length(t_ind)-n_t	
-  
+  n_c = length(t_ind)-n_t
+
   #! Total number of units
   n_tot = n_t*n_c
-  
+
   #! Build parts of the constraint matrix
   #! Part 1
   if (approximate == 1 | n_controls == 1) {
@@ -30,14 +30,14 @@
     col_ind_1 = 1:(n_tot+n_t)
     ones_1 = c(rep(1, n_tot), rep(-1*n_controls, n_t))
   }
-  
+
   #! Part 2
   row_ind_2 = sort(rep(1:n_c, n_t))+n_t
   col_ind_2 = rep(seq(1, n_t*n_c, n_c), n_c)+(sort(rep(1:n_c, n_t))-1)
   ones_2 = rep(1, n_tot)
   #! Current max row index
   row_ind_cur	= max(row_ind_2)
-  
+
   #! Parts 3 and 4: moments and K-S
   mom_ks_covs = NULL
   if ((!is.null(mom_covs) & is.null(mom_targets)) | !is.null(ks_covs)) {
@@ -52,7 +52,7 @@
     if(!is.null(ks_covs)) {
       n_ks_covs = ncol(ks_covs)
     }
-    # Bind moment and K-S covariates		
+    # Bind moment and K-S covariates
     if(!is.null(mom_covs) & is.null(mom_targets) & is.null(ks_covs_aux)) {
       mom_ks_covs = mom_covs
       mom_ks_tols = mom_tols
@@ -71,13 +71,13 @@
       for (i in 1:ncol(ks_covs)) {
         mom_ks_tols = c(mom_ks_tols, rep(ks_tols[i], ks_n_grid[i]), rep(0, max(ks_n_grid)-ks_n_grid[i]))
       }
-    }	
-  }					 
+    }
+  }
   if (!is.null(mom_ks_covs)) {
     n_mom_ks_covs = ncol(mom_ks_covs)
     if ((!is.null(mom_tols) & is.null(mom_targets)) | !is.null(ks_tols)) {
       row_ind_3.4 = sort(rep(1:(2*n_mom_ks_covs)+n_t+n_c, n_tot))
-    }	
+    }
     col_ind_3.4 = NA
     mom_ks_vals_3.4 = NA
     j = 1
@@ -94,7 +94,7 @@
           k = k+1
           if (k >= max(ks_n_grid)) {
             j = j+1
-            k = 0	
+            k = 0
           }
         }
       }
@@ -102,7 +102,7 @@
       if ((!is.null(mom_tols) & is.null(mom_targets)) | !is.null(ks_tols)) {
         temp_mean_2 = temp_mean_1-(mom_ks_tols[i]*rep(1, n_t*n_c))
         temp_mean_3 = -temp_mean_1-(mom_ks_tols[i]*rep(1, n_t*n_c))
-      }	
+      }
       mom_ks_vals_3.4 = c(mom_ks_vals_3.4, temp_mean_2, temp_mean_3)
       if (i == 1) {
         col_ind_3.4 = col_ind_3.4[-1]
@@ -112,7 +112,7 @@
     #! Current max row index
     row_ind_cur	= max(row_ind_3.4)
   }
-  
+
   #! Moment target part
   rows_target = NULL
   cols_target = NULL
@@ -120,7 +120,7 @@
   if (!is.null(mom_covs) & !is.null(mom_targets)) {
     n_mom_covs = ncol(mom_covs)
     rows_target = sort(rep(1:(4*n_mom_covs)+row_ind_cur, n_tot))
-    
+
     for (i in 1:n_mom_covs) {
       cols_target = c(cols_target, rep(1:n_tot, 4))
       temp_treatment_1 = (mom_covs[t_ind==1, i])[sort(rep(1:n_t, n_c))] - (mom_targets[i] + mom_tols[i])
@@ -131,8 +131,8 @@
     }
     row_ind_cur = max(rows_target)
   }
-  
-  
+
+
   #! Part 5: exact
   rows_exact = NULL
   cols_exact = NULL
@@ -145,13 +145,13 @@
       dist_exact_cov = abs(outer(exact_covs[t_ind==1, i], exact_covs[t_ind==0, i], "-"))
       dist_exact_cov = t(dist_exact_cov)
       vals_exact = c(vals_exact, as.vector(dist_exact_cov))
-    }	
+    }
     row_ind_5 = rows_exact
     col_ind_5 = cols_exact
     exact_vals_5 = vals_exact
     row_ind_cur	= max(row_ind_5)
   }
-  
+
   #! Part 6: near-exact
   rows_near_exact = NULL
   cols_near_exact = NULL
@@ -164,13 +164,13 @@
       dist_near_exact_cov = abs(outer(near_exact_covs[t_ind==1, i], near_exact_covs[t_ind==0, i], "-"))
       dist_near_exact_cov = t(dist_near_exact_cov)
       vals_near_exact = c(vals_near_exact, as.vector(dist_near_exact_cov))
-    }	
+    }
     row_ind_6 = rows_near_exact
     col_ind_6 = cols_near_exact
     near_exact_vals_6 = vals_near_exact
     row_ind_cur	= max(row_ind_6)
   }
-  
+
   #! Part 7: fine
   bvec_7 = NULL
   rows_fine = NULL
@@ -181,7 +181,7 @@
     fine_covs_2 = rep(NA, nrow(fine_covs))
     n_fine_covs = ncol(fine_covs)
     j = 1
-    for (i in 1:n_fine_covs) {	
+    for (i in 1:n_fine_covs) {
       aux = factor(fine_covs[, i])
       fine_covs_2 = cbind(fine_covs_2, diag(nlevels(aux))[aux,])
       if (j == 1) {
@@ -196,14 +196,14 @@
       dist_fine_cov = outer(fine_covs_2[t_ind==1, i], fine_covs_2[t_ind==0, i], "-")
       dist_fine_cov = t(dist_fine_cov)
       vals_fine = c(vals_fine, as.vector(dist_fine_cov))
-    }	
+    }
     row_ind_7 = rows_fine
     col_ind_7 = cols_fine
     fine_vals_7 = vals_fine
     bvec_7 = rep(0, n_fine_cats)
     row_ind_cur	= max(row_ind_7)
   }
-  
+
   #! Part 8: near-fine
   bvec_8 = NULL
   rows_near_fine = NULL
@@ -214,7 +214,7 @@
     near_fine_covs_2 = rep(NA, nrow(near_fine_covs))
     n_near_fine_covs = ncol(near_fine_covs)
     j = 1
-    for (i in 1:n_near_fine_covs) {  
+    for (i in 1:n_near_fine_covs) {
       near_aux = factor(near_fine_covs[, i])
       near_fine_covs_2 = cbind(near_fine_covs_2, diag(nlevels(near_aux))[near_aux,])
       if (j == 1) {
@@ -240,7 +240,7 @@
     bvec_8 = rep(0, n_near_fine_cats)
     row_ind_cur	= max(row_ind_8)
   }
-  
+
   #! Part 9: Far
   rows_ind_far_pairs = list()
   if (!is.null(far_covs)) {
@@ -253,10 +253,10 @@
       #! Far on average constraints
       if (!is.null(far_groups)) {
         far_group = far_groups[j]
-        row_ind_far_all = sort(c(rep(row_ind_cur+1, n_tot)))		
-        col_ind_far_all = rep(1:n_tot, 1)			
+        row_ind_far_all = sort(c(rep(row_ind_cur+1, n_tot)))
+        col_ind_far_all = rep(1:n_tot, 1)
         temp_mean_3 = (-rep(far_cov[t_ind==0], n_t)+((far_cov[t_ind==1])[sort(rep(1:n_t, n_c))]))-(far_group*rep(1, n_t*n_c))
-        vals_far_all = c(temp_mean_3)		
+        vals_far_all = c(temp_mean_3)
         row_ind_cur	= max(row_ind_far_all)
       }
       #! Far on all pairs constraints
@@ -275,31 +275,31 @@
           rows_ind_far_pairs[[j]] = -1
           vals_far_pairs = NULL
         }
-      }		
+      }
       #! Put together
       if (!is.null(far_groups) && is.null(far_pairs)) {
         row_ind_9 = c(row_ind_9, row_ind_far_all)
         col_ind_9 = c(col_ind_9, col_ind_far_all)
         far_cov_vals_9 = c(far_cov_vals_9, vals_far_all)
       }
-      if (is.null(far_groups) && !is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (is.null(far_groups) && !is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         row_ind_9 = c(row_ind_9, rows_ind_far_pairs[[j]])
         col_ind_9 = c(col_ind_9, cols_ind_far_pairs)
         far_cov_vals_9 = c(far_cov_vals_9, vals_far_pairs)
       }
-      if (!is.null(far_groups) && !is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (!is.null(far_groups) && !is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         row_ind_9 = c(row_ind_9, row_ind_far_all, rows_ind_far_pairs[[j]])
         col_ind_9 = c(col_ind_9, col_ind_far_all, cols_ind_far_pairs)
         far_cov_vals_9 = c(far_cov_vals_9, vals_far_all, vals_far_pairs)
       }
-      if (!is.null(far_groups) && !is.null(far_pairs) && rows_ind_far_pairs[[j]] == -1) {
+      if (!is.null(far_groups) && !is.null(far_pairs) && all(rows_ind_far_pairs[[j]] == -1)) {
         row_ind_9 = c(row_ind_9, row_ind_far_all)
         col_ind_9 = c(col_ind_9, col_ind_far_all)
         far_cov_vals_9 = c(far_cov_vals_9, vals_far_all)
       }
     }
   }
-  
+
   #! Part 10: Near
   rows_ind_near_pairs = list()
   if (!is.null(near_covs)) {
@@ -312,10 +312,10 @@
       #! Near on average constraints
       if (!is.null(near_groups)) {
         near_group = near_groups[j]
-        row_ind_near_all = sort(c(rep(row_ind_cur+1, n_tot)))  	
-        col_ind_near_all = rep(1:n_tot, 1)			
+        row_ind_near_all = sort(c(rep(row_ind_cur+1, n_tot)))
+        col_ind_near_all = rep(1:n_tot, 1)
         temp_mean_4 = (-rep(near_cov[t_ind==0], n_t)+((near_cov[t_ind==1])[sort(rep(1:n_t, n_c))]))-(near_group*rep(1, n_t*n_c))
-        vals_near_all = c(temp_mean_4)		
+        vals_near_all = c(temp_mean_4)
         row_ind_cur	= max(row_ind_near_all)
       }
       #! Near on all pairs constraints
@@ -334,31 +334,31 @@
           rows_ind_near_pairs[[j]] = -1
           vals_near_pairs = NULL
         }
-      }		
+      }
       #! Put together
       if (!is.null(near_groups) && is.null(near_pairs)) {
         row_ind_10 = c(row_ind_10, row_ind_near_all)
         col_ind_10 = c(col_ind_10, col_ind_near_all)
         near_cov_vals_10 = c(near_cov_vals_10, vals_near_all)
       }
-      if (is.null(near_groups) && !is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (is.null(near_groups) && !is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         row_ind_10 = c(row_ind_10, rows_ind_near_pairs[[j]])
         col_ind_10 = c(col_ind_10, cols_ind_near_pairs)
         near_cov_vals_10 = c(near_cov_vals_10, vals_near_pairs)
       }
-      if (!is.null(near_groups) && !is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (!is.null(near_groups) && !is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         row_ind_10 = c(row_ind_10, row_ind_near_all, rows_ind_near_pairs[[j]])
         col_ind_10 = c(col_ind_10, col_ind_near_all, cols_ind_near_pairs)
         near_cov_vals_10 = c(near_cov_vals_10, vals_near_all, vals_near_pairs)
       }
-      if (!is.null(near_groups) && !is.null(near_pairs) && rows_ind_near_pairs[[j]] == -1) {
+      if (!is.null(near_groups) && !is.null(near_pairs) && all(rows_ind_near_pairs[[j]] == -1)) {
         row_ind_10 = c(row_ind_10, row_ind_near_all)
         col_ind_10 = c(col_ind_10, col_ind_near_all)
         near_cov_vals_10 = c(near_cov_vals_10, vals_near_all)
       }
     }
   }
-  
+
   # Part 11: use controls
   if (!is.null(use_controls)) {
     use_controls = use_controls[(n_t+1):(n_t+n_c)]
@@ -366,19 +366,19 @@
     col_ind_11 = (1:n_tot)[use_controls_aux==1]
     row_ind_11 = rep(row_ind_cur+1, length(col_ind_11))
     use_controls_vals_11 = rep(1, length(col_ind_11))
-    
+
     row_ind_cur = max(row_ind_11)
   }
-  
+
   # Part 12: total_groups
   if (!is.null(total_groups)) {
     row_ind_12 = rep(row_ind_cur+1, n_t*n_c)
     col_ind_12 = 1:(n_t*n_c)
     ones_12 = rep(1, n_t*n_c)
-    
+
     row_ind_cur = max(row_ind_12)
   }
-  
+
   #! Put all the parts of the constraint matrix together
   #! Parts 1 and 2
   row_ind = c(row_ind_1, row_ind_2)
@@ -424,19 +424,19 @@
   if (!is.null(far_covs)) {
     row_ind = c(row_ind, row_ind_9)
     col_ind = c(col_ind, col_ind_9)
-    vals = c(vals, far_cov_vals_9)	
-  }	
+    vals = c(vals, far_cov_vals_9)
+  }
   #! Part 10
   if (!is.null(near_covs)) {
     row_ind = c(row_ind, row_ind_10)
     col_ind = c(col_ind, col_ind_10)
-    vals = c(vals, near_cov_vals_10)  
-  }	
-  #! Part 11				 
+    vals = c(vals, near_cov_vals_10)
+  }
+  #! Part 11
   if (!is.null(use_controls)) {
     row_ind = c(row_ind, row_ind_11)
     col_ind = c(col_ind, col_ind_11)
-    vals = c(vals, use_controls_vals_11)	
+    vals = c(vals, use_controls_vals_11)
   }
   #! Part 12
   if (!is.null(total_groups)) {
@@ -444,13 +444,13 @@
     col_ind = c(col_ind, col_ind_12)
     vals = c(vals, ones_12)
   }
-  
+
   aux = cbind(row_ind, col_ind, vals)[order(col_ind), ]
-  
+
   aux = aux[(aux[, 3] != 0),]
   cnstrn_mat = simple_triplet_matrix(i = aux[, 1], j = aux[, 2], v = aux[, 3])
-  
-  #! Output  
+
+  #! Output
   return(list(cnstrn_mat = cnstrn_mat, bvec_7 = bvec_7, bvec_8 = bvec_8, rows_ind_far_pairs = rows_ind_far_pairs, rows_ind_near_pairs = rows_ind_near_pairs))
-    
-}
\ No newline at end of file
+
+}
diff --git a/R/nmatch.r b/R/nmatch.r
index f2d51ad..3ff1cdd 100644
--- a/R/nmatch.r
+++ b/R/nmatch.r
@@ -188,17 +188,17 @@ nmatch = function(dist_mat, subset_weight = NULL, total_pairs = NULL,
         cols_far = c(cols_far, col_ind_far_all)
         vals_far = c(vals_far, vals_far_all)
       }
-      if (is.null(far_groups) && !is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (is.null(far_groups) && !is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         rows_far = c(rows_far, rows_ind_far_pairs[[j]])
         cols_far = c(cols_far, cols_ind_far_pairs)
         vals_far = c(vals_far, vals_far_pairs)
       }
-      if (!is.null(far_groups) && !is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (!is.null(far_groups) && !is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         rows_far = c(rows_far, row_ind_far_all, rows_ind_far_pairs[[j]])
         cols_far = c(cols_far, col_ind_far_all, cols_ind_far_pairs)
         vals_far = c(vals_far, vals_far_all, vals_far_pairs)
       }
-      if (!is.null(far_groups) && !is.null(far_pairs) && rows_ind_far_pairs[[j]] == -1) {
+      if (!is.null(far_groups) && !is.null(far_pairs) && all(rows_ind_far_pairs[[j]] == -1)) {
         rows_far = c(rows_far, row_ind_far_all)
         cols_far = c(cols_far, col_ind_far_all)
         vals_far = c(vals_far, vals_far_all)
@@ -281,17 +281,17 @@ nmatch = function(dist_mat, subset_weight = NULL, total_pairs = NULL,
         cols_near = c(cols_near, col_ind_near_all)
         vals_near = c(vals_near, vals_near_all)
       }
-      if (is.null(near_groups) && !is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (is.null(near_groups) && !is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         rows_near = c(rows_near, rows_ind_near_pairs[[j]])
         cols_near = c(cols_near, cols_ind_near_pairs)
         vals_near = c(vals_near, vals_near_pairs)
       }
-      if (!is.null(near_groups) && !is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (!is.null(near_groups) && !is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         rows_near = c(rows_near, row_ind_near_all, rows_ind_near_pairs[[j]])
         cols_near = c(cols_near, col_ind_near_all, cols_ind_near_pairs)
         vals_near = c(vals_near, vals_near_all, vals_near_pairs)
       }
-      if (!is.null(near_groups) && !is.null(near_pairs) && rows_ind_near_pairs[[j]] == -1) {
+      if (!is.null(near_groups) && !is.null(near_pairs) && all(rows_ind_near_pairs[[j]] == -1)) {
         rows_near = c(rows_near, row_ind_near_all)
         cols_near = c(cols_near, col_ind_near_all)
         vals_near = c(vals_near, vals_near_all)
@@ -516,7 +516,7 @@ nmatch = function(dist_mat, subset_weight = NULL, total_pairs = NULL,
       if (!is.null(far_groups)) {
         bvec = c(bvec, rep(0, 1))
       }
-      if (!is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (!is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         bvec = c(bvec, rep(0, length(table(rows_ind_far_pairs[[j]]))))
       }
     }
@@ -531,7 +531,7 @@ nmatch = function(dist_mat, subset_weight = NULL, total_pairs = NULL,
       if (!is.null(near_groups)) {
         bvec = c(bvec, rep(0, 1))
       }
-      if (!is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (!is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         bvec = c(bvec, rep(0, length(table(rows_ind_near_pairs[[j]]))))
       }
     }
@@ -575,7 +575,7 @@ nmatch = function(dist_mat, subset_weight = NULL, total_pairs = NULL,
       if (!is.null(far_groups)) { 
         sense = c(sense, rep("G", 1)) 
       }
-      if (!is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (!is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         sense = c(sense, rep("E", length(table(rows_ind_far_pairs[[j]]))))
       }
     }
@@ -590,7 +590,7 @@ nmatch = function(dist_mat, subset_weight = NULL, total_pairs = NULL,
       if (!is.null(near_groups)) { 
         sense = c(sense, rep("L", 1)) 
       }
-      if (!is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (!is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         sense = c(sense, rep("E", length(table(rows_ind_near_pairs[[j]]))))
       }
     }
diff --git a/R/onattach.r b/R/onattach.r
new file mode 100644
index 0000000..35cdcb9
--- /dev/null
+++ b/R/onattach.r
@@ -0,0 +1,6 @@
+.onAttach <- 
+  function(libname, pkgname) {
+    packageStartupMessage("\nThank you for using our package! To acknowledge our work, please cite the package as: \n")
+    packageStartupMessage(" Zubizarreta, Jose R., Kilcioglu, Cinar, Vielma, Juan P., and Cohn, Eric R. (2023). designmatch: Matched Samples that are Balanced and Representative by Design.")
+    packageStartupMessage(" R package version 0.5.2. https://cran.r-project.org/web/packages/designmatch/ \n")
+  }
\ No newline at end of file
diff --git a/R/problemparameters.r b/R/problemparameters.r
index 9ece313..29430fd 100644
--- a/R/problemparameters.r
+++ b/R/problemparameters.r
@@ -10,14 +10,14 @@
                              far_covs, far_pairs, far_groups,
                              use_controls,
                              approximate) {
-  
+
   #! Number of treated units and controls
   n_t = sum(t_ind)
   n_c = length(t_ind)-n_t
-  
+
   #! Number of dec. vars.
   n_dec_vars = n_t*n_c
-  
+
   #! Number of moment covariates
   n_mom_covs = 0
   if(!is.null(mom_covs)) {
@@ -28,11 +28,11 @@
   n_ks_covs = 0
   if(!is.null(ks_covs)) {
     n_ks_covs = ncol(ks_covs)
-    
+
     if ((length(ks_n_grid)==1) && (n_ks_covs > 1)) {
       ks_n_grid = rep(ks_n_grid, n_ks_covs)
     }
-    
+
   }
   #! Parameters used to minimize the K-S statistic
   ks_covs_aux = NULL
@@ -41,14 +41,14 @@
   }
   if (!is.null(ks_covs)) {
     max_ks_n_grid = max(ks_n_grid)
-    #! Grid of values	
+    #! Grid of values
     ks_grid = matrix(0, nrow = max_ks_n_grid, ncol = n_ks_covs)
     for (i in 1:n_ks_covs) {
       ks_covs_t_aux = ks_covs[, i][t_ind==1]
       ks_grid_aux = quantile(ks_covs_t_aux, probs = seq(1/ks_n_grid[i], 1, 1/ks_n_grid[i]))
       ks_grid_aux = c(ks_grid_aux, rep(0, max_ks_n_grid-ks_n_grid[i]))
       ks_grid[, i] = ks_grid_aux
-    }	
+    }
     #! Auxiliary covariates
     ks_covs_aux = matrix(0, nrow = length(t_ind), ncol = max_ks_n_grid*n_ks_covs)
     for (i in 1:n_ks_covs) {
@@ -58,14 +58,14 @@
       }
     }
   }
-  
+
   #! Coeffs. of the obj. fun., cvec
   if (is.null(dist_mat)) {
     if (approximate == 1 | n_controls == 1) {
       cvec = -(1*rep(1, n_t*n_c))
     }
     else {
-      cvec = c(-(1*rep(1, n_t*n_c)), rep(0, n_t)) 
+      cvec = c(-(1*rep(1, n_t*n_c)), rep(0, n_t))
     }
   }
   if (!is.null(dist_mat)) {
@@ -75,9 +75,9 @@
     else {
       cvec = c(as.vector(matrix(t(dist_mat), nrow = 1, byrow = TRUE))-(subset_weight*rep(1, n_t*n_c)), rep(0, n_t))
     }
-    
+
   }
-  
+
   #! Constraint matrix, Amat
   constraintmat_out = .constraintmatrix(t_ind, n_controls, total_groups,
                                        mom_covs, mom_tols, mom_targets,
@@ -90,13 +90,13 @@
                                        far_covs, far_pairs, far_groups,
                                        use_controls,
                                        approximate)
-  
+
   cnstrn_mat = constraintmat_out$cnstrn_mat
   bvec_7 = constraintmat_out$bvec_7
   bvec_8 = constraintmat_out$bvec_8
   rows_ind_far_pairs = constraintmat_out$rows_ind_far_pairs
   rows_ind_near_pairs = constraintmat_out$rows_ind_near_pairs
-  
+
   # Constraint vector, bvec
   #! Parts 1 and 2
   if (approximate == 1 | n_controls == 1) {
@@ -105,37 +105,37 @@
   else {
     bvec = c(rep(0, n_t), rep(1, n_c))
   }
-  
+
   #! Part 3: moments
   if (!is.null(mom_covs) & is.null(mom_targets)) {
-    bvec = c(bvec, rep(0, 2*n_mom_covs))	
-  }	
-  
+    bvec = c(bvec, rep(0, 2*n_mom_covs))
+  }
+
   #! Part 4: K-S
   if (!is.null(ks_covs)) {
     bvec = c(bvec, rep(0, 2*n_ks_covs*max_ks_n_grid))
   }
-  
+
   #! Part 3b: Target
   if (!is.null(mom_covs) & !is.null(mom_targets)) {
     bvec = c(bvec, rep(0, 4*n_mom_covs))
   }
-  
+
   #! Part 5: exact
   if (!is.null(exact_covs)) {
-    bvec = c(bvec, rep(0, ncol(exact_covs))) 
-  }	
-  
+    bvec = c(bvec, rep(0, ncol(exact_covs)))
+  }
+
   #! Part 6: near-exact
   if (!is.null(near_exact_covs)) {
-    bvec = c(bvec, near_exact_devs) 
+    bvec = c(bvec, near_exact_devs)
   }
-  
+
   #! Part 7: fine
   if (!is.null(fine_covs)) {
-    bvec = c(bvec, bvec_7) 
+    bvec = c(bvec, bvec_7)
   }
-  
+
   #! Part 8: near-fine
   if (!is.null(near_fine_covs)) {
     n_near_fine_covs = ncol(near_fine_covs)
@@ -147,9 +147,9 @@
     bvec_8_aux = rep(NA, length(bvec_8)*2)
     bvec_8_aux[1:length(bvec_8)] = -near_fine_devs_aux
     bvec_8_aux[(length(bvec_8)+1):(2*length(bvec_8))] = near_fine_devs_aux
-    bvec = c(bvec, bvec_8_aux) 
+    bvec = c(bvec, bvec_8_aux)
   }
-  
+
   #! Part 9: far
   if (!is.null(far_covs)) {
     n_far_covs = ncol(far_covs)
@@ -157,12 +157,12 @@
       if (!is.null(far_groups)) {
         bvec = c(bvec, rep(0, 1))
       }
-      if (!is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (!is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         bvec = c(bvec, rep(0, length(table(rows_ind_far_pairs[[j]]))))
       }
     }
   }
-  
+
   #! Part 10: near
   if (!is.null(near_covs)) {
     n_near_covs = ncol(near_covs)
@@ -170,17 +170,17 @@
       if (!is.null(near_groups)) {
         bvec = c(bvec, rep(0, 1))
       }
-      if (!is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (!is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         bvec = c(bvec, rep(0, length(table(rows_ind_near_pairs[[j]]))))
       }
     }
   }
-  
+
   #! Part 11: use controls
   if (!is.null(use_controls)) {
-    bvec = c(bvec, sum(use_controls)) 
+    bvec = c(bvec, sum(use_controls))
   }
-  
+
   #! Part 12: total_groups
   if (!is.null(total_groups)) {
     if (!is.null(n_controls)) {
@@ -189,20 +189,20 @@
     else {
       bvec = c(bvec, total_groups)
     }
-    
+
   }
-  
+
   # Upper bounds, ub
   #! Parts 1 and 2
   #! Part 3: moments
   #! Part 4: K-S
   if (approximate == 1 | n_controls == 1) {
     ub = rep(1, n_t*n_c)
-  }	
+  }
   else {
     ub = c(rep(1, n_t*n_c), rep(1, n_t))
   }
-  
+
   # Sense, sense
   #! Parts 1 and 2
   #! Part 3: moments
@@ -213,69 +213,69 @@
   else {
     sense = c(rep("E", n_t), rep("L", n_c), rep("L", 2*n_mom_covs*(is.null(mom_targets))), rep("L", 2*n_ks_covs*max_ks_n_grid))
   }
-  
+
   #! Part 3b: Target
   if (!is.null(mom_covs) & !is.null(mom_targets)) {
     sense = c(sense, rep("L", 4*n_mom_covs))
   }
-  
+
   #! Part 5: exact
   if (!is.null(exact_covs)) {
-    sense = c(sense, rep("E", ncol(exact_covs))) 
+    sense = c(sense, rep("E", ncol(exact_covs)))
   }
-  
+
   #! Part 6: near-exact
   if (!is.null(near_exact_covs)) {
-    sense = c(sense, rep("L", ncol(near_exact_covs))) 
+    sense = c(sense, rep("L", ncol(near_exact_covs)))
   }
-  
+
   #! Part 7: fine
   if (!is.null(fine_covs)) {
-    sense = c(sense, rep("E", length(bvec_7))) 
+    sense = c(sense, rep("E", length(bvec_7)))
   }
-  
+
   #! Part 8: near-fine
   if (!is.null(near_fine_covs)) {
-    #sense = c(sense, rep(c("G", "L"), length(bvec_8))) 
+    #sense = c(sense, rep(c("G", "L"), length(bvec_8)))
     sense = c(sense, rep("G", length(bvec_8)), rep("L", length(bvec_8)))
   }
-  
+
   #! Part 9: far
   if (!is.null(far_covs)) {
     n_far_covs = ncol(far_covs)
     for (j in 1:n_far_covs) {
-      if (!is.null(far_groups)) { 
-        sense = c(sense, rep("G", 1)) 
+      if (!is.null(far_groups)) {
+        sense = c(sense, rep("G", 1))
       }
-      if (!is.null(far_pairs) && rows_ind_far_pairs[[j]] != -1) {
+      if (!is.null(far_pairs) && all(rows_ind_far_pairs[[j]] != -1)) {
         sense = c(sense, rep("E", length(table(rows_ind_far_pairs[[j]]))))
       }
     }
   }
-  
+
   #! Part 10: near
   if (!is.null(near_covs)) {
     n_near_covs = ncol(near_covs)
     for (j in 1:n_near_covs) {
-      if (!is.null(near_groups)) { 
-        sense = c(sense, rep("L", 1)) 
+      if (!is.null(near_groups)) {
+        sense = c(sense, rep("L", 1))
       }
-      if (!is.null(near_pairs) && rows_ind_near_pairs[[j]] != -1) {
+      if (!is.null(near_pairs) && all(rows_ind_near_pairs[[j]] != -1)) {
         sense = c(sense, rep("E", length(table(rows_ind_near_pairs[[j]]))))
       }
     }
   }
-  
+
   #! Part 11: use controls
   if (!is.null(use_controls)) {
-    sense = c(sense, "E") 
+    sense = c(sense, "E")
   }
-  
+
   #! Part 12: total_groups
   if (!is.null(total_groups)) {
     sense = c(sense, "E")
   }
-  
+
   # Variable types, vtype
   #! Parts 1 and 2
   #! Part 3: moments
@@ -289,7 +289,7 @@
   else {
     vtype = c(rep("B", n_t*n_c), rep("B", n_t))
   }
-  
+
   #! Part 5: exact
   #! Part 6: near-exact
   #! Part 7: fine
@@ -298,17 +298,17 @@
   #! Part 10: near
   #! Part 11: use controls
   #! Part 12: total_groups
-  
+
   # c_index
-  c_index = rep(1:n_c, n_t)	
-  
+  c_index = rep(1:n_c, n_t)
+
   # Output
-  return(list(n_t = n_t, n_c = n_c,  
-              cvec = cvec, 
-              Amat = cnstrn_mat, 
-              bvec = bvec, 
-              ub = ub, 
+  return(list(n_t = n_t, n_c = n_c,
+              cvec = cvec,
+              Amat = cnstrn_mat,
+              bvec = bvec,
+              ub = ub,
               sense = sense,
               vtype = vtype,
               c_index = c_index))
-}
\ No newline at end of file
+}