Upload Troubled_Cell_Detector.py (as of March 17 2020)

Troubled_Cell_Detector.py

+# -*- coding: utf-8 -*-
+"""
+@author: Laura C. Kühle
+
+TODO: Check cutoff_factor/whether _is_troubled_cell works correctly \
+    -> Discuss! -> Seems alright; check paper (Siegfried Müller, Nils Gerhard?)
+"""
+
+import numpy as np
+
+class TroubledCellDetector(object):
+
+    def __init__(self, config):
+        pass
+
+    def get_name(self):
+        return self.function_name
+
+    def get_cells(self, multiwavelet_coeffs, projection):
+        pass
+
+
+class NoDetection(TroubledCellDetector):
+
+    def __init__(self, config):
+        # Set name of function
+        self.function_name = 'NoDetection'
+        pass
+
+    def get_cells(self, multiwavelet_coeffs, projection):
+        return []
+
+
+class Boxplot(TroubledCellDetector):
+
+    def __init__(self, config):
+        # Set name of function
+        self.function_name = 'Boxplot'
+
+        # Unpack necessary configurations
+        self.num_grid_cells = config.pop('num_grid_cells')
+        self.fold_len = config.pop('fold_len', 16)
+        self.whisker_len = config.pop('whisker_len', 3)
+
+        # Pop unnecessary parameter
+        config.pop('polynom_degree')
+        config.pop('cell_len')
+
+    def get_cells(self, multiwavelet_coeffs, projection):
+        self.num_grid_cells = len(multiwavelet_coeffs[0])
+        indexed_coeffs = [[multiwavelet_coeffs[0, i], i]
+                          for i in range(self.num_grid_cells)]
+
+        if self.num_grid_cells < self.fold_len:
+            self.fold_len = self.num_grid_cells
+
+        num_folds = int(self.num_grid_cells/self.fold_len)
+        troubled_cells = []
+
+        for fold in range(num_folds):
+            sorted_fold = sorted(indexed_coeffs[fold * self.fold_len:
+                                                (fold+1) * self.fold_len])
+
+            # fold_len/4.0 had comment: int((P+3)/2)/2.0
+            balance_factor = self.fold_len/4.0 % 1  # former: int(...)
+            boundary_index = int(self.fold_len/4.0 - balance_factor)
+
+            first_quartil = (1-balance_factor)\
+                * sorted_fold[boundary_index-1][0]\
+                + balance_factor * sorted_fold[boundary_index][0]
+            third_quartil = (1-balance_factor)\
+                * sorted_fold[3*boundary_index-1][0]\
+                + balance_factor * sorted_fold[3*boundary_index][0]
+
+            lower_bound = first_quartil\
+                - self.whisker_len * (third_quartil-first_quartil)
+            upper_bound = third_quartil\
+                + self.whisker_len * (third_quartil-first_quartil)
+
+            # Check for lower extreme outliers and add respective cells
+            for cell in sorted_fold:
+                if cell[0] < lower_bound:
+                    troubled_cells.append(cell[1])
+                else:
+                    break
+
+            # Check for lower extreme outliers and add respective cells
+            for cell in sorted_fold[::-1][:]:
+                if cell[0] > upper_bound:
+                    troubled_cells.append(cell[1])
+                else:
+                    break
+
+        return sorted(troubled_cells)
+
+
+class Theoretical(TroubledCellDetector):
+
+    def __init__(self, config):
+        # Set name of function
+        self.function_name = 'Theoretical'
+
+        # Unpack necessary configurations
+        self.polynom_degree = config.pop('polynom_degree')
+        self.num_grid_cells = config.pop('num_grid_cells')
+        self.cell_len = config.pop('cell_len')
+        self.cutoff_factor = config.pop('cutoff_factor',
+                                        np.sqrt(2) * self.cell_len)
+        # comment to line above: or 2 or 3
+
+        self.multiwavelet_coeffs = []
+        self.projection = []
+        self.max_avg = None
+
+    def get_cells(self, multiwavelet_coeffs, projection):
+        self.multiwavelet_coeffs = multiwavelet_coeffs
+        self.projection = projection
+        troubled_cells = []
+        self.max_avg = max(1, max(abs(self.projection[0][degree])
+                                  for degree in range(self.polynom_degree+1)))
+
+        for cell in range(int(self.num_grid_cells/2)):
+            if self._is_troubled_cell(cell):
+                troubled_cells.append(cell)
+
+        return troubled_cells
+
+    def _is_troubled_cell(self, cell):
+        max_value = max(abs(self.multiwavelet_coeffs[degree][cell])
+                        for degree in range(self.polynom_degree+1)) \
+            / self.max_avg
+        eps = self.cutoff_factor / (self.cell_len*self.num_grid_cells)
+
+        if max_value > eps:
+            return True
+        return False