add SafeSum C function

Author: paulbkoch

python/interpret-core/interpret/utils/_native.py

@@ -10,6 +10,7 @@
 import sys
 import math
 from contextlib import AbstractContextManager
+from math import prod
 
 import numpy as np
 
@@ -239,6 +240,25 @@ def clean_float(self, val):
         )
         return val_array[0]
 
+    def safe_sum(self, in_tensor, out_tensor, axis):
+        shape = in_tensor.shape
+        n_distant = prod(shape[:axis])
+        n_close = prod(shape[axis + 1 :])
+
+        if prod(out_tensor.shape) != n_distant * n_close:
+            msg = f"in {in_tensor.shape} and out {out_tensor.shape} tensors must have a reducible shape along axis {axis}."
+            raise Exception(msg)
+
+        return_code = self._unsafe.SafeSum(
+            n_distant,
+            shape[axis],
+            n_close,
+            Native._make_pointer(in_tensor, np.float64, None, False),
+            Native._make_pointer(out_tensor, np.float64, None, False),
+        )
+        if return_code:  # pragma: no cover
+            raise Native._get_native_exception(return_code, "SafeSum")
+
     def flat_mean(self, vals, weights=None):
         if weights is not None:
             if vals.shape != weights.shape:
@@ -1034,6 +1054,20 @@ def _initialize(self, is_debug):
         ]
         self._unsafe.CleanFloats.restype = None
 
+        self._unsafe.SafeSum.argtypes = [
+            # int64_t countDistant
+            ct.c_int64,
+            # int64_t countAxis
+            ct.c_int64,
+            # int64_t countClose
+            ct.c_int64,
+            # double * in
+            ct.c_void_p,
+            # double * out
+            ct.c_void_p,
+        ]
+        self._unsafe.SafeSum.restype = ct.c_int32
+
         self._unsafe.SafeMean.argtypes = [
             # int64_t countBags
             ct.c_int64,

shared/libebm/inc/libebm.h

@@ -314,6 +314,9 @@ EBM_API_INCLUDE void EBM_CALLING_CONVENTION SetTraceLevel(TraceEbm traceLevel);
 EBM_API_INCLUDE const char* EBM_CALLING_CONVENTION GetTraceLevelString(TraceEbm traceLevel);
 
 EBM_API_INCLUDE void EBM_CALLING_CONVENTION CleanFloats(IntEbm count, double* valsInOut);
+
+EBM_API_INCLUDE ErrorEbm EBM_CALLING_CONVENTION SafeSum(
+      IntEbm countDistant, IntEbm countAxis, IntEbm countClose, const double* in, double* out);
 EBM_API_INCLUDE ErrorEbm EBM_CALLING_CONVENTION SafeMean(
       IntEbm countBags, IntEbm countTensorBins, const double* vals, const double* weights, double* tensorOut);
 EBM_API_INCLUDE ErrorEbm EBM_CALLING_CONVENTION SafeStandardDeviation(

shared/libebm/interpretable_numerics.cpp

@@ -1480,6 +1480,123 @@ static double Mean(const size_t cSamples,
    return mean;
 }
 
+// we don't care if an extra log message is outputted due to the non-atomic nature of the decrement to this value
+static int g_cLogEnterSafeSumCount = 25;
+static int g_cLogExitSafeSumCount = 25;
+
+EBM_API_BODY ErrorEbm EBM_CALLING_CONVENTION SafeSum(
+      IntEbm countDistant, IntEbm countAxis, IntEbm countClose, const double* in, double* out) {
+
+   LOG_COUNTED_N(&g_cLogEnterSafeSumCount,
+         Trace_Info,
+         Trace_Verbose,
+         "Entered SafeSum: "
+         "countDistant=%" IntEbmPrintf ", "
+         "countAxis=%" IntEbmPrintf ", "
+         "countClose=%" IntEbmPrintf ", "
+         "in=%p, "
+         "out=%p",
+         countDistant,
+         countAxis,
+         countClose,
+         static_cast<const void*>(in),
+         static_cast<const void*>(out));
+
+   if(nullptr == in) {
+      LOG_0(Trace_Error, "ERROR SafeSum nullptr == in");
+      return Error_IllegalParamVal;
+   }
+
+   if(nullptr == out) {
+      LOG_0(Trace_Error, "ERROR SafeSum nullptr == out");
+      return Error_IllegalParamVal;
+   }
+
+   if(countDistant <= IntEbm{0}) {
+      if(countDistant < IntEbm{0}) {
+         LOG_0(Trace_Error, "ERROR SafeSum countDistant < IntEbm{0}");
+         return Error_IllegalParamVal;
+      }
+      return Error_None;
+   }
+   if(IsConvertError<size_t>(countDistant)) {
+      LOG_0(Trace_Error, "ERROR SafeSum IsConvertError<size_t>(countDistant)");
+      return Error_IllegalParamVal;
+   }
+   const size_t cDistant = static_cast<size_t>(countDistant);
+
+   if(countClose <= IntEbm{0}) {
+      if(countClose < IntEbm{0}) {
+         LOG_0(Trace_Error, "ERROR SafeSum countClose < IntEbm{0}");
+         return Error_IllegalParamVal;
+      }
+      return Error_None;
+   }
+   if(IsConvertError<size_t>(countClose)) {
+      LOG_0(Trace_Error, "ERROR SafeSum IsConvertError<size_t>(countClose)");
+      return Error_IllegalParamVal;
+   }
+   const size_t cClose = static_cast<size_t>(countClose);
+
+   if(IsMultiplyError(sizeof(double), cClose)) {
+      LOG_0(Trace_Error, "ERROR SafeSum IsMultiplyError(sizeof(double), cClose)");
+      return Error_IllegalParamVal;
+   }
+   const size_t cCloseBytes = sizeof(double) * cClose;
+
+   if(IsMultiplyError(cCloseBytes, cDistant)) {
+      LOG_0(Trace_Error, "ERROR SafeSum IsMultiplyError(cCloseBytes, cDistant)");
+      return Error_IllegalParamVal;
+   }
+   const size_t cNonAxisBytes = cCloseBytes * cDistant;
+
+   if(countAxis <= IntEbm{1}) {
+      if(IntEbm{1} == countAxis) {
+         memcpy(out, in, cNonAxisBytes);
+      } else if(countAxis < IntEbm{0}) {
+         LOG_0(Trace_Error, "ERROR SafeSum countAxis < IntEbm{0}");
+         return Error_IllegalParamVal;
+      }
+      return Error_None;
+   }
+   if(IsConvertError<size_t>(countAxis)) {
+      LOG_0(Trace_Error, "ERROR SafeSum IsConvertError<size_t>(countAxis)");
+      return Error_IllegalParamVal;
+   }
+   const size_t cAxis = static_cast<size_t>(countAxis);
+
+   if(IsMultiplyError(cNonAxisBytes, cAxis)) {
+      LOG_0(Trace_Error, "ERROR SafeSum IsMultiplyError(cNonAxisBytes, cAxis)");
+      return Error_IllegalParamVal;
+   }
+
+   const double* const pOutEnd = IndexByte(out, cNonAxisBytes);
+   const size_t cAxisBytes = cCloseBytes * cAxis;
+   const size_t cNext = sizeof(*in) - cAxisBytes;
+   const size_t cNextGrouping = cAxisBytes - cCloseBytes;
+   size_t i = 0;
+   do {
+      const double* const pCloseEnd = IndexByte(out, cCloseBytes);
+      do {
+         const size_t iEnd = i + cAxisBytes;
+         double sum = *IndexByte(in, i);
+         i += cCloseBytes;
+         do {
+            sum += *IndexByte(in, i);
+            i += cCloseBytes;
+         } while(iEnd != i);
+         i += cNext;
+         *out = sum;
+         ++out;
+      } while(pCloseEnd != out);
+      i += cNextGrouping;
+   } while(pOutEnd != out);
+
+   LOG_COUNTED_0(&g_cLogExitSafeSumCount, Trace_Info, Trace_Verbose, "Exited SafeSum");
+
+   return Error_None;
+}
+
 // we don't care if an extra log message is outputted due to the non-atomic nature of the decrement to this value
 static int g_cLogEnterSafeMeanCount = 25;
 static int g_cLogExitSafeMeanCount = 25;

shared/libebm/libebm_exports.def

@@ -5,6 +5,7 @@ EXPORTS
   SetTraceLevel
   GetTraceLevelString
   CleanFloats
+  SafeSum
   SafeMean
   SafeStandardDeviation
   MeasureRNG

shared/libebm/libebm_exports.txt

@@ -4,6 +4,7 @@
       SetTraceLevel;
       GetTraceLevelString;
       CleanFloats;
+      SafeSum;
       SafeMean;
       SafeStandardDeviation;
       MeasureRNG;

shared/libebm/tests/SuggestGraphBoundsTest.cpp

@@ -302,6 +302,50 @@ TEST_CASE("SuggestGraphBounds, 2 cuts, overflow diff") {
    CHECK(std::numeric_limits<double>::infinity() == highGraphBound);
 }
 
+TEST_CASE("SafeSum, 4x3x2") {
+   const double in[]{1.0,
+         2.0,
+         3.0,
+         4.0,
+         5.0,
+         6.0,
+         7.0,
+         8.0,
+         9.0,
+         10.0,
+         11.0,
+         12.0,
+         13.0,
+         14.0,
+         15.0,
+         16.0,
+         17.0,
+         18.0,
+         19.0,
+         20.0,
+         21.0,
+         22.0,
+         23.0,
+         24.0};
+
+   double out[sizeof(in) / sizeof(in[0]) / 3];
+
+   const ErrorEbm error = SafeSum(4, 3, 2, in, out);
+   CHECK(Error_None == error);
+
+   CHECK(out[0] == 1.0 + 3.0 + 5.0);
+   CHECK(out[1] == 2.0 + 4.0 + 6.0);
+
+   CHECK(out[2] == 7.0 + 9.0 + 11.0);
+   CHECK(out[3] == 8.0 + 10.0 + 12.0);
+
+   CHECK(out[4] == 13.0 + 15.0 + 17.0);
+   CHECK(out[5] == 14.0 + 16.0 + 18.0);
+
+   CHECK(out[6] == 19.0 + 21.0 + 23.0);
+   CHECK(out[7] == 20.0 + 22.0 + 24.0);
+}
+
 TEST_CASE("SafeMean, 4 values") {
    double vals[]{1.0, 2.5, 10, 100};
    const size_t cVals = sizeof(vals) / sizeof(vals[0]);

shared/libebm/tests/boosting_unusual_inputs.cpp

@@ -2363,7 +2363,7 @@ static double RandomizedTesting(const AccelerationFlags acceleration) {
             }
          } else {
             double modelSum = 0.0;
-            for(IntEbm iTerm = 0; iTerm < static_cast<IntEbm>(terms.size()); ++iTerm) {
+            for(size_t iTerm = 0; iTerm < terms.size(); ++iTerm) {
                const auto term = terms[iTerm];
                size_t cScores = 3 <= classesCount ? static_cast<size_t>(classesCount) : size_t{1};
                for(size_t iDim = 0; iDim < term.size(); ++iDim) {
@@ -2372,12 +2372,12 @@ static double RandomizedTesting(const AccelerationFlags acceleration) {
                }
 
                std::vector<double> model(cScores);
-               test.GetBestTermScoresRaw(iTerm, &model[0]);
+               test.GetBestTermScoresRaw(static_cast<IntEbm>(iTerm), &model[0]);
                for(size_t iScore = 0; iScore < cScores; ++iScore) {
                   modelSum += model[iScore];
                }
 
-               test.GetCurrentTermScoresRaw(iTerm, &model[0]);
+               test.GetCurrentTermScoresRaw(static_cast<IntEbm>(iTerm), &model[0]);
                for(size_t iScore = 0; iScore < cScores; ++iScore) {
                   modelSum += model[iScore];
                }