started refactoring and deletion of two useless lines

Author: guglielmopadula

pygem/bffd.py

@@ -56,23 +56,11 @@ class BFFD(CFFD):
     '''
     def __init__(self, n_control_points=None):
         super().__init__(n_control_points)
+
         def linfun(x):
-            return np.mean(x.reshape(-1,3),axis=0)
+            return np.mean(x.reshape(-1, 3), axis=0)
 
-        self.linconstraint=linfun
+        self.fun = linfun
 
     def __call__(self, src_pts):
         return super().__call__(src_pts)
-
-if __name__ == "__main__":
-        from pygem import BFFD
-        import numpy as np
-        bffd = BFFD()
-        bffd.read_parameters('tests/test_datasets/parameters_test_ffd_sphere.prm')
-        original_mesh_points = np.load('tests/test_datasets/meshpoints_sphere_orig.npy')
-        b=bffd.linconstraint(original_mesh_points)
-        bffd.valconstraint=b
-        bffd.indices=np.arange(np.prod(bffd.n_control_points)*3).tolist()
-        bffd.M=np.eye(len(bffd.indices))
-        new_mesh_points = bffd(original_mesh_points)
-        assert np.isclose(np.linalg.norm(bffd.linconstraint(new_mesh_points)-b),np.array([0.]))

pygem/cffd.py

@@ -34,8 +34,8 @@ class CFFD(FFD):
         y, normalized with the box length y.
     :cvar numpy.ndarray array_mu_z: collects the displacements (weights) along
         z, normalized with the box length z.
-    :cvar callable linconstraint: it defines the F of the constraint F(x)=c.
-    :cvar numpy.ndarray valconstraint: it defines the c of the constraint F(x)=c.
+    :cvar callable fun: it defines the F of the constraint F(x)=c.
+    :cvar numpy.ndarray fixval: it defines the c of the constraint F(x)=c.
     :cvar list indices: it defines the indices of the control points 
         that are moved to enforce the constraint. The control index is obtained by doing:
         all_indices=np.arange(n_x*n_y*n_z*3).reshape(n_x,n_y,n_z,3).tolist().
@@ -53,31 +53,31 @@ class CFFD(FFD):
         >>>     x=x.reshape(-1)
         >>>     return A@x
         >>> b=fun(original_mesh_points)
-        >>> cffd.linconstraint=fun
-        >>> cffd.valconstraint=b
+        >>> cffd.fun=fun
+        >>> cffd.fixval=b
         >>> cffd.indices=np.arange(np.prod(cffd.n_control_points)*3).tolist()
         >>> cffd.M=np.eye(len(cffd.indices))
         >>> new_mesh_points = cffd(original_mesh_points)
         >>> assert np.isclose(np.linalg.norm(fun(new_mesh_points)-b),np.array([0.]))
     """
-
     def __init__(self, n_control_points=None):
         super().__init__(n_control_points)
-        self.linconstraint = None 
-        self.valconstraint = None
+        self.fun = None
+        self.fixval = None
         self.indices = None
         self.M = None
 
     def __call__(self, src_pts):
-        saved_parameters=self._save_parameters()
+        saved_parameters = self._save_parameters()
         A, b = self._compute_linear_map(src_pts, saved_parameters.copy())
         d = A @ saved_parameters[self.indices] + b
-        deltax = np.linalg.inv(self.M) @ A.T @ np.linalg.inv((A @ np.linalg.inv(self.M)@ A.T)) @ (self.valconstraint - d)
+        deltax = np.linalg.inv(self.M) @ A.T @ np.linalg.inv(
+            (A @ np.linalg.inv(self.M) @ A.T)) @ (self.fixval - d)
         saved_parameters[self.indices] = saved_parameters[self.indices] + deltax
         self._load_parameters(saved_parameters)
         return self.ffd(src_pts)
 
-    def ffd(self,src_pts):
+    def ffd(self, src_pts):
         '''
         Performs Classic Free Form Deformation.
     
@@ -87,35 +87,35 @@ def ffd(self,src_pts):
         '''
         return super().__call__(src_pts)
 
-
     def _save_parameters(self):
         '''
         Saves the FFD control points in an array of shape [n_x,ny,nz,3].
 
         :return: the FFD control points in an array of shape [n_x,ny,nz,3].
         :rtype: numpy.ndarray
         '''
-        tmp = np.zeros([*self.n_control_points,3])
+        tmp = np.zeros([*self.n_control_points, 3])
         tmp[:, :, :, 0] = self.array_mu_x
         tmp[:, :, :, 1] = self.array_mu_y
         tmp[:, :, :, 2] = self.array_mu_z
         return tmp.reshape(-1)
-    
-    def _load_parameters(self,tmp):
+
+    def _load_parameters(self, tmp):
         '''
         Loads the FFD control points from an array of shape [n_x,ny,nz,3].
 
         :param np.ndarray tmp: the array of FFD control points.
         :rtype: None
         '''
-        tmp = tmp.reshape(*self.n_control_points,3)
+        tmp = tmp.reshape(*self.n_control_points, 3)
         self.array_mu_x = tmp[:, :, :, 0]
         self.array_mu_y = tmp[:, :, :, 1]
         self.array_mu_z = tmp[:, :, :, 2]
 
 
 # I see that a similar function already exists in pygem.utils, but it does not work for inputs and outputs of different dimensions
-    def _compute_linear_map(self,src_pts,saved_parameters):
+
+    def _compute_linear_map(self, src_pts, saved_parameters):
         '''
         Computes the coefficient and the intercept of the linear map from the control points to the output.
         
@@ -124,21 +124,25 @@ def _compute_linear_map(self,src_pts,saved_parameters):
         :return: a tuple containing the coefficient and the intercept.
         :rtype: tuple(np.ndarray,np.ndarray)
         '''
-        saved_parameters_bak=saved_parameters.copy() #saving ffd parameters
-        n_indices=len(self.indices)
-        inputs=np.zeros([n_indices+1,n_indices+1])
-        outputs=np.zeros([n_indices+1,self.valconstraint.shape[0]])
+        n_indices = len(self.indices)
+        inputs = np.zeros([n_indices + 1, n_indices + 1])
+        outputs = np.zeros([n_indices + 1, self.fixval.shape[0]])
         np.random.seed(0)
-        for i in range(n_indices+1): ##now we generate the interpolation points
-            tmp=np.random.rand(1,n_indices)
-            tmp=tmp.reshape(1,-1)
-            inputs[i]=np.hstack([tmp, np.ones((tmp.shape[0], 1))]) #dependent variable
-            saved_parameters[self.indices]=tmp 
-            self._load_parameters(saved_parameters) #loading the depent variable as a control point
-            def_pts=super().__call__(src_pts) #computing the deformation with the dependent variable
-            outputs[i]=self.linconstraint(def_pts) #computing the independent variable
-        sol=np.linalg.lstsq(inputs,outputs,rcond=None) #computation of the linear map
-        A=sol[0].T[:,:-1] #coefficient
-        b=sol[0].T[:,-1] #intercept
-        self._load_parameters(saved_parameters_bak) #restoring the original FFD parameters
-        return A,b
+        for i in range(n_indices +
+                       1):  ##now we generate the interpolation points
+            tmp = np.random.rand(1, n_indices)
+            tmp = tmp.reshape(1, -1)
+            inputs[i] = np.hstack([tmp, np.ones(
+                (tmp.shape[0], 1))])  #dependent variable
+            saved_parameters[self.indices] = tmp
+            self._load_parameters(
+                saved_parameters
+            )  #loading the depent variable as a control point
+            def_pts = super().__call__(
+                src_pts)  #computing the deformation with the dependent variable
+            outputs[i] = self.fun(def_pts)  #computing the independent variable
+        sol = np.linalg.lstsq(inputs, outputs,
+                              rcond=None)  #computation of the linear map
+        A = sol[0].T[:, :-1]  #coefficient
+        b = sol[0].T[:, -1]  #intercept
+        return A, b

pygem/vffd.py

@@ -12,6 +12,7 @@
 import numpy as np
 from copy import deepcopy
 
+
 class VFFD(CFFD):
     '''
     Class that handles the Volumetric Free Form Deformation on the mesh points.
@@ -58,32 +59,33 @@ class VFFD(CFFD):
         >>> new_mesh_points = vffd(original_mesh_points)
         >>> assert np.isclose(np.linalg.norm(vffd.linconstraint(new_mesh_points)-b),np.array([0.]))
     '''
-
-    def __init__(self,triangles ,n_control_points=None):
+    def __init__(self, triangles, n_control_points=None):
         super().__init__(n_control_points)
-        self.triangles=triangles
-        self.vweight=[1/3,1/3,1/3]
+        self.triangles = triangles
+        self.vweight = [1 / 3, 1 / 3, 1 / 3]
+
         def volume(x):
-            x=x.reshape(-1,3)
-            mesh=x[self.triangles]
-            return np.sum(np.linalg.det(mesh)) 
-        self.linconstraint=volume
+            x = x.reshape(-1, 3)
+            mesh = x[self.triangles]
+            return np.sum(np.linalg.det(mesh))
 
+        self.fun = volume
 
-    def __call__(self,src_pts):
-        self.vweight=np.abs(self.vweight)/np.sum(np.abs(self.vweight))
-        indices_bak=deepcopy(self.indices)
-        self.indices=np.array(self.indices)
-        indices_x=self.indices[self.indices%3==0].tolist()
-        indices_y=self.indices[self.indices%3==1].tolist()
-        indices_z=self.indices[self.indices%3==2].tolist()
-        indexes=[indices_x,indices_y,indices_z]
-        diffvolume=self.valconstraint-self.linconstraint(self.ffd(src_pts))
+    def __call__(self, src_pts):
+        self.vweight = np.abs(self.vweight) / np.sum(np.abs(self.vweight))
+        indices_bak = deepcopy(self.indices)
+        self.indices = np.array(self.indices)
+        indices_x = self.indices[self.indices % 3 == 0].tolist()
+        indices_y = self.indices[self.indices % 3 == 1].tolist()
+        indices_z = self.indices[self.indices % 3 == 2].tolist()
+        indexes = [indices_x, indices_y, indices_z]
+        diffvolume = self.fixval - self.fun(self.ffd(src_pts))
         for i in range(3):
-            self.indices=indexes[i]
-            self.M=np.eye(len(self.indices))
-            self.valconstraint=self.linconstraint(self.ffd(src_pts))+self.vweight[i]*(diffvolume)
-            _=super().__call__(src_pts)
-        tmp=super().__call__(src_pts)
-        self.indices=indices_bak
-        return tmp
+            self.indices = indexes[i]
+            self.M = np.eye(len(self.indices))
+            self.fixval = self.fun(
+                self.ffd(src_pts)) + self.vweight[i] * (diffvolume)
+            _ = super().__call__(src_pts)
+        tmp = super().__call__(src_pts)
+        self.indices = indices_bak
+        return tmp

tests/test_bffd.py

@@ -9,10 +9,11 @@ class TestBFFD(TestCase):
     def test_nothing_happens(self):
         np.random.seed(0)
         cffd = BFFD()
-        original_mesh_points = np.load("tests/test_datasets/meshpoints_sphere_orig.npy")
+        original_mesh_points = np.load(
+            "tests/test_datasets/meshpoints_sphere_orig.npy")
         A = np.random.rand(3, original_mesh_points.reshape(-1).shape[0])
-        b = cffd.linconstraint(original_mesh_points)
-        cffd.valconstraint = b
+        b = cffd.fun(original_mesh_points)
+        cffd.fixval = b
         cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
         cffd.M = np.eye(len(cffd.indices))
         new_mesh_points = cffd(original_mesh_points)
@@ -21,13 +22,14 @@ def test_nothing_happens(self):
     def test_constraint(self):
         np.random.seed(0)
         cffd = BFFD()
-        cffd.read_parameters("tests/test_datasets/parameters_test_ffd_sphere.prm")
-        original_mesh_points = np.load("tests/test_datasets/meshpoints_sphere_orig.npy")
-        b = cffd.linconstraint(original_mesh_points)
-        cffd.valconstraint = b
+        cffd.read_parameters(
+            "tests/test_datasets/parameters_test_ffd_sphere.prm")
+        original_mesh_points = np.load(
+            "tests/test_datasets/meshpoints_sphere_orig.npy")
+        b = cffd.fun(original_mesh_points)
+        cffd.fixval = b
         cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
         cffd.M = np.eye(len(cffd.indices))
         new_mesh_points = cffd(original_mesh_points)
-        assert np.isclose(
-            np.linalg.norm(cffd.linconstraint(new_mesh_points) - b), np.array([0.0])
-        )
+        assert np.isclose(np.linalg.norm(cffd.fun(new_mesh_points) - b),
+                          np.array([0.0]))

tests/test_cffd.py

@@ -9,16 +9,17 @@ class TestCFFD(TestCase):
     def test_nothing_happens(self):
         np.random.seed(0)
         cffd = CFFD()
-        original_mesh_points = np.load("tests/test_datasets/meshpoints_sphere_orig.npy")
+        original_mesh_points = np.load(
+            "tests/test_datasets/meshpoints_sphere_orig.npy")
         A = np.random.rand(3, original_mesh_points.reshape(-1).shape[0])
 
         def fun(x):
             x = x.reshape(-1)
             return A @ x
 
         b = fun(original_mesh_points)
-        cffd.linconstraint = fun
-        cffd.valconstraint = b
+        cffd.fun = fun
+        cffd.fixval = b
         cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
         cffd.M = np.eye(len(cffd.indices))
         new_mesh_points = cffd(original_mesh_points)
@@ -27,34 +28,38 @@ def fun(x):
     def test_constraint(self):
         np.random.seed(0)
         cffd = CFFD()
-        cffd.read_parameters("tests/test_datasets/parameters_test_ffd_sphere.prm")
-        original_mesh_points = np.load("tests/test_datasets/meshpoints_sphere_orig.npy")
+        cffd.read_parameters(
+            "tests/test_datasets/parameters_test_ffd_sphere.prm")
+        original_mesh_points = np.load(
+            "tests/test_datasets/meshpoints_sphere_orig.npy")
         A = np.random.rand(3, original_mesh_points.reshape(-1).shape[0])
 
         def fun(x):
             x = x.reshape(-1)
             return A @ x
 
         b = fun(original_mesh_points)
-        cffd.linconstraint = fun
-        cffd.valconstraint = b
+        cffd.fun = fun
+        cffd.fixval = b
         cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
         cffd.M = np.eye(len(cffd.indices))
         new_mesh_points = cffd(original_mesh_points)
-        assert np.isclose(np.linalg.norm(fun(new_mesh_points) - b), np.array([0.0]))
+        assert np.isclose(np.linalg.norm(fun(new_mesh_points) - b),
+                          np.array([0.0]))
 
     def test_interpolation(self):
         cffd = CFFD()
-        original_mesh_points = np.load("tests/test_datasets/meshpoints_sphere_orig.npy")
+        original_mesh_points = np.load(
+            "tests/test_datasets/meshpoints_sphere_orig.npy")
         A = np.random.rand(3, original_mesh_points.reshape(-1).shape[0])
 
         def fun(x):
             x = x.reshape(-1)
             return A @ x
 
         b = fun(original_mesh_points)
-        cffd.linconstraint = fun
-        cffd.valconstraint = b
+        cffd.fixval = b
+        cffd.fun = fun
         cffd.indices = np.arange(np.prod(cffd.n_control_points) * 3).tolist()
         cffd.M = np.eye(len(cffd.indices))
         save_par = cffd._save_parameters()
@@ -64,8 +69,7 @@ def fun(x):
             save_par[cffd.indices] = tmp
             cffd._load_parameters(tmp)
             assert np.allclose(
-                np.linalg.norm(
-                    C @ tmp + d - cffd.linconstraint(cffd.ffd(original_mesh_points))
-                ),
+                np.linalg.norm(C @ tmp + d -
+                               cffd.fun(cffd.ffd(original_mesh_points))),
                 0.0,
             )