diff --git a/examples/python/kernel_stat_collector.py b/examples/python/kernel_stat_collector.py
new file mode 100644
index 0000000000000000000000000000000000000000..5924b6b30ca7524d12050764b6e43dffe7eb4917
--- /dev/null
+++ b/examples/python/kernel_stat_collector.py
@@ -0,0 +1,88 @@
+import loopy as lp
+import numpy as np
+from loopy.kernel_stat_collector.kernel_stat_collector import KernelStatCollector
+from loopy.kernel_stat_collector.kernel_stat_collector import KernelStatOptions as kso  # noqa
+
+
+knl = lp.make_kernel(
+        "{[i,k,j]: 0<=i<n and 0<=k<m and 0<=j<ell}",
+        [
+            "c[i, j] = sum(k, a[i, k]*b[k, j])"
+        ],
+        name="matmul",
+        assumptions="n,m,ell >= 1",
+        lang_version=(2018, 2),
+        )
+
+knl = lp.add_and_infer_dtypes(knl, dict(a=np.float32, b=np.float32))
+
+lsize = 16
+knl = lp.split_iname(knl, "i", lsize, outer_tag="g.0", inner_tag="l.1")
+knl = lp.split_iname(knl, "j", lsize, outer_tag="g.1", inner_tag="l.0")
+knl = lp.split_iname(knl, "k", lsize)
+knl = lp.add_prefetch(knl, "a", ["k_inner", "i_inner"], default_tag="l.auto")
+knl = lp.add_prefetch(knl, "b", ["j_inner", "k_inner"], default_tag="l.auto")
+
+n = 2**8
+m = 2**9
+ell = 2**7
+param_dict = {'n': n, 'm': m, 'ell': ell}
+
+kp = KernelStatCollector(
+        #"NVIDIA", "GEFORCE",
+        #"NVIDIA", "K40C",
+        evaluate_polys=False,
+        count_madds=False,  # TODO enables this after madd counting branch is merged
+        include_kernel_params_in_ptx_filename=True,
+        )
+stats = kp.collect_stats(
+        knl,
+        [
+            kso.WALL_TIME,
+            kso.MEM_ACCESS_MAP,
+            kso.OP_MAP,
+            kso.SYNC_MAP,
+            kso.GRID_SIZES,
+            kso.FLOP_RATE,
+            kso.MEM_BANDWIDTH,
+            kso.SAVE_PTX,
+            kso.GENERATED_CODE,
+        ],
+        param_dict=param_dict,
+        )
+
+print(stats[kso.GENERATED_CODE])
+
+inv_giga = 10**-9
+
+print("\nWall time:")
+print(stats[kso.WALL_TIME])
+print("\nMem access map:")
+print(lp.stringify_stats_mapping(stats[kso.MEM_ACCESS_MAP]))
+print("\nOp map:")
+print(lp.stringify_stats_mapping(stats[kso.OP_MAP]))
+print("\nSync map:")
+print(lp.stringify_stats_mapping(stats[kso.SYNC_MAP]))
+print("\nGrid sizes:")
+print(stats[kso.GRID_SIZES])
+print("\nFlop rate (GFLOP/s):")
+print(stats[kso.FLOP_RATE]*inv_giga)
+print("\nMem throughput rate (GB/s) (total data accessed, data footprint only):")
+print(stats[kso.MEM_BANDWIDTH][0]*inv_giga,
+    stats[kso.MEM_BANDWIDTH][1]*inv_giga, "\n")
+
+kp2 = kp.copy(evaluate_polys=True)
+
+stats = kp2.collect_stats(
+        knl, stats, param_dict=param_dict)
+
+print("Now change eval_polys to True =========================================")
+print("\nMem access map:")
+print(lp.stringify_stats_mapping(stats[kso.MEM_ACCESS_MAP]))
+print("\nOp map:")
+print(lp.stringify_stats_mapping(stats[kso.OP_MAP]))
+print("\nSync map:")
+print(lp.stringify_stats_mapping(stats[kso.SYNC_MAP]))
+
+print("\nWall time:")
+print(stats[kso.WALL_TIME])
diff --git a/loopy/kernel_stat_collector/kernel_stat_collector.py b/loopy/kernel_stat_collector/kernel_stat_collector.py
new file mode 100644
index 0000000000000000000000000000000000000000..7fa936d217c3ee1e8600051e5df71b5f6d637ab2
--- /dev/null
+++ b/loopy/kernel_stat_collector/kernel_stat_collector.py
@@ -0,0 +1,501 @@
+import pyopencl as cl
+import loopy as lp
+from loopy.preprocess import prepare_for_caching
+
+
+def create_rand_args(ctx, knl, param_dict):
+
+    queue = cl.CommandQueue(ctx)
+    info = lp.generate_code_v2(knl).implemented_data_info
+    args, arg_data = lp.auto_test.make_ref_args(knl, info, queue, param_dict)
+    args.clear()
+    del args
+
+    rand_args = lp.auto_test.make_args(
+        knl, info, queue, arg_data, param_dict)
+
+    del arg_data[:]
+    del arg_data
+    return rand_args
+
+
+def find_cl_device_candidates(platform_name, device_name):
+    candidates = [
+        dev
+        for plat in cl.get_platforms()
+        if platform_name.lower() in plat.name.lower()
+        for dev in plat.get_devices()
+        if device_name.lower() in dev.name.lower()]
+
+    if not candidates:
+        raise ValueError("no device matched the platform/device ID tuple "
+                "'%s' and '%s'" % (platform_name, device_name))
+
+    return candidates
+
+
+def get_cached_cl_context_for_knl(
+            knl, ctx_cache, platform_name=None, device_name=None):
+
+    if knl.target is not None and knl.target.device is not None:
+        # kernel has a device already, see if we can use it
+        knl_platform_name = knl.target.device.platform.name
+        knl_device_name = knl.target.device.name
+
+        # check for mismatch between platforms/devices
+        if (platform_name is not None
+                and platform_name not in knl_platform_name) or (
+                device_name is not None
+                and device_name not in knl_device_name):
+            raise ValueError("kernel target platform %s and/or device %s "
+                "do not match provided platform %s and/or device %s."
+                % (knl_platform_name, knl_device_name,
+                platform_name, device_name))
+
+        cache_key = (knl_platform_name, knl_device_name, "ctx")
+        try:
+            return ctx_cache[cache_key]
+        except KeyError:
+            ctx = cl.Context([find_cl_device_candidates(
+                knl_platform_name, knl_device_name)[-1]]
+                )
+            ctx_cache[cache_key] = ctx
+            return ctx
+
+    elif platform_name is None or device_name is None:
+        # platform+device not provided, and knl did not specify plat+dev
+        ctx = cl.create_some_context()  # interactive mode
+        platform_name = ctx.devices[0].platform.name
+        device_name = ctx.devices[0].name
+        ctx_cache[(platform_name, device_name, "ctx")] = ctx
+        return ctx
+
+    else:
+        # platform + device provided and knl did not specify plat+dev
+        cache_key = (platform_name, device_name, "ctx")
+        try:
+            return ctx_cache[cache_key]
+        except KeyError:
+            ctx = cl.Context([find_cl_device_candidates(
+                platform_name, device_name)[-1]]
+                )
+            ctx_cache[cache_key] = ctx
+            return ctx
+
+
+def get_avg_exec_time(
+            knl,
+            param_dict,
+            ctx,  # TODO make ctx+queue optional params
+            queue,
+            n_time_trials=100,
+            n_warmup_time_trials=10,
+            ):
+
+    if param_dict is None:
+        raise ValueError(
+            "Wall time requires dictionary of kernel parameters.")
+
+    arg_arrays = create_rand_args(ctx, knl, param_dict)
+    knl = lp.set_options(knl, no_numpy=True)
+    compiled = lp.CompiledKernel(ctx, knl)
+
+    wtimes = []
+
+    import time
+    for t in range(n_time_trials + n_warmup_time_trials):
+        queue.finish()
+        tstart = time.time()
+        evt, out = compiled(queue, **arg_arrays)
+        queue.finish()
+        tend = time.time()
+        wtimes.append(tend-tstart)
+
+    import numpy as np
+    return np.average(wtimes[n_warmup_time_trials:])
+
+
+def write_ptx(ctx, knl, filename=None):
+    cl_program = cl.Program(
+        ctx, lp.generate_code_v2(knl).device_code()
+        ).build(options=knl.options.cl_build_options)
+    ptx_src = cl_program.binaries[0]
+    if filename:
+        ptx_src_file = open(filename, 'w')
+        ptx_src_file.write(ptx_src.decode('utf-8', 'ignore'))
+    else:
+        from loopy.diagnostic import warn_with_kernel
+        warn_with_kernel(knl, "write_ptx_no_filename",
+            "No filename passed to write_ptx for kernel %s, "
+            "writing to stdout" % (knl.name))
+        import sys
+        sys.stdout.write(ptx_src.decode('utf-8', 'ignore')+"\n")
+
+
+class KernelStatOptions:
+    WALL_TIME = "wall_time"
+    MEM_ACCESS_MAP = "mem_access_map"
+    OP_MAP = "op_map"
+    SYNC_MAP = "sync_map"
+    GRID_SIZES = "grid_sizes"
+    FLOP_RATE = "flop_rate"
+    MEM_BANDWIDTH = "mem_bandwidth"
+    SAVE_PTX = "save_ptx"
+    GENERATED_CODE = "generated_code"
+    # TODO mem access to footprint ratio
+
+
+class KernelStatCollector(object):
+
+    def __init__(
+                self,
+                platform_name=None,
+                device_name=None,
+                interactive=False,
+                n_warmup_time_trials=10,
+                n_time_trials=100,
+                evaluate_polys=True,
+                subgroup_size=32,
+                count_redundant_work=True,
+                count_madds=True,
+                count_within_subscripts=False,
+                include_kernel_params_in_ptx_filename=False,
+                ptx_filename_suffix="",
+                ):
+
+        # TODO create cache to store kernels (for executing w/different params)
+        self.ctx_cache = {}
+        self.stat_cache = {}
+        self.platform_name = platform_name
+        self.device_name = device_name
+        self.interactive = interactive
+
+        self.n_warmup_time_trials = n_warmup_time_trials
+        self.n_time_trials = n_time_trials
+
+        self.evaluate_polys = evaluate_polys
+        self.subgroup_size = subgroup_size
+        self.count_redundant_work = count_redundant_work
+        self.count_madds = count_madds
+        self.count_within_subscripts = count_within_subscripts
+
+        self.include_kernel_params_in_ptx_filename = \
+                include_kernel_params_in_ptx_filename
+        self.ptx_filename_suffix = ptx_filename_suffix
+
+    def copy(
+                self,
+                platform_name=None,
+                device_name=None,
+                interactive=None,
+                n_warmup_time_trials=None,
+                n_time_trials=None,
+                evaluate_polys=None,
+                subgroup_size=None,
+                count_redundant_work=None,
+                count_madds=None,
+                count_within_subscripts=None,
+                include_kernel_params_in_ptx_filename=None,
+                ptx_filename_suffix=None,
+                ):
+
+        platform_name_new = self.platform_name \
+            if platform_name is None else platform_name
+        device_name_new = self.device_name if device_name is None else device_name
+        interactive_new = self.interactive if interactive is None else interactive
+        n_warmup_time_trials_new = self.n_warmup_time_trials \
+            if n_warmup_time_trials is None else n_warmup_time_trials
+        n_time_trials_new = self.n_time_trials \
+            if n_time_trials is None else n_time_trials
+        evaluate_polys_new = self.evaluate_polys \
+            if evaluate_polys is None else evaluate_polys
+        count_redundant_work_new = self.count_redundant_work \
+            if count_redundant_work is None else count_redundant_work
+        subgroup_size_new = self.subgroup_size \
+            if subgroup_size is None else subgroup_size
+        count_madds_new = self.count_madds if count_madds is None else count_madds
+        count_within_subscripts_new = self.count_within_subscripts \
+            if count_within_subscripts is None else count_within_subscripts
+        include_kernel_params_in_ptx_filename_new = \
+            self.include_kernel_params_in_ptx_filename \
+            if include_kernel_params_in_ptx_filename is None \
+            else include_kernel_params_in_ptx_filename
+        ptx_filename_suffix_new = self.ptx_filename_suffix \
+            if ptx_filename_suffix is None else ptx_filename_suffix
+
+        stat_collector_new = KernelStatCollector(
+            platform_name=platform_name_new,
+            device_name=device_name_new,
+            interactive=interactive_new,
+            n_warmup_time_trials=n_warmup_time_trials_new,
+            n_time_trials=n_time_trials_new,
+            evaluate_polys=evaluate_polys_new,
+            subgroup_size=subgroup_size_new,
+            count_redundant_work=count_redundant_work_new,
+            count_madds=count_madds_new,
+            count_within_subscripts=count_within_subscripts_new,
+            include_kernel_params_in_ptx_filename=  # noqa
+                include_kernel_params_in_ptx_filename_new,
+            ptx_filename_suffix=ptx_filename_suffix_new,
+            )
+        stat_collector_new.ctx_cache = self.ctx_cache
+        stat_collector_new.stat_cache = self.stat_cache
+
+        return stat_collector_new
+
+    def get_cached_cl_context(self, knl):
+        ctx = get_cached_cl_context_for_knl(
+            knl, self.ctx_cache, self.platform_name, self.device_name)
+        if self.platform_name is None:
+            self.platform_name = ctx.devices[0].platform.name
+        if self.device_name is None:
+            self.device_name = ctx.devices[0].name
+        return ctx
+
+    def get_cached_stats_map(
+                self,
+                knl,
+                stat_option,  # KernelStatOptions
+                ):
+
+        cache_key = (prepare_for_caching(knl), stat_option)
+        # TODO avoid multiple calls to prepare_for_caching()?
+
+        try:
+            return self.stat_cache[cache_key]
+        except KeyError:
+            if stat_option == KernelStatOptions.MEM_ACCESS_MAP:
+                from loopy.statistics import get_mem_access_map
+                stats_map = get_mem_access_map(
+                    knl,
+                    count_redundant_work=self.count_redundant_work,
+                    subgroup_size=self.subgroup_size,
+                    )
+            elif stat_option == KernelStatOptions.OP_MAP:
+                from loopy.statistics import get_op_map
+                if self.count_madds:
+                    # TODO once madd counting branch is merged, remove conditional
+                    try:
+                        stats_map = get_op_map(
+                            knl,
+                            count_redundant_work=self.count_redundant_work,
+                            count_within_subscripts=self.count_within_subscripts,
+                            subgroup_size=self.subgroup_size,
+                            count_madds=self.count_madds,
+                            )
+                    except TypeError:
+                        raise NotImplementedError(
+                            "count_madds requires unmerged madd counting branch.")
+                else:
+                    stats_map = get_op_map(
+                        knl,
+                        count_redundant_work=self.count_redundant_work,
+                        count_within_subscripts=self.count_within_subscripts,
+                        subgroup_size=self.subgroup_size,
+                        )
+            elif stat_option == KernelStatOptions.SYNC_MAP:
+                from loopy.statistics import get_synchronization_map
+                stats_map = get_synchronization_map(
+                    knl,
+                    subgroup_size=self.subgroup_size,
+                    )
+            self.stat_cache[cache_key] = stats_map
+            return stats_map
+
+    def get_stats_map_and_evaluate_if_required(
+                self,
+                knl,
+                stat_option,  # KernelStatOptions
+                param_dict=None,
+                ):
+        stats_map = self.get_cached_stats_map(knl, stat_option)
+        if self.evaluate_polys:
+            if param_dict is None:
+                raise ValueError(
+                    "Cannot evaluate polynomials without param_dict.")
+            return stats_map.eval(param_dict)
+        else:
+            return stats_map
+
+    def get_grid_sizes(
+                self,
+                knl,
+                param_dict=None,
+                ):
+
+        cache_key = (prepare_for_caching(knl), KernelStatOptions.GRID_SIZES)
+        # TODO avoid multiple calls to prepare_for_caching()?
+
+        try:
+            grid_sizes = self.stat_cache[cache_key]
+        except KeyError:
+
+            global_size, local_size = knl.get_grid_size_upper_bounds()
+
+            from islpy import PwQPolynomial
+            gsize_pwqs = []
+            lsize_pwqs = []
+            for gsize in global_size:
+                gsize_pwqs.append(PwQPolynomial.from_pw_aff(gsize))
+            for lsize in local_size:
+                lsize_pwqs.append(PwQPolynomial.from_pw_aff(lsize))
+            grid_sizes = [gsize_pwqs, lsize_pwqs]
+            self.stat_cache[cache_key] = grid_sizes
+
+        if self.evaluate_polys:
+            if param_dict is None:
+                raise ValueError("Cannont evaluate polynomials without param_dict.")
+            return [g.eval_with_dict(param_dict) for g in grid_sizes[0]], \
+                   [l.eval_with_dict(param_dict) for l in grid_sizes[1]]
+        else:
+            return grid_sizes
+
+    def collect_stats(
+                self,
+                knl,
+                stat_options=[],
+                param_dict=None,
+                ):
+
+        stats_found = {}
+        kso = KernelStatOptions
+
+        if kso.WALL_TIME in stat_options or \
+                kso.FLOP_RATE in stat_options or \
+                kso.MEM_BANDWIDTH in stat_options:
+
+            if param_dict is None:
+                raise ValueError(
+                    "Wall time requires dictionary of kernel parameters.")
+
+            ctx = self.get_cached_cl_context(knl)
+            queue = cl.CommandQueue(ctx)
+
+            stats_found[kso.WALL_TIME] = get_avg_exec_time(
+                knl, param_dict, ctx, queue,
+                self.n_time_trials,
+                self.n_warmup_time_trials,
+                )
+
+        if kso.MEM_ACCESS_MAP in stat_options or \
+                kso.MEM_BANDWIDTH in stat_options:
+            stats_found[kso.MEM_ACCESS_MAP] = \
+                    self.get_stats_map_and_evaluate_if_required(
+                        knl,
+                        kso.MEM_ACCESS_MAP,
+                        param_dict=param_dict,
+                        )
+
+        if kso.OP_MAP in stat_options or \
+                kso.FLOP_RATE in stat_options:
+            stats_found[kso.OP_MAP] = \
+                self.get_stats_map_and_evaluate_if_required(
+                    knl,
+                    kso.OP_MAP,
+                    param_dict=param_dict,
+                    )
+
+        if kso.SYNC_MAP in stat_options:
+            stats_found[kso.SYNC_MAP] = \
+                self.get_stats_map_and_evaluate_if_required(
+                    knl,
+                    kso.SYNC_MAP,
+                    param_dict=param_dict,
+                    )
+
+        if kso.GRID_SIZES in stat_options:
+            stats_found[kso.GRID_SIZES] = self.get_grid_sizes(
+                knl,
+                param_dict=param_dict,
+                )
+
+        if kso.FLOP_RATE in stat_options:
+            import numpy as np
+            from loopy.statistics import _process_subgroup_size
+            sgs_processed = _process_subgroup_size(knl, self.subgroup_size)
+            # count madds as 2 ops
+            # (count all flops once and then count the madds again)
+
+            # flops counted w/subgroup granularity
+            float_ops = sgs_processed*(
+                stats_found[kso.OP_MAP].filter_by(
+                    dtype=[np.float32, np.float64],
+                    count_granularity=[lp.CountGranularity.SUBGROUP],
+                    ).sum() +
+                stats_found[kso.OP_MAP].filter_by(
+                    dtype=[np.float32, np.float64],
+                    count_granularity=[lp.CountGranularity.SUBGROUP],
+                    name=["madd"]
+                    ).sum())
+
+            # flops counted w/workitem granularity (should be zero)
+            float_ops += stats_found[kso.OP_MAP].filter_by(
+                dtype=[np.float32, np.float64],
+                count_granularity=[lp.CountGranularity.WORKITEM],
+                ).sum() + stats_found[kso.OP_MAP].filter_by(
+                dtype=[np.float32, np.float64],
+                count_granularity=[lp.CountGranularity.WORKITEM],
+                name=["madd"]
+                ).sum()
+            # TODO after ToCountMap gets version of sum that allows
+            # counting w/specified count granularity, update this
+
+            if not self.evaluate_polys:
+                float_ops = float_ops.eval_with_dict(param_dict)
+            stats_found[kso.FLOP_RATE] = float_ops/stats_found[kso.WALL_TIME]
+
+        if kso.MEM_BANDWIDTH in stat_options:
+            # first get footprint of data moved
+            from loopy import gather_access_footprint_bytes
+            from loopy.statistics import _process_subgroup_size
+            sgs_processed = _process_subgroup_size(knl, self.subgroup_size)
+            footsize_bytes = 0
+            for access, count in stats_found[kso.MEM_ACCESS_MAP].filter_by(
+                        mtype=["global"]).items():
+                direction = "write" if access.direction == "store" else "read"
+                footsize_bytes += gather_access_footprint_bytes(knl)[
+                        (access.variable, direction)].eval_with_dict(param_dict)
+
+            # mem access counted w/subgroup granularity
+            data_moved_bytes = stats_found[kso.MEM_ACCESS_MAP].filter_by(
+                mtype=["global"],
+                count_granularity=[lp.CountGranularity.SUBGROUP],
+                ).to_bytes().sum()*sgs_processed
+            # mem access counted w/workitem granularity
+            data_moved_bytes += stats_found[kso.MEM_ACCESS_MAP].filter_by(
+                mtype=["global"],
+                count_granularity=[lp.CountGranularity.WORKITEM],
+                ).to_bytes().sum()
+            # if these polys have not alread been evaluated, evaluate them
+            if not self.evaluate_polys:
+                data_moved_bytes = data_moved_bytes.eval_with_dict(param_dict)
+
+            # TODO decide on better way to handle multiple count granularities here
+            # (uniform access only counted once per warp)
+
+            stats_found[kso.MEM_BANDWIDTH] = (
+                data_moved_bytes/stats_found[kso.WALL_TIME],
+                footsize_bytes/stats_found[kso.WALL_TIME]
+                )
+
+        if kso.SAVE_PTX in stat_options:
+            ctx = self.get_cached_cl_context(knl)
+            if self.include_kernel_params_in_ptx_filename:
+                if param_dict is None:
+                    raise ValueError("Cannot include kernel params "
+                            "in ptx filename, no param dict passed.")
+                write_ptx(
+                    ctx,
+                    knl,
+                    filename="ptx_"+knl.name+"_"+"_".join(
+                        ["%s%d" % (p, v) for p, v in param_dict.items()]
+                        )+self.ptx_filename_suffix+".ptx"
+                    )
+            else:
+                write_ptx(ctx, knl)
+
+        if kso.GENERATED_CODE in stat_options:
+            stats_found[kso.GENERATED_CODE] = lp.generate_code_v2(
+                                              knl).device_code()
+
+        return stats_found
diff --git a/test/test_stat_collector.py b/test/test_stat_collector.py
new file mode 100644
index 0000000000000000000000000000000000000000..585a1bb4408ade0829bc9b34a7cf9be017e6e672
--- /dev/null
+++ b/test/test_stat_collector.py
@@ -0,0 +1,191 @@
+from __future__ import division, print_function
+
+__copyright__ = "Copyright (C) 2019 James Stevens"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import sys
+from pyopencl.tools import (  # noqa
+        pytest_generate_tests_for_pyopencl
+        as pytest_generate_tests)
+import loopy as lp
+import numpy as np
+from loopy.kernel_stat_collector.kernel_stat_collector import KernelStatCollector
+from loopy.kernel_stat_collector.kernel_stat_collector import KernelStatOptions as kso  # noqa
+from loopy.version import LOOPY_USE_LANGUAGE_VERSION_2018_2  # noqa
+from loopy.statistics import CountGranularity as CG
+
+
+def test_flop_rate_calculation(ctx_factory):
+
+    # Set platform and device
+    ctx = ctx_factory()
+    platform_name = ctx.devices[0].platform.name
+    device_name = ctx.devices[0].name
+
+    # Create a kernel
+    bsize = 16
+    knl = lp.make_kernel(
+            "{[i,k,j]: 0<=i<n and 0<=k<m and 0<=j<ell}",
+            [
+                "c[i, j] = sum(k, a[i, k]*b[k, j])"
+            ],
+            name="matmul", assumptions="n,m,ell >= 1")
+    knl = lp.add_and_infer_dtypes(knl, dict(a=np.float32, b=np.float32))
+    knl = lp.split_iname(knl, "i", bsize, outer_tag="g.0", inner_tag="l.1")
+    knl = lp.split_iname(knl, "j", bsize, outer_tag="g.1", inner_tag="l.0")
+    knl = lp.split_iname(knl, "k", bsize)
+    knl = lp.add_prefetch(knl, "a", ["k_inner", "i_inner"], default_tag="l.auto")
+    knl = lp.add_prefetch(knl, "b", ["j_inner", "k_inner"], default_tag="l.auto")
+
+    # Set parameter values
+    n = 512
+    m = 256
+    ell = 128
+    param_dict = {'n': n, 'm': m, 'ell': ell}
+
+    kp = KernelStatCollector(
+            platform_name=platform_name,
+            device_name=device_name,
+            evaluate_polys=False,
+            count_madds=False,  # TODO enable after madd counting branch is merged
+            )
+    stats = kp.collect_stats(
+            knl,
+            [
+                kso.WALL_TIME,
+                kso.OP_MAP,
+                kso.FLOP_RATE,
+            ],
+            param_dict=param_dict,
+            )
+
+    # Measured time + flop rate
+    time_measured = stats[kso.WALL_TIME]
+    flop_rate_measured = stats[kso.FLOP_RATE]
+
+    # Expected flop count and rate
+    flop_ct_expected = 2*m*n*ell
+    flop_rate_expected = flop_ct_expected/time_measured
+
+    # Make sure measured matches expected
+    assert flop_rate_measured == flop_rate_expected
+
+
+# TODO test accuracy of stat outputs:
+# (many of these tests may resemble those found in test_statistics.py)
+
+# kso.MEM_ACCESS_MAP
+
+# kso.OP_MAP
+
+# kso.SYNC_MAP
+
+# kso.GRID_SIZES
+
+# kso.MEM_BANDWIDTH  (both full and footprint only)
+
+
+# TODO test other stat collector functionality:
+
+# make sure evaluate_polys argument yields evaluated polynomials
+
+# If we call collect_stats() without passing a param dict, make sure we don't
+# get a ValueError unless it makes sense:
+# We should be able to collect op/mem/sync counts without a param_dict
+# unless eval_polys=True.
+# We should get a ValueError if we try to collect a stat that requires executing
+# the kernel without passing a param_dict (e.t., FLOP/s).
+
+def test1(ctx_factory):
+    ctx = ctx_factory()
+    platform_name = ctx.devices[0].platform.name
+    device_name = ctx.devices[0].name
+
+    lsize = 16
+    n = (lsize - 2) * 2 ** 4
+
+    dtype = np.float32
+
+    knl = lp.make_kernel(
+        "{[i,j]: 0<=i,j<n}",
+        [
+            "result[i,j] = u[i, j+1] + u[i+1, j] - 4*u[i+1, j+1] + u[i+1, j+2] + u[i+2, j+1]"
+        ],
+        name="fd",
+        assumptions="n >= 1 and n mod %d = 0" % (lsize - 2),
+        lang_version=(2018, 2),
+    )
+
+    knl = lp.add_and_infer_dtypes(knl, dict(u=dtype, result=dtype))
+    knl = lp.split_iname(knl, "i", lsize - 2, outer_tag="g.1", inner_tag="l.1")
+    knl = lp.split_iname(knl, "j", lsize - 2, outer_tag="g.0", inner_tag="l.0")
+    knl = lp.add_prefetch(knl, "u", ["i_inner", "j_inner"],
+                          fetch_bounding_box=True, default_tag=None)
+    knl = lp.tag_inames(knl, "u_dim_0:l.1, u_dim_1:l.0")
+
+    param_dict = {'n': n}
+
+    kp = KernelStatCollector(
+        platform_name=platform_name,
+        device_name=device_name,
+        evaluate_polys=False,
+        count_madds=False,  # TODO enable after madd counting branch is merged
+    )
+
+    stats = kp.collect_stats(
+        knl,
+        [
+            kso.WALL_TIME,
+            kso.OP_MAP,
+            kso.FLOP_RATE,
+            kso.GRID_SIZES,
+            kso.MEM_BANDWIDTH,
+            kso.MEM_ACCESS_MAP,
+            kso.OP_MAP,
+            kso.SYNC_MAP,
+        ],
+        param_dict=param_dict,
+    )
+    time_measured = stats[kso.WALL_TIME]
+    flop_rate_measured = stats[kso.FLOP_RATE]
+
+    # Expected flop count and rate
+    flop_ct_expected = n
+    flop_rate_expected = flop_ct_expected / time_measured
+    assert flop_rate_measured // flop_rate_expected == 1462
+    assert 100 * flop_rate_measured // stats[kso.MEM_BANDWIDTH][0] == 70
+
+    print(str(stats[kso.GRID_SIZES]))
+    assert str(stats[kso.GRID_SIZES]) == '[[PwQPolynomial("[n] -> { 1/14 * n }"), PwQPolynomial("[n] -> { 1/14 * n }")], [PwQPolynomial("[n] -> { 16 }"), PwQPolynomial("[n] -> { 16 }")]]'
+    assert len(stats[kso.SYNC_MAP]) == 2
+    assert stats[kso.SYNC_MAP]["kernel_launch"].eval_with_dict(param_dict) == 1
+    assert stats[kso.OP_MAP][lp.Op(np.float32, 'mul', CG.SUBGROUP)].eval_with_dict(param_dict) == 2048
+    assert len(stats[kso.MEM_ACCESS_MAP]) == 4
+
+
+if __name__ == "__main__":
+    if len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        from pytest import main
+        main([__file__])
+