Slight code tidying - mostly allman bracketting and localising scope

rig_cpp_common/runtime/rig_cpp_common/maths/binomial.cpp

@@ -1,9 +1,9 @@
 #include "binomial.h"
 
 // Common includes
-#include "rig_cpp_common/fixed_point_number.h"
-#include "rig_cpp_common/random/mars_kiss64.h"
-#include "rig_cpp_common/maths/ln.h"
+#include "../fixed_point_number.h"
+#include "../random/mars_kiss64.h"
+#include "ln.h"
 #include "recip.h"
 
 // Namespaces
@@ -24,12 +24,14 @@ namespace
 uint32_t randbin_bg_core(uint32_t n, S1615 ln_1_min_p,
                          MarsKiss64 &rng)
 {
-  uint32_t y = 0, x = 0;
   if (ln_1_min_p >= 0)
-    return x;
-
-  S1615 recip_ln_1_min_p = Reciprocal(ln_1_min_p);
+  {
+    return 0;
+  }
 
+  const S1615 recip_ln_1_min_p = Reciprocal(ln_1_min_p);
+  unsigned int y = 0;
+  unsigned int x = 0;
   while (1)
   {
     // Strip off the sign bit of u
@@ -38,7 +40,9 @@ uint32_t randbin_bg_core(uint32_t n, S1615 ln_1_min_p,
     // random variate in [0,1] with 15 fractional bits
     y += (MulS1615(Ln(u) - 363408, recip_ln_1_min_p) >> 15) + 1;
     if (y > n)
+    {
       break;
+    }
     x += 1;
   }
 
@@ -56,7 +60,7 @@ namespace Maths
 {
 
 // Sample from the binomial distribution with the given parameters
-uint32_t Binomial(uint32_t n, S1615 p, MarsKiss64 &rng)
+unsigned int Binomial(unsigned int n, S1615 p, MarsKiss64 &rng)
 {
   if (p > 16384)
   {
@@ -66,7 +70,7 @@ uint32_t Binomial(uint32_t n, S1615 p, MarsKiss64 &rng)
   }
   else
   {
-    return randbin_bg_core(n, Ln(32768-p), rng);
+    return randbin_bg_core(n, Ln(32768 - p), rng);
   }
 }
 
@@ -75,9 +79,10 @@ uint32_t Binomial(uint32_t n, S1615 p, MarsKiss64 &rng)
 // Since we pass log(1-p) to randbin_bg_core, replace this with
 // log((denom - num)/denom) = log(denom-num) - log(denom) in case we can't represent
 // p given the resolution of our fixed point types.
-uint32_t Binomial(uint32_t n, uint32_t num, uint32_t denom, MarsKiss64 &rng)
+unsigned int Binomial(unsigned int n, unsigned int num, unsigned int denom,
+                      MarsKiss64 &rng)
 {
-  if ((num<<1) > denom)
+  if ((num << 1) > denom)
   {
     // If p > 0.5, sample from the binomial with 1-p, subtract the result
     // from n. This is more efficient, and gives the same distribution.

rig_cpp_common/runtime/rig_cpp_common/maths/binomial.h

@@ -1,7 +1,7 @@
 #pragma once
 
 // Common includes
-#include "rig_cpp_common/fixed_point_number.h"
+#include "../fixed_point_number.h"
 
 // Forward declarations
 namespace Common
@@ -23,8 +23,9 @@ namespace Common
 {
 namespace Maths
 {
-  uint32_t Binomial(uint32_t n, S1615 p, MarsKiss64 &rng);
+  unsigned int Binomial(unsigned int n, S1615 p, MarsKiss64 &rng);
 
-  uint32_t Binomial(uint32_t n, uint32_t numerator, uint32_t denominator, MarsKiss64 &rng);
+  unsigned int Binomial(unsigned int n, unsigned int numerator, unsigned int denominator,
+                        MarsKiss64 &rng);
 } // Maths
 } // Common

rig_cpp_common/runtime/rig_cpp_common/maths/hypergeometric.cpp

@@ -1,11 +1,11 @@
 #include "hypergeometric.h"
 
 // Common includes
-#include "rig_cpp_common/fixed_point_number.h"
-#include "rig_cpp_common/random/mars_kiss64.h"
+#include "../fixed_point_number.h"
+#include "../random/mars_kiss64.h"
+#include "ln.h"
+#include "exp.h"
 #include "logfact.h"
-#include "rig_cpp_common/maths/ln.h"
-#include "rig_cpp_common/maths/exp.h"
 
 // Namespaces
 using namespace Common::Maths;
@@ -22,54 +22,56 @@ namespace
 // Since the number of times through the loop is the value returned,
 // the expected execution time is proportional to
 // nsample * (ngood / (ngood+nbad))
-uint32_t randhg_hin_core(uint32_t ngood, uint32_t nbad, uint32_t nsample,
-                         MarsKiss64 &rng)
+unsigned int randhg_hin_core(unsigned int ngood, unsigned int nbad, unsigned int nsample,
+                             MarsKiss64 &rng)
 {
-  uint32_t n = ngood + nbad, x;
-  S1615 ln_p, p, u;
+  const unsigned int n = ngood + nbad;
 
   // We can never sample a value greater than either ngood or nsample
-  uint32_t maxval = (ngood < nsample) ? ngood : nsample;
+  const unsigned int maxval = (ngood < nsample) ? ngood : nsample;
 
+  unsigned int x;
   do
   {
+    // The log probability of returning the first value
+    S1615 ln_p;
+    if(nsample < nbad)
+    {
+      ln_p = LogFact(nbad) - LogFact(n) + LogFact(n-nsample) - LogFact(nbad-nsample);
+      x = 0;
+    }
+    else
+    {
+      ln_p = LogFact(ngood) - LogFact(n) + LogFact(nsample) - LogFact(nsample-nbad);
+      x = nsample - nbad;
+    }
+    // Offset the log probability so that the format of the probability p
+    // uses 30 fractional bits
+    ln_p += 340695;
 
-  // The log probability of returning the first value
-  if(nsample < nbad)
-  {
-    ln_p = LogFact(nbad) - LogFact(n) + LogFact(n-nsample) - LogFact(nbad-nsample);
-    x = 0;
-  }
-  else
-  {
-    ln_p = LogFact(ngood) - LogFact(n) + LogFact(nsample) - LogFact(nsample-nbad);
-    x = nsample - nbad;
-  }
-  // Offset the log probability so that the format of the probability p
-  // uses 30 fractional bits
-  ln_p += 340695;
-
-  // Sample from a uniform distribution [0, 1]
-  // Note, we're actually using 30 bits for the fractional part
-  u = (S1615)(rng.GetNext() & 0x3fffffff);
+    // Sample from a uniform distribution [0, 1]
+    // Note, we're actually using 30 bits for the fractional part
+    S1615 u = (S1615)(rng.GetNext() & 0x3fffffff);
 
-  // For successive possible values to return, subtract the probability
-  // of returning that value from u. When u <= p, return the current value x
-  p = ExpS1615(ln_p);
-  while(u > p)
-  {
-    u -= p;
-    // The log probability of x calculated in terms of the log probability
-    // of (x-1)
-    ln_p += Ln(ngood-x);
-    ln_p -= Ln(x+1);
-    ln_p += Ln(nsample-x);
-    ln_p -= Ln(nbad-nsample+1+x);
-    p = ExpS1615(ln_p);
-    x++;
-    if(x > maxval)
-      break;
-  }
+    // For successive possible values to return, subtract the probability
+    // of returning that value from u. When u <= p, return the current value x
+    S1615 p = ExpS1615(ln_p);
+    while(u > p)
+    {
+      u -= p;
+      // The log probability of x calculated in terms of the log probability
+      // of (x-1)
+      ln_p += Ln(ngood-x);
+      ln_p -= Ln(x+1);
+      ln_p += Ln(nsample-x);
+      ln_p -= Ln(nbad-nsample+1+x);
+      p = ExpS1615(ln_p);
+      x++;
+      if(x > maxval)
+      {
+        break;
+      }
+    }
 
   // With small probability, x might exceed maxval before u <= p,
   // due to numerical issues. If this happens repeat the procedure
@@ -89,7 +91,8 @@ namespace Maths
 {
 
 // Sample from a hypergeometric distribution with the given parameters
-uint32_t Hypergeom(uint32_t ngood, uint32_t nbad, uint32_t nsample, MarsKiss64 &rng)
+unsigned int Hypergeom(unsigned int ngood, unsigned int nbad, unsigned int nsample,
+                       MarsKiss64 &rng)
 {
   if(ngood <= nbad)
   {

rig_cpp_common/runtime/rig_cpp_common/maths/hypergeometric.h

@@ -1,8 +1,5 @@
 #pragma once
 
-// Common includes
-#include "rig_cpp_common/fixed_point_number.h"
-
 // Forward declarations
 namespace Common
 {
@@ -13,7 +10,6 @@ namespace Common
 }
 
 // Namespaces
-using namespace Common::FixedPointNumber;
 using namespace Common::Random;
 
 //-----------------------------------------------------------------------------
@@ -23,6 +19,7 @@ namespace Common
 {
 namespace Maths
 {
-  uint32_t Hypergeom(uint32_t ngood, uint32_t nbad, uint32_t nsample, MarsKiss64 &rng);
+  unsigned int Hypergeom(unsigned int ngood, unsigned int nbad, unsigned int nsample,
+                         MarsKiss64 &rng);
 } // Maths
 } // Common

rig_cpp_common/runtime/rig_cpp_common/maths/logfact.cpp

@@ -1,9 +1,7 @@
-// Common includes
-#include "rig_cpp_common/fixed_point_number.h"
+#include "logfact.h"
 
 // Maths includes
-#include "rig_cpp_common/maths/ln.h"
-#include "logfact.h"
+#include "ln.h"
 
 // Namespaces
 using namespace Common::Maths;
@@ -16,7 +14,7 @@ namespace
 {
 
 // Log factorial for arguments 0-63
-const S1615 logFact[64] = {
+const S1615 g_LogFact[64] = {
   0,0,22713,58712,104138,156876,215588,279352,347491,419490,494941,573515,
   654941,738989,825465,914203,1005055,1097894,1192605,1289089,1387253,1487016,
   1588303,1691047,1795186,1900662,2007423,2115421,2224611,2334950,2446401,
@@ -36,10 +34,12 @@ namespace Maths
 {
 
 // Log factorial
-S1615 LogFact(uint32_t n)
+S1615 LogFact(unsigned int n)
 {
   if (n < 64) // Lookup table for argument values 0-63
-    return logFact[n];  
+  {
+    return g_LogFact[n];
+  }
   else
   {
     S1615 n1 = (S1615)(n << 15);

rig_cpp_common/runtime/rig_cpp_common/maths/logfact.h

@@ -1,7 +1,7 @@
 #pragma once
 
 // Common includes
-#include "rig_cpp_common/fixed_point_number.h"
+#include "../fixed_point_number.h"
 
 // Namespaces
 using namespace Common::FixedPointNumber;
@@ -13,6 +13,6 @@ namespace Common
 {
 namespace Maths
 {
-  S1615 LogFact(uint32_t n);
+  S1615 LogFact(unsigned int n);
 } // Maths
 } // Common

rig_cpp_common/runtime/rig_cpp_common/maths/recip.cpp

@@ -1,7 +1,3 @@
-// Common includes
-#include "rig_cpp_common/fixed_point_number.h"
-
-// Maths includes
 #include "recip.h"
 
 // Namespaces
@@ -16,7 +12,7 @@ namespace
 
 // Reciprocal lookup table for values between 0.1 and 1.1
 // (2**15/np.linspace(0.1,1.1,65)).astype(int)[:-1]
-const S1615 recips[64] = {
+const S1615 g_Reciprocals[64] = {
   327680, 283398, 249660, 223101, 201649, 183960, 169125, 156503,
   145635, 136178, 127875, 120525, 113975, 108100, 102801,  97997,
   93622,  89621,  85948,  82565,  79437,  76538,  73843,  71331,
@@ -29,12 +25,12 @@ const S1615 recips[64] = {
 
 // Single iteration of Runge-Kutta method, starting at the closest
 // x value in the lookup table
-S1615 reciprocal_core(S1615 x)
+S1615 ReciprocalCore(S1615 x)
 {
   int32_t i0 = (x - 3276) >> 9; // Closest index in table
   S1615 x0 = (i0 << 9) + 3276;  // Corresponding closest x
 
-  S1615 y = recips[i0]; // Closest y
+  S1615 y = g_Reciprocals[i0]; // Closest y
   S1615 h = (x - x0);   // Step size from x0 to x
 
   S1615 k1 = -MulS1615(y, y);
@@ -68,28 +64,28 @@ S1615 Reciprocal(S1615 x)
 {
   // Record the sign, and operate on abs(x)
   int32_t sign = 1;
-  unsigned int left_shift = 0;
-  unsigned int right_shift = 0;
   if (x < 0)
   {
     x = -x;
     sign = -1;
   }
 
   // Shift until x lies in the range [0.1,1.1]
+  unsigned int right_shift = 0;
   while (x >= 36044)
   {
     x = x >> 1;
     right_shift += 1;
   }
+  unsigned int left_shift = 0;
   while (x < 3276)
   {
     x = x << 1;
     left_shift += 1;
   }
 
   // Get the reciprocal, unshift, multiply by the sign
-  return sign * ((reciprocal_core(x) << left_shift) >> right_shift);
+  return sign * ((ReciprocalCore(x) << left_shift) >> right_shift);
 }
 
 } // Maths

rig_cpp_common/runtime/rig_cpp_common/maths/recip.h

@@ -1,7 +1,7 @@
 #pragma once
 
 // Common includes
-#include "rig_cpp_common/fixed_point_number.h"
+#include "../fixed_point_number.h"
 
 // Namespaces
 using namespace Common::FixedPointNumber;