From b7d900afa11995f6e8b8472f698f28cc5575de50 Mon Sep 17 00:00:00 2001
From: David Hammer <dhammer@mailbox.org>
Date: Fri, 17 Sep 2021 16:27:51 +0200
Subject: [PATCH] Fix test, sanitize slopes_pc constant

---
 src/calng/agipd_gpu.py          | 28 +++++++++++++++++-----------
 src/calng/agipd_gpu_kernels.cpp | 17 +++++------------
 src/calng/base_correction.py    |  3 ++-
 src/calng/dssc_gpu_kernels.cpp  |  8 +++++---
 src/tests/test_dssc_kernels.py  | 10 +++++-----
 5 files changed, 34 insertions(+), 32 deletions(-)

diff --git a/src/calng/agipd_gpu.py b/src/calng/agipd_gpu.py
index 56a36ac1..8cfbdd0f 100644
--- a/src/calng/agipd_gpu.py
+++ b/src/calng/agipd_gpu.py
@@ -80,21 +80,27 @@ class AgipdGpuRunner(base_gpu.BaseGpuRunner):
         # pc has funny shape (11, 352, 128, 512) from file
         # this is (fi, memory cell, y, x)
         slopes_pc_map = slopes_pc_map.astype(np.float32)
-        pc_high_m = slopes_pc_map[0]
-        pc_high_I = slopes_pc_map[1]
-        pc_med_m = slopes_pc_map[3]
-        pc_med_I = slopes_pc_map[4]
-        frac_high_med = pc_high_m / pc_med_m
-        # TODO: handle NaN somehow?
-        if np.isnan(frac_high_med).any():
-            ...
-        # TODO: verify formula
-        md_additional_offset = (pc_high_I - pc_med_I * frac_high_med).astype(np.float32)
+        # the following may contain NaNs, though...
+        hg_slope = slopes_pc_map[0]
+        hg_intercept = slopes_pc_map[1]
+        mg_slope = slopes_pc_map[3]
+        mg_intercept = slopes_pc_map[4]
+        # from agipdlib.py: replace NaN with median (per memory cell)
+        # note: suffixes in agipdlib are "_m" and "_l", should probably be "_I"
+        for naughty_array in (hg_slope, hg_intercept, mg_slope, mg_intercept):
+            medians = np.nanmedian(naughty_array, axis=(1, 2))
+            nan_bool = np.isnan(naughty_array)
+            nan_cell, _, _ = np.where(nan_bool)
+            naughty_array[nan_bool] = medians[nan_cell]
+            # TODO: clamp values outside range (0.8 to 1.2) * median to median
+
+        frac_hg_mg = hg_slope / mg_slope
+        md_additional_offset = (hg_intercept - mg_intercept * frac_hg_mg).astype(np.float32)
         rel_gain_map = np.ones(
             (3, self.constant_memory_cells, self.pixels_y, self.pixels_x),
             dtype=np.float32,
         )
-        rel_gain_map[1] = rel_gain_map[0] * frac_high_med
+        rel_gain_map[1] = rel_gain_map[0] * frac_hg_mg
         rel_gain_map[2] = rel_gain_map[1] * 4.48
         # TODO: enable overriding this based on user input
         self.rel_gain_pc_map_gpu.set(np.transpose(rel_gain_map, (1, 3, 2, 0)))
diff --git a/src/calng/agipd_gpu_kernels.cpp b/src/calng/agipd_gpu_kernels.cpp
index e5fc623d..898e0b87 100644
--- a/src/calng/agipd_gpu_kernels.cpp
+++ b/src/calng/agipd_gpu_kernels.cpp
@@ -5,10 +5,11 @@
 
 extern "C" {
 	/*
-	  Perform correction: offset
-	  Take cell_table into account when getting correction values
-	  Converting to float for doing the correction
-	  Converting to output dtype at the end
+	  Perform corrections; see agipd_gpu.CorrectionFlags
+	  Note that THRESHOLD and OFFSET should for any later corrections to make sense
+	  Will take cell_table into account when getting correction values
+	  Will convert from input dtype to float for correction
+	  Will convert to output dtype for output
 	*/
 	__global__ void correct(const {{input_data_dtype}}* data,
 							const unsigned short* cell_table,
@@ -140,12 +141,4 @@ extern "C" {
 			{% endif %}
 		}
 	}
-
-	/*
-	  Same as correction, except don't do any correction
-	*/
-	__global__ void only_cast(const {{input_data_dtype}}* data,
-							  unsigned char* gain_map,
-							  {{output_data_dtype}}* output) {
-	}
 }
diff --git a/src/calng/base_correction.py b/src/calng/base_correction.py
index 6bdc4d37..a7076c48 100644
--- a/src/calng/base_correction.py
+++ b/src/calng/base_correction.py
@@ -553,7 +553,8 @@ class BaseCorrection(calibrationBase.CalibrationReceiverBaseDevice):
             enabled &= self._correction_flag_class.NONE
         self._correction_flag_enabled = enabled
         self._correction_flag_preview = preview
-        self.log.INFO(str(enabled))
+        self.log.INFO(f"Corrections for dataOutput: {str(enabled)}")
+        self.log.INFO(f"Corrections for preview: {str(preview)}")
 
     def _update_output_schema(self, data):
         """Updates the schema of dataOutput based on parameter data (a Hash)
diff --git a/src/calng/dssc_gpu_kernels.cpp b/src/calng/dssc_gpu_kernels.cpp
index b82159c0..a5bf8d8e 100644
--- a/src/calng/dssc_gpu_kernels.cpp
+++ b/src/calng/dssc_gpu_kernels.cpp
@@ -4,12 +4,14 @@
 
 extern "C" {
 	/*
-	  Perform correction: offset
+	  Perform corrections: NONE or OFFSET
 	  Take cell_table into account when getting correction values
-	  Converting to float for doing the correction
+	  Converting to float while correcting
 	  Converting to output dtype at the end
+	  Shape of input data: memory cell, 1, y, x
+	  Shape of offset constant: x, y, memory cell
 	*/
-	__global__ void correct(const {{input_data_dtype}}* data,
+	__global__ void correct(const {{input_data_dtype}}* data, // shape: memory cell, 1, y, x
 							const unsigned short* cell_table,
 	                        const unsigned char corr_flags,
 							const float* offset_map,
diff --git a/src/tests/test_dssc_kernels.py b/src/tests/test_dssc_kernels.py
index 02369c98..1ed78062 100644
--- a/src/tests/test_dssc_kernels.py
+++ b/src/tests/test_dssc_kernels.py
@@ -43,28 +43,28 @@ kernel_runner = dssc_gpu.DsscGpuRunner(
 
 def test_only_cast():
     kernel_runner.load_data(raw_data)
-    kernel_runner.only_cast()
+    kernel_runner.correct(dssc_gpu.CorrectionFlags.NONE)
     assert np.allclose(
         kernel_runner.processed_data_gpu.get(), raw_data.astype(output_dtype)
     )
 
 
 def test_correct():
-    kernel_runner.load_constants(offset_map)
+    kernel_runner.load_offset_map(offset_map)
     kernel_runner.load_data(raw_data)
     kernel_runner.load_cell_table(cell_table)
-    kernel_runner.correct()
+    kernel_runner.correct(dssc_gpu.CorrectionFlags.OFFSET)
     assert np.allclose(kernel_runner.processed_data_gpu.get(), corrected_data)
 
 
 def test_correct_oob_cells():
-    kernel_runner.load_constants(offset_map)
+    kernel_runner.load_offset_map(offset_map)
     kernel_runner.load_data(raw_data)
     # here, half the cell IDs will be out of bounds
     wild_cell_table = cell_table * 2
     kernel_runner.load_cell_table(wild_cell_table)
     # should not crash
-    kernel_runner.correct()
+    kernel_runner.correct(dssc_gpu.CorrectionFlags.OFFSET)
     # should correct as much as possible
     assert np.allclose(
         kernel_runner.processed_data_gpu.get(),
-- 
GitLab