From 42be5ac36fe201663c8d11bb90ca083f3dbc08ce Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Max=20L=C3=BCbke?= <mluebke@uni-potsdam.de>
Date: Fri, 20 Dec 2024 07:33:00 +0100
Subject: [PATCH] fix: handle zero distance case in inverseDistanceWeighting to
 prevent division by zero

---
 src/Chemistry/ChemistryModule.cpp | 59 ++-----------------------------
 1 file changed, 3 insertions(+), 56 deletions(-)

diff --git a/src/Chemistry/ChemistryModule.cpp b/src/Chemistry/ChemistryModule.cpp
index c95796b79..e6572c092 100644
--- a/src/Chemistry/ChemistryModule.cpp
+++ b/src/Chemistry/ChemistryModule.cpp
@@ -66,7 +66,8 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
       distance += std::pow(
           rescaled[key_comp_i][point_i] - rescaled[key_comp_i][data_set_n], 2);
     }
-    weights[point_i] = 1 / std::sqrt(distance);
+
+    weights[point_i] = distance != 0 ? 1 / std::sqrt(distance) : 0;
     assert(!std::isnan(weights[point_i]));
     inv_sum += weights[point_i];
   }
@@ -97,63 +98,9 @@ inverseDistanceWeighting(const std::vector<std::int32_t> &to_calc,
     key_delta /= inv_sum;
 
     results[output_comp_i] = from[output_comp_i] + key_delta;
+    assert(!std::isnan(results[output_comp_i]));
   }
 
-  // if (!has_h) {
-  //   double new_val = 0;
-  //   for (int j = 0; j < data_set_n; j++) {
-  //     new_val += weights[j] * output[j][0];
-  //   }
-  //   results[0] = new_val / inv_sum;
-  // }
-
-  // if (!has_h) {
-  //   double new_val = 0;
-  //   for (int j = 0; j < data_set_n; j++) {
-  //     new_val += weights[j] * output[j][1];
-  //   }
-  //   results[1] = new_val / inv_sum;
-  // }
-
-  // for (std::uint32_t i = 0; i < to_calc.size(); i++) {
-  //   const std::uint32_t interp_i = to_calc[i];
-
-  //   // rescale input between 0 and 1
-  //   for (int j = 0; j < input.size(); j++) {
-  //     buffer[j] = input[j].at(i);
-  //   }
-
-  //   buffer[buffer_size - 1] = from[interp_i];
-
-  //   const double min = *std::min_element(buffer, buffer + buffer_size);
-  //   const double max = *std::max_element(buffer, buffer + buffer_size);
-
-  //   for (int j = 0; j < input.size(); j++) {
-  //     buffer[j] = ((max - min) != 0 ? (buffer[j] - min) / (max - min) : 1);
-  //   }
-  //   from_rescaled =
-  //       ((max - min) != 0 ? (from[interp_i] - min) / (max - min) : 0);
-
-  //   double inv_sum = 0;
-
-  //   // calculate distances for each point
-  //   for (int i = 0; i < input.size(); i++) {
-  //     const double distance = std::pow(buffer[i] - from_rescaled, 2);
-
-  //     buffer[i] = distance > 0 ? (1 / std::sqrt(distance)) : 0;
-  //     inv_sum += buffer[i];
-  //   }
-  //   // calculate new values
-  //   double new_val = 0;
-  //   for (int i = 0; i < output.size(); i++) {
-  //     new_val += buffer[i] * output[i][interp_i];
-  //   }
-  //   results[interp_i] = new_val / inv_sum;
-  //   if (std::isnan(results[interp_i])) {
-  //     std::cout << "nan with new_val = " << output[0][i] << std::endl;
-  //   }
-  // }
-
   return results;
 }