diff --git a/pybench/__init__.py b/pybench/__init__.py
new file mode 100644
index 0000000..9d17fe8
--- /dev/null
+++ b/pybench/__init__.py
@@ -0,0 +1,11 @@
+from os import path
+
+from benchmark import Benchmark, parser, timed  # NOQA: export
+
+# Use README as module documentation
+readme = path.join(path.dirname(__file__), '..', 'README.rst')
+if path.exists(readme):
+    with open(readme) as f:
+        __doc__ = f.read()
+del readme
+del f
diff --git a/pybench.py b/pybench/benchmark.py
similarity index 81%
rename from pybench.py
rename to pybench/benchmark.py
index eedc967..94a4ea9 100644
--- a/pybench.py
+++ b/pybench/benchmark.py
@@ -15,11 +15,10 @@
 import time
 from warnings import warn
 
-# Use README as module documentation
-readme = path.join(path.dirname(__file__), 'README.rst')
-if path.exists(readme):
-    with open(readme) as f:
-        __doc__ = f.read()
+import pandas as pd
+import xray
+
+from utils import value_combinations
 
 # Imports for plot, warn if those fail but do not die
 try:
@@ -95,7 +94,7 @@ def parser(**kwargs):
 
 class Benchmark(object):
     """An abstract base class for benchmarks."""
-    params = []
+    params = {}
     """The parameters to run the benchmark for: a list of pairs, each of which
     is the parameter name and a list of benchmark values."""
     repeats = 3
@@ -121,6 +120,8 @@ class Benchmark(object):
         * node_threshold: eliminate nodes below this threshold
         * edge_threshold: eliminate edges below this threshold
         * format: comma-separated list of output formats (supported by dot)"""
+    regions = ['total']
+    """Regions to time."""
     profileregions = ['total']
     """Regions to create profile graphs for."""
     meta = {}
@@ -130,7 +131,7 @@ class Benchmark(object):
     series = {}
     """Benchmark series created from several invocations of the script e.g.
     for parallel runs on variable number of processors."""
-    suffix = '.dat'
+    suffix = '.nc'
     """Suffix for the result file to write."""
 
     def __init__(self, **kwargs):
@@ -147,7 +148,6 @@ def __init__(self, **kwargs):
                 setattr(self, k, v)
         if isinstance(self.method, str):
             self.method = getattr(self, self.method, self.method)
-        self.regions = defaultdict(float)
         self.profiles = {}
         self.meta['benchmark_version'] = get_git_revision()
         self.meta['pybench_version'] = get_git_revision(path.dirname(__file__))
@@ -157,6 +157,15 @@ def __init__(self, **kwargs):
             self.meta['jobid'] = getenv('PBS_JOBID')
         if getenv('PBS_JOBNAME'):
             self.meta['jobname'] = getenv('PBS_JOBNAME')
+        self.timings = defaultdict(float)
+        self.data = xray.Dataset()
+
+    def _init_data(self, params=None):
+        params = dict(params or self.params)
+        shape = tuple(len(p) for p in params.values())
+        array = np.zeros(shape=shape)
+        array.fill(np.nan)
+        return xray.DataArray(array, coords=params, name=self.benchmark)
 
     @property
     def name(self):
@@ -174,13 +183,13 @@ def timed_region(self, name, normalize=1.0):
             self.profiles[name].enable()
         t_ = self.timer()
         yield
-        self.regions[name] += (self.timer() - t_) * normalize
+        self.timings[name] += (self.timer() - t_) * normalize
         if name in self.profiles:
             self.profiles[name].disable()
 
     def register_timing(self, name, value):
         """Register the timing `value` for the region identified by `name`."""
-        self.regions[name] += value
+        self.timings[name] += value
 
     def _args(self, kwargs):
         """Parse name, params and method from the kwargs dictionary."""
@@ -283,17 +292,13 @@ def profile(self, **kwargs):
         out = path.join(profiledir, name)
         if rank == 0 and not path.exists(profiledir):
             makedirs(profiledir)
-        if params:
-            pkeys, pvals = zip(*params)
-        else:
-            pkeys, pvals = (), ()
-        for pvalues in product(*pvals):
+        for param in value_combinations(params):
             if rank == 0:
-                print 'Profile', name, 'for parameters', ', '.join('%s=%s' % (k, v) for k, v in zip(pkeys, pvalues))
-            kwargs.update(dict(zip(pkeys, pvalues)))
+                print 'Profile', name, 'for parameters', ', '.join('%s=%s' % (k, v) for k, v in sorted(param.items()))
+            kwargs.update(param)
             # Dry run
             method(**kwargs)
-            suff = '_'.join('%s%s' % (k, v) for k, v in zip(pkeys, pvals))
+            suff = '_'.join('%s%s' % (k, v) for k, v in sorted(param.items()))
             for r in regions:
                 self.profiles[r] = Profile()
             if 'total' in regions:
@@ -327,29 +332,23 @@ def run(self, **kwargs):
         warmups = kwargs.pop('warmups', self.warmups)
         average = kwargs.pop('average', self.average)
 
-        timings = self.result.get('timings') or {}
-        self.result = {'name': name,
-                       'description': description,
-                       'params': sorted(params),
-                       'repeats': repeats,
-                       'warmups': warmups,
-                       'average': average.__name__,
-                       'method': method.__name__,
-                       'regions': self.regions.keys(),
-                       'meta': self.meta,
-                       'series': self.series,
-                       'timings': timings}
-        if params:
-            pkeys, pvals = zip(*sorted(params))
-        else:
-            pkeys, pvals = (), ()
+        # Check that all params are defined at the benchmark level
+        for param, vals in params.items():
+            if param not in self.params:
+                raise ValueError('Unknown parameter ' + param)
+            for v in vals:
+                if v not in self.params[param]:
+                    raise ValueError('Invalid value %s for parameter %s (valid: %s)' %
+                                     (v, param, self.params[param]))
+
         self.meta['start_time'] = str(datetime.now())
-        for pvalues in product(*pvals):
+
+        for param in value_combinations(params):
             if rank == 0:
-                pstr = ', '.join('%s=%s' % (k, v) for k, v in zip(pkeys, pvalues))
+                pstr = ', '.join('%s=%s' % (k, v) for k, v in param.items())
                 sstr = ', '.join('%s=%s' % (k, v) for k, v in self.series.items())
                 print 'Benchmark', name, 'for parameters', pstr, 'series', sstr
-            kwargs.update(dict(zip(pkeys, pvalues)))
+            kwargs.update(param)
 
             if rank == 0:
                 print '  Running', warmups, 'warmup runs'
@@ -357,22 +356,27 @@ def run(self, **kwargs):
                 method(**kwargs)
 
             def bench():
-                self.regions = defaultdict(float)
+                self.timings = defaultdict(float)
                 with self.timed_region('total'):
                     method(**kwargs)
-                return self.regions
+                return self.timings
             if rank == 0:
                 print '  Running', repeats, 'benchmark runs'
             times = [bench() for _ in range(repeats)]
-            # Average over all timed regions
-            times = dict((k, average(d[k] for d in times))
-                         for k in self.regions.keys())
-            if pvalues:
-                timings[pvalues] = times
-            else:
-                self.result['timings'] = times
+            # Average over all timings
+            if times:
+                for k in self.timings.keys():
+                    self[k].loc[param] = average(d[k] for d in times)
         self.meta['end_time'] = str(datetime.now())
-        return self.result
+        return self
+
+    def __getitem__(self, region):
+        if region not in self.data:
+            self.data[region] = self._init_data(params=self.params)
+        return self.data[region]
+
+    def __call__(self, region, **kwargs):
+        return self.data[region].loc[kwargs]
 
     def _file(self, filename=None, suffix=None):
         """Return a filepath specified by given `filename` and `suffix`, which
@@ -386,133 +390,62 @@ def _file(self, filename=None, suffix=None):
         return path.join(self.resultsdir, filename + suffix)
 
     def _read(self, filename=None, suffix=None):
-        """Read a file specified by given `filename` and `suffix`, which
-        default to the global name and suffix attributes if not given, and
-        evaluate its contents."""
-        with open(self._file(filename, suffix)) as f:
-            return eval(f.read())
+        """Read a netCDF file specified by given `filename` and `suffix`, which
+        default to the global name and suffix attributes if not given."""
+        return xray.open_dataset(self._file(filename, suffix))
 
     def load(self, filename=None, suffix=None):
         """Load results from a file specified by given `filename` and `suffix`,
         which default to the global name and suffix attributes if not given."""
         try:
-            self.result = self._read(filename)
+            self.data = self._read(filename)
         except IOError:
-            self.result = {}
-        return self.result
+            self.data = xray.Dataset()
+        return self
 
     def save(self, filename=None, suffix=None):
         """Save results to a file specified by given `filename` and `suffix`,
         which default to the global name and suffix attributes if not given."""
         if rank > 0:
             return
-        with open(self._file(filename, suffix), 'w') as f:
-            pprint(self.result, f)
-
-    def combine(self, files):
-        """Combine results given by the dictionary `files`, with file names as
-        keys and prefixes as values. The prefix is prepended to the regions."""
-        result = {'name': self.name, 'series': self.series}
-        timings = defaultdict(dict)
-        regions = set()
-        for name, pref in files.items():
-            res = self._read(name)
-            for key in ['description', 'meta', 'params']:
-                result[key] = res[key]
-            for k, v in res['timings'].items():
-                # Parametrized benchmark
-                if isinstance(v, dict):
-                    for r, t in v.items():
-                        timings[k][pref + ' ' + r] = t
-                        regions.add(pref + ' ' + r)
-                # Non-parametrized benchmark
-                else:
-                    timings[pref + ' ' + k] = v
-                    regions.add(pref + ' ' + k)
-        result['timings'] = timings
-        result['regions'] = list(regions)
-        self.result = result
-        return result
-
-    def combine_series(self, series, filename=None, aggregate={}, merge=False):
+        self.data.to_netcdf(self._file(filename, suffix))
+        return self
+
+    def combine(self, name, labels, files):
+        """Combine results read from the file names in `files` along a new
+        dimension `name` with the labels `labels`."""
+        arrays = [self._read(f) if isinstance(f, str) else f for f in files]
+        self.data = xray.concat(arrays, pd.Index(labels, name=name))
+        return self
+
+    def combine_series(self, series, filename=None, coords=None):
         """Combine the results of one or more series of benchmarks.
 
         :param series: a dictionary with the series names as keys and the list
             of values defing the series as value.
         :param filename: the basename for the files to combine (defaults to the
             global name property if not given)
-        :param aggregate: dictionary of regions to aggregate where the key is
-            the resulting region and the value is a list of regions to sum
-        :param merge: if set to `True`, the given series is merged with the
-            existing parameter values. The default setting of `False` assumes
-            that all series are added as new parameters.
+        :param coords: dictionary of coordinate axes to relabel, where the key
+            is the coordinate and the value is the list of new labels
         """
-        filename = filename or self.name
-        if merge:
-            pkeys, pvals = zip(*sorted(self.params))
-            for k, v in sorted(series):
-                # The key already exists in the params
-                if k in pkeys:
-                    i = pkeys.index(k)
-                    for p in v:
-                        if p not in pvals[i]:
-                            self.params[i][1].append(p)
-                if k not in pkeys:
-                    self.params.append((k, v))
-        else:
-            self.params = self.params + series
-            pkeys, pvals = zip(*sorted(self.params))
-        result = {'name': self.name, 'params': self.params}
-        timings = self.result.get('timings') or {}
-        skeys, svals = zip(*sorted(series))
-        for svalues in product(*svals):
-            suff = '_'.join('%s%s' % (k, v) for k, v in zip(skeys, svalues))
+        # Add the series as additional dimensions to the parameter space
+        self.params.update(series)
+        self.data = xray.Dataset()
+
+        filename = filename or self.benchmark
+
+        # Read file for each combination of the series
+        for s in value_combinations(series):
+            suff = '_'.join('%s%s' % (k, v) for k, v in sorted(s.items()))
             fname = '%s_%s' % (filename, suff)
-            try:
-                res = self._read(fname)
-            except IOError:
-                warn("Series not found: " + str(svalues))
-                continue
-            for key in ['description', 'meta', 'regions']:
-                result[key] = res[key]
-            for target, regions in aggregate.items():
-                # FIXME: this won't currently work with a param series
-                if all([r in res['timings'] for r in regions]):
-                    res['timings'][target] = sum(res['timings'][region]
-                                                 for region in regions)
-            if pkeys == skeys:
-                timings[svalues] = res['timings']
-            else:
-                rkeys = zip(*res['params'])[0]
-                for k, v in res['timings'].items():
-                    key = zip(*sorted(zip(rkeys, k) + zip(skeys, svalues)))[1]
-                    timings[key] = v
-        result['timings'] = timings
-        self.result = result
-        return result
+            for k, v in self._read(fname).data_vars.items():
+                self[k].loc[s] = v.values
 
-    def dataframe(self, **kwargs):
-        """Return results as a pandas DataFrame
+        # Re-label coordinates if requested
+        for k, v in (coords or {}).items():
+            self.data.coords[k] = v
 
-        :param kwargs: keyword arguments override values given in the results
-            * params: benchmark parameters
-            * regions: regions to output
-            * timings: benchmark timings
-            * skip: parameters to skip
-        """
-        import pandas as pd
-        params = kwargs.pop('params', self.result['params'])
-        regions = kwargs.pop('regions', self.result['regions'])
-        timings = kwargs.pop('timings', self.result['timings'])
-        skip = kwargs.pop('skip', [])
-
-        pkeys, pvals = zip(*sorted(params))
-        idx = [pkeys.index(s) for s in skip]
-        df = pd.DataFrame([dict(list((pkeys[i], p) for i, p in enumerate(pv)
-                                     if i not in idx) +
-                                list((r, timings[pv][r]) for r in regions))
-                           for pv in product(*pvals)])
-        return df.set_index([p for p in pkeys if p not in skip])
+        return self
 
     def table(self, **kwargs):
         """Export results as html or latex table (requires pandas).
@@ -521,27 +454,22 @@ def table(self, **kwargs):
             * filename: base name of output file
             * dataframe: pandas DataFrame to export (when given, the keyword
               arguments params, regions, skip and timings are ignored)
-            * params: benchmark parameters
             * tabledir: output directory
             * regions: regions to output
-            * timings: benchmark timings
-            * skip: parameters to skip
             * format: comma-separated list of output formats (html, latex, both)
+            * precision: precision of floating point values
         """
         if rank > 0:
             return
-        filename = kwargs.pop('filename', self.result['name'])
-        df = kwargs.get('dataframe')
-        if df is None:
-            params = kwargs.pop('params', self.result['params'])
-            regions = kwargs.pop('regions', self.result['regions'])
-            skip = kwargs.pop('skip', [])
-            timings = kwargs.pop('timings', self.result['timings'])
-            df = self.dataframe(params=params, regions=regions, skip=skip, timings=timings)
+        filename = kwargs.pop('filename', self.name)
+        df = kwargs.get('dataframe') or self.data.to_dataframe()
+        if kwargs.get('regions'):
+            df = df[kwargs.get('regions')]
         tabledir = kwargs.pop('tabledir', self.tabledir)
         formats = kwargs.pop('format', 'html').split(',')
         if not path.exists(tabledir):
             makedirs(tabledir)
+        pd.set_option('display.precision', kwargs.pop('precision', 4))
 
         # Reset the index only if it is a MultiIndex
         if hasattr(df.index, 'levels'):
@@ -553,6 +481,7 @@ def table(self, **kwargs):
         for fmt in formats:
             with open(path.join(tabledir, "%s.%s" % (filename, fmt)), 'w') as f:
                 f.write(render[fmt](df))
+        return self
 
     def lookup(self, region, params, keyset=()):
         """Retrieve a specific timing from benchmark results
@@ -565,9 +494,7 @@ def lookup(self, region, params, keyset=()):
         if isinstance(params, list):
             params = dict(params)
         params.update(keyset)
-        pvals = zip(*sorted(params.items()))[1]
-        timings = self.result['timings'].get(pvals)
-        return timings[region] if timings is not None else np.nan
+        return self(region=region, **params)
 
     def subplot(self, ax, xaxis, kind='plot', **kwargs):
         """Plot a graph into the given axes
@@ -642,9 +569,8 @@ def subplot(self, ax, xaxis, kind='plot', **kwargs):
         legend = kwargs.get('legend', {'loc': 'best'})
         lines = kwargs.get('lines', [])
         linewidth = kwargs.pop('linewidth', 2)
-        regions = kwargs.pop('regions', self.result['regions'])
+        regions = kwargs.pop('regions', self.regions)
         ticksize = kwargs.get('ticksize')
-        timings = kwargs.pop('timings', self.result['timings'])
         title = kwargs.pop('title', self.name)
         transform = kwargs.get('transform')
         xlabel = kwargs.pop('xlabel', None)
@@ -666,7 +592,6 @@ def subplot(self, ax, xaxis, kind='plot', **kwargs):
         groups = dict(kwargs.pop('groups'))
         groups, gvals = zip(*groups.items()) if groups else ([], [])
         params = dict(kwargs.pop('params'))
-        pvals = zip(*sorted(params))[1]
 
         nregions = len(regions)
         ngroups = int(np.prod([len(g) for g in gvals]))
@@ -860,9 +785,9 @@ def plot(self, xaxis, **kwargs):
         """
         if rank > 0:
             return
-        figname = kwargs.pop('figname', self.result['name'])
+        figname = kwargs.pop('figname', self.name)
         figsize = kwargs.pop('figsize', (9, 6))
-        params = dict(kwargs.pop('params', self.result['params']))
+        params = dict(kwargs.pop('params', self.params))
         groups = kwargs.get('groups', [])
         legend = kwargs.get('legend', {'loc': 'best'})
         format = kwargs.pop('format', 'svg')
@@ -885,9 +810,6 @@ def plot(self, xaxis, **kwargs):
 
         kwargs['xvals'] = kwargs.pop('xvals', params.pop(xaxis))
         kwargs['groups'] = zip(groups, [params.pop(g) for g in groups])
-        pkeys, pvals = zip(*sorted(params.items()))
-
-        nv = len(list(product(*pvals)))
 
         def save(fig, fname, outline, extra_artists=[]):
             if not format:
@@ -903,14 +825,14 @@ def save(fig, fname, outline, extra_artists=[]):
                         outline += ['<td><img src="%s"></td>' % fname]
             plt.close(fig)
 
+        param_combinations = value_combinations(params)
         for kind in kinds.split(','):
             outline = []
             if subplot:
                 axes = []
                 fig = plt.figure(figname + '_' + kind, figsize=figsize, dpi=300)
-            for p, pv in enumerate(product(*pvals), 1):
-                pdict = zip(pkeys, pv)
-                fsuff = '_'.join('%s%s' % (k, str(v).replace('.', '_')) for k, v in pdict)
+            for p, param in enumerate(param_combinations, 1):
+                fsuff = '_'.join('%s%s' % (k, str(v).replace('.', '_')) for k, v in param.items())
                 # Append speedup to file base name if any
                 if speedup:
                     fsuff += '_speedup' + ''.join(speedup)
@@ -926,7 +848,7 @@ def save(fig, fname, outline, extra_artists=[]):
                             kargs['title'] = None
                             kargs['axis'] = 'tight'
                             kargs.update(subplotargs[r, c])
-                            self.subplot(ax[r][c], xaxis, kind, params=pdict, **kargs)
+                            self.subplot(ax[r][c], xaxis, kind, params=param, **kargs)
                     # Adjust space between subplots
                     fig.subplots_adjust(hspace=hspace, wspace=wspace)
                     if title:
@@ -943,7 +865,7 @@ def save(fig, fname, outline, extra_artists=[]):
                     save(fig, '%s_%s_%s' % (figname, kind, fsuff), outline, extra_artists)
                     outline += ['</tr>']
                 elif subplot:
-                    ax = fig.add_subplot(1, nv, p, sharey=(axes[p-2] if p > 1 else None))
+                    ax = fig.add_subplot(1, len(param_combinations), p, sharey=(axes[p-2] if p > 1 else None))
                     axes.append(ax)
                     kargs = copy(kwargs)
                     kargs['legend'] = False
@@ -952,11 +874,11 @@ def save(fig, fname, outline, extra_artists=[]):
                         kargs['ylabel'] = None
                     if subplotargs:
                         kargs.update(subplotargs[p])
-                    self.subplot(ax, xaxis, kind, params=pdict, **kargs)
+                    self.subplot(ax, xaxis, kind, params=param, **kargs)
                 else:
                     fig = plt.figure(figname + '_' + fsuff, figsize=figsize, dpi=300)
                     ax = fig.add_subplot(111)
-                    self.subplot(ax, xaxis, kind, params=pdict, **kwargs)
+                    self.subplot(ax, xaxis, kind, params=param, **kwargs)
                     outline += ['<tr>']
                     save(fig, '%s_%s_%s' % (figname, kind, fsuff), outline)
                     outline += ['</tr>']
diff --git a/pybench/utils.py b/pybench/utils.py
new file mode 100644
index 0000000..fe004ed
--- /dev/null
+++ b/pybench/utils.py
@@ -0,0 +1,18 @@
+from collections import OrderedDict
+from itertools import product
+
+
+def value_combinations(dct):
+    """Return a list of dictionaries with all combinations of values from
+    dictionary `dct`, where each value is a list. The input ::
+
+        {'a': [1, 2], 'b': [3, 4, 5]}
+
+    yields ::
+
+        [{'a': 1, 'b': 3}, {'a': 1, 'b': 4}, {'a': 1, 'b': 5},
+         {'a': 2, 'b': 3}, {'a': 2, 'b': 4}, {'a': 2, 'b': 5}]
+    """
+    dct = dict(dct)
+    keys = sorted(dct)
+    return [OrderedDict(zip(keys, p)) for p in product(*(dct[k] for k in keys))]
diff --git a/setup.py b/setup.py
index 0bf021f..c34cc15 100644
--- a/setup.py
+++ b/setup.py
@@ -22,4 +22,4 @@
           'Programming Language :: Python :: 2.6',
           'Programming Language :: Python :: 2.7',
       ],
-      py_modules=['pybench'])
+      packages=['pybench'])
diff --git a/test_pybench.py b/test_pybench.py
index cb7d899..e263a4f 100644
--- a/test_pybench.py
+++ b/test_pybench.py
@@ -1,45 +1,67 @@
-from itertools import product
-from time import sleep
+import pytest
 
 from pybench import Benchmark
 
 
+class TestBenchmark(Benchmark):
+    params = {'a': range(1, 4), 'b': range(2, 7, 2)}
+    regions = ['total', 'r1', 'r2']
+
+    def test(self, a=None, b=None):
+        with self.timed_region('r1'):
+            pass
+        with self.timed_region('r2'):
+            pass
+
+
 class TimedRegion(Benchmark):
+    regions = ['total', 'stuff']
+
+    def test(self):
+        with self.timed_region('stuff'):
+            pass
+
+
+class InvalidRegion(Benchmark):
     def test(self):
         with self.timed_region('stuff'):
             pass
 
 
 class Parametrized(Benchmark):
-    params = [('a', range(3)), ('b', range(3))]
+    params = {'a': range(3), 'b': range(3)}
 
     def test(self, a=None, b=None):
         pass
 
 
 def test_no_params():
-    assert Benchmark().run(method=lambda: None)['timings']['total'] > 0.0
+    assert Benchmark().run(method=lambda: None)(region='total') > 0.0
 
 
 def test_no_repeats():
-    assert not Benchmark().run(method=lambda: None, repeats=0)['timings']
+    assert all(Benchmark().run(method=lambda: None, repeats=0).data.isnull())
+
+
+def test_invalid_region():
+    with pytest.raises(ValueError):
+        InvalidRegion().run()
+
 
+def test_invalid_param():
+    with pytest.raises(ValueError):
+        Parametrized().run(params={'c': [1]})
 
-def test_sleep():
-    def myfunc(n, duration):
-        for _ in range(n):
-            sleep(duration)
-    times = Benchmark().run(method=myfunc,
-                            params=[('n', range(3)),
-                                    ('duration', (0.001, 0.002))])
-    for (n, d), t in times['timings'].items():
-        assert abs(n*d - t['total']) < 1e-3
+
+def test_invalid_param_value():
+    with pytest.raises(ValueError):
+        Parametrized().run(params={'a': [3]})
 
 
 def test_timed_region():
-    result = TimedRegion().run()
-    assert result['timings']['total'] > 0.0
-    assert result['timings']['stuff'] > 0.0
+    bench = TimedRegion().run()
+    assert bench(region='total') > 0.0
+    assert bench(region='stuff') > 0.0
 
 
 def test_save_load(tmpdir):
@@ -50,32 +72,26 @@ def test_save_load(tmpdir):
     TimedRegion().load(d) == result
 
 
-def test_combine_regions(tmpdir):
-    b = TimedRegion()
-    da = tmpdir.join('a').strpath
-    db = tmpdir.join('b').strpath
-    keys = b.run()['timings'].keys()
-    b.save(da)
-    b.save(db)
-    result = TimedRegion().combine({da: 'a', db: 'b'})
-    assert all('a ' + k in result['timings'] for k in keys)
-    assert all('b ' + k in result['timings'] for k in keys)
+def test_combine_arrays(tmpdir):
+    a = TestBenchmark().run()
+    b = TestBenchmark().run()
+    c = TestBenchmark().combine('c', ['a', 'b'], [a.data, b.data])
+    assert 'c' in c.data.coords
+    assert all(c.data.coords['c'] == ['a', 'b'])
+    assert (c(c='a') == a.data).all()
+    assert (c(c='b') == b.data).all()
 
 
 def test_parametrized():
-    result = Parametrized().run()
-    _, pvalues = zip(*result['params'])
-    assert all(result['timings'][p]['total'] > 0.0
-               for p in product(*pvalues))
+    assert (Parametrized().run().data > 0.0).all()
 
 
 def test_combine_parametrized(tmpdir):
-    b = Parametrized()
     da = tmpdir.join('a').strpath
     db = tmpdir.join('b').strpath
-    params = b.run()['timings'].keys()
-    b.save(da)
-    b.save(db)
-    result = Parametrized().combine({da: 'a', db: 'b'})
-    assert all('a total' in result['timings'][p] for p in params)
-    assert all('b total' in result['timings'][p] for p in params)
+    b = Parametrized().run().save(da).save(db)
+    c = Parametrized().combine('c', ['a', 'b'], [da, db])
+    assert 'c' in c.data.coords
+    assert all(c.data.coords['c'] == ['a', 'b'])
+    assert (c(c='a') == b.data).all()
+    assert (c(c='b') == b.data).all()