initial commit

7 files changed, 1786 insertions, 0 deletions

diff --git a/bin/analyze.py b/bin/analyze.py
new file mode 100755
index 0000000..57803fe
--- /dev/null
+++ b/bin/analyze.py
@@ -0,0 +1,40 @@
+#!/usr/bin/env python3
+
+import json
+import numpy as np
+import os
+from scipy.cluster.vq import kmeans2
+import struct
+import sys
+import tarfile
+from dfatool import running_mean, MIMOSA
+
+voltage = float(sys.argv[1])
+shunt = float(sys.argv[2])
+filename = sys.argv[3]
+
+mim = MIMOSA(voltage, shunt)
+
+charges, triggers = mim.load_data(filename)
+trigidx = mim.trigger_edges(triggers)
+triggers = []
+cal_edges = mim.calibration_edges(running_mean(mim.currents_nocal(charges[0:trigidx[0]]), 10))
+calfunc, caldata = mim.calibration_function(charges, cal_edges)
+vcalfunc = np.vectorize(calfunc, otypes=[np.float64])
+
+json_out = {
+    'triggers' : len(trigidx),
+    'first_trig' : trigidx[0] * 10,
+    'calibration' : caldata,
+    'trace' : mim.analyze_states(charges, trigidx, vcalfunc)
+}
+
+basename, _ = os.path.splitext(filename)
+
+# TODO also look for interesting gradients inside each state
+
+with open(basename + ".json", "w") as f:
+    json.dump(json_out, f)
+    f.close()
+
+#print(kmeans2(charges[:firstidx], np.array([130 * ua_step, 3.6 / 987 * 1000000, 3.6 / 99300 * 1000000])))
diff --git a/bin/cal-hist b/bin/cal-hist
new file mode 100755
index 0000000..96c4469
--- /dev/null
+++ b/bin/cal-hist
@@ -0,0 +1,66 @@
+#!/usr/bin/env python3
+
+import csv
+import numpy as np
+import os
+import struct
+import sys
+import tarfile
+import matplotlib.pyplot as plt
+from dfatool import running_mean, MIMOSA
+
+voltage = float(sys.argv[1])
+shunt = float(sys.argv[2])
+mimfile = sys.argv[3]
+
+mim = MIMOSA(voltage, shunt)
+
+charges, triggers = mim.load_data(mimfile)
+trigidx = mim.trigger_edges(triggers)
+cal_edges = mim.calibration_edges(running_mean(mim.currents_nocal(charges[0:trigidx[0]]), 10))
+
+#charges = charges[charges > 20000]
+#charges = charges[charges < 21000]
+
+def show_hist(data):
+	bins = np.max(data) - np.min(data)
+	if bins == 0:
+		bins = 1
+	if bins > 1000:
+		bins = bins / 10
+	#bins = 500
+	n, bins, patches = plt.hist(data, bins, normed=0, facecolor='green', alpha=0.75)
+	plt.grid(True)
+	plt.show()
+	print(np.histogram(data, bins=bins))
+
+#show_hist(charges[cal_edges[0]:cal_edges[1]])
+#show_hist(charges[cal_edges[4]:cal_edges[5]])
+#show_hist(charges[cal_edges[2]:cal_edges[3]])
+#show_hist(charges[trigidx[7]:trigidx[8]])
+#show_hist(np.array(charges))
+
+#print(charges[cal_edges[0]:cal_edges[1]])
+#print(charges[cal_edges[4]:cal_edges[5]])
+#print(charges[cal_edges[2]:cal_edges[3]])
+
+plt.hist(mim.charge_to_current_nocal(charges[cal_edges[2]:cal_edges[3]]) * 1e-3, 100, normed=0, facecolor='blue', alpha=0.8)
+plt.xlabel('mA MimosaCMD')
+plt.ylabel('#')
+plt.grid(True)
+plt.show()
+plt.hist(mim.charge_to_current_nocal(charges[cal_edges[4]:cal_edges[5]]) * 1e-3, 100, normed=0, facecolor='blue', alpha=0.8)
+plt.xlabel('mA MimosaCMD')
+plt.ylabel('#')
+plt.grid(True)
+plt.show()
+plt.hist(mim.charge_to_current_nocal(charges[cal_edges[0]:cal_edges[1]]) * 1e-3, 100, normed=0, facecolor='blue', alpha=0.8)
+plt.xlabel('mA MimosaCMD')
+plt.ylabel('#')
+plt.grid(True)
+plt.show()
+plt.hist(charges[cal_edges[0]:cal_edges[1]], 100, normed=0, facecolor='blue', alpha=0.8)
+plt.xlabel('Rohwert MimosaCMD')
+plt.ylabel('#')
+plt.grid(True)
+plt.show()
diff --git a/bin/dfatool b/bin/dfatool
new file mode 100755
index 0000000..00b4196
--- /dev/null
+++ b/bin/dfatool
@@ -0,0 +1,366 @@
+#!/usr/bin/env perl
+
+use strict;
+use warnings;
+use 5.020;
+
+use Getopt::Long;
+use IO::Handle;
+use Kratos::DFADriver;
+use List::Util qw(max sum);
+use Time::Progress;
+
+our $VERSION = '0.00';
+my %opt;
+
+GetOptions(
+	\%opt,
+	qw{
+	  exclude-states=s@
+	  ignore-nested-calls
+	  logging!
+	  plot=s
+	  no-cache
+	  no-update
+	  state-duration=i
+	  shunt=f
+	  trace-filter=s@
+	  trace-revisit=i
+	  trigger-pin=i
+	  trigger-port=s
+	  voltage=f
+	  offset=i
+	  zomg-fasta-nao
+	  }
+);
+
+if ( @ARGV < 2 ) {
+	show_usage();
+}
+
+@{ $opt{'exclude-states'} }
+  = split( qr{,}, join( q{,}, @{ $opt{'exclude-states'} // [] } ) );
+
+my ( $command, $xml_file, @data_files ) = @ARGV;
+
+my $driver = Kratos::DFADriver->new(
+	cache => $opt{'no-cache'} ? 0 : 1,
+	data_file       => $data_files[0],
+	excluded_states => $opt{'exclude-states'} // [],
+	fast_analysis   => $opt{'zomg-fasta-nao'} // 0,
+	ignore_nested   => $opt{'ignore-nested-calls'} // 0,
+	logging         => $opt{logging} // 1,
+	state_duration  => $opt{'state-duration'} // 1000,
+	trigger_pin     => $opt{'trigger-pin'},
+	trigger_port    => $opt{'trigger-port'},
+	merge_args      => $opt{'plot'}
+	? [ map { "--$_" } split( qr{,}, $opt{'plot'} ) ]
+	: [],
+	mimosa_offset  => $opt{offset}          // 130,
+	mimosa_shunt   => $opt{shunt}           // 330,
+	mimosa_voltage => $opt{voltage}         // 3.60,
+	trace_filter   => $opt{'trace-filter'}  // [],
+	trace_revisit  => $opt{'trace-revisit'} // 2,
+	xml_file       => $xml_file,
+);
+
+my %action = (
+	clean => sub {
+		$driver->launchpad_log_clean;
+	},
+	enable => sub {
+		$driver->write_acc_files;
+	},
+	disable => sub {
+		$driver->rm_acc_files;
+	},
+	'to-tikz' => sub {
+		say $driver->to_tikz;
+	},
+	'to-dfa' => sub {
+		say join(
+			"\n",
+			map {
+				join( ' -> ', map { "$_->[0]($_->[1],$_->[2])" } @{$_} )
+			} $driver->runs
+		);
+	},
+	flash => sub {
+		say "Compiling Kratos and flashing Launchpad";
+		$driver->launchpad_flash;
+	},
+	maketest => sub {
+		$driver->write_test_files;
+	},
+	rmtest => sub {
+		$driver->rm_test_files;
+	},
+	reset => sub {
+		say "Resetting MSP430";
+		$driver->launchpad_reset;
+	},
+	log => sub {
+		say "Resetting MSP430";
+		$driver->launchpad_reset;
+		say "Connecting to Launchpad";
+		$driver->launchpad_log_init;
+		say "Starting measurement";
+		$driver->mimosa->start;
+		say "Calibrating MIMOSA";
+		$driver->mimosa->calibrate;
+		STDOUT->autoflush(1);
+		print "Waiting for sync";
+
+		while ( not $driver->launchpad_log_is_synced ) {
+			$driver->launchpad_log_read;
+			print q{.};
+		}
+		print "\r\e[2KSynced with DriverEval app\n";
+		my ( $iter, $run, $maxrun ) = $driver->launchpad_log_status;
+		my $timer = Time::Progress->new(
+			min => 0,
+			max => $maxrun
+		);
+		while ( $run < 1 ) {
+			$driver->launchpad_log_read;
+			( $iter, $run, $maxrun ) = $driver->launchpad_log_status;
+		}
+		while ( $driver->launchpad_log_is_synced ) {
+			$driver->launchpad_log_read;
+			( $iter, $run, $maxrun ) = $driver->launchpad_log_status;
+			print $timer->report(
+"\r\e[2K%40b %p (${run}/${maxrun})  %L elapsed  %E remaining  in iteration ${iter}",
+				$run
+			);
+			if ( $run == $maxrun ) {
+				printf( "\r\e[2KIteration %d done after %d seconds\n",
+					$iter - 1, $timer->elapsed($run) );
+				say "Stopping measurement";
+				$driver->mimosa->stop;
+				say "Archiving files";
+				$driver->archive_files;
+				return;
+			}
+
+			if ( my @errors = $driver->launchpad_get_errors ) {
+				say "\r\e[2KErrors in iteration ${iter}, run ${run}";
+				say join( "\n", @errors );
+				say "Aborting measurement. Current run will not be saved.";
+				$driver->mimosa->kill;
+				exit 1;
+			}
+		}
+	},
+	analyze => sub {
+		$driver->analyze(@data_files);
+		$driver->assess_model;
+		if ( not $opt{'no-update'} ) {
+			$driver->update_model;
+		}
+	},
+	'analyze-tex' => sub {
+		$driver->analyze(@data_files);
+		$driver->assess_model_tex;
+	},
+	analyzesingle => sub {
+		if ( $opt{'no-cache'} or not $driver->log->load_cache ) {
+			say "Analyzing DriverEval iterations (this may take a while)";
+			$driver->log->load_archive;
+			$driver->log->preprocess;
+			$driver->log->save_cache;
+		}
+		say "Processing";
+		$driver->log->analyze;
+		$driver->assess_model;
+		if ( not $opt{'no-update'} ) {
+			$driver->update_model;
+		}
+	},
+	validate => sub {
+		$driver->validate_model(@data_files);
+	},
+	crossvalidate => sub {
+		printf("./dfatool crossvalidate %s %s\n",
+			$xml_file,
+			join(q{ }, @data_files)
+		);
+		$driver->crossvalidate_model(@data_files);
+	},
+	ls => sub {
+		for my $file (@data_files) {
+			my $log = $driver->log($file);
+			$log->load_archive;
+			my $setup = $log->setup;
+			say $file;
+			printf(
+				"  %.2fV @ %3dΩ, %dms per state, max revisit %d\n",
+				$setup->{mimosa_voltage}, $setup->{mimosa_shunt},
+				$setup->{state_duration}, $setup->{trace_revisit}
+			);
+			if ( $setup->{excluded_states} and @{ $setup->{excluded_states} } )
+			{
+				printf( "  excluded states: %s\n", $setup->{excluded_states} );
+			}
+			if ( $setup->{trace_filter} and @{ $setup->{trace_filter} } ) {
+				printf( "  trace filter: %s\n",
+					join( q{ | }, @{ $setup->{trace_filter} } ) );
+			}
+		}
+	},
+	list => sub {
+		for my $file (@data_files) {
+			my $log = $driver->log($file);
+			if ( $opt{'no-cache'} or not $driver->log->load_cache ) {
+				$log->load_archive;
+				$log->preprocess;
+				$log->save_cache;
+			}
+			$log->analyze;
+			my $data  = $log->data;
+			my $setup = $log->setup;
+			say $file;
+			printf(
+				"  %.2fV @ %3dΩ, %dms per state, max revisit %d\n",
+				$setup->{mimosa_voltage}, $setup->{mimosa_shunt},
+				$setup->{state_duration}, $setup->{trace_revisit}
+			);
+			if ( $setup->{excluded_states} and @{ $setup->{excluded_states} } )
+			{
+				printf( "  excluded states: %s\n", $setup->{excluded_states} );
+			}
+			if ( $setup->{trace_filter} and @{ $setup->{trace_filter} } ) {
+				printf( "  trace filter: %s\n",
+					join( q{ | }, @{ $setup->{trace_filter} } ) );
+			}
+			printf( "  MIMOSA offset: %5s %5s %5s\n", 'inf', '100k', '1k' );
+			for my $cal ( @{ $data->{calibration} } ) {
+				printf( "                 %5.f %5.f %5.f µW\n",
+					@{$cal}{qw{r0_err_uW r2_err_uW r1_err_uW}},
+				);
+			}
+			for my $state ( sort keys %{ $data->{aggregate}{state} } ) {
+				if ( $state ne 'UNINITIALIZED' ) {
+					my $val = $data->{aggregate}{state}{$state};
+					printf(
+						"  %15s : %.f±%.f = %.f, clip %.f%% ^ %.f%%\n",
+						$state,                   $val->{power}{mean},
+						$val->{power}{std_inner}, $val->{power}{median},
+						$val->{clip}{mean},       $val->{clip}{max}
+					);
+				}
+			}
+		}
+	},
+	show => sub {
+		for my $file (@data_files) {
+			my $log = $driver->log($file);
+			$log->load_archive;
+			if ( $opt{'no-cache'} or not $driver->log->load_cache ) {
+				$log->load_archive;
+				$log->preprocess;
+				$log->save_cache;
+			}
+			$log->analyze;
+			for my $trace ( @{ $log->data->{traces} } ) {
+				my ( @data, @widths );
+				printf( '%3d', $trace->{id} );
+				for my $elem ( @{ $trace->{trace} } ) {
+					if (    $elem->{isa} eq 'state'
+						and $elem->{name} ne 'UNINITIALIZED' )
+					{
+						push( @widths, max (length( $elem->{name} ), 9) );
+						printf(' → %9s', $elem->{name});
+						my @powers
+						  = map { $_->{uW_mean} } @{ $elem->{offline} };
+						push( @data, sum(@powers) / @powers );
+					}
+					elsif ( $elem->{isa} eq 'transition' ) {
+						my $args = join( q{, }, @{ $elem->{args} // [qw[?]] } );
+						my $pstr = "$elem->{name}($args)";
+						push( @widths, max (length($pstr), 9) );
+						printf(' → %9s', $pstr);
+						my @energies
+						  = map { $_->{uW_mean_delta} * ( $_->{us} - 20 ) }
+						  @{ $elem->{offline} };
+						push( @data, sum(@energies) / @energies );
+					}
+				}
+				print "\n   ";
+				for my $i ( 0 .. $#data ) {
+					printf( " → \%$widths[$i]d", $data[$i] );
+				}
+				print "\n";
+			}
+		}
+	},
+	reset => sub {
+		$driver->reset_model;
+	},
+);
+
+$SIG{INT} = $SIG{TERM} = sub {
+	say "\r\e[2KTerminating MIMOSA daemon";
+	$driver->mimosa->kill;
+	say "Goodbye";
+	exit 0;
+};
+
+sub show_usage {
+	say STDERR "Usage: $0 <action> <DFA driver XML file>";
+	say STDERR 'Supported actions: ' . join( q{ }, sort keys %action );
+	exit 1;
+}
+
+if ( exists $action{$command} ) {
+	$action{$command}();
+}
+elsif ( $command eq 'loop' ) {
+	$action{clean}();
+	$action{enable}();
+	$action{maketest}();
+	$action{flash}();
+	while (1) {
+		$action{log}();
+	}
+}
+else {
+	show_usage();
+}
+
+__END__
+
+=head1 NAME
+
+=head1 SYNOPSIS
+
+=head1 VERSION
+
+=head1 DESCRIPTION
+
+=head1 OPTIONS
+
+=over
+
+=back
+
+=head1 EXIT STATUS
+
+=head1 CONFIGURATION
+
+None.
+
+=head1 DEPENDENCIES
+
+=over
+
+=back
+
+=head1 BUGS AND LIMITATIONS
+
+=head1 AUTHOR
+
+Copyright (C) 2016 by Daniel Friesel E<lt>derf@finalrewind.orgE<gt>
+
+=head1 LICENSE
+
+  0. You just DO WHAT THE FUCK YOU WANT TO.
diff --git a/bin/gradient b/bin/gradient
new file mode 100755
index 0000000..d2f43ef
--- /dev/null
+++ b/bin/gradient
@@ -0,0 +1,32 @@
+#!/usr/bin/env python3
+
+import csv
+import numpy as np
+import os
+import struct
+import sys
+import tarfile
+import matplotlib.pyplot as plt
+from dfatool import running_mean, MIMOSA
+
+voltage = float(sys.argv[1])
+shunt = float(sys.argv[2])
+mimfile = sys.argv[3]
+
+mim = MIMOSA(voltage, shunt)
+
+charges, triggers = mim.load_data(mimfile)
+charges = charges[2000000:3000000]
+
+currents = running_mean(mim.charge_to_current_nocal(charges), 10) * 1e-6
+xr = np.arange(len(currents)) * 1e-5
+threshold = 3e-5
+grad = np.gradient(currents, 2)
+tp = np.abs(grad) > threshold
+plt.plot( xr, currents, "r-")
+plt.plot( xr, grad, "y-")
+plt.plot( xr[tp], grad[tp], "bo")
+plt.xlabel('Zeit [s]')
+plt.ylabel('Strom [A]')
+plt.grid(True)
+plt.show()
diff --git a/bin/merge.py b/bin/merge.py
new file mode 100755
index 0000000..1ff7e74
--- /dev/null
+++ b/bin/merge.py
@@ -0,0 +1,1033 @@
+#!/usr/bin/env python3
+
+import getopt
+import json
+import numpy as np
+import os
+import re
+import sys
+import plotter
+from copy import deepcopy
+from dfatool import aggregate_measures, regression_measures, is_numeric, powerset
+from matplotlib.patches import Polygon
+from scipy import optimize
+
+opts = {}
+
+def load_json(filename):
+    with open(filename, "r") as f:
+        return json.load(f)
+
+def save_json(data, filename):
+    with open(filename, "w") as f:
+        return json.dump(data, f)
+
+def print_data(aggregate):
+    for key in sorted(aggregate.keys()):
+        data = aggregate[key]
+        name, params = key
+        print("%s @ %s : ~ = %.f (%.f, %.f)  µ_σ_outer = %.f  n = %d" %
+            (name, params, np.median(data['means']), np.percentile(data['means'], 25),
+                np.percentile(data['means'], 75), np.mean(data['stds']), len(data['means'])))
+
+def flatten(somelist):
+    return [item for sublist in somelist for item in sublist]
+
+def mimosa_data(elem):
+    means = [x['uW_mean'] for x in elem['offline']]
+    durations = [x['us'] - 20 for x in elem['offline']]
+    stds = [x['uW_std'] for x in elem['offline']]
+    energies = [x['uW_mean'] * (x['us'] - 20) for x in elem['offline']]
+    clips = [x['clip_rate'] for x in elem['offline']]
+    substate_thresholds = []
+    substate_data = []
+    timeouts = []
+    rel_energies = []
+    if 'timeout' in elem['offline'][0]:
+        timeouts = [x['timeout'] for x in elem['offline']]
+    if 'uW_mean_delta' in elem['offline'][0]:
+        rel_energies = [x['uW_mean_delta'] * (x['us'] - 20) for x in elem['offline']]
+    for x in elem['offline']:
+        if 'substates' in x:
+            substate_thresholds.append(x['substates']['threshold'])
+            substate_data.append(x['substates']['states'])
+
+    return means, stds, durations, energies, rel_energies, clips, timeouts, substate_thresholds
+
+def online_data(elem):
+    means = [int(x['power']) for x in elem['online']]
+    durations = [int(x['time']) for x in elem['online']]
+
+    return means, durations
+
+# parameters = statistic variables such as txpower, bitrate etc.
+# variables = function variables/parameters set by linear regression
+def str_to_param_function(function_string, parameters, variables):
+    rawfunction = function_string
+    dependson = [False] * len(parameters)
+
+    for i in range(len(parameters)):
+        if rawfunction.find("global(%s)" % (parameters[i])) >= 0:
+            dependson[i] = True
+            rawfunction = rawfunction.replace("global(%s)" % (parameters[i]), "arg[%d]" % (i))
+    for i in range(len(variables)):
+        rawfunction = rawfunction.replace("param(%d)" % (i), "param[%d]" % (i))
+    fitfunc = eval("lambda param, arg: " + rawfunction);
+
+    return fitfunc, dependson
+
+def mk_function_data(name, paramdata, parameters, dependson, datatype):
+    X = [[] for i in range(len(parameters))]
+    Xm = [[] for i in range(len(parameters))]
+    Y = []
+    Ym = []
+
+    num_valid = 0
+    num_total = 0
+
+    for key, val in paramdata.items():
+        if key[0] == name and len(key[1]) == len(parameters):
+            valid = True
+            num_total += 1
+            for i in range(len(parameters)):
+                if dependson[i] and not is_numeric(key[1][i]):
+                    valid = False
+            if valid:
+                num_valid += 1
+                Y.extend(val[datatype])
+                Ym.append(np.median(val[datatype]))
+                for i in range(len(parameters)):
+                    if dependson[i] or is_numeric(key[1][i]):
+                        X[i].extend([int(key[1][i])] * len(val[datatype]))
+                        Xm[i].append(int(key[1][i]))
+                    else:
+                        X[i].extend([0] * len(val[datatype]))
+                        Xm[i].append(0)
+
+    for i in range(len(parameters)):
+        X[i] = np.array(X[i])
+        Xm[i] = np.array(Xm[i])
+    X = tuple(X)
+    Xm = tuple(Xm)
+    Y = np.array(Y)
+    Ym = np.array(Ym)
+
+    return X, Y, Xm, Ym, num_valid, num_total
+
+def raw_num0_8(num):
+    return (8 - num1(num))
+
+def raw_num0_16(num):
+    return (16 - num1(num))
+
+def raw_num1(num):
+    return bin(int(num)).count("1")
+
+num0_8 = np.vectorize(raw_num0_8)
+num0_16 = np.vectorize(raw_num0_16)
+num1 = np.vectorize(raw_num1)
+
+def try_fits(name, datatype, paramidx, paramdata):
+    functions = {
+        'linear' : lambda param, arg: param[0] + param[1] * arg,
+        'logarithmic' : lambda param, arg: param[0] + param[1] * np.log(arg),
+        'logarithmic1' : lambda param, arg: param[0] + param[1] * np.log(arg + 1),
+        'exponential' : lambda param, arg: param[0] + param[1] * np.exp(arg),
+        #'polynomial' : lambda param, arg: param[0] + param[1] * arg + param[2] * arg ** 2,
+        'square' : lambda param, arg: param[0] + param[1] * arg ** 2,
+        'fractional' : lambda param, arg: param[0] + param[1] / arg,
+        'sqrt' : lambda param, arg: param[0] + param[1] * np.sqrt(arg),
+        #'num0_8' : lambda param, arg: param[0] + param[1] * num0_8(arg),
+        #'num0_16' : lambda param, arg: param[0] + param[1] * num0_16(arg),
+        #'num1' : lambda param, arg: param[0] + param[1] * num1(arg),
+    }
+    results = dict([[key, []] for key in functions.keys()])
+    errors = {}
+
+    allvalues = [(*key[1][:paramidx], *key[1][paramidx+1:]) for key in paramdata.keys() if key[0] == name]
+    allvalues = list(set(allvalues))
+
+    for value in allvalues:
+        X = []
+        Xm = []
+        Y = []
+        Ym = []
+        num_valid = 0
+        num_total = 0
+        for key, val in paramdata.items():
+            if key[0] == name and len(key[1]) > paramidx and (*key[1][:paramidx], *key[1][paramidx+1:]) == value:
+                num_total += 1
+                if is_numeric(key[1][paramidx]):
+                    num_valid += 1
+                    Y.extend(val[datatype])
+                    Ym.append(np.median(val[datatype]))
+                    X.extend([int(key[1][paramidx])] * len(val[datatype]))
+                    Xm.append(int(key[1][paramidx]))
+                    if int(key[1][paramidx]) == 0:
+                        functions.pop('fractional', None)
+                    if int(key[1][paramidx]) <= 0:
+                        functions.pop('logarithmic', None)
+                    if int(key[1][paramidx]) < 0:
+                        functions.pop('logarithmic1', None)
+                        functions.pop('sqrt', None)
+                    if int(key[1][paramidx]) > 64:
+                        functions.pop('exponential', None)
+
+        # there should be -at least- two values when fitting
+        if num_valid > 1:
+            Y = np.array(Y)
+            Ym = np.array(Ym)
+            X = np.array(X)
+            Xm = np.array(Xm)
+            for kind, function in functions.items():
+                results[kind] = {}
+                errfunc = lambda P, X, y: function(P, X) - y
+                try:
+                    res = optimize.least_squares(errfunc, [0, 1], args=(X, Y), xtol=2e-15)
+                    measures = regression_measures(function(res.x, X), Y)
+                    for k, v in measures.items():
+                        if not k in results[kind]:
+                            results[kind][k] = []
+                        results[kind][k].append(v)
+                except:
+                    pass
+
+    for function_name, result in results.items():
+        if len(result) > 0 and function_name in functions:
+            errors[function_name] = {}
+            for measure in result.keys():
+                errors[function_name][measure] = np.mean(result[measure])
+
+    return errors
+
+def fit_function(function, name, datatype, parameters, paramdata, xaxis=None, yaxis=None):
+
+    variables = list(map(lambda x: float(x), function['params']))
+    fitfunc, dependson = str_to_param_function(function['raw'], parameters, variables)
+
+    X, Y, Xm, Ym, num_valid, num_total = mk_function_data(name, paramdata, parameters, dependson, datatype)
+
+    if num_valid > 0:
+
+        if num_valid != num_total:
+            num_invalid = num_total - num_valid
+            print("Warning: fit(%s): %d of %d states had incomplete parameter hashes" % (name, num_invalid, len(paramdata)))
+
+        errfunc = lambda P, X, y: fitfunc(P, X) - y
+        try:
+            res = optimize.least_squares(errfunc, variables, args=(X, Y), xtol=2e-15) # loss='cauchy'
+        except ValueError as err:
+            function['error'] = str(err)
+            return
+        #x1 = optimize.curve_fit(lambda param, *arg: fitfunc(param, arg), X, Y, functionparams)
+        measures = regression_measures(fitfunc(res.x, X), Y)
+
+        if res.status <= 0:
+            function['error'] = res.message
+            return
+
+        if 'fit' in opts:
+            for i in range(len(parameters)):
+                plotter.plot_param_fit(function['raw'], name, fitfunc, res.x, parameters, datatype, i, X, Y, xaxis, yaxis)
+
+        function['params'] = list(res.x)
+        function['fit'] = measures
+
+    else:
+        function['error'] = 'log contained no numeric parameters'
+
+def assess_function(function, name, datatype, parameters, paramdata):
+
+    variables = list(map(lambda x: float(x), function['params']))
+    fitfunc, dependson = str_to_param_function(function['raw'], parameters, variables)
+    X, Y, Xm, Ym, num_valid, num_total = mk_function_data(name, paramdata, parameters, dependson, datatype)
+
+    if num_valid > 0:
+        return regression_measures(fitfunc(variables, X), Y)
+    else:
+        return None
+
+def xv_assess_function(name, funbase, what, validation, mae, smape):
+    goodness = assess_function(funbase, name, what, parameters, validation)
+    if goodness != None:
+        if not name in mae:
+            mae[name] = []
+        if not name in smape:
+            smape[name] = []
+        mae[name].append(goodness['mae'])
+        smape[name].append(goodness['smape'])
+
+def xv2_assess_function(name, funbase, what, validation, mae, smape, rmsd):
+    goodness = assess_function(funbase, name, what, parameters, validation)
+    if goodness != None:
+        if goodness['mae'] < 10e9:
+            mae.append(goodness['mae'])
+            rmsd.append(goodness['rmsd'])
+            smape.append(goodness['smape'])
+        else:
+            print('[!] Ignoring MAE of %d (SMAPE %.f)' % (goodness['mae'], goodness['smape']))
+
+# Returns the values used for each parameter in the measurement, e.g.
+# { 'txpower' : [1, 2, 4, 8], 'length' : [16] }
+# non-numeric values such as '' are skipped
+def param_values(parameters, by_param):
+    paramvalues = dict([[param, set()] for param in parameters])
+
+    for _, paramvalue in by_param.keys():
+        for i, param in enumerate(parameters):
+            if is_numeric(paramvalue[i]):
+                paramvalues[param].add(paramvalue[i])
+
+    return paramvalues
+
+def param_key(elem):
+    name = elem['name']
+    paramtuple = ()
+
+    if 'parameter' in elem:
+        paramkeys = sorted(elem['parameter'].keys())
+        paramtuple = tuple([elem['parameter'][x] for x in paramkeys])
+
+    return (name, paramtuple)
+
+#def param_arg_key(elem):
+#    # Argumentbasierte Parametrisierung ist erstmal out of scope
+#    #if 'args' in elem:
+#    #    argtuple = tuple(elem['args'])
+#
+#    return (name, paramtuple, argtuple)
+
+def add_data_to_aggregate(aggregate, key, isa, data):
+    if not key in aggregate:
+        aggregate[key] = {
+            'isa' : isa,
+        }
+        for datakey in data.keys():
+            aggregate[key][datakey] = []
+    for datakey, dataval in data.items():
+        aggregate[key][datakey].extend(dataval)
+
+def fake_add_data_to_aggregate(aggregate, key, isa, database, idx):
+    timeout_val = []
+    if len(database['timeouts']):
+        timeout_val = [database['timeouts'][idx]]
+    rel_energy_val = []
+    if len(database['rel_energies']):
+        rel_energy_val = [database['rel_energies'][idx]]
+    add_data_to_aggregate(aggregate, key, isa, {
+        'means' : [database['means'][idx]],
+        'stds' : [database['stds'][idx]],
+        'durations' : [database['durations'][idx]],
+        'energies' : [database['energies'][idx]],
+        'rel_energies' : rel_energy_val,
+        'clip_rate' : [database['clip_rate'][idx]],
+        'timeouts' : timeout_val,
+    })
+
+def weight_by_name(aggdata):
+    total = {}
+    count = {}
+    weight = {}
+    for key in aggdata.keys():
+        if not key[0] in total:
+            total[key[0]] = 0
+        total[key[0]] += len(aggdata[key]['means'])
+        count[key] = len(aggdata[key]['means'])
+    for key in aggdata.keys():
+        weight[key] = float(count[key]) / total[key[0]]
+    return weight
+
+# returns the mean standard deviation of all measurements of 'what'
+# (e.g. power consumption or timeout) for state/transition 'name' where
+# parameter 'index' is dynamic and all other parameters are fixed.
+# I.e., if parameters are a, b, c ∈ {1,2,3} and 'index' corresponds to b', then
+# this function returns the mean of the standard deviations of (a=1, b=*, c=1),
+# (a=1, b=*, c=2), and so on
+def mean_std_by_param(data, keys, name, what, index):
+    partitions = []
+    for key in keys:
+        partition = []
+        for k, v in data.items():
+            if (*k[1][:index], *k[1][index+1:]) == key and k[0] == name:
+                partition.extend(v[what])
+        partitions.append(partition)
+    return np.mean([np.std(partition) for partition in partitions])
+
+# returns the mean standard deviation of all measurements of 'what'
+# (e.g. power consumption or timeout) for state/transition 'name' where the
+# trace of previous transitions is fixed except for a single transition,
+# whose occurence or absence is silently ignored.
+# this is done separately for each transition (-> returns a dictionary)
+def mean_std_by_trace_part(data, transitions, name, what):
+    ret = {}
+    for transition in transitions:
+        keys = set(map(lambda x : (x[0], x[1], tuple([y for y in x[2] if y != transition])), data.keys()))
+        ret[transition] = {}
+        partitions = []
+        for key in keys:
+            partition = []
+            for k, v in data.items():
+                key_without_transition = (k[0], k[1], tuple([y for y in k[2] if y != transition]))
+                if key[0] == name and key == key_without_transition:
+                    partition.extend(v[what])
+            if len(partition):
+                partitions.append(partition)
+        ret[transition] = np.mean([np.std(partition) for partition in partitions])
+    return ret
+
+
+def load_run_elem(index, element, trace, by_name, by_param, by_trace):
+    means, stds, durations, energies, rel_energies, clips, timeouts, sub_thresholds = mimosa_data(element)
+
+    online_means = []
+    online_durations = []
+    if element['isa'] == 'state':
+        online_means, online_durations = online_data(element)
+
+    key = param_key(element)
+    pre_trace = tuple(map(lambda x : x['name'], trace[1:index:2]))
+    trace_key = (*key, pre_trace)
+    name = element['name']
+
+    add_data_to_aggregate(by_name, name, element['isa'], {
+        'means' : means,
+        'stds' : stds,
+        'durations' : durations,
+        'energies' : energies,
+        'rel_energies' : rel_energies,
+        'clip_rate' : clips,
+        'timeouts' : timeouts,
+        'sub_thresholds' : sub_thresholds,
+        'param' : [key[1]] * len(means),
+        'online_means' : online_means,
+        'online_durations' : online_durations,
+    })
+    add_data_to_aggregate(by_param, key, element['isa'], {
+        'means' : means,
+        'stds' : stds,
+        'durations' : durations,
+        'energies' : energies,
+        'rel_energies' : rel_energies,
+        'clip_rate' : clips,
+        'timeouts' : timeouts,
+        'sub_thresholds' : sub_thresholds,
+        'online_means' : online_means,
+        'online_durations' : online_durations,
+    })
+    add_data_to_aggregate(by_trace, trace_key, element['isa'], {
+        'means' : means,
+        'stds' : stds,
+        'durations' : durations,
+        'energies' : energies,
+        'rel_energies' : rel_energies,
+        'clip_rate' : clips,
+        'timeouts' : timeouts,
+        'sub_thresholds' : sub_thresholds,
+        'online_means' : online_means,
+        'online_durations' : online_durations,
+    })
+
+def fmap(name, funtype):
+    if funtype == 'linear':
+        return "global(%s)" % name
+    if funtype == 'logarithmic':
+        return "np.log(global(%s))" % name
+    if funtype == 'logarithmic1':
+        return "np.log(global(%s) + 1)" % name
+    if funtype == 'exponential':
+        return "np.exp(global(%s))" % name
+    if funtype == 'square':
+        return "global(%s)**2" % name
+    if funtype == 'fractional':
+        return "1 / global(%s)" % name
+    if funtype == 'sqrt':
+        return "np.sqrt(global(%s))" % name
+    if funtype == 'num0_8':
+        return "num0_8(global(%s))" % name
+    if funtype == 'num0_16':
+        return "num0_16(global(%s))" % name
+    if funtype == 'num1':
+        return "num1(global(%s))" % name
+    return "ERROR"
+
+def fguess_to_function(name, datatype, aggdata, parameters, paramdata, yaxis):
+    best_fit = {}
+    fitguess = aggdata['fit_guess']
+    params = list(filter(lambda x : x in fitguess, parameters))
+    if len(params) > 0:
+        for param in params:
+            best_fit_val = np.inf
+            for func_name, fit_val in fitguess[param].items():
+                if fit_val['rmsd'] < best_fit_val:
+                    best_fit_val = fit_val['rmsd']
+                    best_fit[param] = func_name
+        buf = '0'
+        pidx = 0
+        for elem in powerset(best_fit.items()):
+            buf += " + param(%d)" % pidx
+            pidx += 1
+            for fun in elem:
+                buf += " * %s" % fmap(*fun)
+        aggdata['function']['estimate'] = {
+            'raw' : buf,
+            'params' : list(np.ones((pidx))),
+            'base' : [best_fit[param] for param in params]
+        }
+        fit_function(
+            aggdata['function']['estimate'], name, datatype, parameters,
+            paramdata, yaxis=yaxis)
+
+def param_measures(name, paramdata, key, fun):
+    mae = []
+    smape = []
+    rmsd = []
+    for pkey, pval in paramdata.items():
+        if pkey[0] == name:
+            # Median ist besseres Maß für MAE / SMAPE,
+            # Mean ist besseres für SSR. Da least_squares SSR optimiert
+            # nutzen wir hier auch Mean.
+            goodness = aggregate_measures(fun(pval[key]), pval[key])
+            mae.append(goodness['mae'])
+            rmsd.append(goodness['rmsd'])
+            if 'smape' in goodness:
+                smape.append(goodness['smape'])
+    ret = {
+        'mae' : np.mean(mae),
+        'rmsd' : np.mean(rmsd)
+    }
+    if len(smape):
+        ret['smape'] = np.mean(smape)
+
+    return ret
+
+def keydata(name, val, paramdata, tracedata, key):
+    ret = {
+        'count' : len(val[key]),
+        'median' : np.median(val[key]),
+        'mean'   : np.mean(val[key]),
+        'mean_goodness' : aggregate_measures(np.mean(val[key]), val[key]),
+        'median_goodness' : aggregate_measures(np.median(val[key]), val[key]),
+        'param_mean_goodness' : param_measures(name, paramdata, key, np.mean),
+        'param_median_goodness' : param_measures(name, paramdata, key, np.median),
+        'std_inner' : np.std(val[key]),
+        'std_param' : np.mean([np.std(paramdata[x][key]) for x in paramdata.keys() if x[0] == name]),
+        'std_trace' : np.mean([np.std(tracedata[x][key]) for x in tracedata.keys() if x[0] == name]),
+        'std_by_param' : {},
+        'fit_guess' : {},
+        'function' : {},
+    }
+
+    return ret
+
+def splitidx_kfold(length, num_slices):
+    pairs = []
+    indexes = np.arange(length)
+    for i in range(0, num_slices):
+        training = np.delete(indexes, slice(i, None, num_slices))
+        validation = indexes[i::num_slices]
+        pairs.append((training, validation))
+    return pairs
+
+def splitidx_srs(length, num_slices):
+    pairs = []
+    for i in range(0, num_slices):
+        shuffled = np.random.permutation(np.arange(length))
+        border = int(length * float(2) / 3)
+        training = shuffled[:border]
+        validation = shuffled[border:]
+        pairs.append((training, validation))
+    return pairs
+
+def val_run(aggdata, split_fun, count):
+    mae = []
+    smape = []
+    rmsd = []
+    pairs = split_fun(len(aggdata), count)
+    for i in range(0, count):
+        training = aggdata[pairs[i][0]]
+        validation = aggdata[pairs[i][1]]
+        median = np.median(training)
+        goodness = aggregate_measures(median, validation)
+        mae.append(goodness['mae'])
+        rmsd.append(goodness['rmsd'])
+        if 'smape' in goodness:
+            smape.append(goodness['smape'])
+
+    mae_mean = np.mean(mae)
+    rmsd_mean = np.mean(rmsd)
+    if len(smape):
+        smape_mean = np.mean(smape)
+    else:
+        smape_mean = -1
+
+    return mae_mean, smape_mean, rmsd_mean
+
+# by_trace is not part of the cross-validation process
+def val_run_fun(aggdata, by_trace, name, key, funtype1, funtype2, splitfun, count):
+    aggdata = aggdata[name]
+    isa = aggdata['isa']
+    mae = []
+    smape = []
+    rmsd = []
+    estimates = []
+    pairs = splitfun(len(aggdata[key]), count)
+    for i in range(0, count):
+        bpa_training = {}
+        bpa_validation = {}
+
+        for idx in pairs[i][0]:
+            bpa_key = (name, aggdata['param'][idx])
+            fake_add_data_to_aggregate(bpa_training, bpa_key, isa, aggdata, idx)
+        for idx in pairs[i][1]:
+            bpa_key = (name, aggdata['param'][idx])
+            fake_add_data_to_aggregate(bpa_validation, bpa_key, isa, aggdata, idx)
+
+        fake_by_name = { name : aggdata }
+        ares = analyze(fake_by_name, bpa_training, by_trace, parameters)
+        if name in ares[isa] and funtype2 in ares[isa][name][funtype1]['function']:
+            xv2_assess_function(name, ares[isa][name][funtype1]['function'][funtype2], key, bpa_validation, mae, smape, rmsd)
+            if funtype2 == 'estimate':
+                if 'base' in ares[isa][name][funtype1]['function'][funtype2]:
+                    estimates.append(tuple(ares[isa][name][funtype1]['function'][funtype2]['base']))
+                else:
+                    estimates.append(None)
+    return mae, smape, rmsd, estimates
+
+# by_trace is not part of the cross-validation process
+def val_run_fun_p(aggdata, by_trace, name, key, funtype1, funtype2, splitfun, count):
+    aggdata = dict([[x, aggdata[x]] for x in aggdata if x[0] == name])
+    isa = aggdata[list(aggdata.keys())[0]]['isa']
+    mae = []
+    smape = []
+    rmsd = []
+    estimates = []
+    pairs = splitfun(len(aggdata.keys()), count) # pairs are by_param index arrays
+    keys = sorted(aggdata.keys())
+    for i in range(0, count):
+        bpa_training = dict([[keys[x], aggdata[keys[x]]] for x in pairs[i][0]])
+        bpa_validation = dict([[keys[x], aggdata[keys[x]]] for x in pairs[i][1]])
+        bna_training = {}
+        for val in bpa_training.values():
+            for idx in range(0, len(val[key])):
+                fake_add_data_to_aggregate(bna_training, name, isa, val, idx)
+
+        ares = analyze(bna_training, bpa_training, by_trace, parameters)
+        if name in ares[isa] and funtype2 in ares[isa][name][funtype1]['function']:
+            xv2_assess_function(name, ares[isa][name][funtype1]['function'][funtype2], key, bpa_validation, mae, smape, rmsd)
+            if funtype2 == 'estimate':
+                if 'base' in ares[isa][name][funtype1]['function'][funtype2]:
+                    estimates.append(tuple(ares[isa][name][funtype1]['function'][funtype2]['base']))
+                else:
+                    estimates.append(None)
+    return mae, smape, rmsd, estimates
+
+def crossvalidate(by_name, by_param, by_trace, model, parameters):
+    param_mc_count = 200
+    paramv = param_values(parameters, by_param)
+    for name in sorted(by_name.keys()):
+        isa = by_name[name]['isa']
+        by_name[name]['means'] = np.array(by_name[name]['means'])
+        by_name[name]['energies'] = np.array(by_name[name]['energies'])
+        by_name[name]['rel_energies'] = np.array(by_name[name]['rel_energies'])
+        by_name[name]['durations'] = np.array(by_name[name]['durations'])
+
+        if isa == 'state':
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['means'], splitidx_srs, 200)
+            print('%16s,   static      power,             Monte Carlo: MAE %8.f µW,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['means'], splitidx_kfold, 10)
+            print('%16s,   static      power,             10-fold sys: MAE %8.f µW,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+        else:
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['energies'], splitidx_srs, 200)
+            print('%16s,   static     energy,             Monte Carlo: MAE %8.f pJ,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['energies'], splitidx_kfold, 10)
+            print('%16s,   static     energy,             10-fold sys: MAE %8.f pJ,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['rel_energies'], splitidx_srs, 200)
+            print('%16s,   static rel_energy,             Monte Carlo: MAE %8.f pJ,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['rel_energies'], splitidx_kfold, 10)
+            print('%16s,   static rel_energy,             10-fold sys: MAE %8.f pJ,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['durations'], splitidx_srs, 200)
+            print('%16s,   static   duration,             Monte Carlo: MAE %8.f µs,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+            mae_mean, smape_mean, rms_mean = val_run(by_name[name]['durations'], splitidx_kfold, 10)
+            print('%16s,   static   duration,             10-fold sys: MAE %8.f µs,  SMAPE %6.2f%%,  RMS %d' % (name, mae_mean, smape_mean, rms_mean))
+
+        def print_estimates(estimates, total):
+            histogram = {}
+            buf = '    '
+            for estimate in estimates:
+                if not estimate in histogram:
+                    histogram[estimate] = 1
+                else:
+                    histogram[estimate] += 1
+            for estimate, count in sorted(histogram.items(), key=lambda kv: kv[1], reverse=True):
+                buf += '  %.f%% %s' % (count * 100 / total, estimate)
+            if len(estimates):
+                print(buf)
+
+        def val_run_funs(by_name, by_trace, name, key1, key2, key3, unit):
+            mae, smape, rmsd, estimates = val_run_fun(by_name, by_trace, name, key1, key2, key3, splitidx_srs, param_mc_count)
+            print('%16s, %8s %10s,             Monte Carlo: MAE %8.f %s,  SMAPE %6.2f%%,  RMS %d' % (
+                name, key3, key2, np.mean(mae), unit, np.mean(smape), np.mean(rmsd)))
+            print_estimates(estimates, param_mc_count)
+            mae, smape, rmsd, estimates = val_run_fun(by_name, by_trace, name, key1, key2, key3, splitidx_kfold, 10)
+            print('%16s, %8s %10s,             10-fold sys: MAE %8.f %s,  SMAPE %6.2f%%,  RMS %d' % (
+                name, key3, key2, np.mean(mae), unit, np.mean(smape), np.mean(rmsd)))
+            print_estimates(estimates, 10)
+            mae, smape, rmsd, estimates = val_run_fun_p(by_param, by_trace, name, key1, key2, key3, splitidx_srs, param_mc_count)
+            print('%16s, %8s %10s, param-aware Monte Carlo: MAE %8.f %s,  SMAPE %6.2f%%,  RMS %d' % (
+                name, key3, key2, np.mean(mae), unit, np.mean(smape), np.mean(rmsd)))
+            print_estimates(estimates, param_mc_count)
+            mae, smape, rmsd, estimates = val_run_fun_p(by_param, by_trace, name, key1, key2, key3, splitidx_kfold, 10)
+            print('%16s, %8s %10s, param-aware 10-fold sys: MAE %8.f %s,  SMAPE %6.2f%%,  RMS %d' % (
+                name, key3, key2, np.mean(mae), unit, np.mean(smape), np.mean(rmsd)))
+            print_estimates(estimates, 10)
+
+        if 'power' in model[isa][name] and 'function' in model[isa][name]['power']:
+            if 'user' in model[isa][name]['power']['function']:
+                val_run_funs(by_name, by_trace, name, 'means', 'power', 'user', 'µW')
+            if 'estimate' in model[isa][name]['power']['function']:
+                val_run_funs(by_name, by_trace, name, 'means', 'power', 'estimate', 'µW')
+        if 'timeout' in model[isa][name] and 'function' in model[isa][name]['timeout']:
+            if 'user' in model[isa][name]['timeout']['function']:
+                val_run_funs(by_name, by_trace, name, 'timeouts', 'timeout', 'user', 'µs')
+            if 'estimate' in model[isa][name]['timeout']['function']:
+                val_run_funs(by_name, by_trace, name, 'timeouts', 'timeout', 'estimate', 'µs')
+        if 'duration' in model[isa][name] and 'function' in model[isa][name]['duration']:
+            if 'user' in model[isa][name]['duration']['function']:
+                val_run_funs(by_name, by_trace, name, 'durations', 'duration', 'user', 'µs')
+            if 'estimate' in model[isa][name]['duration']['function']:
+                val_run_funs(by_name, by_trace, name, 'durations', 'duration', 'estimate', 'µs')
+        if 'energy' in model[isa][name] and 'function' in model[isa][name]['energy']:
+            if 'user' in model[isa][name]['energy']['function']:
+                val_run_funs(by_name, by_trace, name, 'energies', 'energy', 'user', 'pJ')
+            if 'estimate' in model[isa][name]['energy']['function']:
+                val_run_funs(by_name, by_trace, name, 'energies', 'energy', 'estimate', 'pJ')
+        if 'rel_energy' in model[isa][name] and 'function' in model[isa][name]['rel_energy']:
+            if 'user' in model[isa][name]['rel_energy']['function']:
+                val_run_funs(by_name, by_trace, name, 'rel_energies', 'rel_energy', 'user', 'pJ')
+            if 'estimate' in model[isa][name]['rel_energy']['function']:
+                val_run_funs(by_name, by_trace, name, 'rel_energies', 'rel_energy', 'estimate', 'pJ')
+
+    return
+    for i, param in enumerate(parameters):
+        user_mae = {}
+        user_smape = {}
+        estimate_mae = {}
+        estimate_smape = {}
+        for val in paramv[param]:
+            bpa_training = dict([[x, by_param[x]] for x in by_param if x[1][i] != val])
+            bpa_validation = dict([[x, by_param[x]] for x in by_param if x[1][i] == val])
+            to_pop = []
+            for name in by_name.keys():
+                if not any(map(lambda x : x[0] == name, bpa_training.keys())):
+                    to_pop.append(name)
+            for name in to_pop:
+                by_name.pop(name, None)
+            ares = analyze(by_name, bpa_training, by_trace, parameters)
+            for name in sorted(ares['state'].keys()):
+                state = ares['state'][name]
+                if 'function' in state['power']:
+                    if 'user' in state['power']['function']:
+                        xv_assess_function(name, state['power']['function']['user'], 'means', bpa_validation, user_mae, user_smape)
+                    if 'estimate' in state['power']['function']:
+                        xv_assess_function(name, state['power']['function']['estimate'], 'means', bpa_validation, estimate_mae, estimate_smape)
+            for name in sorted(ares['transition'].keys()):
+                trans = ares['transition'][name]
+                if 'timeout' in trans and 'function' in trans['timeout']:
+                    if 'user' in trans['timeout']['function']:
+                        xv_assess_function(name, trans['timeout']['function']['user'], 'timeouts', bpa_validation, user_mae, user_smape)
+                    if 'estimate' in trans['timeout']['function']:
+                        xv_assess_function(name, trans['timeout']['function']['estimate'], 'timeouts', bpa_validation, estimate_mae, estimate_smape)
+
+        for name in sorted(user_mae.keys()):
+            if by_name[name]['isa'] == 'state':
+                print('user function %s power by %s: MAE %.f µW,  SMAPE %.2f%%' % (
+                    name, param, np.mean(user_mae[name]), np.mean(user_smape[name])))
+            else:
+                print('user function %s timeout by %s: MAE %.f µs,  SMAPE %.2f%%' % (
+                    name, param, np.mean(user_mae[name]), np.mean(user_smape[name])))
+        for name in sorted(estimate_mae.keys()):
+            if by_name[name]['isa'] == 'state':
+                print('estimate function %s power by %s: MAE %.f µW,  SMAPE %.2f%%' % (
+                    name, param, np.mean(estimate_mae[name]), np.mean(estimate_smape[name])))
+            else:
+                print('estimate function %s timeout by %s: MAE %.f µs,  SMAPE %.2f%%' % (
+                    name, param, np.mean(estimate_mae[name]), np.mean(estimate_smape[name])))
+
+def validate(by_name, by_param, parameters):
+    aggdata = {
+        'state' : {},
+        'transition' : {},
+    }
+    for key, val in by_name.items():
+        name = key
+        isa = val['isa']
+        model = data['model'][isa][name]
+
+        if isa == 'state':
+            aggdata[isa][name] = {
+                'power' : {
+                    'goodness' : aggregate_measures(model['power']['static'], val['means']),
+                    'median' : np.median(val['means']),
+                    'mean'   : np.mean(val['means']),
+                    'std_inner' : np.std(val['means']),
+                    'function' : {},
+                },
+                'online_power' : {
+                    'goodness' : regression_measures(np.array(val['online_means']), np.array(val['means'])),
+                    'median' : np.median(val['online_means']),
+                    'mean'   : np.mean(val['online_means']),
+                    'std_inner' : np.std(val['online_means']),
+                    'function' : {},
+                },
+                'online_duration' : {
+                    'goodness' : regression_measures(np.array(val['online_durations']), np.array(val['durations'])),
+                    'median' : np.median(val['online_durations']),
+                    'mean'   : np.mean(val['online_durations']),
+                    'std_inner' : np.std(val['online_durations']),
+                    'function' : {},
+                },
+                'clip' : {
+                    'mean' : np.mean(val['clip_rate']),
+                    'max'  : max(val['clip_rate']),
+                },
+                'timeout' : {},
+            }
+            if 'function' in model['power']:
+                aggdata[isa][name]['power']['function'] = {
+                    'estimate' : {
+                        'fit' : assess_function(model['power']['function']['estimate'],
+                            name, 'means', parameters, by_param),
+                    },
+                    'user': {
+                        'fit' : assess_function(model['power']['function']['user'],
+                            name, 'means', parameters, by_param),
+                    },
+                }
+        if isa == 'transition':
+            aggdata[isa][name] = {
+                'duration' : {
+                    'goodness' : aggregate_measures(model['duration']['static'], val['durations']),
+                    'median' : np.median(val['durations']),
+                    'mean'   : np.mean(val['durations']),
+                    'std_inner' : np.std(val['durations']),
+                    'function' : {},
+                },
+                'energy' : {
+                    'goodness' : aggregate_measures(model['energy']['static'], val['energies']),
+                    'median' : np.median(val['energies']),
+                    'mean'   : np.mean(val['energies']),
+                    'std_inner' : np.std(val['energies']),
+                    'function' : {},
+                },
+                'rel_energy' : {
+                    'goodness' : aggregate_measures(model['rel_energy']['static'], val['rel_energies']),
+                    'median' : np.median(val['rel_energies']),
+                    'mean'   : np.mean(val['rel_energies']),
+                    'std_inner' : np.std(val['rel_energies']),
+                    'function' : {},
+                },
+                'clip' : {
+                    'mean' : np.mean(val['clip_rate']),
+                    'max'  : max(val['clip_rate']),
+                },
+                'timeout' : {},
+            }
+            if 'function' in model['timeout']:
+                aggdata[isa][name]['timeout'] = {
+                    'median' : np.median(val['timeouts']),
+                    'mean' : np.mean(val['timeouts']),
+                    'function': {
+                        'estimate' : {
+                            'fit' : assess_function(model['timeout']['function']['estimate'],
+                                name, 'timeouts', parameters, by_param),
+                        },
+                        'user': {
+                            'fit' : assess_function(model['timeout']['function']['user'],
+                                name, 'timeouts', parameters, by_param),
+                        },
+                    },
+                }
+    return aggdata
+
+def analyze(by_name, by_param, by_trace, parameters):
+    aggdata = {
+        'state' : {},
+        'transition' : {},
+    }
+    transition_names = list(map(lambda x: x[0], filter(lambda x: x[1]['isa'] == 'transition', by_name.items())))
+    for name, val in by_name.items():
+        isa = val['isa']
+        model = data['model'][isa][name]
+
+        aggdata[isa][name] = {
+            'power' : keydata(name, val, by_param, by_trace, 'means'),
+            'duration' : keydata(name, val, by_param, by_trace, 'durations'),
+            'energy' : keydata(name, val, by_param, by_trace, 'energies'),
+            'clip' : {
+                'mean' : np.mean(val['clip_rate']),
+                'max'  : max(val['clip_rate']),
+            },
+            'timeout' : {},
+        }
+
+        aggval = aggdata[isa][name]
+        aggval['power']['std_outer'] = np.mean(val['stds'])
+
+        if isa == 'transition':
+            aggval['rel_energy'] = keydata(name, val, by_param, by_trace, 'rel_energies')
+
+        if isa == 'transition' and 'function' in data['model']['transition'][name]['timeout']:
+            aggval['timeout'] = keydata(name, val, by_param, by_trace, 'timeouts')
+
+        for i, param in enumerate(sorted(data['model']['parameter'].keys())):
+            values = list(set([key[1][i] for key in by_param.keys() if key[0] == name and key[1][i] != '']))
+            allvalues = [(*key[1][:i], *key[1][i+1:]) for key in by_param.keys() if key[0] == name]
+            #allvalues = list(set(allvalues))
+            if len(values) > 1:
+                if isa == 'state':
+                    aggval['power']['std_by_param'][param] = mean_std_by_param(
+                        by_param, allvalues, name, 'means', i)
+                    if aggval['power']['std_by_param'][param] > 0 and aggval['power']['std_param'] / aggval['power']['std_by_param'][param] < 0.6:
+                        aggval['power']['fit_guess'][param] = try_fits(name, 'means', i, by_param)
+                else:
+                    aggval['duration']['std_by_param'][param] = mean_std_by_param(
+                        by_param, allvalues, name, 'durations', i)
+                    if aggval['duration']['std_by_param'][param] > 0 and aggval['duration']['std_param'] / aggval['duration']['std_by_param'][param] < 0.6:
+                        aggval['duration']['fit_guess'][param] = try_fits(name, 'durations', i, by_param)
+                    aggval['energy']['std_by_param'][param] = mean_std_by_param(
+                        by_param, allvalues, name, 'energies', i)
+                    if aggval['energy']['std_by_param'][param] > 0 and aggval['energy']['std_param'] / aggval['energy']['std_by_param'][param] < 0.6:
+                        aggval['energy']['fit_guess'][param] = try_fits(name, 'energies', i, by_param)
+                    aggval['rel_energy']['std_by_param'][param] = mean_std_by_param(
+                        by_param, allvalues, name, 'rel_energies', i)
+                    if aggval['rel_energy']['std_by_param'][param] > 0 and aggval['rel_energy']['std_param'] / aggval['rel_energy']['std_by_param'][param] < 0.6:
+                        aggval['rel_energy']['fit_guess'][param] = try_fits(name, 'rel_energies', i, by_param)
+                if isa == 'transition' and 'function' in data['model']['transition'][name]['timeout']:
+                    aggval['timeout']['std_by_param'][param] = mean_std_by_param(
+                        by_param, allvalues, name, 'timeouts', i)
+                    if aggval['timeout']['std_by_param'][param] > 0 and aggval['timeout']['std_param'] / aggval['timeout']['std_by_param'][param] < 0.6:
+                        aggval['timeout']['fit_guess'][param] = try_fits(name, 'timeouts', i, by_param)
+
+        if isa == 'state':
+            fguess_to_function(name, 'means', aggval['power'], parameters, by_param,
+                'estimated %s power [µW]' % name)
+            if 'function' in model['power'] and 'user' in model['power']['function']:
+                aggval['power']['function']['user'] = {
+                    'raw' : model['power']['function']['user']['raw'],
+                    'params' : model['power']['function']['user']['params'],
+                }
+                fit_function(
+                    aggval['power']['function']['user'], name, 'means', parameters, by_param,
+                    yaxis='%s power [µW]' % name)
+            if aggval['power']['std_param'] > 0 and aggval['power']['std_trace'] / aggval['power']['std_param'] < 0.5:
+                aggval['power']['std_by_trace'] = mean_std_by_trace_part(by_trace, transition_names, name, 'means')
+        else:
+            fguess_to_function(name, 'durations', aggval['duration'], parameters, by_param,
+                'estimated %s duration [µs]' % name)
+            fguess_to_function(name, 'energies', aggval['energy'], parameters, by_param,
+                'estimated %s energy [pJ]' % name)
+            fguess_to_function(name, 'rel_energies', aggval['rel_energy'], parameters, by_param,
+                'estimated relative %s energy [pJ]' % name)
+            if 'function' in model['duration'] and 'user' in model['duration']['function']:
+                aggval['duration']['function']['user'] = {
+                    'raw' : model['duration']['function']['user']['raw'],
+                    'params' : model['duration']['function']['user']['params'],
+                }
+                fit_function(
+                    aggval['duration']['function']['user'], name, 'durations', parameters, by_param,
+                    yaxis='%s duration [µs]' % name)
+            if 'function' in model['energy'] and 'user' in model['energy']['function']:
+                aggval['energy']['function']['user'] = {
+                    'raw' : model['energy']['function']['user']['raw'],
+                    'params' : model['energy']['function']['user']['params'],
+                }
+                fit_function(
+                    aggval['energy']['function']['user'], name, 'energies', parameters, by_param,
+                    yaxis='%s energy [pJ]' % name)
+            if 'function' in model['rel_energy'] and 'user' in model['rel_energy']['function']:
+                aggval['rel_energy']['function']['user'] = {
+                    'raw' : model['rel_energy']['function']['user']['raw'],
+                    'params' : model['rel_energy']['function']['user']['params'],
+                }
+                fit_function(
+                    aggval['rel_energy']['function']['user'], name, 'rel_energies', parameters, by_param,
+                    yaxis='%s rel_energy [pJ]' % name)
+            if 'function' in model['timeout'] and 'user' in model['timeout']['function']:
+                fguess_to_function(name, 'timeouts', aggval['timeout'], parameters, by_param,
+                    'estimated %s timeout [µs]' % name)
+                aggval['timeout']['function']['user'] = {
+                    'raw' : model['timeout']['function']['user']['raw'],
+                    'params' : model['timeout']['function']['user']['params'],
+                }
+                fit_function(
+                    aggval['timeout']['function']['user'], name, 'timeouts', parameters,
+                    by_param, yaxis='%s timeout [µs]' % name)
+                if aggval['timeout']['std_param'] > 0 and aggval['timeout']['std_trace'] / aggval['timeout']['std_param'] < 0.5:
+                    aggval['timeout']['std_by_trace'] = mean_std_by_trace_part(by_trace, transition_names, name, 'timeouts')
+
+    return aggdata
+
+try:
+    raw_opts, args = getopt.getopt(sys.argv[1:], "", [
+        "fit", "states", "transitions", "params", "clipping", "timing",
+        "histogram", "substates", "validate", "crossvalidate", "ignore-trace-idx="])
+    for option, parameter in raw_opts:
+        optname = re.sub(r'^--', '', option)
+        opts[optname] = parameter
+    if 'ignore-trace-idx' in opts:
+        opts['ignore-trace-idx'] = int(opts['ignore-trace-idx'])
+except getopt.GetoptError as err:
+    print(err)
+    sys.exit(2)
+
+data = load_json(args[0])
+by_name = {}
+by_param = {}
+by_trace = {}
+parameters = sorted(data['model']['parameter'].keys())
+
+for arg in args:
+    mdata = load_json(arg)
+    for runidx, run in enumerate(mdata['traces']):
+        if 'ignore-trace-idx' not in opts or opts['ignore-trace-idx'] != runidx:
+            for i, elem in enumerate(run['trace']):
+                if elem['name'] != 'UNINITIALIZED':
+                    load_run_elem(i, elem, run['trace'], by_name, by_param, by_trace)
+
+if 'states' in opts:
+    if 'params' in opts:
+        plotter.plot_states_param(data['model'], by_param)
+    else:
+        plotter.plot_states(data['model'], by_name)
+    if 'timing' in opts:
+        plotter.plot_states_duration(data['model'], by_name)
+        plotter.plot_states_duration(data['model'], by_param)
+    if 'clipping' in opts:
+        plotter.plot_states_clips(data['model'], by_name)
+if 'transitions' in opts:
+    plotter.plot_transitions(data['model'], by_name)
+    if 'timing' in opts:
+        plotter.plot_transitions_duration(data['model'], by_name)
+        plotter.plot_transitions_timeout(data['model'], by_param)
+    if 'clipping' in opts:
+        plotter.plot_transitions_clips(data['model'], by_name)
+if 'histogram' in opts:
+    for key in sorted(by_name.keys()):
+        plotter.plot_histogram(by_name[key]['means'])
+if 'substates' in opts:
+    if 'params' in opts:
+        plotter.plot_substate_thresholds_p(data['model'], by_param)
+    else:
+        plotter.plot_substate_thresholds(data['model'], by_name)
+
+if 'validate' in opts:
+    data['aggregate'] = validate(by_name, by_param, parameters)
+elif 'crossvalidate' in opts:
+    crossvalidate(by_name, by_param, by_trace, data['model'], parameters)
+else:
+    data['aggregate'] = analyze(by_name, by_param, by_trace, parameters)
+
+# TODO optionally also plot data points for states/transitions which do not have
+# a function, but may depend on a parameter (visualization is always good!)
+
+save_json(data, args[0])
diff --git a/bin/mim-vs-keysight.py b/bin/mim-vs-keysight.py
new file mode 100755
index 0000000..2cfa976
--- /dev/null
+++ b/bin/mim-vs-keysight.py
@@ -0,0 +1,225 @@
+#!/usr/bin/env python3
+
+import csv
+import numpy as np
+import os
+import struct
+import sys
+import tarfile
+import matplotlib.pyplot as plt
+from dfatool import running_mean, MIMOSA, Keysight
+
+voltage = float(sys.argv[1])
+shunt = float(sys.argv[2])
+
+mimfile = "../data/20161114_arb_%d.mim" % shunt
+csvfile = "../data/20161114_%d_arb.csv" % shunt
+
+mim = MIMOSA(voltage, shunt)
+ks = Keysight()
+
+charges, triggers = mim.load_data(mimfile)
+timestamps, currents = ks.load_data(csvfile)
+
+def calfunc330(charge):
+    if charge < 140.210488888889:
+        return 0
+    if charge <= 526.507377777778:
+        return float(charge) * 0.0941215500652876 + -13.196828549634
+    else:
+        return float(charge) * 0.0897304193584184 + -47.2437278033012 + 36.358862
+
+def calfunc82(charge):
+    if charge < 126.993600:
+        return 0
+    if charge <= 245.464889:
+        return charge * 0.306900 + -38.974361
+    else:
+        return charge * 0.356383 + -87.479495 + 36.358862
+
+
+def calfunc33(charge):
+    if charge < 127.000000:
+        return 0
+    if charge <= 127.211911:
+        return charge * 171.576006 + -21790.152700
+    else:
+        return charge * 0.884357 + -112.500777 + 36.358862
+
+calfuncs = {
+    33 : calfunc33,
+    82 : calfunc82,
+    330 : calfunc330,
+}
+
+vcalfunc = np.vectorize(calfuncs[int(shunt)], otypes=[np.float64])
+
+#plt.plot(np.arange(0, 1000, 0.01), vcalfunc(np.arange(0, 1000, 0.01)))
+#plt.xlabel('Rohdatenwert')
+#plt.ylabel('Strom [µA]')
+#plt.show()
+#sys.exit(0)
+
+mim_x = np.arange(len(charges)-199) * 1e-5
+mim_y = running_mean(mim.charge_to_current_nocal(charges), 200) * 1e-6
+cal_y = running_mean(vcalfunc(charges), 200) * 1e-6
+ks_x = timestamps[:len(timestamps)-9]
+ks_y = running_mean(currents, 10)
+
+# look for synchronization opportunity in first 5 seconds
+mim_sync_idx = 0
+ks_sync_idx = 0
+for i in range(0, 500000):
+    if mim_sync_idx == 0 and mim_y[i] > 0.001:
+        mim_sync_idx = i
+    if ks_sync_idx == 0 and ks_y[i] > 0.001:
+        ks_sync_idx = i
+
+mim_x = mim_x - mim_x[mim_sync_idx]
+ks_x = ks_x - ks_x[ks_sync_idx]
+
+mim_max_start = int(len(mim_y) * 0.4)
+mim_max_end = int(len(mim_y) * 0.6)
+mim_start_end = int(len(mim_y) * 0.1)
+mim_end_start = int(len(mim_y) * 0.9)
+mim_max = np.max(mim_y[mim_max_start:mim_max_end])
+mim_min1 = np.min(mim_y[:mim_start_end])
+mim_min2 = np.min(mim_y[mim_end_start:])
+mim_center = 0
+mim_start = 0
+mim_end = 0
+for i, y in enumerate(mim_y):
+    if y == mim_max and i / len(mim_y) > 0.4 and i / len(mim_y) < 0.6:
+        mim_center = i
+    elif y == mim_min1 and i / len(mim_y) < 0.1:
+        mim_start = i
+    elif y == mim_min2 and i / len(mim_y) > 0.9:
+        mim_end = i
+
+plt.plot([mim_x[mim_center]], [mim_y[mim_center]], "yo")
+plt.plot([mim_x[mim_start]], [mim_y[mim_start]], "yo")
+plt.plot([mim_x[mim_end]], [mim_y[mim_end]], "yo")
+#
+mimhandle, = plt.plot(mim_x, mim_y, "r-", label='MIMOSA')
+#calhandle, = plt.plot(mim_x, cal_y, "g-", label='MIMOSA (autocal)')
+kshandle, = plt.plot(ks_x, ks_y, "b-", label='Keysight')
+#plt.legend(handles=[mimhandle, calhandle, kshandle])
+plt.xlabel('Zeit [s]')
+plt.ylabel('Strom [A]')
+plt.grid(True)
+
+ks_steps_up = []
+ks_steps_down = []
+mim_steps_up = []
+mim_steps_down = []
+
+skip = 0
+for i, gradient in enumerate(np.gradient(ks_y, 10000)):
+    if gradient > 0.5e-9 and i - skip > 200 and ks_x[i] < mim_x[mim_center] and ks_x[i] > 5:
+        plt.plot([ks_x[i]], [ks_y[i]], "go")
+        ks_steps_up.append(i)
+        skip = i
+    elif gradient < -0.5e-9 and i - skip > 200 and ks_x[i] > mim_x[mim_center] and ks_x[i] < mim_x[mim_end]:
+        plt.plot([ks_x[i]], [ks_y[i]], "g*")
+        ks_steps_down.append(i)
+        skip = i
+
+j = 0
+for i, ts in enumerate(mim_x):
+    if j < len(ks_steps_up) and ts > ks_x[ks_steps_up[j]]:
+        mim_steps_up.append(i)
+        j += 1
+
+j = 0
+for i, ts in enumerate(mim_x):
+    if j < len(ks_steps_down) and ts > ks_x[ks_steps_down[j]]:
+        mim_steps_down.append(i)
+        j += 1
+
+print(ks_steps_up)
+print(mim_steps_up)
+
+mim_values = []
+cal_values = []
+ks_values = []
+
+for i in range(1, len(ks_steps_up)):
+    mim_values.append(np.mean(mim_y[mim_steps_up[i-1]:mim_steps_up[i]]))
+    cal_values.append(np.mean(cal_y[mim_steps_up[i-1]:mim_steps_up[i]]))
+    ks_values.append(np.mean(ks_y[ks_steps_up[i-1]:ks_steps_up[i]]))
+    print("step %d avg %5.3f vs %5.3f vs %5.3f mA" %
+        (i, np.mean(ks_y[ks_steps_up[i-1]:ks_steps_up[i]]) * 1e3,
+            np.mean(mim_y[mim_steps_up[i-1]:mim_steps_up[i]]) * 1e3,
+            np.mean(cal_y[mim_steps_up[i-1]:mim_steps_up[i]]) * 1e3))
+for i in range(1, len(ks_steps_down)):
+    mim_values.append(np.mean(mim_y[mim_steps_down[i-1]:mim_steps_down[i]]))
+    cal_values.append(np.mean(cal_y[mim_steps_down[i-1]:mim_steps_down[i]]))
+    ks_values.append(np.mean(ks_y[ks_steps_down[i-1]:ks_steps_down[i]]))
+    print("step %d avg %5.3f vs %5.3f vs %5.3f mA" %
+        (i, np.mean(ks_y[ks_steps_down[i-1]:ks_steps_down[i]]) * 1e3,
+            np.mean(mim_y[mim_steps_down[i-1]:mim_steps_down[i]]) * 1e3,
+            np.mean(cal_y[mim_steps_down[i-1]:mim_steps_down[i]]) * 1e3))
+
+mim_values = np.array(mim_values)
+cal_values = np.array(cal_values)
+ks_values = np.array(ks_values)
+
+plt.show()
+
+plt.hist(ks_y[ks_steps_up[48]:ks_steps_up[49]] * 1e3, 100, normed=0, facecolor='blue', alpha=0.8)
+plt.xlabel('mA Keysight')
+plt.ylabel('#')
+plt.grid(True)
+plt.show()
+plt.hist(mim_y[mim_steps_up[48]:mim_steps_up[49]] * 1e3, 100, normed=0, facecolor='blue', alpha=0.8)
+plt.xlabel('mA MimosaGUI')
+plt.ylabel('#')
+plt.grid(True)
+plt.show()
+
+mimhandle, = plt.plot(ks_values * 1e3, mim_values * 1e3, "ro", label='Unkalibriert', markersize=4)
+calhandle, = plt.plot(ks_values * 1e3, cal_values * 1e3, "bs", label='Kalibriert', markersize=4)
+plt.legend(handles=[mimhandle, calhandle])
+plt.xlabel('mA Keysight')
+plt.ylabel('mA MIMOSA')
+plt.grid(True)
+plt.show()
+
+mimhandle, = plt.plot(ks_values * 1e3, (mim_values - ks_values) * 1e3, "ro", label='Unkalibriert', markersize=4)
+calhandle, = plt.plot(ks_values * 1e3, (cal_values - ks_values) * 1e3, "bs", label='Kalibriert', markersize=4)
+plt.legend(handles=[mimhandle, calhandle])
+plt.xlabel('Sollstrom [mA]')
+plt.ylabel('Messfehler MIMOSA [mA]')
+plt.grid(True)
+plt.show()
+
+mimhandle, = plt.plot(ks_values * 1e3, (mim_values - ks_values) * 1e3, "r--", label='Unkalibriert')
+calhandle, = plt.plot(ks_values * 1e3, (cal_values - ks_values) * 1e3, "b-", label='Kalibriert')
+plt.legend(handles=[mimhandle, calhandle])
+plt.xlabel('Sollstrom [mA]')
+plt.ylabel('Messfehler MIMOSA [mA]')
+plt.grid(True)
+plt.show()
+
+mimhandle, = plt.plot(ks_values * 1e3, (mim_values - ks_values) / ks_values * 100, "ro", label='Unkalibriert', markersize=4)
+calhandle, = plt.plot(ks_values * 1e3, (cal_values - ks_values) / ks_values * 100, "bs", label='Kalibriert', markersize=4)
+plt.legend(handles=[mimhandle, calhandle])
+plt.xlabel('Sollstrom [mA]')
+plt.ylabel('Messfehler MIMOSA [%]')
+plt.grid(True)
+plt.show()
+
+mimhandle, = plt.plot(ks_values * 1e3, (mim_values - ks_values) / ks_values * 100, "r--", label='Unkalibriert')
+calhandle, = plt.plot(ks_values * 1e3, (cal_values - ks_values) / ks_values * 100, "b-", label='Kalibriert')
+plt.legend(handles=[mimhandle, calhandle])
+plt.xlabel('Sollstrom [mA]')
+plt.ylabel('Messfehler MIMOSA [%]')
+plt.grid(True)
+plt.show()
+
+#mimhandle, = plt.plot(mim_x, np.gradient(mim_y, 10000), "r-", label='MIMOSA')
+#kshandle, = plt.plot(ks_x, np.gradient(ks_y, 10000), "b-", label='Keysight')
+#plt.legend(handles=[mimhandle, kshandle])
+#plt.xlabel('Zeit [s]')
+#plt.ylabel('Strom [A]')
+#plt.show()
diff --git a/bin/mimosawatch b/bin/mimosawatch
new file mode 100755
index 0000000..88d34a3
--- /dev/null
+++ b/bin/mimosawatch
@@ -0,0 +1,24 @@
+#!/usr/bin/env perl
+
+use strict;
+use warnings;
+use 5.020;
+
+use File::Slurp qw(slurp);
+use IO::Handle;
+
+our $VERSION = '0.00';
+
+STDOUT->autoflush(1);
+
+while (sleep(1)) {
+	if ( -e '/tmp/mimosa/mimosa_scale_100000.tmp' ) {
+		my $raw_data = slurp('/tmp/mimosa/mimosa_scale_100000.tmp');
+		my @data = map { $_ >> 8 } unpack('L*', $raw_data);
+		my @buzzer = map { $_ & 0x0f } unpack('L*', $raw_data);
+
+		if (@data > 8) {
+			printf("\r\e[2K" . ('%8d ' x 8), @data[-8 .. -1]);
+		}
+	}
+}