__copyright__ = """ Copyright (C) 2013 Andreas Kloeckner Copyright (C) 2018 Matt Wala Copyright (C) 2018 Hao Gao """ __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 logging import sys import time import numpy as np import pytest from arraycontext import pytest_generate_tests_for_array_contexts from boxtree.array_context import ( PytestPyOpenCLArrayContextFactory, _acf, # noqa: F401 ) from boxtree.cost import ( FMMCostModel, _PythonFMMCostModel, make_pde_aware_translation_cost_model, ) logger = logging.getLogger(__name__) pytest_generate_tests = pytest_generate_tests_for_array_contexts([ PytestPyOpenCLArrayContextFactory, ]) # {{{ test_compare_cl_and_py_cost_model @pytest.mark.opencl @pytest.mark.parametrize( ("nsources", "ntargets", "dims", "dtype"), [ (50000, 50000, 3, np.float64) ] ) def test_compare_cl_and_py_cost_model(actx_factory, nsources, ntargets, dims, dtype): actx = actx_factory() # {{{ Generate sources, targets and target_radii from boxtree.tools import make_normal_particle_array as p_normal sources = p_normal(actx.queue, nsources, dims, dtype, seed=15) targets = p_normal(actx.queue, ntargets, dims, dtype, seed=18) rng = np.random.default_rng(22) target_radii = rng.uniform(0.0, 0.05, (ntargets,)).astype(dtype) # }}} # {{{ Generate tree and traversal from boxtree import TreeBuilder tb = TreeBuilder(actx.context) tree, _ = tb( actx.queue, sources, targets=targets, target_radii=target_radii, stick_out_factor=0.15, max_particles_in_box=30, debug=True ) from boxtree.traversal import FMMTraversalBuilder tg = FMMTraversalBuilder(actx.context, well_sep_is_n_away=2) trav_dev, _ = tg(actx.queue, tree, debug=True) trav = trav_dev.get(queue=actx.queue) # }}} # {{{ Construct cost models cl_cost_model = FMMCostModel(None) python_cost_model = _PythonFMMCostModel(None) constant_one_params = cl_cost_model.get_unit_calibration_params().copy() for ilevel in range(trav.tree.nlevels): constant_one_params["p_fmm_lev%d" % ilevel] = 10 xlat_cost = make_pde_aware_translation_cost_model(dims, trav.tree.nlevels) # }}} # {{{ Test process_form_multipoles from pymbolic import evaluate nlevels = trav.tree.nlevels p2m_cost = np.zeros(nlevels, dtype=np.float64) for ilevel in range(nlevels): p2m_cost[ilevel] = evaluate( xlat_cost.p2m(ilevel), context=constant_one_params ) p2m_cost_dev = actx.from_numpy(p2m_cost) actx.queue.finish() start_time = time.time() cl_form_multipoles = cl_cost_model.process_form_multipoles( actx.queue, trav_dev, p2m_cost_dev ) actx.queue.finish() logger.info("OpenCL time for process_form_multipoles: %gs", time.time() - start_time) start_time = time.time() python_form_multipoles = python_cost_model.process_form_multipoles( actx.queue, trav, p2m_cost ) logger.info("Python time for process_form_multipoles: %gs", time.time() - start_time) assert np.array_equal(actx.to_numpy(cl_form_multipoles), python_form_multipoles) # }}} # {{{ Test process_coarsen_multipoles m2m_cost = np.zeros(nlevels - 1, dtype=np.float64) for target_level in range(nlevels - 1): m2m_cost[target_level] = evaluate( xlat_cost.m2m(target_level + 1, target_level), context=constant_one_params ) m2m_cost_dev = actx.from_numpy(m2m_cost) actx.queue.finish() start_time = time.time() cl_coarsen_multipoles = cl_cost_model.process_coarsen_multipoles( actx.queue, trav_dev, m2m_cost_dev ) actx.queue.finish() logger.info("OpenCL time for coarsen_multipoles: %gs", time.time() - start_time) start_time = time.time() python_coarsen_multipoles = python_cost_model.process_coarsen_multipoles( actx.queue, trav, m2m_cost ) logger.info("Python time for coarsen_multipoles: %gs", time.time() - start_time) assert cl_coarsen_multipoles == python_coarsen_multipoles # }}} # {{{ Test process_direct actx.queue.finish() start_time = time.time() cl_ndirect_sources_per_target_box = \ cl_cost_model.get_ndirect_sources_per_target_box(actx.queue, trav_dev) cl_direct = cl_cost_model.process_direct( actx.queue, trav_dev, cl_ndirect_sources_per_target_box, 5.0 ) actx.queue.finish() logger.info("OpenCL time for process_direct: %gs", time.time() - start_time) start_time = time.time() python_ndirect_sources_per_target_box = \ python_cost_model.get_ndirect_sources_per_target_box(actx.queue, trav) python_direct = python_cost_model.process_direct( actx.queue, trav, python_ndirect_sources_per_target_box, 5.0 ) logger.info("Python time for process_direct: %gs", time.time() - start_time) assert np.array_equal(actx.to_numpy(cl_direct), python_direct) # }}} # {{{ Test aggregate_over_boxes start_time = time.time() cl_direct_aggregate = cl_cost_model.aggregate_over_boxes(cl_direct) actx.queue.finish() logger.info("OpenCL time for aggregate_over_boxes: %gs", time.time() - start_time) start_time = time.time() python_direct_aggregate = python_cost_model.aggregate_over_boxes(python_direct) logger.info("Python time for aggregate_over_boxes: %gs", time.time() - start_time) assert cl_direct_aggregate == python_direct_aggregate # }}} # {{{ Test process_list2 nlevels = trav.tree.nlevels m2l_cost = np.zeros(nlevels, dtype=np.float64) for ilevel in range(nlevels): m2l_cost[ilevel] = evaluate( xlat_cost.m2l(ilevel, ilevel), context=constant_one_params ) m2l_cost_dev = actx.from_numpy(m2l_cost) actx.queue.finish() start_time = time.time() cl_m2l_cost = cl_cost_model.process_list2(actx.queue, trav_dev, m2l_cost_dev) actx.queue.finish() logger.info("OpenCL time for process_list2: %gs", time.time() - start_time) start_time = time.time() python_m2l_cost = python_cost_model.process_list2(actx.queue, trav, m2l_cost) logger.info("Python time for process_list2: %gs", time.time() - start_time) assert np.array_equal(actx.to_numpy(cl_m2l_cost), python_m2l_cost) # }}} # {{{ Test process_list 3 m2p_cost = np.zeros(nlevels, dtype=np.float64) for ilevel in range(nlevels): m2p_cost[ilevel] = evaluate( xlat_cost.m2p(ilevel), context=constant_one_params ) m2p_cost_dev = actx.from_numpy(m2p_cost) actx.queue.finish() start_time = time.time() cl_m2p_cost = cl_cost_model.process_list3(actx.queue, trav_dev, m2p_cost_dev) actx.queue.finish() logger.info("OpenCL time for process_list3: %gs", time.time() - start_time) start_time = time.time() python_m2p_cost = python_cost_model.process_list3(actx.queue, trav, m2p_cost) logger.info("Python time for process_list3: %gs", time.time() - start_time) assert np.array_equal(actx.to_numpy(cl_m2p_cost), python_m2p_cost) # }}} # {{{ Test process_list4 p2l_cost = np.zeros(nlevels, dtype=np.float64) for ilevel in range(nlevels): p2l_cost[ilevel] = evaluate( xlat_cost.p2l(ilevel), context=constant_one_params ) p2l_cost_dev = actx.from_numpy(p2l_cost) actx.queue.finish() start_time = time.time() cl_p2l_cost = cl_cost_model.process_list4(actx.queue, trav_dev, p2l_cost_dev) actx.queue.finish() logger.info("OpenCL time for process_list4: %gs", time.time() - start_time) start_time = time.time() python_p2l_cost = python_cost_model.process_list4(actx.queue, trav, p2l_cost) logger.info("Python time for process_list4: %gs", time.time() - start_time) assert np.array_equal(actx.to_numpy(cl_p2l_cost), python_p2l_cost) # }}} # {{{ Test process_refine_locals l2l_cost = np.zeros(nlevels - 1, dtype=np.float64) for ilevel in range(nlevels - 1): l2l_cost[ilevel] = evaluate( xlat_cost.l2l(ilevel, ilevel + 1), context=constant_one_params ) l2l_cost_dev = actx.from_numpy(l2l_cost) actx.queue.finish() start_time = time.time() cl_refine_locals_cost = cl_cost_model.process_refine_locals( actx.queue, trav_dev, l2l_cost_dev ) actx.queue.finish() logger.info("OpenCL time for refine_locals: %gs", time.time() - start_time) start_time = time.time() python_refine_locals_cost = python_cost_model.process_refine_locals( actx.queue, trav, l2l_cost ) logger.info("Python time for refine_locals: %gs", time.time() - start_time) assert cl_refine_locals_cost == python_refine_locals_cost # }}} # {{{ Test process_eval_locals l2p_cost = np.zeros(nlevels, dtype=np.float64) for ilevel in range(nlevels): l2p_cost[ilevel] = evaluate( xlat_cost.l2p(ilevel), context=constant_one_params ) l2p_cost_dev = actx.from_numpy(l2p_cost) actx.queue.finish() start_time = time.time() cl_l2p_cost = cl_cost_model.process_eval_locals( actx.queue, trav_dev, l2p_cost_dev) actx.queue.finish() logger.info("OpenCL time for process_eval_locals: %gs", time.time() - start_time) start_time = time.time() python_l2p_cost = python_cost_model.process_eval_locals( actx.queue, trav, l2p_cost) logger.info("Python time for process_eval_locals: %gs", time.time() - start_time) assert np.array_equal(actx.to_numpy(cl_l2p_cost), python_l2p_cost) # }}} # }}} # {{{ test_estimate_calibration_params @pytest.mark.opencl def test_estimate_calibration_params(actx_factory): from boxtree.pyfmmlib_integration import ( FMMLibExpansionWrangler, FMMLibTreeIndependentDataForWrangler, Kernel, ) rng = np.random.default_rng(seed=42) nsources_list = [1000, 2000, 3000, 4000] ntargets_list = [1000, 2000, 3000, 4000] dims = 3 dtype = np.float64 actx = actx_factory() traversals = [] traversals_dev = [] level_to_orders = [] timing_results = [] def fmm_level_to_order(tree, ilevel): return 10 for nsources, ntargets in zip(nsources_list, ntargets_list): # {{{ Generate sources, targets and target_radii from boxtree.tools import make_normal_particle_array as p_normal sources = p_normal(actx.queue, nsources, dims, dtype, seed=15) targets = p_normal(actx.queue, ntargets, dims, dtype, seed=18) rng = np.random.default_rng(22) target_radii = rng.uniform(0.0, 0.05, (ntargets,)).astype(dtype) # }}} # {{{ Generate tree and traversal from boxtree import TreeBuilder tb = TreeBuilder(actx.context) tree, _ = tb( actx.queue, sources, targets=targets, target_radii=target_radii, stick_out_factor=0.15, max_particles_in_box=30, debug=True ) from boxtree.traversal import FMMTraversalBuilder tg = FMMTraversalBuilder(actx.context, well_sep_is_n_away=2) trav_dev, _ = tg(actx.queue, tree, debug=True) trav = trav_dev.get(queue=actx.queue) traversals.append(trav) traversals_dev.append(trav_dev) # }}} tree_indep = FMMLibTreeIndependentDataForWrangler( trav.tree.dimensions, Kernel.LAPLACE) wrangler = FMMLibExpansionWrangler(tree_indep, trav, fmm_level_to_order=fmm_level_to_order) level_to_orders.append(wrangler.level_orders) timing_data = {} from boxtree.fmm import drive_fmm src_weights = rng.random(size=tree.nsources, dtype=tree.coord_dtype) drive_fmm(wrangler, (src_weights,), timing_data=timing_data) timing_results.append(timing_data) time_field_name = "process_elapsed" def test_params_sanity(test_params): param_names = ["c_p2m", "c_m2m", "c_p2p", "c_m2l", "c_m2p", "c_p2l", "c_l2l", "c_l2p"] for name in param_names: assert isinstance(test_params[name], np.float64) def test_params_equal(test_params1, test_params2): param_names = ["c_p2m", "c_m2m", "c_p2p", "c_m2l", "c_m2p", "c_p2l", "c_l2l", "c_l2p"] for name in param_names: assert test_params1[name] == test_params2[name] python_cost_model = _PythonFMMCostModel(make_pde_aware_translation_cost_model) python_model_results = [] for icase in range(len(traversals)-1): traversal = traversals[icase] level_to_order = level_to_orders[icase] python_model_results.append(python_cost_model.cost_per_stage( actx.queue, traversal, level_to_order, _PythonFMMCostModel.get_unit_calibration_params(), )) python_params = python_cost_model.estimate_calibration_params( python_model_results, timing_results[:-1], time_field_name=time_field_name ) test_params_sanity(python_params) cl_cost_model = FMMCostModel(make_pde_aware_translation_cost_model) cl_model_results = [] for icase in range(len(traversals_dev)-1): traversal = traversals_dev[icase] level_to_order = level_to_orders[icase] cl_model_results.append(cl_cost_model.cost_per_stage( actx.queue, traversal, level_to_order, FMMCostModel.get_unit_calibration_params(), )) cl_params = cl_cost_model.estimate_calibration_params( cl_model_results, timing_results[:-1], time_field_name=time_field_name ) test_params_sanity(cl_params) test_params_equal(cl_params, python_params) # }}} # {{{ test_cost_model_op_counts_agree_with_constantone_wrangler class OpCountingTranslationCostModel: """A translation cost model which assigns at cost of 1 to each operation.""" def __init__(self, dim, nlevels): pass @staticmethod def direct(): return 1 @staticmethod def p2l(level): return 1 l2p = p2l p2m = p2l m2p = p2l @staticmethod def m2m(src_level, tgt_level): return 1 l2l = m2m m2l = m2m @pytest.mark.opencl @pytest.mark.parametrize( ("nsources", "ntargets", "dims", "dtype"), [ (5000, 5000, 3, np.float64) ] ) def test_cost_model_op_counts_agree_with_constantone_wrangler( actx_factory, nsources, ntargets, dims, dtype): actx = actx_factory() from boxtree.tools import make_normal_particle_array as p_normal sources = p_normal(actx.queue, nsources, dims, dtype, seed=16) targets = p_normal(actx.queue, ntargets, dims, dtype, seed=19) rng = np.random.default_rng(20) target_radii = rng.uniform(0, 0.04, (ntargets,)).astype(dtype) from boxtree import TreeBuilder tb = TreeBuilder(actx.context) tree, _ = tb( actx.queue, sources, targets=targets, target_radii=target_radii, stick_out_factor=0.15, max_particles_in_box=30, debug=True ) from boxtree.traversal import FMMTraversalBuilder tg = FMMTraversalBuilder(actx.context, well_sep_is_n_away=2) trav_dev, _ = tg(actx.queue, tree, debug=True) trav = trav_dev.get(queue=actx.queue) from boxtree.constant_one import ( ConstantOneExpansionWrangler, ConstantOneTreeIndependentDataForWrangler, ) tree_indep = ConstantOneTreeIndependentDataForWrangler() wrangler = ConstantOneExpansionWrangler(tree_indep, trav) timing_data = {} from boxtree.fmm import drive_fmm src_weights = rng.random(size=tree.nsources, dtype=tree.coord_dtype) drive_fmm(wrangler, (src_weights,), timing_data=timing_data) cost_model = FMMCostModel( translation_cost_model_factory=OpCountingTranslationCostModel ) level_to_order = np.array([1 for _ in range(tree.nlevels)]) modeled_time = cost_model.cost_per_stage( actx.queue, trav_dev, level_to_order, FMMCostModel.get_unit_calibration_params(), ) mismatches = [] for stage in timing_data: if timing_data[stage]["ops_elapsed"] != modeled_time[stage]: mismatches.append( (stage, timing_data[stage]["ops_elapsed"], modeled_time[stage])) assert not mismatches, "\n".join(str(s) for s in mismatches) # {{{ Test per-box cost total_cost = 0.0 for stage in timing_data: total_cost += timing_data[stage]["ops_elapsed"] per_box_cost = cost_model.cost_per_box( actx.queue, trav_dev, level_to_order, FMMCostModel.get_unit_calibration_params(), ) total_aggregate_cost = cost_model.aggregate_over_boxes(per_box_cost) assert total_cost == ( total_aggregate_cost + modeled_time["coarsen_multipoles"] + modeled_time["refine_locals"] ) # }}} # }}} # You can test individual routines by typing # $ python test_cost_model.py 'test_routine(_acf)' if __name__ == "__main__": if len(sys.argv) > 1: exec(sys.argv[1]) else: from pytest import main main([__file__]) # vim: foldmethod=marker