Fix zero total population bug in segregationcalc

src/segregation.cpp

@@ -5,76 +5,52 @@ using namespace Rcpp;
 NumericVector segregationcalc(NumericMatrix distmat,
                               NumericVector grouppop,
                               NumericVector fullpop) {
-
   // Vector to hold dissimilarity indices
   NumericVector diVec(distmat.ncol());
 
   // Population parameters
   int T = sum(fullpop);
-  double P = (double)sum(grouppop) / T;
+  double P = sum(grouppop) / T;
 
   // Calculate denominators
-  double d = (double)1 / (2 * T * P * (1 - P));
+  double d = 1.0 / (2.0 * T * P * (1 - P));
 
   // Get the number of unique plans
   NumericVector cd1 = distmat(_, 0);
   arma::vec cdVec1 = as<arma::vec> (cd1);
   arma::vec cdLabs = arma::unique(cdVec1);
 
-  // Range to look over for cd's
-  int start;
+  // Range to look over for CDs
   int end = max(cd1) + 1;
-  if(min(cd1) == 1){
-    start = 1;
-  }else{
-    start = 0;
-  }
+  int start = min(cd1) == 1 ? 1 : 0;
 
   // Loop over possible plans
-  for(int i = 0; i < distmat.ncol(); i++){
-
-    // Create dissimilarity objects
-    double dissim = 0;
+  for (int i = 0; i < distmat.ncol(); i++) {
+    double dissim = 0.0;
+    arma::vec cds = as<arma::vec> (distmat(_, i));
 
-    // Get a plan
-    NumericVector cdvec = distmat(_,i);
-    arma::vec cds = as<arma::vec> (cdvec);
-
-    // Loop over congressional districts
-    for(int j = start; j < end; j++){
-
-      // Initialize counts of groups
-      int tpop = 0;
-      int gpop = 0;
+    for (int j = start; j < end; j++) {
+      double tpop = 0.0;
+      double gpop = 0.0;
 
       // Which precincts in the plan are in this cd?
       arma::uvec findCds = find(cds == j);
 
       // Loop over precincts
-      for(int k = 0; k < findCds.size(); k++){
-
-        // Add population counts
+      for (int k = 0; k < findCds.size(); k++) {
         tpop += fullpop(findCds(k));
         gpop += grouppop(findCds(k));
-
       }
 
-      // Get district proportions
-      double p = (double)gpop / tpop;
-
-      // Add to dissimilarity index
-      dissim += (double)d * tpop * std::abs(p - P);
-
+      if (tpop > 0.0) {
+        dissim += d * tpop * std::abs(gpop / tpop - P);
+      }
     }
 
-    // Add to vector
     diVec(i) = dissim;
-
   }
 
-  // Return vector
   return diVec;
-
 }