From ec79a05759952c7a730b6b5fc5db27d58dffe10b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andreas=20Kl=C3=B6ckner?= Date: Tue, 23 Jul 2019 23:41:08 +0200 Subject: [PATCH 01/16] Conda CI: Upgrade symengine to 0.4, Python to 3.7 --- .test-conda-env-py3.yml | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index a696dc8d..1b7ee47c 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -10,7 +10,9 @@ dependencies: - pocl - islpy - pyopencl -- python=3.6 -- symengine=0.3.0 -- python-symengine=0.3.0 +- python=3.7 +- symengine=0.4.0 +- python-symengine=0.4.0 + # things not in here: loopy boxtree pymbolic pyfmmlib +# (see .test-conda-env-py3-requirements.txt) \ No newline at end of file -- GitLab From 20bff213b6a92fa5ddca03510a1f00643ea31ecd Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Wed, 24 Jul 2019 02:22:36 +0200 Subject: [PATCH 02/16] Use sympy 1.4 for python 3.7 --- asv.conf.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/asv.conf.json b/asv.conf.json index 9632c4fb..6d6fd675 100644 --- a/asv.conf.json +++ b/asv.conf.json @@ -68,7 +68,7 @@ // }, "matrix": { "numpy" : [""], - "sympy" : ["1.0"], + "sympy" : ["1.4"], "pyopencl" : [""], "islpy" : [""], "pocl" : [""], -- GitLab From 43f60d6f86b2e6364281185b33dd4010ca435117 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Wed, 11 Sep 2019 06:27:29 +0200 Subject: [PATCH 03/16] Update symengine to 0.5 --- .test-conda-env-py3.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 1b7ee47c..89c055eb 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -11,8 +11,8 @@ dependencies: - islpy - pyopencl - python=3.7 -- symengine=0.4.0 -- python-symengine=0.4.0 +- symengine=0.5.0 +- python-symengine=0.5.0 # things not in here: loopy boxtree pymbolic pyfmmlib # (see .test-conda-env-py3-requirements.txt) \ No newline at end of file -- GitLab From b65a4ab5592e3365dc4f1f5eadd0f10cd23cd011 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Thu, 12 Sep 2019 07:28:05 +0200 Subject: [PATCH 04/16] Fix bad merge --- .test-conda-env-py3.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index e0434a18..0c9031a0 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -12,8 +12,8 @@ dependencies: - islpy - pyopencl - python=3 -- symengine=0.3.0 -- python-symengine=0.3.0 +- symengine=0.5.0 +- python-symengine=0.5.0 - pyfmmlib - pip -- GitLab From 7e76990df983c967a9e1ad08b0f3b574bc40cd38 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 21 Jan 2020 16:25:12 +0100 Subject: [PATCH 05/16] Check symengine 0.5.1 --- .test-conda-env-py3.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 0c9031a0..609b7d95 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -13,7 +13,7 @@ dependencies: - pyopencl - python=3 - symengine=0.5.0 -- python-symengine=0.5.0 +- python-symengine=0.5.1 - pyfmmlib - pip -- GitLab From a9af19b2ded203bbbaa354648638d229a392debc Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Sun, 9 Feb 2020 21:03:04 -0600 Subject: [PATCH 06/16] Distribute an integer multiplication Reduces number of operations of M2M by 10% for Laplace 2D order 16 --- sumpy/expansion/multipole.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 068fea23..e3759d43 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -170,8 +170,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): src_coeff_exprs = list(src_coeff_exprs) for i, mi in enumerate(src_expansion.get_coefficient_identifiers()): - src_coeff_exprs[i] *= mi_factorial(mi) * \ - sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(mi) + src_coeff_exprs[i] *= sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(mi) result = [0] * len(self.get_full_coefficient_identifiers()) @@ -219,21 +218,15 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): continue contrib = dim_coeffs_to_translate[src_index] - for idim in range(self.dim): n = tgt_mi[idim] k = src_mi[idim] assert n >= k - from sympy import binomial - contrib *= (binomial(n, k) - * sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k)) + contrib /= mi_factorial((n-k,)) + contrib *= sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k) result[i] += contrib - # Defer division by target factorial until the very end - if d == self.dim-1: - result[i] /= mi_factorial(tgt_mi) - dim_coeffs_to_translate = result[:] mi_to_index = tgt_mi_to_index -- GitLab From 647c95a2dd6ba35497285c527488118847042275 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 11 Feb 2020 06:02:40 +0100 Subject: [PATCH 07/16] Update to symengine-0.6.0 --- .test-conda-env-py3.yml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml index 609b7d95..63214622 100644 --- a/.test-conda-env-py3.yml +++ b/.test-conda-env-py3.yml @@ -12,8 +12,7 @@ dependencies: - islpy - pyopencl - python=3 -- symengine=0.5.0 -- python-symengine=0.5.1 +- python-symengine=0.6.0 - pyfmmlib - pip -- GitLab From b99da0248fb7c5f1f0c2d7e338e3c5b2de97bd5f Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Fri, 21 Feb 2020 14:39:37 -0600 Subject: [PATCH 08/16] O(p^d) algorithm for L2L --- sumpy/expansion/__init__.py | 40 ++++++++++++++++++++++- sumpy/expansion/local.py | 63 +++++++++++++++++++++++++++++++++++-- 2 files changed, 100 insertions(+), 3 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 85c93850..6c87e252 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -218,6 +218,45 @@ class ExpansionTermsWrangler(object): return type(self)(**new_kwargs) + @memoize_method + def _get_coeff_identifier_split(self): + """ + This splits the coefficients into a O(p) number of sets + so that for each set, all the identifiers have the form, + (m_1, m_2, ..., m_{j-1}, c, m_{j+1}, ... , m_d) + where c is a constant. The intersection of two sets might + not be empty. + + If this is an instance of LinearPDEBasedExpansionTermsWrangler, + then the number of sets will be O(1). + + Returns a List[Tuple[j, List[identifiers]]] + """ + res = [] + mis = self.get_full_coefficient_identifiers() + coeff_ident_enumerate_dict = dict((tuple(mi), i) for + (i, mi) in enumerate(mis)) + + max_mi = None + if isinstance(self, LinearPDEBasedExpansionTermsWrangler): + pde_dict = self.get_pde().eq + for ident in pde_dict.keys(): + if ident not in coeff_ident_enumerate_dict: + break + else: + max_mi_idx = max(coeff_ident_enumerate_dict[ident] for + ident in pde_dict.keys()) + max_mi = mis[max_mi_idx] + + if max_mi is None: + max_mi = [max(mi[d] for mi in mis) for d in range(self.dim)] + + for d in range(self.dim): + filtered_mis = [mi for mi in mis if mi[d] < max_mi[d]] + for i in range(max_mi[d]): + res.append((d, [mi for mi in filtered_mis if mi[d] == i])) + return res + class FullExpansionTermsWrangler(ExpansionTermsWrangler): @@ -231,7 +270,6 @@ class FullExpansionTermsWrangler(ExpansionTermsWrangler): get_stored_mpole_coefficients_from_full = ( get_full_kernel_derivatives_from_stored) - # {{{ sparse matrix-vector multiplication def _spmv(spmat, x, sparse_vectors): diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index a58c5981..8efad734 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -295,8 +295,67 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): lexp_mi_terms.append( coeff * kernel_deriv * tgt_rscale**sum(lexp_mi)) result.append(sym.Add(*lexp_mi_terms)) + logger.info("building translation operator: done") + return result - else: + rscale_ratio = sym.UnevaluatedExpr(tgt_rscale/src_rscale) + + from sumpy.tools import MiDerivativeTaker + from math import factorial + src_coeffs = ( + src_expansion.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored( + src_coeff_exprs, src_rscale)) + src_mis = \ + src_expansion.expansion_terms_wrangler.get_full_coefficient_identifiers() + + src_mi_to_index = dict((mi, i) for i, mi in enumerate(src_mis)) + + tgt_mis_all = \ + self.expansion_terms_wrangler.get_coefficient_identifiers() + tgt_mi_to_index = dict((mi, i) for i, mi in enumerate(tgt_mis_all)) + + tgt_split = self.expansion_terms_wrangler._get_coeff_identifier_split() + + p = max(sum(mi) for mi in src_mis) + Y = src_coeffs + result = [0] * len(tgt_mis_all) + + # O(1) iterations + for const_dim in set(d for d, _ in tgt_split): + # Use the const_dim as the first dimension to vary so that the below + # algorithm is O(p^{d+1}) for full and O(p^{d}) for compressed + dims = [const_dim] + list(range(const_dim)) + \ + list(range(const_dim+1, self.dim)) + # Start with source coefficients + Y = src_coeffs + # O(1) iterations + for d in dims: + C = Y + Y = [0] * len(src_mis) + # Only O(p^{d-1}) operations are used in compressed + # All of them are used in full. + for i, s in enumerate(src_mis): + # O(p) iterations + for q in range(p+1-sum(s)): + src_mi = list(s) + src_mi[d] += q + src_mi = tuple(src_mi) + if src_mi in src_mi_to_index: + Y[i] += (dvec[d]/src_rscale) ** q * \ + C[src_mi_to_index[src_mi]] / factorial(q) + + # This is O(p) in full and O(1) in compressed + for d, tgt_mis in tgt_split: + if d != const_dim: + continue + # O(p^{d-1}) iterations + for mi in tgt_mis: + if mi not in src_mi_to_index: + continue + result[tgt_mi_to_index[mi]] = Y[src_mi_to_index[mi]] * rscale_ratio ** sum(mi) + + # {{{ simpler, functionally equivalent code + if 0: 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 @@ -314,7 +373,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): result = [ (taker.diff(mi).xreplace(replace_dict) * rscale_ratio**sum(mi)) for mi in self.get_coefficient_identifiers()] - + # }}} logger.info("building translation operator: done") return result -- GitLab From b7b780d0eff2a1915fe552c11976a659a3534663 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Fri, 21 Feb 2020 14:40:44 -0600 Subject: [PATCH 09/16] Remove previous code --- sumpy/expansion/local.py | 85 ++++------------------------------------ 1 file changed, 8 insertions(+), 77 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 8efad734..e12b513a 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -134,83 +134,14 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): 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())) - - # Sort multi-indices so that last dimension varies fastest - sorted_target_mis = sorted(self.get_full_coefficient_identifiers()) - dim = self.dim - - # Start with an invalid "seen" multi-index - seen_mi = [-1]*dim - # 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 preceding dimension. - local_multiplier = [0]*dim - - # 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 - # (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: - - # {{{ 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. - - 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] - - # {{{ 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] + 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())) def translate_from(self, src_expansion, src_coeff_exprs, src_rscale, dvec, tgt_rscale): -- GitLab From 3a625ac4341459ad1231a08b1def08e820aaef36 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 25 Feb 2020 09:49:34 -0600 Subject: [PATCH 10/16] Improved M2M --- sumpy/expansion/__init__.py | 11 ++-- sumpy/expansion/multipole.py | 110 +++++++++++++++++++++-------------- 2 files changed, 73 insertions(+), 48 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index 6c87e252..a1bf7a27 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -221,11 +221,10 @@ class ExpansionTermsWrangler(object): @memoize_method def _get_coeff_identifier_split(self): """ - This splits the coefficients into a O(p) number of sets + This splits the coefficients into a O(p) number of disjoint sets so that for each set, all the identifiers have the form, (m_1, m_2, ..., m_{j-1}, c, m_{j+1}, ... , m_d) - where c is a constant. The intersection of two sets might - not be empty. + where c is a constant. If this is an instance of LinearPDEBasedExpansionTermsWrangler, then the number of sets will be O(1). @@ -251,10 +250,14 @@ class ExpansionTermsWrangler(object): if max_mi is None: max_mi = [max(mi[d] for mi in mis) for d in range(self.dim)] + seen_mis = set() for d in range(self.dim): filtered_mis = [mi for mi in mis if mi[d] < max_mi[d]] for i in range(max_mi[d]): - res.append((d, [mi for mi in filtered_mis if mi[d] == i])) + new_mis = [mi for mi in filtered_mis if mi[d] == i and \ + mi not in seen_mis] + seen_mis.update(new_mis) + res.append((d, new_mis)) return res diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index bdd66ac7..4587233a 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -164,14 +164,17 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): from sumpy.tools import mi_factorial src_mi_to_index = dict((mi, i) for i, mi in enumerate( - src_expansion.get_coefficient_identifiers())) + src_expansion.get_full_coefficient_identifiers())) tgt_mi_to_index = dict((mi, i) for i, mi in enumerate( self.get_full_coefficient_identifiers())) - src_coeff_exprs = list(src_coeff_exprs) + src_coeff_exprs_full = \ + [0]*len(src_expansion.get_full_coefficient_identifiers()) + for i, mi in enumerate(src_expansion.get_coefficient_identifiers()): - src_coeff_exprs[i] *= sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(mi) + src_coeff_exprs_full[src_mi_to_index[mi]] = src_coeff_exprs[i] * \ + sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(mi) result = [0] * len(self.get_full_coefficient_identifiers()) @@ -182,54 +185,73 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): # d_x^i d_y^j \binom{m}{i} \binom{n}{j}$ # and can be rewritten as follows. # - # Let $T_{m, n} = \sum_{i\le m} A_{i, n} d_x^i \binom{m}{i}$. + # Let $Y_{m, n} = \sum_{i\le m} A_{i, n} d_x^i \binom{m}{i}$. # - # Then, $B_{m, n} = \sum_{j\le n} T_{m, j} d_y^j \binom{n}{j}$. + # Then, $B_{m, n} = \sum_{j\le n} Y_{m, j} d_y^j \binom{n}{j}$. # - # $T_{m, n}$ are $p^2$ number of temporary variables that are + # $Y_{m, n}$ are $p^2$ number of temporary variables that are # reused for different M2M coefficients and costs $p$ per variable. - # Total cost for calculating $T_{m, n}$ is $p^3$ and similar + # Total cost for calculating $Y_{m, n}$ is $p^3$ and similar # for $B_{m, n}$ # In other words, we're better off computing the translation - # one dimension at a time. This is realized here by using - # the output from one dimension as the input to the next. - # per_dim_coeffs_to_translate serves as the array of input - # coefficients for each dimension's translation. - - dim_coeffs_to_translate = src_coeff_exprs - - mi_to_index = src_mi_to_index - for d in range(self.dim): - result = [0] * len(self.get_full_coefficient_identifiers()) - for i, tgt_mi in enumerate( - self.get_full_coefficient_identifiers()): - - src_mis_per_dim = [] - for mi_i in range(tgt_mi[d]+1): - new_mi = list(tgt_mi) - new_mi[d] = mi_i - src_mis_per_dim.append(tuple(new_mi)) - - for src_mi in src_mis_per_dim: - try: - src_index = mi_to_index[src_mi] - except KeyError: - # Omitted coefficients: not life-threatening - continue - - contrib = dim_coeffs_to_translate[src_index] - for idim in range(self.dim): - n = tgt_mi[idim] - k = src_mi[idim] - assert n >= k - contrib /= mi_factorial((n-k,)) - contrib *= sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k) - - result[i] += contrib + # one dimension at a time. If the coefficients in the source + # expansion have the form (0, m) and (m, 0), then we calculate + # the output from (0, m) with the second dimension as the fastest + # varying dimension and then calculate the output from (m, 0) + # with the first dimension as the fastest varying dimension. + # However, the contribution from (0, 0) is counted twice and we + # need to subtract that. + + src_split = \ + src_expansion.expansion_terms_wrangler._get_coeff_identifier_split() + result = [0] * len(self.get_full_coefficient_identifiers()) - dim_coeffs_to_translate = result[:] - mi_to_index = tgt_mi_to_index + non_zero_coeffs_per_dim = [[] for d in range(self.dim)] + for d, src_mi in src_split: + non_zero_coeffs_per_dim[d] += src_mi + + for const_dim in set(d for d, _ in src_split): + dim_coeffs_to_translate = \ + [0] * len(src_expansion.get_full_coefficient_identifiers()) + for mi in non_zero_coeffs_per_dim[const_dim]: + idx = src_mi_to_index[mi] + dim_coeffs_to_translate[idx] = src_coeff_exprs_full[idx] + + # Use the const_dim as the last dimension to vary + dims = list(range(const_dim)) + \ + list(range(const_dim+1, self.dim)) + [const_dim] + for d in dims: + temp = [0] * len(src_expansion.get_full_coefficient_identifiers()) + for i, tgt_mi in enumerate( + src_expansion.get_full_coefficient_identifiers()): + src_mis_per_dim = [] + for mi_i in range(tgt_mi[d]+1): + new_mi = list(tgt_mi) + new_mi[d] = mi_i + src_mis_per_dim.append(tuple(new_mi)) + + for src_mi in src_mis_per_dim: + try: + src_index = src_mi_to_index[src_mi] + except KeyError: + # Omitted coefficients: not life-threatening + continue + + contrib = dim_coeffs_to_translate[src_index] + for idim in range(self.dim): + n = tgt_mi[idim] + k = src_mi[idim] + assert n >= k + contrib /= mi_factorial((n-k,)) + contrib *= sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k) + + temp[i] += contrib + + dim_coeffs_to_translate = temp[:] + + for i, mi in enumerate(src_expansion.get_full_coefficient_identifiers()): + result[tgt_mi_to_index[mi]] = dim_coeffs_to_translate[i] logger.info("building translation operator: done") return ( -- GitLab From 7de8a6f1e5fa0b5a2a631becab0f558a9549e61d Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 25 Feb 2020 12:35:45 -0600 Subject: [PATCH 11/16] Missed a plus --- sumpy/expansion/multipole.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 4587233a..fbba6acb 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -176,8 +176,6 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): src_coeff_exprs_full[src_mi_to_index[mi]] = src_coeff_exprs[i] * \ sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(mi) - result = [0] * len(self.get_full_coefficient_identifiers()) - # This algorithm uses the observation that M2M coefficients # have the following form in 2D # @@ -251,7 +249,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): dim_coeffs_to_translate = temp[:] for i, mi in enumerate(src_expansion.get_full_coefficient_identifiers()): - result[tgt_mi_to_index[mi]] = dim_coeffs_to_translate[i] + result[tgt_mi_to_index[mi]] += dim_coeffs_to_translate[i] logger.info("building translation operator: done") return ( -- GitLab From d7bfbb7bedc36661c6aa2b1c1df643ac8c770e0c Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 25 Feb 2020 12:36:18 -0600 Subject: [PATCH 12/16] Keep the old code --- sumpy/expansion/multipole.py | 38 ++++++++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index fbba6acb..13411e3d 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -251,6 +251,44 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): for i, mi in enumerate(src_expansion.get_full_coefficient_identifiers()): result[tgt_mi_to_index[mi]] += dim_coeffs_to_translate[i] + if 0: # NOQA + src_mi_to_index = dict((mi, i) for i, mi in enumerate( + src_expansion.get_coefficient_identifiers())) + result = [0] * len(self.get_full_coefficient_identifiers()) + + for i, mi in enumerate(src_expansion.get_coefficient_identifiers()): + src_coeff_exprs[i] *= mi_factorial(mi) + + from pytools import generate_nonnegative_integer_tuples_below as gnitb + + for i, tgt_mi in enumerate( + self.get_full_coefficient_identifiers()): + + tgt_mi_plus_one = tuple(mi_i + 1 for mi_i in tgt_mi) + + for src_mi in gnitb(tgt_mi_plus_one): + try: + src_index = src_mi_to_index[src_mi] + except KeyError: + # Omitted coefficients: not life-threatening + continue + + contrib = src_coeff_exprs2[src_index] + + for idim in range(self.dim): + n = tgt_mi[idim] + k = src_mi[idim] + assert n >= k + from sympy import binomial + contrib *= (binomial(n, k) + * sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k)) + + result[i] += ( + contrib + * sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(src_mi)) + + result[i] /= mi_factorial(tgt_mi) + logger.info("building translation operator: done") return ( self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full( -- GitLab From efeb06d99162cf0ff083491e544d67fccb4b8e6d Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 25 Feb 2020 12:42:03 -0600 Subject: [PATCH 13/16] vim fold --- sumpy/expansion/multipole.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 13411e3d..95255ef1 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -251,7 +251,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): for i, mi in enumerate(src_expansion.get_full_coefficient_identifiers()): result[tgt_mi_to_index[mi]] += dim_coeffs_to_translate[i] - if 0: # NOQA + # {{{ simpler, functionally equivalent code + if 0: src_mi_to_index = dict((mi, i) for i, mi in enumerate( src_expansion.get_coefficient_identifiers())) result = [0] * len(self.get_full_coefficient_identifiers()) @@ -288,6 +289,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): * sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(src_mi)) result[i] /= mi_factorial(tgt_mi) + # }}} logger.info("building translation operator: done") return ( -- GitLab From 915e8aab634ca2561e6ad24167416808414c530e Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 25 Feb 2020 14:00:26 -0600 Subject: [PATCH 14/16] Formatting fixes --- sumpy/expansion/__init__.py | 5 +++-- sumpy/expansion/local.py | 9 ++++----- sumpy/expansion/multipole.py | 10 +++++----- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index a1bf7a27..c0122bcc 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -254,8 +254,8 @@ class ExpansionTermsWrangler(object): for d in range(self.dim): filtered_mis = [mi for mi in mis if mi[d] < max_mi[d]] for i in range(max_mi[d]): - new_mis = [mi for mi in filtered_mis if mi[d] == i and \ - mi not in seen_mis] + new_mis = [mi for mi in filtered_mis if mi[d] == i + and mi not in seen_mis] seen_mis.update(new_mis) res.append((d, new_mis)) return res @@ -273,6 +273,7 @@ class FullExpansionTermsWrangler(ExpansionTermsWrangler): get_stored_mpole_coefficients_from_full = ( get_full_kernel_derivatives_from_stored) + # {{{ sparse matrix-vector multiplication def _spmv(spmat, x, sparse_vectors): diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index e12b513a..00490fb4 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -127,8 +127,6 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): for mi in self.get_coefficient_identifiers()] def evaluate(self, coeffs, bvec, rscale, sac=None): - from pytools import factorial - evaluated_coeffs = ( self.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored( coeffs, rscale)) @@ -233,8 +231,9 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): from sumpy.tools import MiDerivativeTaker from math import factorial + src_wrangler = src_expansion.expansion_terms_wrangler src_coeffs = ( - src_expansion.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored( + src_wrangler.get_full_kernel_derivatives_from_stored( src_coeff_exprs, src_rscale)) src_mis = \ src_expansion.expansion_terms_wrangler.get_full_coefficient_identifiers() @@ -283,11 +282,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): for mi in tgt_mis: if mi not in src_mi_to_index: continue - result[tgt_mi_to_index[mi]] = Y[src_mi_to_index[mi]] * rscale_ratio ** sum(mi) + result[tgt_mi_to_index[mi]] = Y[src_mi_to_index[mi]] \ + * rscale_ratio ** sum(mi) # {{{ simpler, functionally equivalent code if 0: - 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. Scaling dvec, then differentiating, diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 95255ef1..36b08e15 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -242,7 +242,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): k = src_mi[idim] assert n >= k contrib /= mi_factorial((n-k,)) - contrib *= sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k) + contrib *= \ + sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k) temp[i] += contrib @@ -274,7 +275,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): # Omitted coefficients: not life-threatening continue - contrib = src_coeff_exprs2[src_index] + contrib = src_coeff_exprs[src_index] for idim in range(self.dim): n = tgt_mi[idim] @@ -284,9 +285,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): contrib *= (binomial(n, k) * sym.UnevaluatedExpr(dvec[idim]/tgt_rscale)**(n-k)) - result[i] += ( - contrib - * sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(src_mi)) + result[i] += (contrib + * sym.UnevaluatedExpr(src_rscale/tgt_rscale)**sum(src_mi)) result[i] /= mi_factorial(tgt_mi) # }}} -- GitLab From 3dbca61bd36a8fedb01cba2443bb11df50c38488 Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Tue, 25 Feb 2020 15:02:53 -0600 Subject: [PATCH 15/16] Ignore naming --- sumpy/expansion/local.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 00490fb4..51687487 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -247,7 +247,6 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): tgt_split = self.expansion_terms_wrangler._get_coeff_identifier_split() p = max(sum(mi) for mi in src_mis) - Y = src_coeffs result = [0] * len(tgt_mis_all) # O(1) iterations @@ -257,11 +256,11 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): dims = [const_dim] + list(range(const_dim)) + \ list(range(const_dim+1, self.dim)) # Start with source coefficients - Y = src_coeffs + Y = src_coeffs # noqa: N806 # O(1) iterations for d in dims: - C = Y - Y = [0] * len(src_mis) + C = Y # noqa: N806 + Y = [0] * len(src_mis) # noqa: N806 # Only O(p^{d-1}) operations are used in compressed # All of them are used in full. for i, s in enumerate(src_mis): -- GitLab From 1bcf666f074bf0f9a8964b4abb068360b7101a3a Mon Sep 17 00:00:00 2001 From: Isuru Fernando Date: Sat, 7 Mar 2020 23:37:26 -0600 Subject: [PATCH 16/16] Remove passing a sac around --- sumpy/e2p.py | 2 +- sumpy/expansion/__init__.py | 31 ++++++++++++++----------------- sumpy/expansion/local.py | 12 ++++++------ sumpy/expansion/multipole.py | 10 +++++----- sumpy/p2e.py | 2 +- 5 files changed, 27 insertions(+), 30 deletions(-) diff --git a/sumpy/e2p.py b/sumpy/e2p.py index 408b7364..7b0072ad 100644 --- a/sumpy/e2p.py +++ b/sumpy/e2p.py @@ -91,7 +91,7 @@ class E2PBase(KernelCacheWrapper): coeff_exprs = [sym.Symbol("coeff%d" % i) for i in range(len(self.expansion.get_coefficient_identifiers()))] - value = self.expansion.evaluate(coeff_exprs, bvec, rscale, sac=sac) + value = self.expansion.evaluate(coeff_exprs, bvec, rscale) result_names = [ sac.assign_unique("result_%d_p" % i, diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py index c0122bcc..9d3a917a 100644 --- a/sumpy/expansion/__init__.py +++ b/sumpy/expansion/__init__.py @@ -114,21 +114,19 @@ class ExpansionBase(object): """ raise NotImplementedError - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): """Form an expansion from a source point. :arg avec: vector from source to center. :arg bvec: vector from center to target. Not usually necessary, except for line-Taylor expansion. - :arg sac: a symbolic assignment collection where temporary - expressions are stored. :returns: a list of :mod:`sympy` expressions representing the coefficients of the expansion. """ raise NotImplementedError - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): """ :return: a :mod:`sympy` expression corresponding to the evaluated expansion with the coefficients @@ -175,11 +173,11 @@ class ExpansionTermsWrangler(object): raise NotImplementedError def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, - rscale, sac=None): + rscale): raise NotImplementedError def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients, - rscale, sac=None): + rscale): raise NotImplementedError @memoize_method @@ -267,7 +265,7 @@ class FullExpansionTermsWrangler(ExpansionTermsWrangler): ExpansionTermsWrangler.get_full_coefficient_identifiers) def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, - rscale, sac=None): + rscale): return stored_kernel_derivatives get_stored_mpole_coefficients_from_full = ( @@ -309,7 +307,7 @@ def _spmv(spmat, x, sparse_vectors): # }}} -def _fast_spmv(reconstruct_matrix, vec, sac, transpose=False): +def _fast_spmv(reconstruct_matrix, vec, transpose=False): if not transpose: res = [0] * len(reconstruct_matrix) stored_idx = 0 @@ -320,8 +318,7 @@ def _fast_spmv(reconstruct_matrix, vec, sac, transpose=False): else: for k, v in deps: res[row] += res[k] * v - new_sym = sym.Symbol(sac.assign_unique("expr", res[row])) - res[row] = new_sym + res[row] = sym.UnevaluatedExpr(res[row]) return res else: res = [] @@ -330,7 +327,7 @@ def _fast_spmv(reconstruct_matrix, vec, sac, transpose=False): if len(deps) == 0: res.append(expr_all[row]) continue - new_sym = sym.Symbol(sac.assign_unique("expr", expr_all[row])) + new_sym = sym.UnevaluatedExpr(expr_all[row]) for k, v in deps: expr_all[k] += new_sym * v res.reverse() @@ -358,26 +355,26 @@ class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler): return self.stored_identifiers def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives, - rscale, sac=None): + rscale): coeff_matrix, reconstruct_matrix, use_reconstruct = \ self.get_coefficient_matrix(rscale) - if not use_reconstruct or sac is None: + if not use_reconstruct: return _spmv(coeff_matrix, stored_kernel_derivatives, sparse_vectors=False) - return _fast_spmv(reconstruct_matrix, stored_kernel_derivatives, sac) + return _fast_spmv(reconstruct_matrix, stored_kernel_derivatives) def get_stored_mpole_coefficients_from_full(self, full_mpole_coefficients, - rscale, sac=None): + rscale): # = M^T x, where M = coeff matrix coeff_matrix, reconstruct_matrix, use_reconstruct = \ self.get_coefficient_matrix(rscale) - if not use_reconstruct or sac is None: + if not use_reconstruct: result = [0] * len(self.stored_identifiers) for row, coeff in enumerate(full_mpole_coefficients): for col, val in coeff_matrix[row]: result[col] += coeff * val return result - return _fast_spmv(reconstruct_matrix, full_mpole_coefficients, sac, + return _fast_spmv(reconstruct_matrix, full_mpole_coefficients, transpose=True) @property diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py index 51687487..c32778ae 100644 --- a/sumpy/expansion/local.py +++ b/sumpy/expansion/local.py @@ -58,7 +58,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase): def get_coefficient_identifiers(self): return list(range(self.order+1)) - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): # no point in heeding rscale here--just ignore it if bvec is None: raise RuntimeError("cannot use line-Taylor expansions in a setting " @@ -99,7 +99,7 @@ class LineTaylorLocalExpansion(LocalExpansionBase): .subs("tau", 0) for i in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): # no point in heeding rscale here--just ignore it from pytools import factorial return sym.Add(*( @@ -116,7 +116,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): Coefficients represent derivative values of the kernel. """ - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): from sumpy.tools import MiDerivativeTaker ppkernel = self.kernel.postprocess_at_source( self.kernel.get_expression(avec), avec) @@ -126,7 +126,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase): taker.diff(mi) * rscale ** sum(mi) for mi in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): evaluated_coeffs = ( self.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored( coeffs, rscale)) @@ -357,7 +357,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): def get_coefficient_identifiers(self): return list(range(-self.order, self.order+1)) - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): if not self.use_rscale: rscale = 1 @@ -375,7 +375,7 @@ class _FourierBesselLocalExpansion(LocalExpansionBase): * sym.exp(sym.I * l * source_angle_rel_center), avec) for l in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): if not self.use_rscale: rscale = 1 diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py index 36b08e15..b7e2769b 100644 --- a/sumpy/expansion/multipole.py +++ b/sumpy/expansion/multipole.py @@ -56,7 +56,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): Coefficients represent the terms in front of the kernel derivatives. """ - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): from sumpy.kernel import DirectionalSourceDerivative kernel = self.kernel @@ -109,7 +109,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): for mi in self.get_full_coefficient_identifiers()] return ( self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full( - result, rscale, sac=sac)) + result, rscale)) def get_scaled_multipole(self, expr, bvec, rscale, nderivatives, nderivatives_for_scaling=None): @@ -127,7 +127,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase): else: return (rscale**nderivatives_for_scaling * expr) - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): if not self.use_rscale: rscale = 1 @@ -347,7 +347,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): def get_coefficient_identifiers(self): return list(range(-self.order, self.order+1)) - def coefficients_from_source(self, avec, bvec, rscale, sac=None): + def coefficients_from_source(self, avec, bvec, rscale): if not self.use_rscale: rscale = 1 @@ -367,7 +367,7 @@ class _HankelBased2DMultipoleExpansion(MultipoleExpansionBase): avec) for l in self.get_coefficient_identifiers()] - def evaluate(self, coeffs, bvec, rscale, sac=None): + def evaluate(self, coeffs, bvec, rscale): if not self.use_rscale: rscale = 1 diff --git a/sumpy/p2e.py b/sumpy/p2e.py index f5472fc4..2aa64804 100644 --- a/sumpy/p2e.py +++ b/sumpy/p2e.py @@ -88,7 +88,7 @@ class P2EBase(KernelCacheWrapper): coeff_names = [ sac.assign_unique("coeff%d" % i, coeff_i) for i, coeff_i in enumerate( - self.expansion.coefficients_from_source(avec, None, rscale, sac))] + self.expansion.coefficients_from_source(avec, None, rscale))] sac.run_global_cse() -- GitLab