From 6fce6af60e7c99e77fcc24084db001367a7813be Mon Sep 17 00:00:00 2001 From: Daniel Friesel Date: Thu, 8 Oct 2020 16:16:12 +0200 Subject: BA Janis import: working detection of optimal number of changepoints --- bin/analyze-archive.py | 100 ++++--------------- lib/loader.py | 17 +++- lib/model.py | 25 ++++- lib/pelt.py | 261 +++++++++++++++++++++++++++++++++++++++++++++++++ lib/utils.py | 2 + 5 files changed, 321 insertions(+), 84 deletions(-) create mode 100644 lib/pelt.py 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 -- cgit v1.2.3