#!/usr/bin/env python3
# vim:tabstop=4:softtabstop=4:shiftwidth=4:textwidth=160:smarttab:expandtab
import getopt
import itertools
import matplotlib.pyplot as plt
import numpy as np
import os
import re
import sys
from dfatool import aggregate_measures, running_mean, MIMOSA
opt = dict()
def show_help():
print('''mimosa-etv - MIMOSA Analyzer and Visualizer
USAGE
mimosa-etv [--skip <count>] [--threshold <power>] [--plot] [--stat] <voltage> <shunt> <file>
DESCRIPTION
mimosa-etv analyzes measurements taken via MIMOSA. Data can be plotted or aggregated on stdout.
OPTIONS
--skip <count>
Skip the first <count> data samples.
--threshold <watts>|mean
Partition data into points with mean power >= <watts> and points with
mean power < <watts>, and print some statistics. higher power is handled
as peaks, whereas low-power measurements constitute the baseline.
If the threshold is set to "mean", the mean power of all measurements
will be used
--threshold-peakcount <num>
Automatically determine threshold so that there are exactly <num> peaks.
A peaks is a group of consecutive measurements with mean power >= threshold
--plot
Show power/time plot
--stat
Show mean voltage, current, and power as well as total energy consumption.
''')
def peak_search(data, lower, upper, direction_function):
while upper - lower > 1e-6:
bs_test = np.mean([lower, upper])
peakcount = itertools.groupby(data, lambda x: x >= bs_test)
peakcount = filter(lambda x: x[0] == True, peakcount)
peakcount = sum(1 for i in peakcount)
direction = direction_function(peakcount, bs_test)
if direction == 0:
return bs_test
elif direction == 1:
lower = bs_test
else:
upper = bs_test
return None
def peak_search2(data, lower, upper, check_function):
for power in np.arange(lower, upper, 1e-6):
peakcount = itertools.groupby(data, lambda x: x >= power)
peakcount = filter(lambda x: x[0] == True, peakcount)
peakcount = sum(1 for i in peakcount)
if check_function(peakcount, power) == 0:
return power
return None
if __name__ == '__main__':
try:
optspec = ('help skip= threshold= threshold-peakcount= plot stat')
raw_opts, args = getopt.getopt(sys.argv[1:], "", optspec.split(' '))
for option, parameter in raw_opts:
optname = re.sub(r'^--', '', option)
opt[optname] = parameter
if 'help' in opt:
show_help()
sys.exit(0)
if 'skip' in opt:
opt['skip'] = int(opt['skip'])
else:
opt['skip'] = 0
if 'threshold' in opt and opt['threshold'] != 'mean':
opt['threshold'] = float(opt['threshold'])
if 'threshold-peakcount' in opt:
opt['threshold-peakcount'] = int(opt['threshold-peakcount'])
except getopt.GetoptError as err:
print(err)
sys.exit(2)
except IndexError:
print('Usage: mimosa-etv <duration>')
sys.exit(2)
except ValueError:
print('Error: duration or skip is not a number')
sys.exit(2)
voltage, shunt, inputfile = args
voltage = float(voltage)
shunt = int(shunt)
mim = MIMOSA(voltage, shunt)
charges, triggers = mim.load_file(inputfile)
currents = mim.charge_to_current_nocal(charges) * 1e-6
powers = currents * voltage
if 'threshold-peakcount' in opt:
bs_mean = np.mean(powers)
# Finding the correct threshold is tricky. If #peaks < peakcont, our
# current threshold may be too low (extreme case: a single peak
# containing all measurements), but it may also be too high (extreme
# case: a single peak containing just one data point). Similarly,
# #peaks > peakcount may be due to baseline noise causing lots of
# small peaks, or due to peak noise (if the threshold is already rather
# high).
# For now, we first try a simple binary search:
# The threshold is probably somewhere around the mean, so if
# #peaks != peakcount and threshold < mean, we go up, and if
# #peaks != peakcount and threshold >= mean, we go down.
# If that doesn't work, we fall back to a linear search in 1 µW steps
def direction_function(peakcount, power):
if peakcount == opt['threshold-peakcount']:
return 0
if power < bs_mean:
return 1
return -1
threshold = peak_search(power, np.min(power), np.max(power), direction_function)
if threshold == None:
threshold = peak_search2(power, np.min(power), np.max(power), direction_function)
if threshold != None:
print('Threshold set to {:.0f} µW : {:.9f}'.format(threshold * 1e6, threshold))
opt['threshold'] = threshold
else:
print('Found no working threshold')
if 'threshold' in opt:
if opt['threshold'] == 'mean':
opt['threshold'] = np.mean(powers)
print('Threshold set to {:.0f} µW : {:.9f}'.format(opt['threshold'] * 1e6, opt['threshold']))
baseline_mean = 0
if np.any(powers < opt['threshold']):
baseline_mean = np.mean(powers[powers < opt['threshold']])
print('Baseline mean: {:.0f} µW : {:.9f}'.format(
baseline_mean * 1e6, baseline_mean))
if np.any(powers >= opt['threshold']):
print('Peak mean: {:.0f} µW : {:.9f}'.format(
np.mean(powers[powers >= opt['threshold']]) * 1e6,
np.mean(powers[powers >= opt['threshold']])))
peaks = []
peak_start = -1
for i, dp in enumerate(powers):
if dp >= opt['threshold'] and peak_start == -1:
peak_start = i
elif dp < opt['threshold'] and peak_start != -1:
peaks.append((peak_start, i))
peak_start = -1
total_energy = 0
delta_energy = 0
for peak in peaks:
duration = (peak[1] - peak[0]) * 1e-5
total_energy += np.mean(powers[peak[0] : peak[1]]) * duration
delta_energy += (np.mean(powers[peak[0] : peak[1]]) - baseline_mean) * duration
delta_powers = powers[peak[0] : peak[1]] - baseline_mean
print('{:.2f}ms peak ({:f} -> {:f})'.format(duration * 1000,
peak[0], peak[1]))
print(' {:f} µJ / mean {:f} µW'.format(
np.mean(powers[peak[0] : peak[1]]) * duration * 1e6,
np.mean(powers[peak[0] : peak[1]]) * 1e6 ))
measures = aggregate_measures(np.mean(delta_powers), delta_powers)
print(' {:f} µW delta mean = {:0.1f}% / {:f} µW error'.format(np.mean(delta_powers) * 1e6, measures['smape'], measures['rmsd'] * 1e6 ))
print('Peak energy mean: {:.0f} µJ : {:.9f}'.format(
total_energy * 1e6 / len(peaks), total_energy / len(peaks)))
print('Average per-peak energy (delta over baseline): {:.0f} µJ : {:.9f}'.format(
delta_energy * 1e6 / len(peaks), delta_energy / len(peaks)))
if 'stat' in opt:
mean_current = np.mean(currents)
mean_power = np.mean(powers)
print('Mean current: {:.0f} µA : {:.9f}'.format(mean_current * 1e6, mean_current))
print('Mean power: {:.0f} µW : {:.9f}'.format(mean_power * 1e6, mean_power))
if 'plot' in opt:
timestamps = np.arange(len(powers)) * 1e-5
pwrhandle, = plt.plot(timestamps, powers, 'b-', label='U*I', markersize=1)
plt.legend(handles=[pwrhandle])
plt.xlabel('Time [s]')
plt.ylabel('Power [W]')
plt.grid(True)
plt.show()