[pyfr] 49/88: Replace @chop with a more powerful @clean decorator.

[Ghislain Vaillant]

This is an automated email from the git hooks/post-receive script.

ghisvail-guest pushed a commit to branch master
in repository pyfr.

commit f562ed7a7069e50f54b9afc274391db426d072ae
Author: Freddie Witherden <freddie at witherden.org>
Date:   Tue Jun 7 10:08:05 2016 -0700

    Replace @chop with a more powerful @clean decorator.
    
    Clean uses a similar approach to GiMMiK that also works to
    coalesce similar values in a matrix in addition to just truncating
    small values.
---
 pyfr/nputil.py | 36 ++++++++++++++++++++++++++++--------
 pyfr/polys.py  | 10 +++++-----
 pyfr/shapes.py | 12 ++++++------
 3 files changed, 39 insertions(+), 19 deletions(-)

diff --git a/pyfr/nputil.py b/pyfr/nputil.py
index 66e3146..64cbefd 100644
--- a/pyfr/nputil.py
+++ b/pyfr/nputil.py
@@ -21,16 +21,36 @@ def block_diag(arrs):
     return out
 
 
-def chop(fn):
-    @ft.wraps(fn)
-    def newfn(*args, **kwargs):
-        arr = fn(*args, **kwargs)
+def clean(origfn=None, tol=1e-10):
+    def cleanfn(fn):
+        @ft.wraps(fn)
+        def newfn(*args, **kwargs):
+            arr = fn(*args, **kwargs).copy()
 
-        # Determine a tolerance and flush
-        arr[abs(arr) < 10*np.finfo(arr.dtype).resolution] = 0
+            # Flush small elements to zero
+            arr[np.abs(arr) < tol] = 0
 
-        return arr
-    return newfn
+            # Coalesce similar elements
+            amfl = np.abs(arr.flat)
+            amix = np.argsort(amfl)
+
+            i, ix = 0, amix[0]
+            for j, jx in enumerate(amix[1:], start=1):
+                if amfl[jx] - amfl[ix] >= tol:
+                    if j - i > 1:
+                        amfl[amix[i:j]] = np.median(amfl[amix[i:j]])
+                    i, ix = j, jx
+
+            if i != j:
+                amfl[amix[i:]] = np.median(amfl[amix[i:]])
+
+            # Fix up the signs and assign
+            arr.flat = np.copysign(amfl, arr.flat)
+
+            return arr
+        return newfn
+
+    return cleanfn(origfn) if origfn else cleanfn
 
 
 _npeval_syms = {
diff --git a/pyfr/polys.py b/pyfr/polys.py
index 2f2b6b0..a08b8b6 100644
--- a/pyfr/polys.py
+++ b/pyfr/polys.py
@@ -5,7 +5,7 @@ from math import sqrt
 
 import numpy as np
 
-from pyfr.nputil import chop
+from pyfr.nputil import clean
 from pyfr.util import lazyprop, subclass_where
 
 
@@ -52,14 +52,14 @@ class BasePolyBasis(object):
         self.order = order
         self.pts = pts
 
-    @chop
+    @clean
     def ortho_basis_at(self, pts):
         if len(pts) and not isinstance(pts[0], Iterable):
             pts = [(p,) for p in pts]
 
         return np.array([self.ortho_basis_at_py(*p) for p in pts]).T
 
-    @chop
+    @clean
     def jac_ortho_basis_at(self, pts):
         if len(pts) and not isinstance(pts[0], Iterable):
             pts = [(p,) for p in pts]
@@ -68,11 +68,11 @@ class BasePolyBasis(object):
 
         return np.array(J).swapaxes(0, 2)
 
-    @chop
+    @clean
     def nodal_basis_at(self, epts):
         return np.linalg.solve(self.vdm, self.ortho_basis_at(epts)).T
 
-    @chop
+    @clean
     def jac_nodal_basis_at(self, epts):
         return np.linalg.solve(self.vdm, self.jac_ortho_basis_at(epts))
 
diff --git a/pyfr/shapes.py b/pyfr/shapes.py
index 0c8013b..77f8b48 100644
--- a/pyfr/shapes.py
+++ b/pyfr/shapes.py
@@ -6,7 +6,7 @@ import re
 
 import numpy as np
 
-from pyfr.nputil import block_diag, chop
+from pyfr.nputil import block_diag, clean
 from pyfr.polys import get_polybasis
 from pyfr.quadrules import get_quadrule
 from pyfr.util import lazyprop
@@ -17,7 +17,7 @@ def _proj_pts(projector, pts):
     return np.vstack(np.broadcast_arrays(*projector(*pts))).T
 
 
- at chop
+ at clean
 def _proj_l2(qrule, basis):
     return np.dot(basis.vdm.T, qrule.wts*basis.ortho_basis_at(qrule.pts))
 
@@ -78,7 +78,7 @@ class BaseShape(object):
         else:
             raise ValueError('Invalid number of shape points')
 
-    @chop
+    @clean
     def opmat(self, expr):
         if not re.match(r'[M0-9\-+*() ]+$', expr):
             raise ValueError('Invalid operator matrix expression')
@@ -142,7 +142,7 @@ class BaseShape(object):
         return block_diag([self.m9]*self.ndims)
 
     @lazyprop
-    @chop
+    @clean
     def m11(self):
         ub = self.ubasis
 
@@ -259,7 +259,7 @@ class BaseShape(object):
 
         return np.vstack(coeffs)
 
-    @chop
+    @clean
     def gbasis_at(self, pts):
         return np.dot(self.gbasis_coeffs, self.ubasis.ortho_basis_at(pts)).T
 
@@ -311,7 +311,7 @@ class BaseShape(object):
     @lazyprop
     def mpts(self):
         return self.std_ele(self.order)
-    
+
     @lazyprop
     def nmpts(self):
         return len(self.mpts)

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/pyfr.git