Bugfixes and refactoring

epde/evaluators.py

@@ -84,17 +84,8 @@ def __call__(self, factor, structural: bool = False, func_args: List[Union[torch
                 self.indexes_vect[tensor_idx] = tuple([subarg[tensor_idx]
                                                        for subarg in func_args])
         value = grid_function(self.indexes_vect)
-        if len(global_var.grid_cache.initial_shape) > 1:
-            value = value.reshape(*global_var.grid_cache.initial_shape)
-            if isinstance(global_var.grid_cache.boundary_width, int):
-                for dim in range(value.ndim):
-                    value[dim] = value[global_var.grid_cache.boundary_width:-global_var.grid_cache.boundary_width]
-            elif isinstance(global_var.grid_cache.boundary_width, (list, tuple)):
-                for dim in range(value.ndim):
-                    value[dim] = value[global_var.grid_cache.boundary_width[dim]:-global_var.grid_cache.boundary_width[dim]]
-            value = value.reshape(-1)
-        else:
-            value = value[global_var.grid_cache.boundary_width:-global_var.grid_cache.boundary_width]
+        value = value[global_var.grid_cache.g_func != 0]
+        value = value.reshape(-1)
         return value
 
 

epde/operators/common/sparsity.py

@@ -52,13 +52,11 @@ def fit(self, X, y, sample_weights):
         norm_sq_features = np.sum(X ** 2, axis=0)
         residual = y - (X @ self.coef_ + self.intercept_)
 
-        # 2. Coordinate Descent Loop
-        for iteration in range(self.max_iter):
-            max_change = self.tol
-
-            if all(self.coef_ == 0):
-                break
+        iteration = 0
+        max_change = np.inf
 
+        # 2. Coordinate Descent Loop
+        while iteration < self.max_iter and not all(self.coef_ == 0):
             # Sort features by instability (highest CV first)
             for j in np.argsort(cv)[::-1]:
                 old_coef = self.coef_[j]
@@ -89,6 +87,7 @@ def fit(self, X, y, sample_weights):
                     self.coef_ = np.array([next(new_coef) if _ else 0 for _ in self.coef_ != 0])
                     self.intercept_ = weights.mean(axis=0)[-1]
                     residual = y - (X @ self.coef_ + self.intercept_)
+                    iteration = 0
                     break
 
                 residual -= (new_coef - old_coef) * X[:, j]
@@ -98,6 +97,7 @@ def fit(self, X, y, sample_weights):
 
             if max_change < self.tol:
                 break
+            iteration += 1
         # print(iteration)
         return self
 

epde/supplementary.py

@@ -455,6 +455,8 @@ def calculate_weights(X, y, sample_weights, grid_shape):
         XTW = X_batch.transpose(0, 2, 1) * weights_batch.transpose(0, 2, 1)
         XTWX = XTW @ X_batch
         XTWy = XTW @ y_batch[..., None]
+        ridge = 1e-6 * np.eye(n_features_aug)
+        XTWX += ridge
 
         # 2. Solve (Fast CPU Vectorized Solver)
         # np.linalg.solve supports batch dimensions!