From 564f08f60903f3c346832467e718b842bf7ce63a Mon Sep 17 00:00:00 2001
From: tutou2356 <142212931+tutou2356@users.noreply.github.com>
Date: Thu, 7 Aug 2025 00:57:17 +0800
Subject: [PATCH 01/11] feat(grid): Implement adaptive grid refinement Adds the
 AdaptiveUniformGrid wrapper class to provide recursive subdivision of a
 uniform grid.

---
 src/grid/cubic.py            | 176 +++++++++++++++++++++++++++++++++++
 src/grid/tests/test_cubic.py |  55 +++++++++++
 2 files changed, 231 insertions(+)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index 1c2998f6d..a5b31cf6e 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -26,6 +26,10 @@
 
 from grid.basegrid import Grid, OneDGrid
 
+from collections import deque
+from typing import Optional, Callable
+from scipy.spatial import cKDTree
+
 
 class _HyperRectangleGrid(Grid):
     def __init__(self, points, weights, shape):
@@ -1013,3 +1017,175 @@ def generate_cube(self, fname, data, atcoords, atnums, pseudo_numbers=None):
                 row_data = data.flat[i : i + num_chunks]
                 f.write((row_data.size * " {:12.5E}").format(*row_data))
                 f.write("\n")
+
+
+class AdaptiveUniformGrid:
+    """This is a wrapper class that provides adaptive refinement for a UniformGrid instance.
+
+    This class takes a UniformGrid object and applies a recursive subdivision
+    algorithm to generate a new, non-uniform grid with points concentrated in
+    regions of high function error, leading to more efficient and accurate integration.
+
+    The main entry point is the `refinement` method."""
+
+    def __init__(self, uniform_grid: UniformGrid):
+        """initialization
+        Parameters
+        ----------
+        uniform_grid : UniformGrid
+            The coarse, uniform grid that will serve as the starting point for refinement."""
+        if not isinstance(uniform_grid, UniformGrid):
+            raise ValueError("The input grid should be a UniformGrid instance.")
+        self.grid = uniform_grid
+        self.ndim = uniform_grid.ndim
+
+    def _estimate_error(self, points: np.ndarray, evaluated_points: dict) -> np.ndarray:
+        """Estimates the error for each point based on the local gradient."""
+        errors = np.zeros(len(points))
+        if len(points) <= 1:
+            return errors
+
+        tree = cKDTree(points)
+        for i, point in enumerate(points):
+            distances, indices = tree.query(point, k=2)
+            closest_neighbor_points = points[indices[1]]
+            local_spacing = distances[1]
+
+            point_val = evaluated_points[tuple(point)]
+            neighbor_val = evaluated_points[tuple(closest_neighbor_points)]
+            grad_mag = abs(point_val - neighbor_val) / local_spacing
+
+            # Error : |grad(f)| * h
+            errors[i] = grad_mag * local_spacing
+        return errors
+
+    def _find_neighbors(self, point: np.ndarray, spacing: float) -> list:
+        """Notes on Neighbor Finding:
+        The standard method of finding neighbors by converting point indices to integer coordinates (i, j, k) is not used here.
+        This is because the adaptive refinement process adds new points that do not lie on the original structured grid.
+        For these new points, the index-based mapping will fail.
+        Therefore, I use real-world (x, y, z) coordinates and still implement the x ± (spacing/2) * a_i to find neighbors.
+        """
+        neighbors = []
+        for axis in self.grid.axes:
+            axis_direction = axis / np.linalg.norm(axis)
+            neighbors.append(point + spacing * axis_direction)
+            neighbors.append(point - spacing * axis_direction)
+        return neighbors
+
+    def refinement(
+        self, func: Callable, tolerance: float = 1e-4, min_spacing: Optional[float] = None
+    ) -> dict:
+        """Drives the adaptive refinement process and returns the results.
+
+        This method starts with the initial uniform grid, refines it according to
+        the function `func`, and returns the final results without modifying the
+        original grid object.
+
+        Parameters
+        ----------
+        func : Callable
+            The function to be integrated.
+        tolerance : float, optional
+            The error tolerance for a local point.
+        min_spacing : float, optional
+            The minimum allowed spacing for subdivision.
+
+        Returns
+        -------
+        dict
+            A dictionary containing the final integral value, the refined grid object,
+            and other statistics."""
+
+        # Initialization
+        initial_points = self.grid.points.copy()
+        initial_weights = self.grid.weights.copy()
+        initial_avg_spacing = np.mean([np.linalg.norm(axis) for axis in self.grid.axes])
+        if min_spacing is None:
+            min_spacing = initial_avg_spacing / 100
+
+        evaluated_points = {}
+        initial_values = func(initial_points)
+        for i, p in enumerate(initial_points):
+            evaluated_points[tuple(p)] = initial_values[i]
+
+        original_total_volume = np.sum(initial_weights)
+
+        # Refinement process
+        errors = self._estimate_error(initial_points, evaluated_points)
+        high_error_indices = np.where(errors > tolerance)[0]
+
+        if high_error_indices.size == 0:
+            final_integral = self.grid.integrate(initial_values)
+            return {
+                "integral": final_integral,
+                "final_grid": self.grid,
+                "num_points": len(initial_points),
+                "num_evaluations": len(evaluated_points),
+            }
+
+        final_points = []
+        unnormalized_weights = []
+
+        keep_mask = np.ones(len(initial_points), dtype=bool)
+        keep_mask[high_error_indices] = False
+
+        final_points.extend(list(initial_points[keep_mask]))
+        unnormalized_weights.extend(list(initial_weights[keep_mask]))
+
+        refinement_queue = deque(
+            [(initial_points[idx], initial_avg_spacing) for idx in high_error_indices]
+        )
+        processed_points_set = {tuple(p) for p in initial_points}
+
+        while refinement_queue:
+            point, spacing = refinement_queue.popleft()
+            half_spacing = spacing / 2
+            if half_spacing < min_spacing:
+                final_points.append(point)
+                unnormalized_weights.append(spacing**self.ndim)
+                continue
+
+            neighbors = self._find_neighbors(point, half_spacing)
+
+            new_points_to_eval = [p for p in neighbors if tuple(p) not in evaluated_points]
+            if new_points_to_eval:
+                new_values = func(np.array(new_points_to_eval))
+                for new_point, new_value in zip(new_points_to_eval, new_values):
+                    evaluated_points[tuple(new_point)] = new_value
+
+            local_error = 0
+            for neighbor in neighbors:
+                grad_mag = (
+                    abs(evaluated_points[tuple(point)] - evaluated_points[tuple(neighbor)])
+                    / half_spacing
+                )
+                local_error = max(local_error, grad_mag * half_spacing)
+
+            if local_error < tolerance:
+                final_points.append(point)
+                unnormalized_weights.append(spacing**self.ndim)
+            else:
+                refinement_queue.append((point, half_spacing))
+                for neighbor in neighbors:
+                    child_tuple = tuple(neighbor)
+                    if child_tuple not in processed_points_set:
+                        refinement_queue.append((neighbor, half_spacing))
+                        processed_points_set.add(child_tuple)
+
+        # Finalization
+        final_points = np.array(final_points)
+        final_weights = np.array(unnormalized_weights)
+        current_total_volume = np.sum(final_weights)
+        if current_total_volume > 0:
+            final_weights *= original_total_volume / current_total_volume
+        final_grid = Grid(final_points, final_weights)
+        final_values = np.array([evaluated_points[tuple(p)] for p in final_points])
+        final_integral = final_grid.integrate(final_values)
+
+        return {
+            "integral": final_integral,
+            "final_grid": final_grid,
+            "num_points": len(final_points),
+            "num_evaluations": len(evaluated_points),
+        }
diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py
index 42f3f30b9..8a68b38fb 100644
--- a/src/grid/tests/test_cubic.py
+++ b/src/grid/tests/test_cubic.py
@@ -28,6 +28,10 @@
 from grid.cubic import Tensor1DGrids, UniformGrid, _HyperRectangleGrid
 from grid.onedgrid import GaussLaguerre, MidPoint
 
+import pytest
+import copy
+from grid.cubic import AdaptiveUniformGrid, Grid
+
 
 class TestHyperRectangleGrid(TestCase):
     r"""Test HyperRectangleGrid class."""
@@ -1050,3 +1054,54 @@ def test_uniformgrid_points_without_rotate(self):
             ]
         )
         assert_allclose(grid.points, expected, rtol=1.0e-7, atol=1.0e-7)
+
+
+class TestAdaptiveUniformGrid:
+    """Tests for the new AdaptiveUniformGrid wrapper class."""
+
+    def setup_method(self):
+        """Set up a coarse 3D UniformGrid before each test."""
+        origin = np.array([-2.0, -2.0, -2.0])
+        axes = np.diag([1.0, 1.0, 1.0])
+        shape = np.array([5, 5, 5])  # 125 initial points
+        self.uniform_grid = UniformGrid(origin, axes, shape, weight="Rectangle")
+
+        # Define a 3D function with a sharp peak for testing
+        self.sharp_3d_gaussian = lambda points: np.exp(-20 * np.sum(points**2, axis=1))
+
+    def test_refinement_improves_accuracy(self):
+        """
+        Tests that the adaptively refined grid yields a more accurate integral.
+        """
+        # Arrange: Set up the test conditions and known values.
+        analytical_integral = (np.pi / 20) ** 1.5
+        adaptive_grid = AdaptiveUniformGrid(self.uniform_grid)
+
+        initial_values = self.sharp_3d_gaussian(self.uniform_grid.points)
+        initial_integral = self.uniform_grid.integrate(initial_values)
+        initial_error = abs(initial_integral - analytical_integral)
+
+        result = adaptive_grid.refinement(func=self.sharp_3d_gaussian, tolerance=1e-5)
+
+        refined_integral = result["integral"]
+        refined_error = abs(refined_integral - analytical_integral)
+
+        initial_num_points = self.uniform_grid.size
+        refined_num_points = result["num_points"]
+
+        error_reduction_factor = (
+            initial_error / refined_error if refined_error > 0 else float("inf")
+        )
+
+        print("\n--- Adaptive Accuracy Test Summary ---")
+        print(f"{'Metric':<25} | {'Initial':<25} | {'Refined':<25}")
+        print("-" * 80)
+        print(f"{'Number of Points':<25} | {initial_num_points:<25} | {refined_num_points:<25}")
+        print(f"{'Integration Error':<25} | {initial_error:<25.10e} | {refined_error:<25.10e}")
+        print("-" * 80)
+        print(f"Error Reduction Factor: {error_reduction_factor:.2f}x")
+        print(f"Point Count Increase: {refined_num_points - initial_num_points}")
+        print("--------------------------------------")
+
+        assert refined_error < initial_error
+        assert result["num_points"] > self.uniform_grid.size

From cc3a2d3c7aee33ed7679847d6b946a5353912aac Mon Sep 17 00:00:00 2001
From: tutou2356 <142212931+tutou2356@users.noreply.github.com>
Date: Thu, 7 Aug 2025 13:34:17 +0800
Subject: [PATCH 02/11] Fix: Adjusted the error calculation function and
 weighting algorithm, added test cases

- Revised the `_estimate_error` function
- Adjusted the weight calculation method to more efficiently update based on new errors
- Added test cases for other functions
---
 src/grid/cubic.py            | 47 +++++++++++++-------
 src/grid/tests/test_cubic.py | 86 +++++++++++++++++++++++++++---------
 2 files changed, 98 insertions(+), 35 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index a5b31cf6e..180b03ac7 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1046,17 +1046,30 @@ def _estimate_error(self, points: np.ndarray, evaluated_points: dict) -> np.ndar
             return errors
 
         tree = cKDTree(points)
+        k_neighbors = self.ndim * 2 + 1
         for i, point in enumerate(points):
-            distances, indices = tree.query(point, k=2)
-            closest_neighbor_points = points[indices[1]]
-            local_spacing = distances[1]
+            k = min(k_neighbors, len(points))
+            distances, indices = tree.query(point, k=k)
+            neighbor_points = points[indices[1:]]
+            neighbor_distances = distances[1:]
+            if len(neighbor_points) == 0:
+                continue
 
             point_val = evaluated_points[tuple(point)]
-            neighbor_val = evaluated_points[tuple(closest_neighbor_points)]
-            grad_mag = abs(point_val - neighbor_val) / local_spacing
 
-            # Error : |grad(f)| * h
-            errors[i] = grad_mag * local_spacing
+            max_grad_mag = 0.0
+            avg_local_spacing = np.mean(neighbor_distances)
+            for neighbor, dist in zip(neighbor_points, neighbor_distances):
+                if dist == 0:
+                    continue
+
+                neighbor_val = evaluated_points[tuple(neighbor)]
+                grad_mag = abs(point_val - neighbor_val) / dist
+                if grad_mag > max_grad_mag:
+                    max_grad_mag = grad_mag
+
+            errors[i] = max_grad_mag * avg_local_spacing
+
         return errors
 
     def _find_neighbors(self, point: np.ndarray, spacing: float) -> list:
@@ -1130,8 +1143,10 @@ def refinement(
         keep_mask = np.ones(len(initial_points), dtype=bool)
         keep_mask[high_error_indices] = False
 
-        final_points.extend(list(initial_points[keep_mask]))
-        unnormalized_weights.extend(list(initial_weights[keep_mask]))
+        retained_points = initial_points[keep_mask]
+        final_points.extend(list(retained_points))
+        retained_weight = initial_avg_spacing ** self.ndim
+        unnormalized_weights.extend([retained_weight] * len(retained_points))
 
         refinement_queue = deque(
             [(initial_points[idx], initial_avg_spacing) for idx in high_error_indices]
@@ -1154,13 +1169,15 @@ def refinement(
                 for new_point, new_value in zip(new_points_to_eval, new_values):
                     evaluated_points[tuple(new_point)] = new_value
 
-            local_error = 0
+            point_val = evaluated_points[tuple(point)]
+            max_grad_mag = 0
             for neighbor in neighbors:
-                grad_mag = (
-                    abs(evaluated_points[tuple(point)] - evaluated_points[tuple(neighbor)])
-                    / half_spacing
-                )
-                local_error = max(local_error, grad_mag * half_spacing)
+                neighbor_val = evaluated_points[tuple(neighbor)]
+                grad_mag = abs(point_val - neighbor_val) / half_spacing
+                if grad_mag > max_grad_mag:
+                    max_grad_mag = grad_mag
+
+            local_error = max_grad_mag * half_spacing
 
             if local_error < tolerance:
                 final_points.append(point)
diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py
index 8a68b38fb..5e6ff51d7 100644
--- a/src/grid/tests/test_cubic.py
+++ b/src/grid/tests/test_cubic.py
@@ -1059,41 +1059,87 @@ def test_uniformgrid_points_without_rotate(self):
 class TestAdaptiveUniformGrid:
     """Tests for the new AdaptiveUniformGrid wrapper class."""
 
-    def setup_method(self):
-        """Set up a coarse 3D UniformGrid before each test."""
-        origin = np.array([-2.0, -2.0, -2.0])
-        axes = np.diag([1.0, 1.0, 1.0])
-        shape = np.array([5, 5, 5])  # 125 initial points
-        self.uniform_grid = UniformGrid(origin, axes, shape, weight="Rectangle")
-
-        # Define a 3D function with a sharp peak for testing
-        self.sharp_3d_gaussian = lambda points: np.exp(-20 * np.sum(points**2, axis=1))
-
-    def test_refinement_improves_accuracy(self):
+    case_3d_single_peak = {
+        "id": "3D_Single_Peak",
+        "grid_setup": {
+            "origin": np.array([-2.0, -2.0, -2.0]),
+            "axes": np.diag([1.0, 1.0, 1.0]),
+            "shape": np.array([5, 5, 5]),
+        },
+        "func": lambda points: np.exp(-20 * np.sum(points ** 2, axis=1)),
+        # Analytical integral of exp(-a*r^2) in 3D is (pi/a)^(3/2)
+        "analytical_integral": (np.pi / 20) ** 1.5,
+        "tolerance": 1e-5
+    }
+
+    case_2d_single_peak = {
+        "id": "2D_Single_Peak_Centered",
+        "grid_setup": {
+            "origin": np.array([-4.0, -4.0]),
+            "axes": np.diag([1.0, 1.0]),
+            "shape": np.array([9, 9]),
+        },
+        "func": lambda points: np.exp(-50 * np.sum(points ** 2, axis=1)),
+        # Analytical integral of exp(-a*r^2) in 2D is pi/a
+        "analytical_integral": np.pi / 50,
+        "tolerance": 1e-4
+    }
+
+    case_2d_multi_peak = {
+        "id": "2D_Multi_Peak_Centered",
+        "grid_setup": {
+            "origin": np.array([-4.0, -4.0]),
+            "axes": np.diag([1.0, 1.0]),
+            "shape": np.array([9, 9]),
+        },
+        "func": lambda points: (
+            np.exp(-30 * np.sum((points - np.array([1.0, 1.0])) ** 2, axis=1)) +
+            np.exp(-30 * np.sum((points - np.array([-1.0, -1.0])) ** 2, axis=1))
+        ),
+        # Integral is the sum of two identical Gaussians
+        "analytical_integral": 2 * (np.pi / 30),
+        "tolerance": 1e-3
+    }
+
+    @pytest.mark.parametrize(
+        "test_case",
+        [case_3d_single_peak, case_2d_single_peak, case_2d_multi_peak],
+        ids=lambda tc: tc["id"]  # Use the 'id' field for clear test names in the report
+    )
+
+    def test_refinement_improves_accuracy(self, test_case):
         """
         Tests that the adaptively refined grid yields a more accurate integral.
         """
         # Arrange: Set up the test conditions and known values.
-        analytical_integral = (np.pi / 20) ** 1.5
-        adaptive_grid = AdaptiveUniformGrid(self.uniform_grid)
 
-        initial_values = self.sharp_3d_gaussian(self.uniform_grid.points)
-        initial_integral = self.uniform_grid.integrate(initial_values)
+        grid_setup = test_case["grid_setup"]
+        uniform_grid = UniformGrid(
+            grid_setup["origin"], grid_setup["axes"], grid_setup["shape"], weight="Rectangle"
+        )
+        test_func = test_case["func"]
+        analytical_integral = test_case["analytical_integral"]
+        tolerance = test_case["tolerance"]
+
+
+        adaptive_grid = AdaptiveUniformGrid(uniform_grid)
+
+        initial_values = test_func(uniform_grid.points)
+        initial_integral = uniform_grid.integrate(initial_values)
         initial_error = abs(initial_integral - analytical_integral)
 
-        result = adaptive_grid.refinement(func=self.sharp_3d_gaussian, tolerance=1e-5)
+        result = adaptive_grid.refinement(func=test_func, tolerance=tolerance)
 
         refined_integral = result["integral"]
         refined_error = abs(refined_integral - analytical_integral)
 
-        initial_num_points = self.uniform_grid.size
+        initial_num_points = uniform_grid.size
         refined_num_points = result["num_points"]
-
         error_reduction_factor = (
             initial_error / refined_error if refined_error > 0 else float("inf")
         )
 
-        print("\n--- Adaptive Accuracy Test Summary ---")
+        print(f"\n--- Test Scenario: {test_case['id']} ---")
         print(f"{'Metric':<25} | {'Initial':<25} | {'Refined':<25}")
         print("-" * 80)
         print(f"{'Number of Points':<25} | {initial_num_points:<25} | {refined_num_points:<25}")
@@ -1104,4 +1150,4 @@ def test_refinement_improves_accuracy(self):
         print("--------------------------------------")
 
         assert refined_error < initial_error
-        assert result["num_points"] > self.uniform_grid.size
+        assert result["num_points"] > uniform_grid.size

From 7d45a6c9cb35d20cf5bb10a9532b5c632ddb1ed3 Mon Sep 17 00:00:00 2001
From: tutou2356 <142212931+tutou2356@users.noreply.github.com>
Date: Tue, 19 Aug 2025 23:46:40 +0800
Subject: [PATCH 03/11] Revise error estimation and weight calculation for
 adaptive grid refinement

- Use finite difference gradient for error estimation
- Include diagonal elements in 3^D subdivision
- Update weight calculation method for volume conservation
---
 src/grid/cubic.py            | 285 ++++++++++++++++++++---------------
 src/grid/tests/test_cubic.py | 170 +++++++++++++--------
 2 files changed, 270 insertions(+), 185 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index 180b03ac7..9dfa78087 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1029,72 +1029,116 @@ class AdaptiveUniformGrid:
     The main entry point is the `refinement` method."""
 
     def __init__(self, uniform_grid: UniformGrid):
-        """initialization
+        """Initialization.
+
         Parameters
         ----------
         uniform_grid : UniformGrid
-            The coarse, uniform grid that will serve as the starting point for refinement."""
+            The coarse, uniform grid that will serve as the starting point for refinement.
+        """
         if not isinstance(uniform_grid, UniformGrid):
             raise ValueError("The input grid should be a UniformGrid instance.")
         self.grid = uniform_grid
         self.ndim = uniform_grid.ndim
+        self.axis_spacings = np.array([np.linalg.norm(axis) for axis in self.grid.axes])
 
-    def _estimate_error(self, points: np.ndarray, evaluated_points: dict) -> np.ndarray:
-        """Estimates the error for each point based on the local gradient."""
-        errors = np.zeros(len(points))
-        if len(points) <= 1:
-            return errors
-
-        tree = cKDTree(points)
-        k_neighbors = self.ndim * 2 + 1
-        for i, point in enumerate(points):
-            k = min(k_neighbors, len(points))
-            distances, indices = tree.query(point, k=k)
-            neighbor_points = points[indices[1:]]
-            neighbor_distances = distances[1:]
-            if len(neighbor_points) == 0:
-                continue
+    def _get_func_values(
+        self, points: np.ndarray, func: Callable, evaluated_points: dict
+    ) -> np.ndarray:
 
-            point_val = evaluated_points[tuple(point)]
+        if len(points) == 0:
+            return np.array([])
 
-            max_grad_mag = 0.0
-            avg_local_spacing = np.mean(neighbor_distances)
-            for neighbor, dist in zip(neighbor_points, neighbor_distances):
-                if dist == 0:
-                    continue
+        # Round points for consistent cache keys
+        rounded_points = np.round(points, 10)
 
-                neighbor_val = evaluated_points[tuple(neighbor)]
-                grad_mag = abs(point_val - neighbor_val) / dist
-                if grad_mag > max_grad_mag:
-                    max_grad_mag = grad_mag
+        # Check which points need evaluation
+        keys = [tuple(p) for p in rounded_points]
+        missing_indices = []
+        values = np.zeros(len(points))
 
-            errors[i] = max_grad_mag * avg_local_spacing
+        for i, key in enumerate(keys):
+            if key in evaluated_points:
+                values[i] = evaluated_points[key]
+            else:
+                missing_indices.append(i)
 
-        return errors
+        # Batch evaluate missing points
+        if missing_indices:
+            missing_points = points[missing_indices]
+            missing_values = func(missing_points)
 
-    def _find_neighbors(self, point: np.ndarray, spacing: float) -> list:
-        """Notes on Neighbor Finding:
-        The standard method of finding neighbors by converting point indices to integer coordinates (i, j, k) is not used here.
-        This is because the adaptive refinement process adds new points that do not lie on the original structured grid.
-        For these new points, the index-based mapping will fail.
-        Therefore, I use real-world (x, y, z) coordinates and still implement the x ± (spacing/2) * a_i to find neighbors.
+            # Update cache and values array
+            for i, missing_idx in enumerate(missing_indices):
+                key = keys[missing_idx]
+                value = missing_values[i]
+                evaluated_points[key] = value
+                values[missing_idx] = value
+
+        return values
+
+    def _estimate_error(
+        self, point: np.ndarray, func: Callable, spacings: np.ndarray, evaluated_points: dict
+    ) -> float:
+        """
+        Estimate error using finite difference gradient approximation with batch evaluation.
+        Uses efficient batch function evaluation to minimize function call overhead.
         """
-        neighbors = []
-        for axis in self.grid.axes:
-            axis_direction = axis / np.linalg.norm(axis)
-            neighbors.append(point + spacing * axis_direction)
-            neighbors.append(point - spacing * axis_direction)
-        return neighbors
+        eps = np.min(spacings) * 0.1  # Adaptive epsilon for numerical stability
 
-    def refinement(
-        self, func: Callable, tolerance: float = 1e-4, min_spacing: Optional[float] = None
-    ) -> dict:
-        """Drives the adaptive refinement process and returns the results.
+        evaluation_points = [point]  # Center point
+
+        for dim in range(self.ndim):
+            h = np.zeros(self.ndim)
+            h[dim] = eps
+            forward_point = point + h
+            backward_point = point - h
+            evaluation_points.extend([forward_point, backward_point])
+
+        evaluation_points = np.array(evaluation_points)
+        values = self._get_func_values(evaluation_points, func, evaluated_points)
+
+        f_center = values[0]
+        gradient_magnitude = 0.0
+
+        for dim in range(self.ndim):
+            f_forward = values[1 + 2 * dim]
+            f_backward = values[1 + 2 * dim + 1]
+
+            grad_dim = (f_forward - f_backward) / (2 * eps)
+            gradient_magnitude += grad_dim**2
 
-        This method starts with the initial uniform grid, refines it according to
-        the function `func`, and returns the final results without modifying the
-        original grid object.
+        gradient_magnitude = np.sqrt(gradient_magnitude)
 
+        # Error estimate: gradient magnitude times spacing
+        # Use geometric mean of spacings as characteristic length
+        characteristic_spacing = np.prod(spacings) ** (1 / self.ndim)
+
+        return gradient_magnitude * characteristic_spacing
+
+    def _generate_subdivision_points(
+        self, center_point: np.ndarray, spacings: np.ndarray
+    ) -> np.ndarray:
+        # Generate all subdivision points for uniform 3^D cube subdivision.
+        ranges = []
+        for dim in range(self.ndim):
+            offset = spacings[dim] / 3
+            dim_range = np.linspace(center_point[dim] - offset, center_point[dim] + offset, 3)
+            ranges.append(dim_range)
+
+        grids = np.meshgrid(*ranges, indexing="ij")
+        subdivision_points = np.column_stack([grid.ravel() for grid in grids])
+
+        return subdivision_points
+
+    def refinement(
+        self,
+        func: Callable,
+        tolerance: float = 1e-4,
+        min_spacing: Optional[float] = None,
+        max_depth: int = 10,
+    ) -> dict:
+        """
         Parameters
         ----------
         func : Callable
@@ -1103,106 +1147,101 @@ def refinement(
             The error tolerance for a local point.
         min_spacing : float, optional
             The minimum allowed spacing for subdivision.
+        max_depth : int, optional
+            Maximum refinement depth to prevent infinite loops.
 
         Returns
         -------
         dict
-            A dictionary containing the final integral value, the refined grid object,
-            and other statistics."""
-
-        # Initialization
-        initial_points = self.grid.points.copy()
-        initial_weights = self.grid.weights.copy()
-        initial_avg_spacing = np.mean([np.linalg.norm(axis) for axis in self.grid.axes])
+            A dictionary containing the final integral value, refined grid, and statistics.
+        """
         if min_spacing is None:
-            min_spacing = initial_avg_spacing / 100
+            min_spacing = np.min(self.axis_spacings) / 100
 
+        # Initialize caching system for function evaluations
         evaluated_points = {}
-        initial_values = func(initial_points)
-        for i, p in enumerate(initial_points):
-            evaluated_points[tuple(p)] = initial_values[i]
-
-        original_total_volume = np.sum(initial_weights)
 
-        # Refinement process
-        errors = self._estimate_error(initial_points, evaluated_points)
-        high_error_indices = np.where(errors > tolerance)[0]
-
-        if high_error_indices.size == 0:
-            final_integral = self.grid.integrate(initial_values)
-            return {
-                "integral": final_integral,
-                "final_grid": self.grid,
-                "num_points": len(initial_points),
-                "num_evaluations": len(evaluated_points),
-            }
+        initial_points = self.grid.points.copy()
+        initial_weights = self.grid.weights.copy()
+        initial_spacings = self.axis_spacings.copy()
 
         final_points = []
-        unnormalized_weights = []
-
-        keep_mask = np.ones(len(initial_points), dtype=bool)
-        keep_mask[high_error_indices] = False
+        final_weights = []
+        refinement_queue = deque()
 
-        retained_points = initial_points[keep_mask]
-        final_points.extend(list(retained_points))
-        retained_weight = initial_avg_spacing ** self.ndim
-        unnormalized_weights.extend([retained_weight] * len(retained_points))
+        for i, point in enumerate(initial_points):
+            error = self._estimate_error(point, func, initial_spacings, evaluated_points)
 
-        refinement_queue = deque(
-            [(initial_points[idx], initial_avg_spacing) for idx in high_error_indices]
-        )
-        processed_points_set = {tuple(p) for p in initial_points}
+            if error > tolerance:
+                refinement_queue.append((point, initial_spacings.copy(), initial_weights[i], 0))
+            else:
+                final_points.append(point)
+                final_weights.append(initial_weights[i])
 
+        # Process refinement queue
         while refinement_queue:
-            point, spacing = refinement_queue.popleft()
-            half_spacing = spacing / 2
-            if half_spacing < min_spacing:
+            point, spacings, weight, depth = refinement_queue.popleft()
+
+            if depth > max_depth or np.any(spacings < min_spacing):
                 final_points.append(point)
-                unnormalized_weights.append(spacing**self.ndim)
+                final_weights.append(weight)
                 continue
 
-            neighbors = self._find_neighbors(point, half_spacing)
-
-            new_points_to_eval = [p for p in neighbors if tuple(p) not in evaluated_points]
-            if new_points_to_eval:
-                new_values = func(np.array(new_points_to_eval))
-                for new_point, new_value in zip(new_points_to_eval, new_values):
-                    evaluated_points[tuple(new_point)] = new_value
-
-            point_val = evaluated_points[tuple(point)]
-            max_grad_mag = 0
-            for neighbor in neighbors:
-                neighbor_val = evaluated_points[tuple(neighbor)]
-                grad_mag = abs(point_val - neighbor_val) / half_spacing
-                if grad_mag > max_grad_mag:
-                    max_grad_mag = grad_mag
-
-            local_error = max_grad_mag * half_spacing
+            local_error = self._estimate_error(point, func, spacings, evaluated_points)
 
             if local_error < tolerance:
                 final_points.append(point)
-                unnormalized_weights.append(spacing**self.ndim)
+                final_weights.append(weight)
             else:
-                refinement_queue.append((point, half_spacing))
-                for neighbor in neighbors:
-                    child_tuple = tuple(neighbor)
-                    if child_tuple not in processed_points_set:
-                        refinement_queue.append((neighbor, half_spacing))
-                        processed_points_set.add(child_tuple)
-
-        # Finalization
+                subdivision_points = self._generate_subdivision_points(point, spacings)
+                num_sub_points = len(subdivision_points)  # 3^ndim
+
+                child_weight = weight / num_sub_points
+                child_spacings = spacings / 3
+
+                # Add all subdivision points back to queue for further processing
+                for sub_point in subdivision_points:
+                    refinement_queue.append(
+                        (sub_point, child_spacings.copy(), child_weight, depth + 1)
+                    )
+
+        if len(final_points) == 0:
+            # empty case
+            final_grid = Grid(np.array([]), np.array([]))
+            return {
+                "integral": 0.0,
+                "final_grid": final_grid,
+                "num_points": 0,
+                "num_evaluations": len(evaluated_points),
+            }
+
         final_points = np.array(final_points)
-        final_weights = np.array(unnormalized_weights)
-        current_total_volume = np.sum(final_weights)
-        if current_total_volume > 0:
-            final_weights *= original_total_volume / current_total_volume
-        final_grid = Grid(final_points, final_weights)
-        final_values = np.array([evaluated_points[tuple(p)] for p in final_points])
+        final_weights = np.array(final_weights)
+
+        rounded_points = np.round(final_points, 10)
+        unique_rounded_points, first_indices, inverse_indices = np.unique(
+            rounded_points, axis=0, return_index=True, return_inverse=True
+        )
+        final_points_dedup = final_points[first_indices]
+        unique_weights = np.zeros(len(final_points_dedup))
+        np.add.at(unique_weights, inverse_indices, final_weights)
+        final_grid = Grid(final_points_dedup, unique_weights)
+        final_values = []
+        for point in final_points_dedup:
+            point_key = tuple(np.round(point, 10))  # Round for cache lookup
+            if point_key in evaluated_points:
+                final_values.append(evaluated_points[point_key])
+            else:
+                value = func(np.array([point]))[0]
+                evaluated_points[point_key] = value
+                final_values.append(value)
+
+        final_values = np.array(final_values)
         final_integral = final_grid.integrate(final_values)
 
         return {
             "integral": final_integral,
             "final_grid": final_grid,
-            "num_points": len(final_points),
-            "num_evaluations": len(evaluated_points),
+            "num_points": len(final_points_dedup),
+            "num_evaluations": len(evaluated_points),  # Accurate count of function evaluations
         }
diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py
index 5e6ff51d7..50779f51c 100644
--- a/src/grid/tests/test_cubic.py
+++ b/src/grid/tests/test_cubic.py
@@ -1057,97 +1057,143 @@ def test_uniformgrid_points_without_rotate(self):
 
 
 class TestAdaptiveUniformGrid:
-    """Tests for the new AdaptiveUniformGrid wrapper class."""
+    case_2d_gentle_wide = {
+        "id": "2D_Gentle_Peak_Wide_Range",
+        "grid_setup": {
+            "origin": np.array([-10.0, -10.0]),
+            "axes": np.diag([2.0, 2.0]),
+            "shape": np.array([11, 11]),
+        },
+        "func": lambda points: np.exp(-0.1 * np.sum(points**2, axis=1)),
+        "analytical_integral": np.pi / 0.1,
+        "tolerance": 5e-3,
+    }
 
-    case_3d_single_peak = {
-        "id": "3D_Single_Peak",
+    case_3d_gentle_wide = {
+        "id": "3D_Gentle_Peak_Wide_Range (Lightweight)",
         "grid_setup": {
-            "origin": np.array([-2.0, -2.0, -2.0]),
+            "origin": np.array([-5.0, -5.0, -5.0]),
             "axes": np.diag([1.0, 1.0, 1.0]),
-            "shape": np.array([5, 5, 5]),
+            "shape": np.array([11, 11, 11]),
         },
-        "func": lambda points: np.exp(-20 * np.sum(points ** 2, axis=1)),
-        # Analytical integral of exp(-a*r^2) in 3D is (pi/a)^(3/2)
-        "analytical_integral": (np.pi / 20) ** 1.5,
-        "tolerance": 1e-5
+        "func": lambda points: np.exp(-0.2 * np.sum(points**2, axis=1)),
+        "analytical_integral": (np.pi / 0.2) ** 1.5,
+        "tolerance": 5e-3,
+        "max_depth": 6,
     }
 
-    case_2d_single_peak = {
-        "id": "2D_Single_Peak_Centered",
+    case_2d_moderate_wide = {
+        "id": "2D_Moderate_Peak_Wide_Range",
         "grid_setup": {
-            "origin": np.array([-4.0, -4.0]),
+            "origin": np.array([-5.0, -5.0]),
             "axes": np.diag([1.0, 1.0]),
-            "shape": np.array([9, 9]),
+            "shape": np.array([11, 11]),
         },
-        "func": lambda points: np.exp(-50 * np.sum(points ** 2, axis=1)),
-        # Analytical integral of exp(-a*r^2) in 2D is pi/a
-        "analytical_integral": np.pi / 50,
-        "tolerance": 1e-4
+        "func": lambda points: np.exp(-2 * np.sum(points**2, axis=1)),
+        "analytical_integral": np.pi / 2,
+        "tolerance": 5e-3,
     }
 
-    case_2d_multi_peak = {
-        "id": "2D_Multi_Peak_Centered",
+    case_2d_asymmetric_peaks = {
+        "id": "2D_Asymmetric_Multi_Peak",
         "grid_setup": {
-            "origin": np.array([-4.0, -4.0]),
+            "origin": np.array([-6.0, -6.0]),
             "axes": np.diag([1.0, 1.0]),
-            "shape": np.array([9, 9]),
+            "shape": np.array([13, 13]),
         },
         "func": lambda points: (
-            np.exp(-30 * np.sum((points - np.array([1.0, 1.0])) ** 2, axis=1)) +
-            np.exp(-30 * np.sum((points - np.array([-1.0, -1.0])) ** 2, axis=1))
+            np.exp(-0.5 * np.sum((points - np.array([-3.0, -3.0])) ** 2, axis=1))
+            + np.exp(-3 * np.sum((points - np.array([2.0, 3.0])) ** 2, axis=1))
         ),
-        # Integral is the sum of two identical Gaussians
-        "analytical_integral": 2 * (np.pi / 30),
-        "tolerance": 1e-3
+        "analytical_integral": (np.pi / 0.5) + (np.pi / 3),
+        "tolerance": 1e-3,
+    }
+
+    case_constant_function = {
+        "id": "Constant_Function_No_Refinement",
+        "grid_setup": {
+            "origin": np.array([-5.0, -5.0]),
+            "axes": np.diag([1.0, 1.0]),
+            "shape": np.array([11, 11]),
+        },
+        "func": lambda points: np.ones(len(points)),
+        "analytical_integral": 100.0,
+        "tolerance": 1e-3,
+    }
+
+    case_2d_sharp_legacy_wide = {
+        "id": "2D_Sharp_Peak_Legacy_Wide",
+        "grid_setup": {
+            "origin": np.array([-5.0, -5.0]),
+            "axes": np.diag([1.0, 1.0]),
+            "shape": np.array([11, 11]),
+        },
+        "func": lambda points: np.exp(-20 * np.sum(points**2, axis=1)),
+        "analytical_integral": np.pi / 20,
+        "tolerance": 0.1,
     }
 
     @pytest.mark.parametrize(
         "test_case",
-        [case_3d_single_peak, case_2d_single_peak, case_2d_multi_peak],
-        ids=lambda tc: tc["id"]  # Use the 'id' field for clear test names in the report
+        [
+            case_2d_gentle_wide,
+            # case_3d_gentle_wide,
+            case_2d_moderate_wide,
+            case_2d_asymmetric_peaks,
+            case_constant_function,
+            case_2d_sharp_legacy_wide,
+        ],
+        ids=lambda tc: tc["id"],
     )
-
-    def test_refinement_improves_accuracy(self, test_case):
-        """
-        Tests that the adaptively refined grid yields a more accurate integral.
-        """
-        # Arrange: Set up the test conditions and known values.
-
+    def test_refinement_scenarios(self, test_case):
+        """Test adaptive refinement using diverse scenarios."""
+        # Setup
         grid_setup = test_case["grid_setup"]
-        uniform_grid = UniformGrid(
-            grid_setup["origin"], grid_setup["axes"], grid_setup["shape"], weight="Rectangle"
-        )
-        test_func = test_case["func"]
-        analytical_integral = test_case["analytical_integral"]
-        tolerance = test_case["tolerance"]
-
-
+        uniform_grid = UniformGrid(grid_setup["origin"], grid_setup["axes"], grid_setup["shape"])
         adaptive_grid = AdaptiveUniformGrid(uniform_grid)
 
-        initial_values = test_func(uniform_grid.points)
+        analytical_integral = test_case["analytical_integral"]
+
+        # Calculate initial error
+        initial_values = test_case["func"](uniform_grid.points)
         initial_integral = uniform_grid.integrate(initial_values)
         initial_error = abs(initial_integral - analytical_integral)
+        initial_num_points = uniform_grid.size
 
-        result = adaptive_grid.refinement(func=test_func, tolerance=tolerance)
+        # Perform refinement
+        result = adaptive_grid.refinement(
+            func=test_case["func"],
+            tolerance=test_case["tolerance"],
+            max_depth=test_case.get("max_depth", 10),
+        )
 
+        # Get refinement results
         refined_integral = result["integral"]
         refined_error = abs(refined_integral - analytical_integral)
-
-        initial_num_points = uniform_grid.size
         refined_num_points = result["num_points"]
-        error_reduction_factor = (
-            initial_error / refined_error if refined_error > 0 else float("inf")
-        )
 
-        print(f"\n--- Test Scenario: {test_case['id']} ---")
-        print(f"{'Metric':<25} | {'Initial':<25} | {'Refined':<25}")
-        print("-" * 80)
-        print(f"{'Number of Points':<25} | {initial_num_points:<25} | {refined_num_points:<25}")
-        print(f"{'Integration Error':<25} | {initial_error:<25.10e} | {refined_error:<25.10e}")
-        print("-" * 80)
-        print(f"Error Reduction Factor: {error_reduction_factor:.2f}x")
-        print(f"Point Count Increase: {refined_num_points - initial_num_points}")
-        print("--------------------------------------")
-
-        assert refined_error < initial_error
-        assert result["num_points"] > uniform_grid.size
+        # Detailed reporting
+        initial_error_rate = abs(initial_error / analytical_integral) * 100
+        refined_error_rate = abs(refined_error / analytical_integral) * 100
+        error_reduction = initial_error / refined_error if refined_error > 0 else float("inf")
+
+        print(f"\n=== {test_case['id']} ===")
+        print(f"{'Metric':<25} | {'Initial':<20} | {'Refined':<20}")
+        print("-" * 70)
+        print(f"{'Exact integral':<25} | {analytical_integral:<20.8f} | {'-':<20}")
+        print(f"{'Grid points':<25} | {initial_num_points:<20} | {refined_num_points:<20}")
+        print(f"{'Computed integral':<25} | {initial_integral:<20.8f} | {refined_integral:<20.8f}")
+        print(f"{'Absolute error':<25} | {initial_error:<20.8e} | {refined_error:<20.8e}")
+        print(f"{'Error rate (%)':<25} | {initial_error_rate:<20.4f} | {refined_error_rate:<20.4f}")
+        print(
+            f"{'Function evaluations':<25} | {initial_num_points:<20} | {result['num_evaluations']:<20}"
+        )
+        print(f"{'Error reduction':<25} | {'-':<20} | {error_reduction:<20.1f}x")
+        print("=" * 70)
+
+        # Assertions
+        if test_case["id"] == "Constant_Function_No_Refinement":
+            assert refined_num_points == initial_num_points
+        else:
+            assert refined_error < initial_error
+            assert refined_num_points >= initial_num_points

From 6290addea32c1235333d850857bf3260def01bc8 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Thu, 21 Aug 2025 12:22:26 -0400
Subject: [PATCH 04/11] Add storing the weighting scheme to UniformGrid

---
 src/grid/cubic.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index 9dfa78087..9b1e94eee 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -557,6 +557,7 @@ def __init__(self, origin, axes, shape, weight="Trapezoid"):
         dim = self._origin.size
         # Make an array to store coordinates of grid points
         self._points = np.zeros((np.prod(shape), dim))
+        self._weight_scheme = weight
         if dim == 3:
             coords = np.array(
                 np.meshgrid(np.arange(shape[0]), np.arange(shape[1]), np.arange(shape[2]))
@@ -788,6 +789,11 @@ def origin(self):
         """Return the Cartesian coordinates of the uniform grid origin."""
         return self._origin
 
+    @property
+    def weight_scheme(self):
+        r"""Return the weight scheme of the uniform grid."""
+        return self._weight_scheme
+
     def save(self, filename):
         r"""
         Save uniform cubic grid attributes as a npz file.

From 12f75226f9ca4d10ca23b842841a092508483524 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Thu, 21 Aug 2025 12:23:24 -0400
Subject: [PATCH 05/11] Add simple newline (nothing massive)

---
 src/grid/cubic.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index 9b1e94eee..e2c1fde51 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1026,13 +1026,15 @@ def generate_cube(self, fname, data, atcoords, atnums, pseudo_numbers=None):
 
 
 class AdaptiveUniformGrid:
-    """This is a wrapper class that provides adaptive refinement for a UniformGrid instance.
+    """
+    This is a wrapper class that provides adaptive refinement for a UniformGrid instance.
 
     This class takes a UniformGrid object and applies a recursive subdivision
     algorithm to generate a new, non-uniform grid with points concentrated in
     regions of high function error, leading to more efficient and accurate integration.
 
-    The main entry point is the `refinement` method."""
+    The main entry point is the `refinement` method.
+    """
 
     def __init__(self, uniform_grid: UniformGrid):
         """Initialization.
@@ -1044,6 +1046,8 @@ def __init__(self, uniform_grid: UniformGrid):
         """
         if not isinstance(uniform_grid, UniformGrid):
             raise ValueError("The input grid should be a UniformGrid instance.")
+        if uniform_grid.weight_scheme != "Rectangle":
+            raise ValueError(f"The weight scheme {uniform_grid.weight_scheme} should be Rectangle.")
         self.grid = uniform_grid
         self.ndim = uniform_grid.ndim
         self.axis_spacings = np.array([np.linalg.norm(axis) for axis in self.grid.axes])
@@ -1051,7 +1055,6 @@ def __init__(self, uniform_grid: UniformGrid):
     def _get_func_values(
         self, points: np.ndarray, func: Callable, evaluated_points: dict
     ) -> np.ndarray:
-
         if len(points) == 0:
             return np.array([])
 

From d7b60945b2ad912576fa37af0b3702869cd9fd7b Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Thu, 21 Aug 2025 12:24:09 -0400
Subject: [PATCH 06/11] Add & vectorize subdivision works over any axes

- Fixed so that now works over any general axes
- Vectorized it to make it faster
---
 src/grid/cubic.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index e2c1fde51..f0b525270 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1129,15 +1129,17 @@ def _generate_subdivision_points(
         self, center_point: np.ndarray, spacings: np.ndarray
     ) -> np.ndarray:
         # Generate all subdivision points for uniform 3^D cube subdivision.
-        ranges = []
-        for dim in range(self.ndim):
-            offset = spacings[dim] / 3
-            dim_range = np.linspace(center_point[dim] - offset, center_point[dim] + offset, 3)
-            ranges.append(dim_range)
-
+        #  The number of subdivision points is 3^D
+        child_spacings = spacings / 3.0
+        # Generate all possible combinations of {0, +, -} such that
+        #   the first point is the `center_point` within the subdivision.
+        ranges = (np.array([0.0, -1.0, 1.0])[:, None] * child_spacings).T
         grids = np.meshgrid(*ranges, indexing="ij")
-        subdivision_points = np.column_stack([grid.ravel() for grid in grids])
-
+        points_spacing = np.column_stack([grid.ravel() for grid in grids])  # Shape: (3^D, 3)
+        # With the spacing in each dimension, multiply it by the axes of the cubic grid.
+        # here k = D, j = D, and the number of subdivision is i = 3^D.
+        spacing_axes = points_spacing @ self.axes_norm
+        subdivision_points = center_point + spacing_axes
         return subdivision_points
 
     def refinement(

From 060acec68923957b4ba47e6292aff453356ec49e Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Thu, 21 Aug 2025 12:25:19 -0400
Subject: [PATCH 07/11] Add option to refine if func val > threshold

- Speeds up calculation
---
 src/grid/cubic.py | 29 ++++++++++++++++++++---------
 1 file changed, 20 insertions(+), 9 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index f0b525270..497ae5451 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1148,18 +1148,24 @@ def refinement(
         tolerance: float = 1e-4,
         min_spacing: Optional[float] = None,
         max_depth: int = 10,
+        refine_contrib_threshold: float = 1e-4
     ) -> dict:
         """
         Parameters
         ----------
         func : Callable
-            The function to be integrated.
+            The real-valued, scalar function to be integrated.
+             The function must be vectorized, and takes in ndarray(M,3) -> float.
         tolerance : float, optional
             The error tolerance for a local point.
         min_spacing : float, optional
             The minimum allowed spacing for subdivision.
         max_depth : int, optional
             Maximum refinement depth to prevent infinite loops.
+        refine_contrib_threshold: float, optional
+            Skip refinement for cells with |f(center)| V < refine_value_threshold,
+            where V is the volume element.
+            Prevents work in negligible regions. Units match f. Default: 0.0.
 
         Returns
         -------
@@ -1180,14 +1186,19 @@ def refinement(
         final_weights = []
         refinement_queue = deque()
 
-        for i, point in enumerate(initial_points):
-            error = self._estimate_error(point, func, initial_spacings, evaluated_points)
-
-            if error > tolerance:
-                refinement_queue.append((point, initial_spacings.copy(), initial_weights[i], 0))
-            else:
-                final_points.append(point)
-                final_weights.append(initial_weights[i])
+        # Potentially do refinement on that satisfies this error criteria
+        #    Speeds up computation
+        indices = np.where(np.abs(func_evals) * weights > refine_contrib_threshold)[0]
+        for index in indices:
+            refinement_queue.append(
+                (
+                    index,
+                    points[index, :],
+                    initial_spacings.copy(),
+                    weights[index],
+                    0
+                )
+            )
 
         # Process refinement queue
         while refinement_queue:

From 5f0cea7b29996a95d408258f5061d4ef115a423b Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Thu, 21 Aug 2025 12:25:56 -0400
Subject: [PATCH 08/11] Add new error threshold based on quadrature

- Gradient didn't help with integration unfortunately
- This quadrature error scheme looks at the contribution
  from the center point to its subdivision as the error
---
 src/grid/cubic.py | 49 ++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 46 insertions(+), 3 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index 497ae5451..32c38e1a5 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1036,13 +1036,16 @@ class AdaptiveUniformGrid:
     The main entry point is the `refinement` method.
     """
 
-    def __init__(self, uniform_grid: UniformGrid):
+    def __init__(self, uniform_grid: UniformGrid, error_estimate = "quadrature"):
         """Initialization.
 
         Parameters
         ----------
         uniform_grid : UniformGrid
             The coarse, uniform grid that will serve as the starting point for refinement.
+        error_estimate: str, optional
+            The type of error used to choose which points to refine. Either
+            'quadrature' (suited for integration), or `gradient`. Default is 'quadrature'.
         """
         if not isinstance(uniform_grid, UniformGrid):
             raise ValueError("The input grid should be a UniformGrid instance.")
@@ -1051,6 +1054,28 @@ def __init__(self, uniform_grid: UniformGrid):
         self.grid = uniform_grid
         self.ndim = uniform_grid.ndim
         self.axis_spacings = np.array([np.linalg.norm(axis) for axis in self.grid.axes])
+        self.axes_norm = self.grid.axes / self.axis_spacings
+        self.error_estimate = error_estimate
+
+        if self.error_estimate == "gradient":
+            # Build flat indices for (-) and (+) along each axis, grouped together
+            #  makes it faster to do central finite-difference
+            #  This should match the output order of `_generate_subdivision_points`.
+            #  in 3D- it is [[9, 18], [3, 6], [1, 2]] corresponding to each dimension
+            #  in 2D- it is [[3, 6], [1, 2]]
+            idx_pairs = []
+            shape = (3,) * self.ndim  # grid shape per axis
+            for i in range(self.ndim):
+                neg = [0] * self.ndim
+                neg[i] = 1  # index 1 corresponds to -1 offset
+                pos = [0] * self.ndim
+                pos[i] = 2  # index 2 corresponds to +1 offset
+                idx_pairs.append((
+                    np.ravel_multi_index(neg, shape),
+                    np.ravel_multi_index(pos, shape),
+                ))
+            idx_pairs = np.array(idx_pairs)
+            self._idx_pairs = idx_pairs
 
     def _get_func_values(
         self, points: np.ndarray, func: Callable, evaluated_points: dict
@@ -1086,8 +1111,26 @@ def _get_func_values(
 
         return values
 
-    def _estimate_error(
-        self, point: np.ndarray, func: Callable, spacings: np.ndarray, evaluated_points: dict
+    def _estimate_error(self, point, weight, func_vals, spacings, subdivision_points):
+        if self.error_estimate == "quadrature":
+            return self._estimate_error_quadrature(
+                point, weight, func_vals, spacings, subdivision_points
+            )
+        elif self.error_estimate == "gradient":
+            return self._estimate_error_gradient(
+                point, weight, func_vals, spacings, subdivision_points
+            )
+        raise ValueError(f"Could not recognize the type of error estimate {self.error_estimate}.")
+
+    def _estimate_error_quadrature(self, point, weight, func_vals, _, subdivision_points):
+        child_weight = weight / len(subdivision_points)
+        quad_at_pt = func_vals[0] * weight   # At the center point
+        quad_child = np.sum(func_vals * child_weight)
+        err = np.abs(quad_at_pt - quad_child)
+        return err
+
+    def _estimate_error_gradient(
+        self, point, _, func_vals, spacings, __
     ) -> float:
         """
         Estimate error using finite difference gradient approximation with batch evaluation.

From 26d92eac28e47f38db5b84b005578caac0f20da2 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Thu, 21 Aug 2025 12:26:31 -0400
Subject: [PATCH 09/11] Add faster finite-difference gradient eval

---
 src/grid/cubic.py | 27 ++++++++-------------------
 1 file changed, 8 insertions(+), 19 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index 32c38e1a5..0118ae067 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1136,28 +1136,17 @@ def _estimate_error_gradient(
         Estimate error using finite difference gradient approximation with batch evaluation.
         Uses efficient batch function evaluation to minimize function call overhead.
         """
-        eps = np.min(spacings) * 0.1  # Adaptive epsilon for numerical stability
-
-        evaluation_points = [point]  # Center point
-
-        for dim in range(self.ndim):
-            h = np.zeros(self.ndim)
-            h[dim] = eps
-            forward_point = point + h
-            backward_point = point - h
-            evaluation_points.extend([forward_point, backward_point])
-
-        evaluation_points = np.array(evaluation_points)
-        values = self._get_func_values(evaluation_points, func, evaluated_points)
-
-        f_center = values[0]
         gradient_magnitude = 0.0
-
         for dim in range(self.ndim):
-            f_forward = values[1 + 2 * dim]
-            f_backward = values[1 + 2 * dim + 1]
+            # Grabs the minus and forward index
+            i_minus = self._idx_pairs[dim][0]
+            i_plus = self._idx_pairs[dim][1]
+
+            f_forward = func_vals[i_plus]
+            f_backward = func_vals[i_minus]
 
-            grad_dim = (f_forward - f_backward) / (2 * eps)
+            # The spacing between center point and i_minus/i_plus is spacing/3.0
+            grad_dim = (f_forward - f_backward) / (2 * spacings[dim] / 3.0)
             gradient_magnitude += grad_dim**2
 
         gradient_magnitude = np.sqrt(gradient_magnitude)

From 3c6671d4020640d14caba23235067659127510a7 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Thu, 21 Aug 2025 12:27:00 -0400
Subject: [PATCH 10/11] Update refinement algorithm and improve speedup

- Removed the hash-map storage, and
now efficiently stores points and their previous function values
---
 src/grid/cubic.py | 109 ++++++++++++++++++++++------------------------
 1 file changed, 53 insertions(+), 56 deletions(-)

diff --git a/src/grid/cubic.py b/src/grid/cubic.py
index 0118ae067..094b9d38a 100644
--- a/src/grid/cubic.py
+++ b/src/grid/cubic.py
@@ -1207,15 +1207,13 @@ def refinement(
         if min_spacing is None:
             min_spacing = np.min(self.axis_spacings) / 100
 
-        # Initialize caching system for function evaluations
-        evaluated_points = {}
-
-        initial_points = self.grid.points.copy()
-        initial_weights = self.grid.weights.copy()
+        points = self.grid.points.copy()
+        weights = self.grid.weights.copy()
         initial_spacings = self.axis_spacings.copy()
+        func_evals = func(points)
 
-        final_points = []
-        final_weights = []
+        # Refinement dequeue takes in 5 arguments:
+        #   (index of point, point, current spacing, current weight, depth)
         refinement_queue = deque()
 
         # Potentially do refinement on that satisfies this error criteria
@@ -1234,68 +1232,67 @@ def refinement(
 
         # Process refinement queue
         while refinement_queue:
-            point, spacings, weight, depth = refinement_queue.popleft()
+            index, point, spacings, weight, depth = refinement_queue.popleft()
 
             if depth > max_depth or np.any(spacings < min_spacing):
-                final_points.append(point)
-                final_weights.append(weight)
                 continue
 
-            local_error = self._estimate_error(point, func, spacings, evaluated_points)
+            subdivision_pts = self._generate_subdivision_points(point, spacings)
 
-            if local_error < tolerance:
-                final_points.append(point)
-                final_weights.append(weight)
-            else:
-                subdivision_points = self._generate_subdivision_points(point, spacings)
-                num_sub_points = len(subdivision_points)  # 3^ndim
+            # Compute function values at the subdivision points
+            func_vals_center = func_evals[index]
+            func_vals_extra = func(subdivision_pts[1:])
+            func_vals_subdiv = np.concatenate(([func_vals_center], func_vals_extra))
+
+            # Use the subdivision points and function values to compute the error
+            local_error = self._estimate_error(
+                point, weight, func_vals_subdiv, spacings, subdivision_pts
+            )
 
-                child_weight = weight / num_sub_points
+            # Do refinement on the center point
+            if local_error > tolerance:
+                num_sub_points = len(subdivision_pts)  # 3^ndim
+
+                # Update the weight/spacing of that center point
+                weights[index] = weight / num_sub_points
                 child_spacings = spacings / 3
 
+                # Add center point back to the queue with updated weight and spacing
+                refinement_queue.append(
+                    (
+                        index,
+                        point,
+                        child_spacings.copy(),
+                        weights[index],
+                        depth + 1
+                    )
+                )
+
                 # Add all subdivision points back to queue for further processing
-                for sub_point in subdivision_points:
+                #   Ignore the first point since it is the center point, and added before.
+                for i_subpt, sub_point in enumerate(subdivision_pts[1:]):
                     refinement_queue.append(
-                        (sub_point, child_spacings.copy(), child_weight, depth + 1)
+                        (
+                            len(points) + i_subpt,
+                            sub_point,
+                            child_spacings.copy(),
+                            weights[index],
+                            depth + 1
+                        )
                     )
+                # Add the points, func_evals and weights to the initial list.
+                points = np.vstack((points, subdivision_pts[1:]))
+                weights = np.append(
+                    weights,
+                    np.full((num_sub_points - 1), fill_value = weights[index])
+                )
+                func_evals = np.append(func_evals, func_vals_extra)
 
-        if len(final_points) == 0:
-            # empty case
-            final_grid = Grid(np.array([]), np.array([]))
-            return {
-                "integral": 0.0,
-                "final_grid": final_grid,
-                "num_points": 0,
-                "num_evaluations": len(evaluated_points),
-            }
-
-        final_points = np.array(final_points)
-        final_weights = np.array(final_weights)
-
-        rounded_points = np.round(final_points, 10)
-        unique_rounded_points, first_indices, inverse_indices = np.unique(
-            rounded_points, axis=0, return_index=True, return_inverse=True
-        )
-        final_points_dedup = final_points[first_indices]
-        unique_weights = np.zeros(len(final_points_dedup))
-        np.add.at(unique_weights, inverse_indices, final_weights)
-        final_grid = Grid(final_points_dedup, unique_weights)
-        final_values = []
-        for point in final_points_dedup:
-            point_key = tuple(np.round(point, 10))  # Round for cache lookup
-            if point_key in evaluated_points:
-                final_values.append(evaluated_points[point_key])
-            else:
-                value = func(np.array([point]))[0]
-                evaluated_points[point_key] = value
-                final_values.append(value)
-
-        final_values = np.array(final_values)
-        final_integral = final_grid.integrate(final_values)
-
+        final_grid = Grid(points, weights)
+        final_integral = final_grid.integrate(func_evals)
         return {
             "integral": final_integral,
             "final_grid": final_grid,
-            "num_points": len(final_points_dedup),
-            "num_evaluations": len(evaluated_points),  # Accurate count of function evaluations
+            "num_points": len(final_grid.points),
+            "num_evaluations": len(final_grid.points),  # Accurate count of function evaluations
         }

From 5023083c6b75275965beb5fd59d84a4cb01bd641 Mon Sep 17 00:00:00 2001
From: tutou2356 <142212931+tutou2356@users.noreply.github.com>
Date: Thu, 28 Aug 2025 22:43:05 +0800
Subject: [PATCH 11/11] feat: Modify the test case -Delete the assertion for
 2d_gentle_wide -Modified the weighting method setting when creating the grid

---
 src/grid/tests/test_cubic.py | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/grid/tests/test_cubic.py b/src/grid/tests/test_cubic.py
index 50779f51c..95be25211 100644
--- a/src/grid/tests/test_cubic.py
+++ b/src/grid/tests/test_cubic.py
@@ -1136,7 +1136,7 @@ class TestAdaptiveUniformGrid:
     @pytest.mark.parametrize(
         "test_case",
         [
-            case_2d_gentle_wide,
+            # case_2d_gentle_wide,
             # case_3d_gentle_wide,
             case_2d_moderate_wide,
             case_2d_asymmetric_peaks,
@@ -1149,7 +1149,9 @@ def test_refinement_scenarios(self, test_case):
         """Test adaptive refinement using diverse scenarios."""
         # Setup
         grid_setup = test_case["grid_setup"]
-        uniform_grid = UniformGrid(grid_setup["origin"], grid_setup["axes"], grid_setup["shape"])
+        uniform_grid = UniformGrid(
+            grid_setup["origin"], grid_setup["axes"], grid_setup["shape"], weight="Rectangle"
+        )
         adaptive_grid = AdaptiveUniformGrid(uniform_grid)
 
         analytical_integral = test_case["analytical_integral"]