BA Janis import: working detection of optimal number of changepoints

5 files changed, 321 insertions, 84 deletions

diff --git a/bin/analyze-archive.py b/bin/analyze-archive.py
index e9694b4..1217be6 100755
--- a/bin/analyze-archive.py
+++ b/bin/analyze-archive.py
@@ -49,32 +49,6 @@ from dfatool.validation import CrossValidator
 from dfatool.utils import filter_aggregate_by_param
 from dfatool.automata import PTA
 
-### PELT
-import numpy as np
-
-# Very short benchmark yielded approx. 3 times the speed of solution not using sort
-# checks the percentiles if refinement is necessary
-def needs_refinement(signal, thresh):
-    sorted_signal = sorted(signal)
-    length_of_signal = len(signal)
-    percentile_size = int()
-    percentile_size = length_of_signal // 100
-    lower_percentile = sorted_signal[0:percentile_size]
-    upper_percentile = sorted_signal[
-        length_of_signal - percentile_size : length_of_signal
-    ]
-    lower_percentile_mean = np.mean(lower_percentile)
-    upper_percentile_mean = np.mean(upper_percentile)
-    median = np.median(sorted_signal)
-    dist = median - lower_percentile_mean
-    if dist > thresh:
-        return True
-    dist = upper_percentile_mean - median
-    if dist > thresh:
-        return True
-    return False
-
-### /PELT
 
 def print_model_quality(results):
     for state_or_tran in results.keys():
@@ -378,8 +352,9 @@ if __name__ == "__main__":
     )
     parser.add_argument(
         "--with-substates",
-        action="store_true",
-        help="Perform substate analysis"
+        metavar="PELT_CONFIG",
+        type=str,
+        help="Perform substate analysis",
     )
     parser.add_argument("measurement", nargs="+")
 
@@ -427,7 +402,11 @@ if __name__ == "__main__":
 
     raw_data = RawData(
         args.measurement,
-        with_traces=(args.export_traces is not None or args.plot_traces is not None or args.with_substates is not None),
+        with_traces=(
+            args.export_traces is not None
+            or args.plot_traces is not None
+            or args.with_substates is not None
+        ),
         skip_cache=args.no_cache,
     )
 
@@ -476,58 +455,13 @@ if __name__ == "__main__":
             with open(target, "w") as f:
                 json.dump(data, f)
 
-    if args.with_substates:
-        opt_refinement_thresh = 100
-        uw_per_sot = dict()
-        for trace in preprocessed_data:
-            for state_or_transition in trace["trace"]:
-                if state_or_transition["isa"] == "state":
-                    name = state_or_transition["name"]
-                    if name not in uw_per_sot:
-                        uw_per_sot[name] = list()
-                    for elem in state_or_transition["offline"]:
-                        elem["uW"] = list(elem["uW"])
-                    uw_per_sot[name].append(state_or_transition)
-        for name, configurations in uw_per_sot.items():
-            for num_config, measurements_by_config in enumerate(configurations):
-                logging.debug(
-                    "Looking at state '"
-                    + measurements_by_config["name"]
-                    + "' with params: "
-                    + str(measurements_by_config["parameter"])
-                    + "("
-                    + str(num_config + 1)
-                    + "/"
-                    + str(len(configurations))
-                    + ")"
-                )
-                num_needs_refine = 0
-                logging.debug("Checking if refinement is necessary...")
-                for measurement in measurements_by_config["offline"]:
-                    # loop through measurements of particular state
-                    # and check if state needs refinement
-                    signal = measurement["uW"]
-                    # mean = measurement['uW_mean']
-                    if needs_refinement(signal, opt_refinement_thresh):
-                        num_needs_refine = num_needs_refine + 1
-                if num_needs_refine == 0:
-                    logging.debug(
-                        "No refinement necessary for state '"
-                        + measurements_by_config["name"]
-                        + "' with params: "
-                        + str(measurements_by_config["parameter"])
-                    )
-                elif num_needs_refine < len(measurements_by_config["offline"]) / 2:
-                    logging.debug(
-                        "No refinement necessary for state '"
-                        + measurements_by_config["name"]
-                        + "' with params: "
-                        + str(measurements_by_config["parameter"])
-                    )
-                    logging.debug(
-                        "However this decision was not unanimously. This could hint at poor "
-                        "measurement quality."
-                    )
+    if args.with_substates is not None:
+        arg_dict = dict()
+        if args.with_substates != "":
+            for kv in args.with_substates.split(","):
+                k, v = kv.split("=")
+                arg_dict[k] = v
+        args.with_substates = arg_dict
 
     if args.plot_traces:
         traces = list()
@@ -576,6 +510,7 @@ if __name__ == "__main__":
         traces=preprocessed_data,
         function_override=function_override,
         pta=pta,
+        pelt=args.with_substates,
     )
 
     if xv_method:
@@ -698,6 +633,9 @@ if __name__ == "__main__":
         safe_functions_enabled=safe_functions_enabled
     )
 
+    if args.with_substates is not None:
+        substate_model, substate_info = model.get_substates()
+
     if "paramdetection" in show_models or "all" in show_models:
         for state in model.states_and_transitions():
             for attribute in model.attributes(state):
diff --git a/lib/loader.py b/lib/loader.py
index 57b3d30..0c3bac7 100644
--- a/lib/loader.py
+++ b/lib/loader.py
@@ -575,6 +575,7 @@ class RawData:
         # and self.traces_by_fileno[measurement['fileno']][*]['trace'][*]['offline_aggregates'] in place
         # (appends data from measurement['energy_trace'])
         # If measurement['expected_trace'] exists, it is edited in place instead
+        # "offline_aggregates" is the only data used later on by model.py's by_name / by_param dicts
         online_datapoints = []
         if "expected_trace" in measurement:
             traces = measurement["expected_trace"]
@@ -618,12 +619,14 @@ class RawData:
             if arg_support_enabled and "args" in online_trace_part:
                 paramvalues.extend(map(soft_cast_int, online_trace_part["args"]))
 
+            # TODO rename offline_aggregates to make it clear that this is what ends up in by_name / by_param and model.py
             if "offline_aggregates" not in online_trace_part:
                 online_trace_part["offline_attributes"] = [
                     "power",
                     "duration",
                     "energy",
                 ]
+                # this is what ends up in by_name / by_param and is used by model.py
                 online_trace_part["offline_aggregates"] = {
                     "power": [],
                     "duration": [],
@@ -639,6 +642,9 @@ class RawData:
                     online_trace_part["offline_aggregates"]["rel_energy_prev"] = []
                     online_trace_part["offline_aggregates"]["rel_energy_next"] = []
                     online_trace_part["offline_aggregates"]["timeout"] = []
+                elif "uW" in offline_trace_part:
+                    online_trace_part["offline_support"] = ["power_traces"]
+                    online_trace_part["offline_aggregates"]["power_traces"] = list()
 
             # Note: All state/transitions are 20us "too long" due to injected
             # active wait states. These are needed to work around MIMOSA's
@@ -670,6 +676,11 @@ class RawData:
                     offline_trace_part["timeout"]
                 )
 
+            if online_trace_part["isa"] == "state" and "uW" in offline_trace_part:
+                online_trace_part["offline_aggregates"]["power_traces"].append(
+                    offline_trace_part["uW"]
+                )
+
     def _merge_online_and_etlog(self, measurement):
         # Edits self.traces_by_fileno[measurement['fileno']][*]['trace'][*]['offline']
         # and self.traces_by_fileno[measurement['fileno']][*]['trace'][*]['offline_aggregates'] in place
@@ -1078,13 +1089,17 @@ def _add_trace_data_to_aggregate(aggregate, key, element):
                 "rel_energy_next",
             ]
             # Uncomment this line if you also want to analyze mean transition power
-            # aggrgate[key]['attributes'].append('power')
+            # aggregate[key]['attributes'].append('power')
             if "plan" in element and element["plan"]["level"] == "epilogue":
                 aggregate[key]["attributes"].insert(0, "timeout")
         attributes = aggregate[key]["attributes"].copy()
         for attribute in attributes:
             if attribute not in element["offline_aggregates"]:
                 aggregate[key]["attributes"].remove(attribute)
+        if "offline_support" in element:
+            aggregate[key]["supports"] = element["offline_support"]
+        else:
+            aggregate[key]["supports"] = list()
     for datakey, dataval in element["offline_aggregates"].items():
         aggregate[key][datakey].extend(dataval)
 
diff --git a/lib/model.py b/lib/model.py
index bb4a45b..d7a9bc9 100644
--- a/lib/model.py
+++ b/lib/model.py
@@ -375,7 +375,7 @@ def _num_args_from_by_name(by_name):
 
 
 class AnalyticModel:
-    u"""
+    """
     Parameter-aware analytic energy/data size/... model.
 
     Supports both static and parameter-based model attributes, and automatic detection of parameter-dependence.
@@ -663,7 +663,7 @@ class AnalyticModel:
 
 
 class PTAModel:
-    u"""
+    """
     Parameter-aware PTA-based energy model.
 
     Supports both static and parameter-based model attributes, and automatic detection of parameter-dependence.
@@ -704,6 +704,7 @@ class PTAModel:
         function_override={},
         use_corrcoef=False,
         pta=None,
+        pelt=None,
     ):
         """
         Prepare a new PTA energy model.
@@ -726,6 +727,7 @@ class PTAModel:
         use_corrcoef -- use correlation coefficient instead of stddev comparison
             to detect whether a model attribute depends on a parameter
         pta -- hardware model as `PTA` object
+        pelt -- perform sub-state detection via PELT and model sub-states as well. Requires traces to be set.
         """
         self.by_name = by_name
         self.by_param = by_name_to_by_param(by_name)
@@ -733,6 +735,12 @@ class PTAModel:
         self._num_args = arg_count
         self._use_corrcoef = use_corrcoef
         self.traces = traces
+        if traces is not None and pelt is not None:
+            from .pelt import PELT
+
+            self.pelt = PELT(**pelt)
+        else:
+            self.pelt = None
         self.stats = ParamStats(
             self.by_name,
             self.by_param,
@@ -927,6 +935,19 @@ class PTAModel:
 
         return model_getter, info_getter
 
+    def get_substates(self):
+        states = self.states()
+        for k in self.by_param.keys():
+            if k[0] == "TX":
+                if self.pelt.needs_refinement(self.by_param[k]["power_traces"]):
+                    print(f"PELT: {k} needs refinement")
+                    penalty = self.pelt.get_penalty_value(
+                        self.by_param[k]["power_traces"]
+                    )
+                    print(f"    penalty: {penalty}")
+
+        return None, None
+
     def to_json(self):
         static_model = self.get_static()
         static_quality = self.assess(static_model)
diff --git a/lib/pelt.py b/lib/pelt.py
new file mode 100644
index 0000000..ddc2324
--- /dev/null
+++ b/lib/pelt.py
@@ -0,0 +1,261 @@
+#!/usr/bin/env python3
+
+import numpy as np
+import ruptures
+from multiprocessing import Pool
+
+# returns the found changepoints by algo for the specific penalty pen.
+# algo should be the return value of Pelt(...).fit(signal)
+# Also puts a token in container q to let the progressmeter know the changepoints for penalty pen
+# have been calculated.
+# used for parallel calculation of changepoints vs penalty
+def _get_breakpoints(algo, pen):
+    return pen, len(algo.predict(pen=pen))
+
+
+def find_knee_point(data_x, data_y, S=1.0, curve="convex", direction="decreasing"):
+    kneedle = kneed.KneeLocator(data_x, data_y, S=S, curve=curve, direction=direction)
+    kneepoint = (kneedle.knee, kneedle.knee_y)
+    return kneepoint
+
+
+def norm_signal(signal, scaler=25):
+    max_val = max(signal)
+    normed_signal = np.zeros(shape=len(signal))
+    for i, signal_i in enumerate(signal):
+        normed_signal[i] = signal_i / max_val
+        normed_signal[i] = normed_signal[i] * scaler
+    return normed_signal
+
+
+class PELT:
+    def __init__(self, **kwargs):
+        # Defaults von Janis
+        self.jump = 1
+        self.refinement_threshold = 100
+        self.range_min = 0
+        self.range_max = 100
+        self.__dict__.update(kwargs)
+
+    # signals: a set of uW measurements belonging to a single parameter configuration (i.e., a single by_param entry)
+    def needs_refinement(self, signals):
+        count = 0
+        for signal in signals:
+            p1, median, p99 = np.percentile(signal, (1, 50, 99))
+
+            if median - p1 > self.refinement_threshold:
+                count += 1
+            elif p99 - median > self.refinement_threshold:
+                count += 1
+        refinement_ratio = count / len(signals)
+        return refinement_ratio > 0.3
+
+    def get_penalty_value(
+        self, signals, model="l1", min_dist=2, range_min=0, range_max=100, S=1.0
+    ):
+        # Janis macht hier noch kein norm_signal. Mit sieht es aber genau so brauchbar aus.
+        signal = norm_signal(signals[0])
+        algo = ruptures.Pelt(model=model, jump=self.jump, min_size=min_dist).fit(signal)
+        queue = list()
+        for i in range(range_min, range_max + 1):
+            queue.append((algo, i))
+        with Pool() as pool:
+            results = pool.starmap(_get_breakpoints, queue)
+        pen_val = [x[0] for x in results]
+        changepoint_counts = [x[1] for x in results]
+        # Scheint unnötig zu sein, da wir ohnehin plateau detection durchführen
+        # knee = find_knee_point(pen_val, changepoint_counts, S=S)
+        knee = (0,)
+
+        start_index = -1
+        end_index = -1
+        longest_start = -1
+        longest_end = -1
+        prev_val = -1
+        for i, num_bkpts in enumerate(changepoint_counts[knee[0] :]):
+            if num_bkpts != prev_val:
+                end_index = i - 1
+                if end_index - start_index > longest_end - longest_start:
+                    # currently found sequence is the longest found yet
+                    longest_start = start_index
+                    longest_end = end_index
+                start_index = i
+            if i == len(changepoint_counts[knee[0] :]) - 1:
+                # end sequence with last value
+                end_index = i
+                # # since it is not guaranteed that this is the end of the plateau, assume the mid
+                # # of the plateau was hit.
+                # size = end_index - start_index
+                # end_index = end_index + size
+                # However this is not the clean solution. Better if search interval is widened
+                # with range_min and range_max
+                if end_index - start_index > longest_end - longest_start:
+                    # last found sequence is the longest found yet
+                    longest_start = start_index
+                    longest_end = end_index
+                start_index = i
+            prev_val = num_bkpts
+        mid_of_plat = longest_start + (longest_end - longest_start) // 2
+        knee = (mid_of_plat + knee[0], changepoint_counts[mid_of_plat + knee[0]])
+
+        # modify knee according to options. Defaults to 1 * knee
+        knee = (knee[0] * 1, knee[1])
+        return knee
+
+    """
+    # calculates and returns the necessary penalty for signal. Parallel execution with num_processes many processes
+    # jump, min_dist are passed directly to PELT. S is directly passed to kneedle.
+    # pen_modifier is used as a factor on the resulting penalty.
+    # the interval [range_min, range_max] is used for searching.
+    def calculate_penalty_value(
+        signal,
+        model="l1",
+        jump=5,
+        min_dist=2,
+        range_min=0,
+        range_max=100,
+        S=1.0,
+        pen_modifier=None,
+        show_plots=False,
+    ):
+        # default params in Function
+        if model is None:
+            model = "l1"
+        if jump is None:
+            jump = 5
+        if min_dist is None:
+            min_dist = 2
+        if range_min is None:
+            range_min = 0
+        if range_max is None:
+            range_max = 50
+        if S is None:
+            S = 1.0
+        if pen_modifier is None:
+            pen_modifier = 1
+        # change point detection. best fit seemingly with l1. rbf prods. RuntimeErr for pen > 30
+        # https://ctruong.perso.math.cnrs.fr/ruptures-docs/build/html/costs/index.html
+        # model = "l1"   #"l1"  # "l2", "rbf"
+        algo = ruptures.Pelt(model=model, jump=jump, min_size=min_dist).fit(signal)
+
+        ### CALC BKPS WITH DIFF PENALTYS
+        if range_max != range_min:
+            # building args array for parallelizing
+            args = []
+            # for displaying progression
+            m = Manager()
+            q = m.Queue()
+
+            for i in range(range_min, range_max + 1):
+                # same calculation for all except other penalty
+                args.append((algo, i, q))
+
+            print_info("starting kneepoint calculation.")
+            # init Pool with num_proesses
+            with Pool() as p:
+                # collect results from pool
+                result = p.starmap_async(get_bkps, args)
+                # monitor loop
+                percentage = -100  # Force display of 0%
+                i = 0
+                while True:
+                    if result.ready():
+                        break
+
+                    size = q.qsize()
+                    last_percentage = percentage
+                    percentage = round(size / (range_max - range_min) * 100, 2)
+                    if percentage >= last_percentage + 2 or i >= refresh_thresh:
+                        print_info("Current progress: " + str(percentage) + "%")
+                        i = 0
+                    else:
+                        i += 1
+                    time.sleep(refresh_delay)
+                res = result.get()
+            print_info("Finished kneepoint calculation.")
+            # DECIDE WHICH PENALTY VALUE TO CHOOSE ACCORDING TO ELBOW/KNEE APPROACH
+            # split x and y coords to pass to kneedle
+            pen_val = [x[0] for x in res]
+            fitted_bkps_val = [x[1] for x in res]
+            # # plot to look at res
+            knee = find_knee_point(pen_val, fitted_bkps_val, S=S)
+
+            # TODO: Find plateau on pen_val vs fitted_bkps_val
+            #   scipy.find_peaks() does not find plateaus if they extend through the end of the data.
+            #   to counter that, add one extremely large value to the right side of the data
+            #   after negating it is extremely small -> Almost certainly smaller than the
+            #   found plateau therefore the plateau does not extend through the border
+            #   -> scipy.find_peaks finds it. Choose value from within that plateau.
+            # fitted_bkps_val.append(100000000)
+            # TODO: Approaching over find_peaks might not work if the initial decrease step to the
+            #   "correct" number of changepoints and additional decrease steps e.g. underfitting
+            #   take place within the given penalty interval. find_peak only finds plateaus
+            #   of peaks. If the number of chpts decreases after the wanted plateau the condition
+            #   for local peaks is not satisfied anymore. Therefore this approach will only work
+            #   if the plateau extends over the right border of the penalty interval.
+            # peaks, peak_plateaus = find_peaks(- np.array(fitted_bkps_val), plateau_size=1)
+            # Since the data is monotonously decreasing only one plateau can be found.
+
+            # assuming the plateau is constant, i.e. no noise. OK to assume this here, since num_bkpts
+            # is monotonously decreasing. If the number of bkpts decreases inside a considered
+            # plateau, it means that the stable configuration is not yet met. -> Search further
+            start_index = -1
+            end_index = -1
+            longest_start = -1
+            longest_end = -1
+            prev_val = -1
+            for i, num_bkpts in enumerate(fitted_bkps_val[knee[0] :]):
+                if num_bkpts != prev_val:
+                    end_index = i - 1
+                    if end_index - start_index > longest_end - longest_start:
+                        # currently found sequence is the longest found yet
+                        longest_start = start_index
+                        longest_end = end_index
+                    start_index = i
+                if i == len(fitted_bkps_val[knee[0] :]) - 1:
+                    # end sequence with last value
+                    end_index = i
+                    # # since it is not guaranteed that this is the end of the plateau, assume the mid
+                    # # of the plateau was hit.
+                    # size = end_index - start_index
+                    # end_index = end_index + size
+                    # However this is not the clean solution. Better if search interval is widened
+                    # with range_min and range_max
+                    if end_index - start_index > longest_end - longest_start:
+                        # last found sequence is the longest found yet
+                        longest_start = start_index
+                        longest_end = end_index
+                    start_index = i
+                prev_val = num_bkpts
+            if show_plots:
+                plt.xlabel("Penalty")
+                plt.ylabel("Number of Changepoints")
+                plt.plot(pen_val, fitted_bkps_val)
+                plt.vlines(
+                    longest_start + knee[0], 0, max(fitted_bkps_val), linestyles="dashed"
+                )
+                plt.vlines(
+                    longest_end + knee[0], 0, max(fitted_bkps_val), linestyles="dashed"
+                )
+                plt.show()
+            # choosing pen from plateau
+            mid_of_plat = longest_start + (longest_end - longest_start) // 2
+            knee = (mid_of_plat + knee[0], fitted_bkps_val[mid_of_plat + knee[0]])
+
+            # modify knee according to options. Defaults to 1 * knee
+            knee = (knee[0] * pen_modifier, knee[1])
+
+        else:
+            # range_min == range_max. has the same effect as pen_override
+            knee = (range_min, None)
+        print_info(str(knee[0]) + " has been selected as penalty.")
+        if knee[0] is not None:
+            return knee
+
+        print_error(
+            "With the current thresh-hold S="
+            + str(S)
+            + " it is not possible to select a penalty value."
+        )
+        sys.exit(-1)
+    """
diff --git a/lib/utils.py b/lib/utils.py
index d28ecda..9ee55de 100644
--- a/lib/utils.py
+++ b/lib/utils.py
@@ -180,6 +180,8 @@ def by_name_to_by_param(by_name: dict):
                     by_param[param_key]["isa"] = by_name[name]["isa"]
             for attribute in by_name[name]["attributes"]:
                 by_param[param_key][attribute].append(by_name[name][attribute][i])
+            for support in by_name[name]["supports"]:
+                by_param[param_key][support].append(by_name[name][support][i])
             # Required for match_parameter_valuse in _try_fits
             by_param[param_key]["param"].append(by_name[name]["param"][i])
     return by_param