From dc1f84381db9fa52b6e6d0242aea1a262e65d6a3 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 9 Dec 2019 15:30:23 -0600 Subject: [PATCH 1/8] Remove expand --- sumpy/expansion/local.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 80912077..08d2c017 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -130,7 +130,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): evaluated_coeffs = ( self.derivative_wrangler.get_full_kernel_derivatives_from_stored( coeffs, rscale)) - bvec = bvec * rscale**-1 + bvec = [b*rscale**-1 for b in bvec] + result = sum( coeff * mi_power(bvec, mi, evaluate=False) @@ -186,16 +187,17 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): result.append(sym.Add(*local_result)) else: from sumpy.tools import MiDerivativeTaker - expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=src_rscale) - taker = MiDerivativeTaker(expr, dvec) - # Rscale/operand magnitude is fairly sensitive to the order of # operations--which is something we don't have fantastic control - # over at the symbolic level. The '.expand()' below moves the two - # canceling "rscales" closer to each other in the hope of helping + # over at the symbolic level. Using rscale=1 when evaluating with + # dvec divided by rscale below moves the two cancelling "rscales" + # closer to each other at the end in the hope of helping # with that. + dvec = [d/src_rscale for d in dvec] + expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=sym.Integer(1)) + taker = MiDerivativeTaker(expr, dvec) result = [ - (taker.diff(mi) * tgt_rscale**sum(mi)).expand() + (taker.diff(mi) * sym.UnevaluatedExpr(tgt_rscale/src_rscale)**sum(mi)) for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") -- GitLab From cad5811c314db3710ecffe6eb8dd3e55f05f6fad Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 9 Dec 2019 16:19:29 -0600 Subject: [PATCH 2/8] p^(2d-2) algorithm for L2L --- sumpy/expansion/local.py | 57 ++++++++++++++++++++++++++++++++-------- 1 file changed, 46 insertions(+), 11 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 08d2c017..06edced7 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -126,19 +126,52 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): for mi in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale): - from sumpy.tools import mi_power, mi_factorial + from pytools import factorial evaluated_coeffs = ( self.derivative_wrangler.get_full_kernel_derivatives_from_stored( coeffs, rscale)) bvec = [b*rscale**-1 for b in bvec] - - result = sum( - coeff - * mi_power(bvec, mi, evaluate=False) - / mi_factorial(mi) - for coeff, mi in zip( - evaluated_coeffs, self.get_full_coefficient_identifiers())) - return result + mi_to_index = dict((mi, i) for i, mi in + enumerate(self.get_full_coefficient_identifiers())) + sorted_mis = sorted(self.get_full_coefficient_identifiers()) + dim = self.dim + + # Start with a previously unseen multi-index + seen_mi = [-1]*dim + # Local sum keep the sum of the terms that differ by that dimension + local_sum = [0]*dim + # Local multiplier keep the scalar that the sum has to be multiplied by + # when adding to the sum of the previous dimension. + local_multiplier = [0]*dim + + for mi in sorted_mis: + # Iterate in reverse order as the mis are sorted in row-major order + for d in reversed(range(dim-1)): + if seen_mi[d] != mi[d]: + # If the dimension d of mi changed from the previous value + # then the sum for dimension d+1 is complete, add it to + # dimension d after multiplying and restart + local_sum[d] += local_sum[d+1]*local_multiplier[d+1] + local_sum[d+1] = 0 + local_multiplier[d+1] = bvec[d]**mi[d] / factorial(mi[d]) + local_sum[dim-1] += evaluated_coeffs[mi_to_index[mi]] * \ + bvec[dim-1]**mi[dim-1] / factorial(mi[dim-1]) + seen_mi = mi + + for d in reversed(range(dim-1)): + local_sum[d] += local_sum[d+1]*local_multiplier[d+1] + + # Above code is functionally equivalent to the following + # + # from sumpy.tools import mi_power, mi_factorial + # result = sum( + # coeff + # * mi_power(bvec, mi, evaluate=False) + # / mi_factorial(mi) + # for coeff, mi in zip( + # evaluated_coeffs, self.get_full_coefficient_identifiers())) + + return local_sum[0] def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, dvec, tgt_rscale): @@ -194,10 +227,12 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # closer to each other at the end in the hope of helping # with that. dvec = [d/src_rscale for d in dvec] - expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=sym.Integer(1)) + expr = src_expansion.evaluate(src_coeff_exprs, dvec, + rscale=sym.Integer(1)) taker = MiDerivativeTaker(expr, dvec) + rscale_ratio = sym.UnevaluatedExpr(tgt_rscale/src_rscale) result = [ - (taker.diff(mi) * sym.UnevaluatedExpr(tgt_rscale/src_rscale)**sum(mi)) + (taker.diff(mi) * rscale_ratio**sum(mi)) for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") -- GitLab From ca5db32e3416e28be2213849c6514f4c85ed3dc9 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 9 Dec 2019 17:49:05 -0600 Subject: [PATCH 3/8] fix tests --- sumpy/expansion/local.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 06edced7..e1482dde 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -226,13 +226,13 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # dvec divided by rscale below moves the two cancelling "rscales" # closer to each other at the end in the hope of helping # with that. - dvec = [d/src_rscale for d in dvec] expr = src_expansion.evaluate(src_coeff_exprs, dvec, rscale=sym.Integer(1)) + replace_dict = dict((d, d/src_rscale) for d in dvec) taker = MiDerivativeTaker(expr, dvec) rscale_ratio = sym.UnevaluatedExpr(tgt_rscale/src_rscale) result = [ - (taker.diff(mi) * rscale_ratio**sum(mi)) + (taker.diff(mi).xreplace(replace_dict) * rscale_ratio**sum(mi)) for mi in self.get_coefficient_identifiers()] logger.info("building translation operator: done") -- GitLab From 438619fe98a685f8d38e4ef2b117b195d00ef738 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Mon, 9 Dec 2019 19:13:41 -0600 Subject: [PATCH 4/8] Fix rescaling --- sumpy/expansion/local.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index e1482dde..842ac677 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -222,12 +222,14 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): from sumpy.tools import MiDerivativeTaker # Rscale/operand magnitude is fairly sensitive to the order of # operations--which is something we don't have fantastic control - # over at the symbolic level. Using rscale=1 when evaluating with - # dvec divided by rscale below moves the two cancelling "rscales" - # closer to each other at the end in the hope of helping - # with that. - expr = src_expansion.evaluate(src_coeff_exprs, dvec, - rscale=sym.Integer(1)) + # over at the symbolic level. Scaling dvec, then differentiating, + # and finally rescaling dvec leaves the expression needing a scaling + # to compensate for differentiating which is done at the end. + # This moves the two cancelling "rscales" closer to each other at + # the end in the hope of helping rscale magnitude. + dvec_scaled = [d*src_rscale for d in dvec] + expr = src_expansion.evaluate(src_coeff_exprs, dvec_scaled, + rscale=src_rscale) replace_dict = dict((d, d/src_rscale) for d in dvec) taker = MiDerivativeTaker(expr, dvec) rscale_ratio = sym.UnevaluatedExpr(tgt_rscale/src_rscale) -- GitLab From 2dbfb32d8e40ae417c7c9b7c1b4ea1eade07224f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Tue, 10 Dec 2019 18:34:15 +0100 Subject: [PATCH 5/8] Initial readability improvements to new L2L --- sumpy/expansion/local.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 842ac677..57e0cf51 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -133,7 +133,9 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): bvec = [b*rscale**-1 for b in bvec] mi_to_index = dict((mi, i) for i, mi in enumerate(self.get_full_coefficient_identifiers())) - sorted_mis = sorted(self.get_full_coefficient_identifiers()) + + # Sort multi-indices so that last dimension varies fastest + sorted_target_mis = sorted(self.get_full_coefficient_identifiers()) dim = self.dim # Start with a previously unseen multi-index @@ -144,8 +146,8 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # when adding to the sum of the previous dimension. local_multiplier = [0]*dim - for mi in sorted_mis: - # Iterate in reverse order as the mis are sorted in row-major order + for mi in sorted_target_mis: + # Iterate in reverse order as the mis are sorted last-varies-fastest for d in reversed(range(dim-1)): if seen_mi[d] != mi[d]: # If the dimension d of mi changed from the previous value -- GitLab From 2233c7990aba1a7c8e3bd0ddfa8c0ba3094a7351 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 10 Dec 2019 20:25:39 +0100 Subject: [PATCH 6/8] Remove empty line with spaces --- sumpy/expansion/local.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 57e0cf51..f07a0fc3 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -133,7 +133,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): bvec = [b*rscale**-1 for b in bvec] mi_to_index = dict((mi, i) for i, mi in enumerate(self.get_full_coefficient_identifiers())) - + # Sort multi-indices so that last dimension varies fastest sorted_target_mis = sorted(self.get_full_coefficient_identifiers()) dim = self.dim -- GitLab From b92cab9ef6658c5bdb8f5da82d32bc60bb70a5dd Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 10 Dec 2019 13:44:13 -0600 Subject: [PATCH 7/8] Add comment about local_sum --- sumpy/expansion/local.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index f07a0fc3..74cb3444 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -146,6 +146,25 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # when adding to the sum of the previous dimension. local_multiplier = [0]*dim + # For the multi-indices for 3D, local_sum, looks like below + # + # Multi-index | coef | local_sum | local_mult + # (0, 0, 0) | c0 | 0, 0, c0 | 0, 1, 1 + # (0, 0, 1) | c1 | 0, 0, c0+c1*dz | 0, 1, 1 + # (0, 0, 2) | c2 | 0, 0, c0+c1*dz+c2*dz^2 | 0, 1, 1 + # (0, 1, 0) | c3 | 0, c0+c1*dz+c2*dz^2, c3 | 0, 1, dy + # (0, 1, 1) | c4 | 0, c0+c1*dz+c2*dz^2, c3+c4*dz | 0, 1, dy + # (0, 1, 2) | c5 | 0, c0+c1*dz+c2*dz^2, c3+c4*dz+c5*dz^2 | 0, 1, dy + # (0, 2, 0) | c6 | 0, c0+c1*dz+c2*dz^2, c6 | 0, 1, dy^2 + # | | +dy*(c3+c4*dz+c5*dz^2) | + # (0, 2, 1) | c7 | 0, c0+c1*dz+c2*dz^2, c6+c7*dz | 0, 1, dy^2 + # | | +dy*(c3+c4*dz+c5*dz^2) | + # (0, 2, 2) | c8 | 0, c0+c1*dz+c2*dz^2, c6+c7*dz+x8*dz^2 | 0, 1, dy^2 + # | | +dy*(c3+c4*dz+c5*dz^2) | + # (1, 0, 0) | c8 | c0+c1*dz+c2*dz^2, 0, 0 | 0, dx, 1 + # | | +dy*(c3+c4*dz+c5*dz^2) | + # | | +dy^2*(c6+c7*dz+c8*dz^2) | + for mi in sorted_target_mis: # Iterate in reverse order as the mis are sorted last-varies-fastest for d in reversed(range(dim-1)): -- GitLab From 3d5051ee514653cf39562f6649c10b858b2c76f9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Fri, 10 Jan 2020 23:10:24 +0100 Subject: [PATCH 8/8] More comment changes to make new L2L easier to understand --- sumpy/expansion/local.py | 47 +++++++++++++++++++++++++++------------- 1 file changed, 32 insertions(+), 15 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 74cb3444..b4b435de 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -127,9 +127,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): def evaluate(self, coeffs, bvec, rscale): from pytools import factorial + evaluated_coeffs = ( self.derivative_wrangler.get_full_kernel_derivatives_from_stored( coeffs, rscale)) + bvec = [b*rscale**-1 for b in bvec] mi_to_index = dict((mi, i) for i, mi in enumerate(self.get_full_coefficient_identifiers())) @@ -138,15 +140,15 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): sorted_target_mis = sorted(self.get_full_coefficient_identifiers()) dim = self.dim - # Start with a previously unseen multi-index + # Start with an invalid "seen" multi-index seen_mi = [-1]*dim - # Local sum keep the sum of the terms that differ by that dimension + # Local sum keep the sum of the terms that differ by each dimension local_sum = [0]*dim # Local multiplier keep the scalar that the sum has to be multiplied by - # when adding to the sum of the previous dimension. + # when adding to the sum of the preceding dimension. local_multiplier = [0]*dim - # For the multi-indices for 3D, local_sum, looks like below + # For the multi-indices in 3D, local_sum looks like this: # # Multi-index | coef | local_sum | local_mult # (0, 0, 0) | c0 | 0, 0, c0 | 0, 1, 1 @@ -166,15 +168,26 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): # | | +dy^2*(c6+c7*dz+c8*dz^2) | for mi in sorted_target_mis: - # Iterate in reverse order as the mis are sorted last-varies-fastest + + # {{{ handle the case where a not-last dimension "clicked over" + + # (where d will be that not-last dimension) + + # Iterate in reverse order of dimensions to properly handle a + # "double click-over". + for d in reversed(range(dim-1)): if seen_mi[d] != mi[d]: # If the dimension d of mi changed from the previous value # then the sum for dimension d+1 is complete, add it to - # dimension d after multiplying and restart + # dimension d after multiplying and restart. + local_sum[d] += local_sum[d+1]*local_multiplier[d+1] local_sum[d+1] = 0 local_multiplier[d+1] = bvec[d]**mi[d] / factorial(mi[d]) + + # }}} + local_sum[dim-1] += evaluated_coeffs[mi_to_index[mi]] * \ bvec[dim-1]**mi[dim-1] / factorial(mi[dim-1]) seen_mi = mi @@ -182,15 +195,19 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): for d in reversed(range(dim-1)): local_sum[d] += local_sum[d+1]*local_multiplier[d+1] - # Above code is functionally equivalent to the following - # - # from sumpy.tools import mi_power, mi_factorial - # result = sum( - # coeff - # * mi_power(bvec, mi, evaluate=False) - # / mi_factorial(mi) - # for coeff, mi in zip( - # evaluated_coeffs, self.get_full_coefficient_identifiers())) + # {{{ simpler, functionally equivalent code + + if 0: + from sumpy.tools import mi_power, mi_factorial + + return sum( + coeff + * mi_power(bvec, mi, evaluate=False) + / mi_factorial(mi) + for coeff, mi in zip( + evaluated_coeffs, self.get_full_coefficient_identifiers())) + + # }}} return local_sum[0] -- GitLab