From 5b7a0fd626b0eb5e2653655e323d450d40334c50 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Fri, 1 Sep 2017 16:30:06 -0500 Subject: [PATCH 01/13] First pass at implementation/tests for MR with least squares --- leap/multistep/__init__.py | 70 ++++++++++++++++++++++++---- leap/multistep/multirate/__init__.py | 54 ++++++++++++++------- test/test_ab.py | 25 +++++++--- test/test_codegen_fortran.py | 23 +++++---- test/test_multirate.py | 23 +++++---- 5 files changed, 144 insertions(+), 51 deletions(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index fc73fe0..07d3363 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -74,6 +74,35 @@ class ABMonomialIntegrationFunctionFamily(ABIntegrationFunctionFamily): return 1/(func_idx+1) * x**(func_idx+1) +class ABTrigMonomialIntegrationFunctionFamily(ABIntegrationFunctionFamily): + # FIXME: Doesn't yet work for on-the-fly coefficients + + def __init__(self, order, alpha): + self.order = order + self.alpha = alpha + + def __len__(self): + return self.order + + def evaluate(self, func_idx, x): + func_idx += 1 + n = func_idx // 2 + if func_idx % 2 == 0: + return np.sin(self.alpha*n*x) + else: + return np.cos(self.alpha*n*x) + + def antiderivative(self, func_idx, x): + func_idx += 1 + n = func_idx // 2 + if func_idx == 1: + return x + elif func_idx % 2 == 0: + return -1/(n*self.alpha) * np.cos(self.alpha*n*x) + else: + return 1/(n*self.alpha) * np.sin(self.alpha*n*x) + + def _emit_func_family_operation(cb, name_gen, function_family, time_values, hist_vars, rhs_func): if isinstance(time_values, var): @@ -84,6 +113,7 @@ def _emit_func_family_operation(cb, name_gen, array = var("array") linear_solve = var("linear_solve") + lstsq = var("lstsq") # use: # Vandermonde^T * ab_coeffs = integrate(t_start, t_end, monomials) @@ -92,7 +122,8 @@ def _emit_func_family_operation(cb, name_gen, cb(vdmt, array(nfunctions*hist_len)) coeff_rhs = var(name_gen("coeff_rhs")) - cb(coeff_rhs, array(hist_len)) + + cb(coeff_rhs, array(nfunctions)) j = var(name_gen("vdm_j")) @@ -104,7 +135,11 @@ def _emit_func_family_operation(cb, name_gen, cb(coeff_rhs[i], rhs_func(i)) ab_coeffs = var(name_gen("ab_coeffs")) - cb(ab_coeffs, linear_solve(vdmt, coeff_rhs, nfunctions, 1)) + + if hist_len == nfunctions: + cb(ab_coeffs, linear_solve(vdmt, coeff_rhs, nfunctions, 1)) + else: + cb(ab_coeffs, lstsq(vdmt, coeff_rhs, hist_len, nfunctions, 1)) return _linear_comb( [ab_coeffs[ii] for ii in range(hist_len)], @@ -127,7 +162,10 @@ def _emit_func_family_operation(cb, name_gen, coeff_rhs[i] = rhs_func(i) - ab_coeffs = la.solve(vdm_t, coeff_rhs) + if hist_len == nfunctions: + ab_coeffs = la.solve(vdm_t, coeff_rhs) + else: + ab_coeffs = la.lstsq(vdm_t, coeff_rhs)[0] return _linear_comb(ab_coeffs, hist_vars) @@ -164,18 +202,32 @@ class AdamsBashforthMethod(Method): .. automethod:: generate """ - def __init__(self, component_id, order, state_filter_name=None, - hist_length=None, static_dt=False): + def __init__(self, component_id, function_family=None, state_filter_name=None, + hist_length=None, static_dt=False, order=None): """ + :arg function_family: Accepts an instance of + :class:`ABIntegrationFunctionFamily` + or an integer, in which case the classical monomial function family + with the order given by the integer is used. :arg static_dt: If *True*, changing the timestep during time integration is not allowed. """ + if function_family is not None and order is not None: + raise ValueError("may not specify both function_family and order") + + if function_family is None: + function_family = order + del order + + if isinstance(function_family, int): + function_family = ABMonomialIntegrationFunctionFamily(function_family) + super(AdamsBashforthMethod, self).__init__() - self.order = order + self.function_family = function_family if hist_length is None: - hist_length = order + hist_length = len(function_family) self.hist_length = hist_length self.static_dt = static_dt @@ -243,7 +295,7 @@ class AdamsBashforthMethod(Method): ab_sum = emit_ab_integration( cb_primary, name_gen, - ABMonomialIntegrationFunctionFamily(self.order), + self.function_family, time_hist, history, 0, t_end) @@ -319,7 +371,7 @@ class AdamsBashforthMethod(Method): cb(self.time_history[i], self.t) from leap.rk import ORDER_TO_RK_METHOD - rk_method = ORDER_TO_RK_METHOD[self.order] + rk_method = ORDER_TO_RK_METHOD[self.function_family.order] rk_tableau = tuple(zip(rk_method.c, rk_method.a_explicit)) rk_coeffs = rk_method.output_coeffs diff --git a/leap/multistep/multirate/__init__.py b/leap/multistep/multirate/__init__.py index cbd79da..4539fe3 100644 --- a/leap/multistep/multirate/__init__.py +++ b/leap/multistep/multirate/__init__.py @@ -66,7 +66,8 @@ class MultiRateHistory(Record): .. automethod:: __init__ """ def __init__(self, interval, func_name, arguments, order=None, - rhs_policy=rhs_policy.late, invalidate_computed_state=False): + rhs_policy=rhs_policy.late, invalidate_computed_state=False, + hist_length=None): """ :arg interval: An integer indicating the interval (relative to the smallest available timestep) at which this history is to be @@ -81,6 +82,9 @@ class MultiRateHistory(Record): :arg invalidate_dependent_state: Whether evaluating this right-hand side should force a recomputation of any state that depended upon now-superseded state. + :arg hist_length: history length. If greater than order, we use a + least-squares solve rather than a linear solve to obtain the AB + coefficients for this history """ super(MultiRateHistory, self).__init__( interval=interval, @@ -88,11 +92,12 @@ class MultiRateHistory(Record): arguments=arguments, order=order, rhs_policy=rhs_policy, - invalidate_computed_state=invalidate_computed_state) + invalidate_computed_state=invalidate_computed_state, + hist_length=hist_length) @property def history_length(self): - return self.order + return self.hist_length class RHS(MultiRateHistory): @@ -303,17 +308,20 @@ class MultiRateMultiStepMethod(Method): # }}} - # {{{ plug default order into rhss + # {{{ plug default order, history length into rhss new_rhss = [] for component_rhss in rhss: new_component_rhss = [] for rhs in component_rhss: order = rhs.order + hist_length = rhs.hist_length if order is None: order = default_order + if hist_length is None: + hist_length = order - new_component_rhss.append(rhs.copy(order=order)) + new_component_rhss.append(rhs.copy(order=order, hist_length=hist_length)) new_rhss.append(tuple(new_component_rhss)) @@ -336,6 +344,10 @@ class MultiRateMultiStepMethod(Method): for component_rhss in self.rhss for rhs in component_rhss) + self.max_hist_length = max(rhs.hist_length + for component_rhss in self.rhss + for rhs in component_rhss) + # {{{ process intervals intervals = sorted(rhs.interval @@ -587,7 +599,7 @@ class MultiRateMultiStepMethod(Method): """Initialize the stepper with an RK method. Return the code that computes the startup history.""" - bootstrap_steps = self.max_order - 1 + bootstrap_steps = self.max_hist_length - 1 final_iglobal_substep = bootstrap_steps * self.nsubsteps @@ -682,12 +694,12 @@ class MultiRateMultiStepMethod(Method): for irhs, rhs in enumerate(component_rhss): if (substeps_from_start % rhs.interval == 0 and (substeps_from_start // rhs.interval - < rhs.order)): + < rhs.history_length)): intervals_from_start = ( substeps_from_start // rhs.interval) - i = rhs.order - 1 - intervals_from_start + i = rhs.history_length - 1 - intervals_from_start assert i >= 0 with cb.if_(self.bootstrap_step, "==", test_step): @@ -740,7 +752,7 @@ class MultiRateMultiStepMethod(Method): key = comp_name, irhs temp_hist_substeps[key] = list(range( - -rhs.interval*(rhs.order-1), 1, rhs.interval)) + -rhs.interval*(rhs.history_length-1), 1, rhs.interval)) if self.static_dt: temp_time_vars[key] = list( @@ -1026,13 +1038,13 @@ class MultiRateMultiStepMethod(Method): if not self.static_dt: for time_var, time_expr in zip( self.time_vars[key], - temp_time_vars[comp_name, irhs][-rhs.order:]): + temp_time_vars[comp_name, irhs][-rhs.history_length:]): cb(time_var, time_expr) cb.fence() for hist_var, hist_expr in zip( self.history_vars[key], - temp_hist_vars[comp_name, irhs][-rhs.order:]): + temp_hist_vars[comp_name, irhs][-rhs.history_length:]): cb(hist_var, hist_expr) cb.fence() @@ -1126,7 +1138,9 @@ class TwoRateAdamsBashforthMethod(MultiRateMultiStepMethod): def __init__(self, method, order, step_ratio, slow_state_filter_name=None, fast_state_filter_name=None, - static_dt=False): + static_dt=False, + hist_length_slow=None, + hist_length_fast=None): from warnings import warn warn("TwoRateAdamsBashforthMethod is a compatibility shim that should no " "longer be used. Use the fully general " @@ -1158,21 +1172,29 @@ class TwoRateAdamsBashforthMethod(MultiRateMultiStepMethod): else: s2f_policy = rhs_policy.late + if hist_length_slow is None: + hist_length_slow = order + + if hist_length_fast is None: + hist_length_slow = order + super(TwoRateAdamsBashforthMethod, self).__init__( order, ( ( "dt", "fast", "=", - MultiRateHistory(1, "f2f", ("fast", "slow",)), + MultiRateHistory(1, "f2f", ("fast", "slow",), + hist_length=hist_length_fast), MultiRateHistory(s2f_interval, "s2f", - ("fast", "slow",), rhs_policy=s2f_policy), + ("fast", "slow",), rhs_policy=s2f_policy, + hist_length=hist_length_fast), ), ( "dt", "slow", "=", MultiRateHistory(step_ratio, "f2s", ("fast", "slow",), - rhs_policy=f2s_policy), + rhs_policy=f2s_policy, hist_length=hist_length_slow), MultiRateHistory(step_ratio, "s2s", ("fast", "slow",), - rhs_policy=s2s_policy), + rhs_policy=s2s_policy, hist_length=hist_length_slow), ),), state_filter_names={ diff --git a/test/test_ab.py b/test/test_ab.py index 7da554e..9e284d4 100755 --- a/test/test_ab.py +++ b/test/test_ab.py @@ -27,21 +27,32 @@ THE SOFTWARE. import sys import pytest from leap.multistep import AdamsBashforthMethod +import leap.multistep as multistep from utils import ( # noqa python_method_impl_interpreter as pmi_int, python_method_impl_codegen as pmi_cg) -@pytest.mark.parametrize("order", [1, 2, 3, 4, 5]) -@pytest.mark.parametrize("static_dt", [True, False]) -def test_ab_accuracy(python_method_impl, order, static_dt, show_dag=False, - plot_solution=False): +@pytest.mark.parametrize(("method", "expected_order"), [ + (AdamsBashforthMethod("y", order, static_dt=static_dt), order) + for order in [1, 3, 5] + for static_dt in [True, False] + ] + [ + (AdamsBashforthMethod("y", + multistep.ABTrigMonomialIntegrationFunctionFamily(4, 0.1), + static_dt=True), 3) + ] + [ + (AdamsBashforthMethod("y", order, hist_length=order+1, + static_dt=static_dt), order) + for order in [1, 3, 5] + for static_dt in [True, False] + ]) +def test_ab_accuracy(python_method_impl, method, expected_order, + show_dag=False, plot_solution=False): from utils import check_simple_convergence - method = AdamsBashforthMethod("y", order, static_dt=static_dt) - #method = AdamsBashforthMethod("y", order, hist_length=order+1) check_simple_convergence(method=method, method_impl=python_method_impl, - expected_order=order, show_dag=show_dag, + expected_order=expected_order, show_dag=show_dag, plot_solution=plot_solution) diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index d04a517..c8294be 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -215,15 +215,18 @@ def test_rk_codegen_fancy(): # }}} -@pytest.mark.parametrize("min_order", [2, 3, 4, 5]) +@pytest.mark.parametrize(("min_order", "hist_length"), [(2,2), (2,3), (3,3), + (3,4), (4,4), (4,5), (5,5), (5,6)]) @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) -def test_multirate_codegen(min_order, method_name): +def test_multirate_codegen(min_order, hist_length, method_name): from leap.multistep.multirate import TwoRateAdamsBashforthMethod stepper = TwoRateAdamsBashforthMethod( method_name, min_order, 4, slow_state_filter_name="slow_filt", - fast_state_filter_name="fast_filt") + fast_state_filter_name="fast_filt", + hist_length_slow=hist_length, + hist_length_fast=hist_length) code = stepper.generate() @@ -403,14 +406,15 @@ def test_adaptive_rk_codegen_error(): ]) -@pytest.mark.parametrize("min_order", [2, 3, 4, 5]) -def test_singlerate_squarewave(min_order): +@pytest.mark.parametrize(("min_order", "hist_length"), + [(2,2), (2,3), (3,3), (3,4), (4,4), (4,5), (5,5), (5,6),]) +def test_singlerate_squarewave(min_order, hist_length): from leap.multistep import AdamsBashforthMethod component_id = 'y' rhs_function = 'y' - stepper = AdamsBashforthMethod("y", min_order) + stepper = AdamsBashforthMethod("y", min_order, hist_length=hist_length) from dagrt.function_registry import ( base_function_registry, register_ode_rhs) @@ -448,9 +452,10 @@ def test_singlerate_squarewave(min_order): @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) -@pytest.mark.parametrize("min_order", [4, 3, 2]) -def test_multirate_squarewave(min_order, method_name): - stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4) +@pytest.mark.parametrize(("min_order","hist_length"), [(4,4), (4,5), (3,3), (3,4), (2,2), (2,3)]) +def test_multirate_squarewave(min_order, hist_length, method_name): + stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4, + hist_length_slow=hist_length, hist_length_fast=hist_length) code = stepper.generate() diff --git a/test/test_multirate.py b/test/test_multirate.py index 07af8e7..383f970 100644 --- a/test/test_multirate.py +++ b/test/test_multirate.py @@ -46,10 +46,11 @@ from utils import ( # noqa class MultirateTimestepperAccuracyChecker(object): """Check that the multirate timestepper has the advertised accuracy.""" - def __init__(self, method, order, step_ratio, static_dt, ode, method_impl, + def __init__(self, method, order, hist_length, step_ratio, static_dt, ode, method_impl, display_dag=False, display_solution=False): self.method = method self.order = order + self.hist_length = hist_length self.step_ratio = step_ratio self.static_dt = static_dt self.ode = ode @@ -61,7 +62,9 @@ class MultirateTimestepperAccuracyChecker(object): def get_code(self): method = TwoRateAdamsBashforthMethod( self.method, self.order, self.step_ratio, - static_dt=self.static_dt) + static_dt=self.static_dt, + hist_length_slow = self.hist_length, + hist_length_fast = self.hist_length) return method.generate() @@ -166,7 +169,7 @@ class MultirateTimestepperAccuracyChecker(object): @pytest.mark.slowtest -@pytest.mark.parametrize("order", [1, 3]) +@pytest.mark.parametrize(("order", "hist_length"), [(1,1),(1,2),(3,3),(3,4)]) @pytest.mark.parametrize("system", [ #"Basic", "Full", @@ -175,7 +178,7 @@ class MultirateTimestepperAccuracyChecker(object): ]) @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) @pytest.mark.parametrize("static_dt", [True, False]) -def test_multirate_accuracy(method_name, order, system, static_dt, step_ratio=2): +def test_multirate_accuracy(method_name, order, hist_length, system, static_dt, step_ratio=2): """Check that the multirate timestepper has the advertised accuracy""" import multirate_test_systems @@ -187,12 +190,12 @@ def test_multirate_accuracy(method_name, order, system, static_dt, step_ratio=2) print("------------------------------------------------------") MultirateTimestepperAccuracyChecker( - method_name, order, step_ratio, static_dt=static_dt, + method_name, order, hist_length, step_ratio, static_dt=static_dt, ode=system(), method_impl=pmi_cg)() -def test_single_rate_identical(order=3): +def test_single_rate_identical(order=3, hist_length=4): from leap.multistep import AdamsBashforthMethod from dagrt.exec_numpy import NumpyInterpreter @@ -204,7 +207,7 @@ def test_single_rate_identical(order=3): # {{{ single rate - single_rate_method = AdamsBashforthMethod("y", order=order) + single_rate_method = AdamsBashforthMethod("y", order=order, hist_length=hist_length) single_rate_code = single_rate_method.generate() def single_rate_rhs(t, y): @@ -244,12 +247,12 @@ def test_single_rate_identical(order=3): ( ( 'dt', 'fast', '=', - MRHistory(1, "f", ("fast", "slow",)), + MRHistory(1, "f", ("fast", "slow",), hist_length=hist_length), ), ( 'dt', 'slow', '=', MRHistory(1, "s", ("fast", "slow",), - rhs_policy=rhs_policy.late), + rhs_policy=rhs_policy.late, hist_length=hist_length), ),) ) multi_rate_code = multi_rate_method.generate() @@ -277,7 +280,7 @@ def test_single_rate_identical(order=3): idx = {"fast": 0, "slow": 1}[event.component_id] if event.t not in multi_rate_values: multi_rate_values[event.t] = [None, None] - + multi_rate_values[event.t][idx] = event.state_component if len(multi_rate_values) > nsteps: -- GitLab From 4e6bc6b2e6cfca03f34d1113565c065cb5804665 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Fri, 1 Sep 2017 16:38:59 -0500 Subject: [PATCH 02/13] Forgot to modify requirements for modded dagrt --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index f8db6ac..b4ec38b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ git+git://github.com/inducer/pytools git+git://github.com/inducer/pymbolic -git+https://gitlab.tiker.net/inducer/dagrt.git +git+https://gitlab.tiker.net/inducer/dagrt.git@cory-lls -- GitLab From db24f5712ad1911d20334c5630f8f331ccc5af95 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Fri, 1 Sep 2017 16:59:12 -0500 Subject: [PATCH 03/13] Flake8 fixes --- leap/multistep/__init__.py | 2 +- leap/multistep/multirate/__init__.py | 14 ++++++++------ test/test_codegen_fortran.py | 11 ++++++----- test/test_multirate.py | 24 ++++++++++++++---------- 4 files changed, 29 insertions(+), 22 deletions(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index 07d3363..8c6e74c 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -162,7 +162,7 @@ def _emit_func_family_operation(cb, name_gen, coeff_rhs[i] = rhs_func(i) - if hist_len == nfunctions: + if hist_len == nfunctions: ab_coeffs = la.solve(vdm_t, coeff_rhs) else: ab_coeffs = la.lstsq(vdm_t, coeff_rhs)[0] diff --git a/leap/multistep/multirate/__init__.py b/leap/multistep/multirate/__init__.py index 4539fe3..855bf57 100644 --- a/leap/multistep/multirate/__init__.py +++ b/leap/multistep/multirate/__init__.py @@ -82,7 +82,7 @@ class MultiRateHistory(Record): :arg invalidate_dependent_state: Whether evaluating this right-hand side should force a recomputation of any state that depended upon now-superseded state. - :arg hist_length: history length. If greater than order, we use a + :arg hist_length: history length. If greater than order, we use a least-squares solve rather than a linear solve to obtain the AB coefficients for this history """ @@ -321,7 +321,8 @@ class MultiRateMultiStepMethod(Method): if hist_length is None: hist_length = order - new_component_rhss.append(rhs.copy(order=order, hist_length=hist_length)) + new_component_rhss.append(rhs.copy(order=order, + hist_length=hist_length)) new_rhss.append(tuple(new_component_rhss)) @@ -1038,7 +1039,8 @@ class MultiRateMultiStepMethod(Method): if not self.static_dt: for time_var, time_expr in zip( self.time_vars[key], - temp_time_vars[comp_name, irhs][-rhs.history_length:]): + temp_time_vars[comp_name, + irhs][-rhs.history_length:]): cb(time_var, time_expr) cb.fence() @@ -1139,7 +1141,7 @@ class TwoRateAdamsBashforthMethod(MultiRateMultiStepMethod): slow_state_filter_name=None, fast_state_filter_name=None, static_dt=False, - hist_length_slow=None, + hist_length_slow=None, hist_length_fast=None): from warnings import warn warn("TwoRateAdamsBashforthMethod is a compatibility shim that should no " @@ -1183,10 +1185,10 @@ class TwoRateAdamsBashforthMethod(MultiRateMultiStepMethod): ( ( "dt", "fast", "=", - MultiRateHistory(1, "f2f", ("fast", "slow",), + MultiRateHistory(1, "f2f", ("fast", "slow",), hist_length=hist_length_fast), MultiRateHistory(s2f_interval, "s2f", - ("fast", "slow",), rhs_policy=s2f_policy, + ("fast", "slow",), rhs_policy=s2f_policy, hist_length=hist_length_fast), ), ( diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index c8294be..9e8caa0 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -215,8 +215,8 @@ def test_rk_codegen_fancy(): # }}} -@pytest.mark.parametrize(("min_order", "hist_length"), [(2,2), (2,3), (3,3), - (3,4), (4,4), (4,5), (5,5), (5,6)]) +@pytest.mark.parametrize(("min_order", "hist_length"), [(2, 2), (2, 3), (3, 3), + (3, 4), (4, 4), (4, 5), (5, 5), (5, 6)]) @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) def test_multirate_codegen(min_order, hist_length, method_name): from leap.multistep.multirate import TwoRateAdamsBashforthMethod @@ -406,8 +406,8 @@ def test_adaptive_rk_codegen_error(): ]) -@pytest.mark.parametrize(("min_order", "hist_length"), - [(2,2), (2,3), (3,3), (3,4), (4,4), (4,5), (5,5), (5,6),]) +@pytest.mark.parametrize(("min_order", "hist_length"), + [(2, 2), (2, 3), (3, 3), (3, 4), (4, 4), (4, 5), (5, 5), (5, 6), ]) def test_singlerate_squarewave(min_order, hist_length): from leap.multistep import AdamsBashforthMethod @@ -452,7 +452,8 @@ def test_singlerate_squarewave(min_order, hist_length): @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) -@pytest.mark.parametrize(("min_order","hist_length"), [(4,4), (4,5), (3,3), (3,4), (2,2), (2,3)]) +@pytest.mark.parametrize(("min_order", "hist_length"), [(4, 4), (4, 5), + (3, 3), (3, 4), (2, 2), (2, 3), ]) def test_multirate_squarewave(min_order, hist_length, method_name): stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4, hist_length_slow=hist_length, hist_length_fast=hist_length) diff --git a/test/test_multirate.py b/test/test_multirate.py index 383f970..c5e1d17 100644 --- a/test/test_multirate.py +++ b/test/test_multirate.py @@ -46,8 +46,8 @@ from utils import ( # noqa class MultirateTimestepperAccuracyChecker(object): """Check that the multirate timestepper has the advertised accuracy.""" - def __init__(self, method, order, hist_length, step_ratio, static_dt, ode, method_impl, - display_dag=False, display_solution=False): + def __init__(self, method, order, hist_length, step_ratio, static_dt, ode, + method_impl, display_dag=False, display_solution=False): self.method = method self.order = order self.hist_length = hist_length @@ -62,9 +62,9 @@ class MultirateTimestepperAccuracyChecker(object): def get_code(self): method = TwoRateAdamsBashforthMethod( self.method, self.order, self.step_ratio, - static_dt=self.static_dt, - hist_length_slow = self.hist_length, - hist_length_fast = self.hist_length) + static_dt=self.static_dt, + hist_length_slow=self.hist_length, + hist_length_fast=self.hist_length) return method.generate() @@ -169,7 +169,8 @@ class MultirateTimestepperAccuracyChecker(object): @pytest.mark.slowtest -@pytest.mark.parametrize(("order", "hist_length"), [(1,1),(1,2),(3,3),(3,4)]) +@pytest.mark.parametrize(("order", "hist_length"), [(1, 1), (1, 2), + (3, 3), (3, 4), ]) @pytest.mark.parametrize("system", [ #"Basic", "Full", @@ -178,7 +179,8 @@ class MultirateTimestepperAccuracyChecker(object): ]) @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) @pytest.mark.parametrize("static_dt", [True, False]) -def test_multirate_accuracy(method_name, order, hist_length, system, static_dt, step_ratio=2): +def test_multirate_accuracy(method_name, order, hist_length, + system, static_dt, step_ratio=2): """Check that the multirate timestepper has the advertised accuracy""" import multirate_test_systems @@ -207,7 +209,8 @@ def test_single_rate_identical(order=3, hist_length=4): # {{{ single rate - single_rate_method = AdamsBashforthMethod("y", order=order, hist_length=hist_length) + single_rate_method = AdamsBashforthMethod("y", order=order, + hist_length=hist_length) single_rate_code = single_rate_method.generate() def single_rate_rhs(t, y): @@ -247,7 +250,8 @@ def test_single_rate_identical(order=3, hist_length=4): ( ( 'dt', 'fast', '=', - MRHistory(1, "f", ("fast", "slow",), hist_length=hist_length), + MRHistory(1, "f", ("fast", "slow",), + hist_length=hist_length), ), ( 'dt', 'slow', '=', @@ -280,7 +284,7 @@ def test_single_rate_identical(order=3, hist_length=4): idx = {"fast": 0, "slow": 1}[event.component_id] if event.t not in multi_rate_values: multi_rate_values[event.t] = [None, None] - + multi_rate_values[event.t][idx] = event.state_component if len(multi_rate_values) > nsteps: -- GitLab From 0cef0bb5f009c4a2d369cc484a713a5b7ebea95d Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 19 Sep 2017 10:37:59 -0500 Subject: [PATCH 04/13] Towards SVD builtin in Fortran for least squares solve --- leap/multistep/__init__.py | 45 ++++++++++++++++++++++++++++++++++---- 1 file changed, 41 insertions(+), 4 deletions(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index 8c6e74c..6ed4405 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -113,7 +113,9 @@ def _emit_func_family_operation(cb, name_gen, array = var("array") linear_solve = var("linear_solve") - lstsq = var("lstsq") + svd = var("svd") + matmul = var("matmul") + transpose = var("transpose") # use: # Vandermonde^T * ab_coeffs = integrate(t_start, t_end, monomials) @@ -122,7 +124,6 @@ def _emit_func_family_operation(cb, name_gen, cb(vdmt, array(nfunctions*hist_len)) coeff_rhs = var(name_gen("coeff_rhs")) - cb(coeff_rhs, array(nfunctions)) j = var(name_gen("vdm_j")) @@ -139,7 +140,40 @@ def _emit_func_family_operation(cb, name_gen, if hist_len == nfunctions: cb(ab_coeffs, linear_solve(vdmt, coeff_rhs, nfunctions, 1)) else: - cb(ab_coeffs, lstsq(vdmt, coeff_rhs, hist_len, nfunctions, 1)) + # Least squares with SVD builtin + u = var(name_gen("u")) + ut = var(name_gen("ut")) + intermed = var(name_gen("intermed")) + Ainv = var(name_gen("Ainv")) + sigma = var(name_gen("sigma")) + sig_array = var(name_gen("sig_array")) + v = var(name_gen("v")) + vt = var(name_gen("vt")) + history_length = var(name_gen("history_length")) + + cb(history_length, hist_len) + + cb(Ainv, array(nfunctions*hist_len)) + cb(intermed, array(nfunctions*hist_len)) + + cb("u, sigma, vt", "`svd`(vdm_transpose, history_length)") + cb(ut, transpose(u, nfunctions)) + cb(v, transpose(vt, hist_len)) + + # Make singular value array + cb(sig_array, array(nfunctions*nfunctions)) + + for j in range(len(function_family)*len(function_family)): + cb(sig_array[j], 0) + + cb.fence() + + for i in range(len(function_family)): + cb(sig_array[i*(nfunctions+1)], sigma[i]**-1) + + cb(intermed, matmul(v, sig_array, nfunctions, nfunctions)) + cb(Ainv, matmul(intermed, ut, nfunctions, nfunctions)) + cb(ab_coeffs, matmul(Ainv, coeff_rhs, nfunctions, 1)) return _linear_comb( [ab_coeffs[ii] for ii in range(hist_len)], @@ -165,7 +199,10 @@ def _emit_func_family_operation(cb, name_gen, if hist_len == nfunctions: ab_coeffs = la.solve(vdm_t, coeff_rhs) else: - ab_coeffs = la.lstsq(vdm_t, coeff_rhs)[0] + # SVD-based least squares solve + u, sigma, v = la.svd(vdm_t, full_matrices=False) + Ainv = np.dot(v.transpose(), np.dot(la.inv(np.diag(sigma)),u.transpose())) + ab_coeffs = np.dot(Ainv, coeff_rhs) return _linear_comb(ab_coeffs, hist_vars) -- GitLab From 36477ee7cd6fc34351d84b9d9726364a3f4579de Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 19 Sep 2017 10:43:28 -0500 Subject: [PATCH 05/13] Flake8 fixes, removing errant debugging variables/attempts --- leap/multistep/__init__.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index 6ed4405..7539791 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -140,7 +140,7 @@ def _emit_func_family_operation(cb, name_gen, if hist_len == nfunctions: cb(ab_coeffs, linear_solve(vdmt, coeff_rhs, nfunctions, 1)) else: - # Least squares with SVD builtin + # Least squares with SVD builtin u = var(name_gen("u")) ut = var(name_gen("ut")) intermed = var(name_gen("intermed")) @@ -149,14 +149,11 @@ def _emit_func_family_operation(cb, name_gen, sig_array = var(name_gen("sig_array")) v = var(name_gen("v")) vt = var(name_gen("vt")) - history_length = var(name_gen("history_length")) - - cb(history_length, hist_len) cb(Ainv, array(nfunctions*hist_len)) cb(intermed, array(nfunctions*hist_len)) - cb("u, sigma, vt", "`svd`(vdm_transpose, history_length)") + cb("u, sigma, vt", svd(vdmt, hist_len)) cb(ut, transpose(u, nfunctions)) cb(v, transpose(vt, hist_len)) @@ -164,12 +161,12 @@ def _emit_func_family_operation(cb, name_gen, cb(sig_array, array(nfunctions*nfunctions)) for j in range(len(function_family)*len(function_family)): - cb(sig_array[j], 0) - + cb(sig_array[j], 0) + cb.fence() for i in range(len(function_family)): - cb(sig_array[i*(nfunctions+1)], sigma[i]**-1) + cb(sig_array[i*(nfunctions+1)], sigma[i]**-1) cb(intermed, matmul(v, sig_array, nfunctions, nfunctions)) cb(Ainv, matmul(intermed, ut, nfunctions, nfunctions)) @@ -201,7 +198,8 @@ def _emit_func_family_operation(cb, name_gen, else: # SVD-based least squares solve u, sigma, v = la.svd(vdm_t, full_matrices=False) - Ainv = np.dot(v.transpose(), np.dot(la.inv(np.diag(sigma)),u.transpose())) + Ainv = np.dot(v.transpose(), np.dot(la.inv(np.diag(sigma)), + u.transpose())) ab_coeffs = np.dot(Ainv, coeff_rhs) return _linear_comb(ab_coeffs, hist_vars) -- GitLab From 3ff7e21ec84ec333e0570e3d6eabd3829a3a2cfe Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 21 Sep 2017 10:04:46 -0500 Subject: [PATCH 06/13] Finishing SVD implementation of least-squares for AB3-4, etc. --- leap/multistep/__init__.py | 3 ++- test/test_codegen_fortran.py | 5 +++++ 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index 7539791..c618185 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -153,7 +153,8 @@ def _emit_func_family_operation(cb, name_gen, cb(Ainv, array(nfunctions*hist_len)) cb(intermed, array(nfunctions*hist_len)) - cb("u, sigma, vt", svd(vdmt, hist_len)) + cb((u, sigma, vt), svd(vdmt, hist_len)) + cb(ut, transpose(u, nfunctions)) cb(v, transpose(vt, hist_len)) diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index 9e8caa0..706d2d4 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -317,6 +317,9 @@ def test_multirate_codegen(min_order, hist_length, method_name): else: fac = 5 + if (hist_length != min_order): + fac = 12 + num_trips_one = 10*fac num_trips_two = 100*fac @@ -442,6 +445,8 @@ def test_singlerate_squarewave(min_order, hist_length): """) code_str = codegen(code) + with open("abmethod_test.f90", "wt") as outf: + outf.write(code_str) run_fortran([ ("abmethod.f90", code_str), -- GitLab From 1abdc38356b67c946da5ba3ad1353798284a671c Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 21 Sep 2017 10:41:25 -0500 Subject: [PATCH 07/13] Fixes some Flake8 issue --- leap/multistep/__init__.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index c618185..ffe2d6f 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -144,13 +144,13 @@ def _emit_func_family_operation(cb, name_gen, u = var(name_gen("u")) ut = var(name_gen("ut")) intermed = var(name_gen("intermed")) - Ainv = var(name_gen("Ainv")) + ainv = var(name_gen("ainv")) sigma = var(name_gen("sigma")) sig_array = var(name_gen("sig_array")) v = var(name_gen("v")) vt = var(name_gen("vt")) - cb(Ainv, array(nfunctions*hist_len)) + cb(ainv, array(nfunctions*hist_len)) cb(intermed, array(nfunctions*hist_len)) cb((u, sigma, vt), svd(vdmt, hist_len)) @@ -170,8 +170,8 @@ def _emit_func_family_operation(cb, name_gen, cb(sig_array[i*(nfunctions+1)], sigma[i]**-1) cb(intermed, matmul(v, sig_array, nfunctions, nfunctions)) - cb(Ainv, matmul(intermed, ut, nfunctions, nfunctions)) - cb(ab_coeffs, matmul(Ainv, coeff_rhs, nfunctions, 1)) + cb(ainv, matmul(intermed, ut, nfunctions, nfunctions)) + cb(ab_coeffs, matmul(ainv, coeff_rhs, nfunctions, 1)) return _linear_comb( [ab_coeffs[ii] for ii in range(hist_len)], @@ -199,9 +199,9 @@ def _emit_func_family_operation(cb, name_gen, else: # SVD-based least squares solve u, sigma, v = la.svd(vdm_t, full_matrices=False) - Ainv = np.dot(v.transpose(), np.dot(la.inv(np.diag(sigma)), + ainv = np.dot(v.transpose(), np.dot(la.inv(np.diag(sigma)), u.transpose())) - ab_coeffs = np.dot(Ainv, coeff_rhs) + ab_coeffs = np.dot(ainv, coeff_rhs) return _linear_comb(ab_coeffs, hist_vars) -- GitLab From 5643f0b8890969571b539630f42d20bba8dc8a00 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Thu, 21 Sep 2017 11:29:52 -0500 Subject: [PATCH 08/13] Placate Flake8 yet again --- test/test_ab.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_ab.py b/test/test_ab.py index 9e284d4..d881284 100755 --- a/test/test_ab.py +++ b/test/test_ab.py @@ -43,7 +43,7 @@ from utils import ( # noqa multistep.ABTrigMonomialIntegrationFunctionFamily(4, 0.1), static_dt=True), 3) ] + [ - (AdamsBashforthMethod("y", order, hist_length=order+1, + (AdamsBashforthMethod("y", order, hist_length=order+1, static_dt=static_dt), order) for order in [1, 3, 5] for static_dt in [True, False] -- GitLab From fea45ce933749098097d0da428fc299031afac6f Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 30 Jan 2018 21:39:44 -0600 Subject: [PATCH 09/13] Remove trig stuff from this branch (whoops) --- leap/multistep/__init__.py | 29 ----------------------------- 1 file changed, 29 deletions(-) diff --git a/leap/multistep/__init__.py b/leap/multistep/__init__.py index ffe2d6f..0787df3 100644 --- a/leap/multistep/__init__.py +++ b/leap/multistep/__init__.py @@ -74,35 +74,6 @@ class ABMonomialIntegrationFunctionFamily(ABIntegrationFunctionFamily): return 1/(func_idx+1) * x**(func_idx+1) -class ABTrigMonomialIntegrationFunctionFamily(ABIntegrationFunctionFamily): - # FIXME: Doesn't yet work for on-the-fly coefficients - - def __init__(self, order, alpha): - self.order = order - self.alpha = alpha - - def __len__(self): - return self.order - - def evaluate(self, func_idx, x): - func_idx += 1 - n = func_idx // 2 - if func_idx % 2 == 0: - return np.sin(self.alpha*n*x) - else: - return np.cos(self.alpha*n*x) - - def antiderivative(self, func_idx, x): - func_idx += 1 - n = func_idx // 2 - if func_idx == 1: - return x - elif func_idx % 2 == 0: - return -1/(n*self.alpha) * np.cos(self.alpha*n*x) - else: - return 1/(n*self.alpha) * np.sin(self.alpha*n*x) - - def _emit_func_family_operation(cb, name_gen, function_family, time_values, hist_vars, rhs_func): if isinstance(time_values, var): -- GitLab From 23c62add40d3803ba33fe0317db164628d5371ec Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 30 Jan 2018 21:53:23 -0600 Subject: [PATCH 10/13] Remove outdated dagrt branch requirement --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index b4ec38b..f8db6ac 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,3 @@ git+git://github.com/inducer/pytools git+git://github.com/inducer/pymbolic -git+https://gitlab.tiker.net/inducer/dagrt.git@cory-lls +git+https://gitlab.tiker.net/inducer/dagrt.git -- GitLab From 91f1168533c0bdfe6d7f6e6dde5dde94058e19b3 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 30 Jan 2018 22:08:36 -0600 Subject: [PATCH 11/13] Remove errant trig test (must have added this) --- test/test_ab.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/test_ab.py b/test/test_ab.py index d881284..69aac39 100755 --- a/test/test_ab.py +++ b/test/test_ab.py @@ -39,10 +39,6 @@ from utils import ( # noqa for order in [1, 3, 5] for static_dt in [True, False] ] + [ - (AdamsBashforthMethod("y", - multistep.ABTrigMonomialIntegrationFunctionFamily(4, 0.1), - static_dt=True), 3) - ] + [ (AdamsBashforthMethod("y", order, hist_length=order+1, static_dt=static_dt), order) for order in [1, 3, 5] -- GitLab From 8c2cf807f8a342fbb13e31c8c05a484d6cb600d7 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Tue, 30 Jan 2018 22:59:26 -0600 Subject: [PATCH 12/13] Yet another Flake8 failure --- test/test_ab.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_ab.py b/test/test_ab.py index 69aac39..f3db445 100755 --- a/test/test_ab.py +++ b/test/test_ab.py @@ -27,7 +27,6 @@ THE SOFTWARE. import sys import pytest from leap.multistep import AdamsBashforthMethod -import leap.multistep as multistep from utils import ( # noqa python_method_impl_interpreter as pmi_int, -- GitLab From 6b2ca4dd181c378e3ec8a44961569e5473781498 Mon Sep 17 00:00:00 2001 From: Cory Mikida Date: Fri, 2 Feb 2018 16:47:58 -0600 Subject: [PATCH 13/13] Pares down testing space a bit --- test/test_codegen_fortran.py | 18 ++++++------------ test/test_multirate.py | 4 ++-- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/test/test_codegen_fortran.py b/test/test_codegen_fortran.py index 706d2d4..781b723 100755 --- a/test/test_codegen_fortran.py +++ b/test/test_codegen_fortran.py @@ -215,18 +215,15 @@ def test_rk_codegen_fancy(): # }}} -@pytest.mark.parametrize(("min_order", "hist_length"), [(2, 2), (2, 3), (3, 3), - (3, 4), (4, 4), (4, 5), (5, 5), (5, 6)]) +@pytest.mark.parametrize("min_order", [2, 3, 4, 5]) @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) -def test_multirate_codegen(min_order, hist_length, method_name): +def test_multirate_codegen(min_order, method_name): from leap.multistep.multirate import TwoRateAdamsBashforthMethod stepper = TwoRateAdamsBashforthMethod( method_name, min_order, 4, slow_state_filter_name="slow_filt", - fast_state_filter_name="fast_filt", - hist_length_slow=hist_length, - hist_length_fast=hist_length) + fast_state_filter_name="fast_filt") code = stepper.generate() @@ -317,9 +314,6 @@ def test_multirate_codegen(min_order, hist_length, method_name): else: fac = 5 - if (hist_length != min_order): - fac = 12 - num_trips_one = 10*fac num_trips_two = 100*fac @@ -409,8 +403,8 @@ def test_adaptive_rk_codegen_error(): ]) -@pytest.mark.parametrize(("min_order", "hist_length"), - [(2, 2), (2, 3), (3, 3), (3, 4), (4, 4), (4, 5), (5, 5), (5, 6), ]) +@pytest.mark.parametrize(("min_order", "hist_length"), [(5, 5), (4, 4), (4, 5), + (3, 3), (3, 4), (2, 2), ]) def test_singlerate_squarewave(min_order, hist_length): from leap.multistep import AdamsBashforthMethod @@ -458,7 +452,7 @@ def test_singlerate_squarewave(min_order, hist_length): @pytest.mark.parametrize("method_name", TwoRateAdamsBashforthMethod.methods) @pytest.mark.parametrize(("min_order", "hist_length"), [(4, 4), (4, 5), - (3, 3), (3, 4), (2, 2), (2, 3), ]) + (3, 3), (3, 4), (2, 2), ]) def test_multirate_squarewave(min_order, hist_length, method_name): stepper = TwoRateAdamsBashforthMethod(method_name, min_order, 4, hist_length_slow=hist_length, hist_length_fast=hist_length) diff --git a/test/test_multirate.py b/test/test_multirate.py index c5e1d17..c701620 100644 --- a/test/test_multirate.py +++ b/test/test_multirate.py @@ -169,7 +169,7 @@ class MultirateTimestepperAccuracyChecker(object): @pytest.mark.slowtest -@pytest.mark.parametrize(("order", "hist_length"), [(1, 1), (1, 2), +@pytest.mark.parametrize(("order", "hist_length"), [(1, 1), (3, 3), (3, 4), ]) @pytest.mark.parametrize("system", [ #"Basic", @@ -197,7 +197,7 @@ def test_multirate_accuracy(method_name, order, hist_length, method_impl=pmi_cg)() -def test_single_rate_identical(order=3, hist_length=4): +def test_single_rate_identical(order=3, hist_length=3): from leap.multistep import AdamsBashforthMethod from dagrt.exec_numpy import NumpyInterpreter -- GitLab