diff --git a/ipynb/Diss/Figures_Diss.ipynb b/ipynb/Diss/Figures_Diss.ipynb index 5fd54ef..d5edc81 100644 --- a/ipynb/Diss/Figures_Diss.ipynb +++ b/ipynb/Diss/Figures_Diss.ipynb @@ -1684,4 +1684,4 @@ }, "nbformat": 4, "nbformat_minor": 1 -} +} \ No newline at end of file diff --git a/ipynb/Extrapolation/Extrapolation_Test.ipynb b/ipynb/Extrapolation/Extrapolation_Test.ipynb index 851a29b..542bd03 100644 --- a/ipynb/Extrapolation/Extrapolation_Test.ipynb +++ b/ipynb/Extrapolation/Extrapolation_Test.ipynb @@ -31,11 +31,10 @@ "# shortcut for romberg SASD\n", "def get_single_dim_romberg(a, b, f, reference_solution, sasd_version,\n", " slice_grouping, slice_version, container_version,\n", - " force_balanced_refinement_tree):\n", + " force_balanced_refinement_tree, grid_version = GridVersion.DEFAULT):\n", " grid_romberg = GlobalRombergGrid(a=a, b=b, modified_basis=False, boundary=True,\n", " slice_grouping=slice_grouping,\n", - " slice_version=slice_version,\n", - " container_version=container_version)\n", + " slice_version=slice_version)\n", "\n", " operation_romberg = Integration(f=f, grid=grid_romberg, dim=dim, reference_solution=reference_solution)\n", "\n", @@ -90,11 +89,11 @@ "\n", "# Manually override functions array:\n", "functions = [\n", - " # f_genz_corner,\n", + " f_genz_corner,\n", " # f_genz_product,\n", " # f_genz_cont,\n", " # f_genz_gaussian,\n", - " f_exp_var,\n", + " # f_exp_var,\n", " # f_genz_osz,\n", " # f_genz_disc\n", "]\n", @@ -189,19 +188,47 @@ " force_balanced_refinement_tree=True)\n", " \n", " # Interpolating Romberg\n", + " r_interpolating_old = get_single_dim_romberg(a, b, f, reference_solution,\n", + " sasd_version=sasd_version,\n", + " slice_grouping=SliceGrouping.GROUPED_OPTIMIZED,\n", + " slice_version=SliceVersion.ROMBERG_DEFAULT, # Makes no difference\n", + " container_version=SliceContainerVersion.LAGRANGE_ROMBERG_OLD,\n", + " force_balanced_refinement_tree=False,\n", + " grid_version=GridVersion.INTERPOLATE_SUB_GRIDS) # Interpolating Romberg\n", + "\n", " r_interpolating = get_single_dim_romberg(a, b, f, reference_solution,\n", " sasd_version=sasd_version,\n", " slice_grouping=SliceGrouping.GROUPED_OPTIMIZED,\n", " slice_version=SliceVersion.ROMBERG_DEFAULT, # Makes no difference\n", " container_version=SliceContainerVersion.LAGRANGE_ROMBERG,\n", - " force_balanced_refinement_tree=False)\n", + " force_balanced_refinement_tree=False,\n", + " grid_version=GridVersion.INTERPOLATE_SUB_GRIDS)\n", " # Interpolating Full grid\n", " r_interpolating_full_grid = get_single_dim_romberg(a, b, f, reference_solution,\n", " sasd_version=sasd_version,\n", " slice_grouping=SliceGrouping.GROUPED_OPTIMIZED,\n", " slice_version=SliceVersion.ROMBERG_DEFAULT, # Makes no difference\n", - " container_version=SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG,\n", - " force_balanced_refinement_tree=False)\n", + " container_version=SliceContainerVersion.LAGRANGE_ROMBERG,\n", + " force_balanced_refinement_tree=False,\n", + " grid_version=GridVersion.INTERPOLATE_FULL_GRID)\n", + " \n", + " # Interpolating Bspline\n", + " r_interpolating_bspline = get_single_dim_romberg(a, b, f, reference_solution,\n", + " sasd_version=sasd_version,\n", + " slice_grouping=SliceGrouping.GROUPED_OPTIMIZED,\n", + " slice_version=SliceVersion.ROMBERG_DEFAULT, # Makes no difference\n", + " container_version=SliceContainerVersion.BSPLINE_ROMBERG,\n", + " force_balanced_refinement_tree=False,\n", + " grid_version=GridVersion.INTERPOLATE_SUB_GRIDS)\n", + " \n", + " # Interpolating Hierarchical Langrange\n", + " r_interpolating_hierarchical_lagrange = get_single_dim_romberg(a, b, f, reference_solution,\n", + " sasd_version=sasd_version,\n", + " slice_grouping=SliceGrouping.GROUPED_OPTIMIZED,\n", + " slice_version=SliceVersion.ROMBERG_DEFAULT, # Makes no difference\n", + " container_version=SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG,\n", + " force_balanced_refinement_tree=False,\n", + " grid_version=GridVersion.INTERPOLATE_SUB_GRIDS)\n", " \n", " # Romberg Constant Subtraction\n", " r_constant_subtraction = get_single_dim_romberg(a, b, f, reference_solution,\n", @@ -259,21 +286,26 @@ " # (r_grouped_optimized, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Default Romberg)'),\n", " # (r_full_grouped_optimized, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Default Romberg, Balanced)'),\n", " # \n", - " # (r_interpolating, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Lagrange Romberg)'),\n", + " (r_interpolating_old, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Lagrange Romberg Old)'),\n", " # (r_interpolating_full_grid, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Lagrange Full Romberg)'),\n", + "\n", + " (r_interpolating, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Lagrange Romberg)'),\n", + " (r_interpolating_bspline, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Bspline Romberg)'),\n", + " (r_interpolating_hierarchical_lagrange, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Hierarchical Lagrange Romberg)'),\n", + "\n", " # \n", " # (r_simson_full_grouped_optimized, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Simpson Romberg, Balanced)'),\n", "\n", " # (r_constant_subtraction, 1, 2, errorOperator, 'Extrapolation Grid (Unit, DefaultConstantSubtraction)'), \n", " \n", - " (r_balanced, 1, 2, errorOperator, 'Balanced Extrapolation Grid'),\n", + " # (r_balanced, 1, 2, errorOperator, 'Balanced Extrapolation Grid'),\n", " ]\n", " \n", " function_name = f.__class__.__name__\n", " filename = \"error_comparison_{}_{}d\".format(function_name, dim)\n", "\n", " performTestcaseArbitraryDim(f, a, b, algorithmArray,\n", - " max_tol, dim, 6, grids=standard_combi_grids, evaluation_points=evaluation_points,\n", + " max_tol, dim, 6, grids=standard_combi_grids, evaluation_points=None,\n", " max_evaluations=max_evaluations,\n", " calc_standard_schemes=True, # enable for standard schemes\n", " minLmin=1,\n", diff --git a/sparseSpACE/BasisFunctions.py b/sparseSpACE/BasisFunctions.py index 04b2f4b..be6043a 100644 --- a/sparseSpACE/BasisFunctions.py +++ b/sparseSpACE/BasisFunctions.py @@ -94,12 +94,15 @@ def get_integral(self, a: float, b: float, coordsD: np.array, weightsD: np.array return result +# Data structure for one lagrange basis class LagrangeBasis(BasisFunction): def __init__(self, p: int, index: int, knots: np.array): self.p = p self.knots = knots #assert p == len(knots) - 1 self.index = index + + # Construct factor of lagrange polynomial for this basis point self.factor = 1 for i, knot in enumerate(self.knots): assert not isinf(self.knots[i]) @@ -109,10 +112,15 @@ def __init__(self, p: int, index: int, knots: np.array): def __call__(self, x: float) -> float: result = 1 + + # Evaluate basis function at point x for i, knot in enumerate(self.knots): if self.index != i: result *= (x - self.knots[i]) - return result * self.factor + + # Complete construction of lagrange polynomial for this basis point evaluated at x + evaluated_basis = result * self.factor + return evaluated_basis def get_first_derivative(self, x: float) -> float: return self.derivative_for_index(x, [self.index]) diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index b96cebe..0e4b7f8 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -517,6 +517,12 @@ def get_grid_levels(self) -> Sequence[int]: # ----------------------------------------------------------------------------------------------------------------- # --- Extrapolation Grid +class GridVersion(Enum): + DEFAULT = 1 # No interpolation are inserted + INTERPOLATE_SUB_GRIDS = 2 # Important missing grid points are inserted as interpolation points + INTERPOLATE_FULL_GRID = 3 # All missing points needed for a full grid are inserted as interpolation points + + class SliceGrouping(Enum): """ This enum specifies the available options for slice grouping. @@ -526,23 +532,30 @@ class SliceGrouping(Enum): GROUPED_OPTIMIZED = 3 # Group slices until container size is the max possible 2^k -class SliceVersion(Enum): +class SliceContainerVersion(Enum): """ This enum specifies the available slice types. """ - ROMBERG_DEFAULT = 1 # Sliced Romberg with default Romberg extrapolation - TRAPEZOID = 2 # Default trapezoidal rule without extrapolation - ROMBERG_DEFAULT_CONST_SUBTRACTION = 3 # Sliced Romberg with default Romberg extrapolation and constant subtraction + # Only extrapolation + ROMBERG_DEFAULT = 1 # Default Romberg extrapolation in non-unit containers + SIMPSON_ROMBERG = 2 # Romberg extrapolation with simpson rule as base quadrature rule + # Interpolation before extrapolation + LAGRANGE_ROMBERG_OLD = 3 # Interpolated grid points are interpolated using Lagrange -class SliceContainerVersion(Enum): + # Interpolation before extrapolation using sparseSpACE grids + LAGRANGE_ROMBERG = 4 # Interpolated grid points are interpolated using Lagrange + HIERARCHICAL_LAGRANGE_ROMBERG = 5 # Interpolated grid points are interpolated using hierarchical Lagrange + BSPLINE_ROMBERG = 6 # Interpolated grid points are interpolated using BSplines + + +class SliceVersion(Enum): """ This enum specifies the available slice types. """ - ROMBERG_DEFAULT = 1 # Default Romberg extrapolation in non-unit containers - LAGRANGE_ROMBERG = 2 # Important missing grid points are interpolated with Lagrange - LAGRANGE_FULL_GRID_ROMBERG = 3 # All missing points needed for a full grid are interpolated - SIMPSON_ROMBERG = 4 # Romberg extrapolation with simpson rule as base quadrature rule + ROMBERG_DEFAULT = 1 # Sliced Romberg with default Romberg extrapolation + TRAPEZOID = 2 # Default trapezoidal rule without extrapolation + ROMBERG_DEFAULT_CONST_SUBTRACTION = 3 # Sliced Romberg with default Romberg extrapolation and constant subtraction class ExtrapolationGrid: @@ -559,9 +572,10 @@ class ExtrapolationGrid: """ def __init__(self, + grid_version: GridVersion = GridVersion.DEFAULT, + container_version: SliceContainerVersion = SliceContainerVersion.ROMBERG_DEFAULT, slice_grouping: SliceGrouping = SliceGrouping.UNIT, slice_version: SliceVersion = SliceVersion.ROMBERG_DEFAULT, - container_version: SliceContainerVersion = SliceContainerVersion.ROMBERG_DEFAULT, force_balanced_refinement_tree: bool = False, print_debug: bool = False): self.print_debug = print_debug @@ -576,24 +590,29 @@ def __init__(self, self.integral_approximation = None # Different grid versions + self.grid_version = grid_version + self.container_version = container_version self.slice_grouping = slice_grouping self.slice_version = slice_version - self.container_version = container_version self.force_balanced_refinement_tree = force_balanced_refinement_tree # Factories self.slice_factory = ExtrapolationGridSliceFactory(self.slice_version) self.container_factory = ExtrapolationGridSliceContainerFactory(self.container_version) - # Interpolation - if self.container_version == SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG: + # Should interpolation grid points be inserted? + if self.grid_version == GridVersion.INTERPOLATE_FULL_GRID: self.max_interpolation_step_width_delta = np.infty - else: + elif self.grid_version == GridVersion.INTERPOLATE_SUB_GRIDS: self.max_interpolation_step_width_delta = 2 ** 1 + else: + self.max_interpolation_step_width_delta = 0 + + # TODO Assert right combinations of grid_version and container_version def interpolation_is_enabled(self) -> bool: - return self.container_version == SliceContainerVersion.LAGRANGE_ROMBERG or \ - self.container_version == SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG + return self.grid_version == GridVersion.INTERPOLATE_SUB_GRIDS or \ + self.grid_version == GridVersion.INTERPOLATE_FULL_GRID # Integration def integrate(self, function: Function = None) -> float: @@ -2254,7 +2273,7 @@ def get_support_points_for_interpolation_geometrically(self, interp_point: float if adaptive and len(support_points) > 2: for support_point in support_points: - weight = self.get_interpolation_weights(support_points, support_point, interp_point) + weight = self.get_interpolation_weight(support_points, support_point, interp_point) if weight < 0: negative_sum += weight @@ -2276,6 +2295,54 @@ def get_support_points_for_interpolation_geometrically(self, interp_point: float return support_points + @staticmethod + def get_levels_for_support_point_grid(support_points: Sequence[float]): + assert len(support_points) >= 2 + + inner_levels = None + + if len(support_points) == 2: + inner_levels = [] + + else: + inner_levels = InterpolationGridSliceContainer.__get_levels_for_support_points_rec( + support_points, + start_index=1, + stop_index=(len(support_points)-1)-1, # -1 for last index and -1 for 2nd last index + level=1 + ) + + assert len(inner_levels) + 2 == len(support_points) + + return [0] + inner_levels + [0] + + @staticmethod + def __get_levels_for_support_points_rec(support_points: Sequence[float], + start_index: float, stop_index: float, + level: int = 1) -> Sequence[int]: + middle_index = int((start_index + stop_index) / 2) + # middle = support_points[middle_index] + + if start_index > stop_index: + return [] + elif start_index == stop_index: + return [level] + + left_levels = InterpolationGridSliceContainer.__get_levels_for_support_points_rec( + support_points, + start_index, + middle_index - 1, + level + 1 + ) + right_levels = InterpolationGridSliceContainer.__get_levels_for_support_points_rec( + support_points, + middle_index + 1, + stop_index, + level + 1 + ) + + return left_levels + [level] + right_levels + def insert_interpolating_slices(self) -> None: """ This method inserts interpolation slices into the container. @@ -2364,56 +2431,6 @@ def get_interpolation_points(self) -> Sequence[float]: return interpolation_points - @abstractmethod - def get_final_weights(self) -> Dict[float, Sequence[float]]: - pass - - def get_interpolation_weight(self, support_points: Sequence[float], basis_point: float, evaluation_point: float)\ - -> float: - pass - - -class LagrangeRombergGridSliceContainer(InterpolationGridSliceContainer): - """ - This class interpolates as many missing points as possible inside the container and executes a default - Romberg method on this interpolated partial full grid. - - For the interpolation we use the closest points to the left and right of this container. - """ - - def __init__(self, function: Function = None): - """ - Constructor of this class. - - :param function: for error computation. - """ - super(LagrangeRombergGridSliceContainer, self).__init__(function) - self.max_interpolation_support_points = 7 - - def get_interpolation_weights(self, support_points, basis_point, evaluation_point) -> float: - return self.get_langrange_basis(support_points, basis_point, evaluation_point) - - @staticmethod - def get_langrange_basis(support_points: Sequence[float], basis_point: float, evaluation_point: float) -> float: - """ - This method computes a lambda expression for the lagrange basis function at basis_point. - Specifically: Let support_points = [x_0, x_1, ... , x_n] - L = product( (evaluation_point - x_k) / (basis_point - x_k) ) - - :param support_points: an array of support points for the interpolation - :param basis_point: determines basis_point of this lagrange basis function - :param evaluation_point: determines at which point the basis is evaluated - :return: the evaluated basis function L - """ - - evaluated_basis = 1 - - for point in support_points: - if point != basis_point: - evaluated_basis *= (evaluation_point - point) / (basis_point - point) - - return evaluated_basis - def get_final_weights(self) -> Dict[float, Sequence[float]]: assert len(self.slices) > 0 @@ -2472,11 +2489,258 @@ def populate_interpolated_point_weights(self, factory: RombergWeights, weight = factory.get_inner_point_weight(normalized_grid_levels[i], normalized_max_level) + print() for support_point in support_points: - interpolated_extrapolated_weight = weight * self.get_langrange_basis(support_points, support_point, - interp_point) + interpolation_weight = self.get_interpolation_weight(support_points, support_point, interp_point) + # print(interpolation_weight) + interpolated_extrapolated_weight = weight * interpolation_weight weight_dictionary[support_point].append(interpolated_extrapolated_weight) + # TODO add a flag for plotting this + # Plot basis functions + # self.plot_interpolation_basis_functions(support_points, interp_point) + # self.plot_interpolation_basis_functions_with_evaluations(support_points, interp_point) + + def get_interpolation_weight(self, support_points: Sequence[float], + basis_point: float, evaluation_point: float) -> float: + pass + + def plot_interpolation_basis_functions(self, support_points: Sequence[float], interpolation_point: float): + pass + + def plot_interpolation_basis_functions_with_evaluations(self, support_points: Sequence[float], + evaluation_point: float): + pass + + +# This class was replaced with LagrangeRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter) +class LagrangeRombergGridSliceContainerOld(InterpolationGridSliceContainer): + """ + This class interpolates as many missing points as possible inside the container and executes a default + Romberg method on this interpolated partial full grid. + + For the interpolation we use the closest points to the left and right of this container. + """ + + def __init__(self, function: Function = None): + """ + Constructor of this class. + + :param function: for error computation. + """ + super(LagrangeRombergGridSliceContainerOld, self).__init__(function) + self.max_interpolation_support_points = 7 + + def get_interpolation_weight(self, support_points, basis_point, evaluation_point) -> float: + return self.get_langrange_basis(support_points, basis_point, evaluation_point) + + @staticmethod + def get_langrange_basis(support_points: Sequence[float], basis_point: float, evaluation_point: float) -> float: + """ + This method computes a lambda expression for the lagrange basis function at basis_point. + Specifically: Let support_points = [x_0, x_1, ... , x_n] + L = product( (evaluation_point - x_k) / (basis_point - x_k) ) + + :param support_points: an array of support points for the interpolation + :param basis_point: determines basis_point of this lagrange basis function + :param evaluation_point: determines at which point the basis is evaluated + :return: the evaluated basis function L + """ + + evaluated_basis = 1 + + for point in support_points: + if point != basis_point: + evaluated_basis *= (evaluation_point - point) / (basis_point - point) + + return evaluated_basis + + def plot_interpolation_basis_functions_with_evaluations(self, support_points, evaluation_point): + # Iterate over all dimensions and plot each basis + # plt.figure() + + columns = 2 + fig, axes = plt.subplots(figsize=(10, 15), nrows=math.ceil(len(support_points) / columns), ncols=columns) + fig.subplots_adjust(hspace=0.5) + fig.suptitle("Basis functions for the interpolation of the point {}".format(evaluation_point)) + + axes = axes.ravel() + + a = support_points[0] + b = support_points[-1] + + for i, basis_point in enumerate(support_points): + if self.function is not None: + evaluation_points = np.linspace(a, b, len(support_points) * 100) + axes[i].plot(evaluation_points, [self.function([point]) for point in evaluation_points], + linestyle=':', color="black") + + # Plot basis functions + for support_point in support_points: + evaluation_points = np.linspace(a, b, len(support_points) * 100) + axes[i].plot(evaluation_points, [self.get_langrange_basis(support_points, support_point, point) + for point in evaluation_points]) + + # Print interpolation grid + axes[i].plot(support_points, np.zeros(len(support_points)), 'bo', markersize=4, color="black") + + # Plot evaluation point + if evaluation_point is not None: + evaluation_point_value = self.get_langrange_basis(support_points, basis_point, evaluation_point) + axes[i].plot(evaluation_point, evaluation_point_value, 'bo', markersize=4, color="red") + axes[i].plot([evaluation_point, evaluation_point], [0, evaluation_point_value], ':', + color="red") + + axes[i].set_title("Evaluation using basis function of point {}".format(basis_point)) + + if len(support_points) < 15: + axes[i].set_xticks(support_points) + axes[i].set_xticklabels(support_points, fontsize=9) + + # Clear unused axes + for i in range(len(support_points), len(axes)): + fig.delaxes(axes[i]) + + fig.show() + + +class GlobalBasisGridRombergGridSliceContainerAdapter(InterpolationGridSliceContainer): + + """ + This class provides an interfaces for interpolating as many missing points as possible inside the container. + It uses the interpolation functionality from the GlobalBasisGrid classes. + + It executes a default Romberg method on this interpolated partial full grid. + For the interpolation we use the closest points to the left and right of this container. + """ + + def __init__(self, function: Function = None): + """ + Constructor of this class. + + :param function: for error computation. + """ + super(GlobalBasisGridRombergGridSliceContainerAdapter, self).__init__(function) + self.max_interpolation_support_points = 7 + self.interpolation_grid = None + + @abstractmethod + def construct_interpolation_grid(self, support_points): + """ + + Parameters + ---------- + support_points + + Returns GlobalBasisGrid + ------- + + """ + pass + + def get_interpolation_weight(self, support_points, basis_point, evaluation_point) -> float: + # Interpolation grid is computed beforehand + interpolation_grid = self.construct_interpolation_grid(support_points) + self.interpolation_grid = interpolation_grid + + support_point_grid_levels = InterpolationGridSliceContainer.get_levels_for_support_point_grid(support_points) + interpolation_grid.set_grid([support_points], [support_point_grid_levels]) + + # Compute weight + assert len(interpolation_grid.basis) == 1 + basis = interpolation_grid.basis[0] + assert len(basis) == len(support_points) + + i = support_points.index(basis_point) # Element must be found + weight = basis[i](evaluation_point) + # TODO which basis function should be chosen? + + return weight + + def plot_interpolation_basis_functions(self, support_points, interpolation_point): + self.interpolation_grid.plot_basis_functions(support_points, interpolation_point, self.function) + + def plot_interpolation_basis_functions_with_evaluations(self, support_points, evaluation_point): + self.interpolation_grid.plot_basis_functions_with_evaluations(support_points, evaluation_point, self.function) + + +class BSplineRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter): + + """ + This class interpolates (BSplines) as many missing points as possible inside the container + It executes a default Romberg method on this interpolated partial full grid. + + For the interpolation we use the closest points to the left and right of this container. + """ + + def __init__(self, function: Function = None): + """ + Constructor of this class. + + :param function: for error computation. + """ + super(BSplineRombergGridSliceContainer, self).__init__(function) + self.max_interpolation_support_points = 7 + + def construct_interpolation_grid(self, support_points): + # Prevent circular dependencies + from sparseSpACE.Grid import GlobalBSplineGrid + + # Interpolation grid is computed beforehand + return GlobalBSplineGrid([support_points[0]], [support_points[-1]], boundary=True) + + +class LagrangeRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter): + + """ + This class interpolates (Hierarchical Lagrange) as many missing points as possible inside the container + It executes a default Romberg method on this interpolated partial full grid. + + For the interpolation we use the closest points to the left and right of this container. + """ + + def __init__(self, function: Function = None): + """ + Constructor of this class. + + :param function: for error computation. + """ + super(LagrangeRombergGridSliceContainer, self).__init__(function) + self.max_interpolation_support_points = 7 + + def construct_interpolation_grid(self, support_points): + # Prevent circular dependencies + from sparseSpACE.Grid import GlobalLagrangeGrid + + # Interpolation grid is computed beforehand + return GlobalLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) + + +class HierarchicalLagrangeRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter): + + """ + This class interpolates (Hierarchical Lagrange) as many missing points as possible inside the container + It executes a default Romberg method on this interpolated partial full grid. + + For the interpolation we use the closest points to the left and right of this container. + """ + + def __init__(self, function: Function = None): + """ + Constructor of this class. + + :param function: for error computation. + """ + super(HierarchicalLagrangeRombergGridSliceContainer, self).__init__(function) + self.max_interpolation_support_points = 7 + + def construct_interpolation_grid(self, support_points): + # Prevent circular dependencies + from sparseSpACE.Grid import GlobalHierarchicalLagrangeGrid + + # Interpolation grid is computed beforehand + return GlobalHierarchicalLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) + class TrapezoidalGridSlice(ExtrapolationGridSlice): """ @@ -2592,11 +2856,16 @@ def __init__(self, slice_container_version: SliceContainerVersion): def get_grid_slice_container(self, function: Function = None) -> ExtrapolationGridSliceContainer: if self.slice_container_version == SliceContainerVersion.ROMBERG_DEFAULT: return RombergGridSliceContainer(function) - elif self.slice_container_version == SliceContainerVersion.LAGRANGE_ROMBERG or \ - self.slice_container_version == SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG: + elif self.slice_container_version == SliceContainerVersion.LAGRANGE_ROMBERG_OLD: + return LagrangeRombergGridSliceContainerOld(function) + elif self.slice_container_version == SliceContainerVersion.LAGRANGE_ROMBERG: return LagrangeRombergGridSliceContainer(function) elif self.slice_container_version == SliceContainerVersion.SIMPSON_ROMBERG: return SimpsonRombergGridSliceContainer(function) + elif self.slice_container_version == SliceContainerVersion.BSPLINE_ROMBERG: + return BSplineRombergGridSliceContainer(function) + elif self.slice_container_version == SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG: + return HierarchicalLagrangeRombergGridSliceContainer(function) else: raise RuntimeError("Wrong ContainerVersion provided.") diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index 9e589fa..d5836a6 100644 --- a/sparseSpACE/Grid.py +++ b/sparseSpACE/Grid.py @@ -10,7 +10,8 @@ from sparseSpACE.ComponentGridInfo import * from typing import Callable, Tuple, Sequence from sparseSpACE.Extrapolation import ExtrapolationGrid, SliceGrouping, SliceVersion, SliceContainerVersion, \ - BalancedExtrapolationGrid + BalancedExtrapolationGrid, GridVersion + # the grid class provides basic functionalities for an abstract grid class Grid(object): @@ -348,7 +349,7 @@ def get_basis(self, d: int, index: int): class LagrangeGrid(BasisGrid): - def __init__(self, a: float, b: float, boundary: bool=True, p: int=3, modified_basis: bool=False): + def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, p: int=3, modified_basis: bool=False): self.boundary = boundary self.a = a self.b = b @@ -366,7 +367,7 @@ def is_high_order_grid(self) -> bool: class LagrangeGrid1D(Grid1d): - def __init__(self, a: float, b: float, boundary: bool=True, p: int=3, modified_basis: bool=False): + def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, p: int=3, modified_basis: bool=False): super().__init__(a=a, b=b, boundary=boundary) self.p = p # max order of lagrange polynomials self.coords_gauss, self.weights_gauss = legendre.leggauss(int(self.p / 2) + 1) @@ -482,7 +483,7 @@ def get_parent(self, p: int, grid_1D, grid_levels): class BSplineGrid(BasisGrid): - def __init__(self, a: float, b: float, boundary: bool=True, p: int=3, modified_basis: bool=False): + def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, p: int=3, modified_basis: bool=False): self.boundary = boundary self.a = a self.b = b @@ -499,7 +500,7 @@ def is_high_order_grid(self) -> bool: return self.p > 1 class BSplineGrid1D(Grid1d): - def __init__(self, a: float, b: float, boundary: bool=True, p: int=3, modified_basis: bool=False): + def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, p: int=3, modified_basis: bool=False): super().__init__(a=a, b=b, boundary=boundary) self.p = p #spline order assert p % 2 == 1 @@ -1189,9 +1190,10 @@ def compute_1D_quad_weights(self, grid_1D: Sequence[float], a: float, b: float, class GlobalRombergGrid(GlobalGrid): def __init__(self, a, b, boundary=True, modified_basis=False, do_cache=True, + grid_version=GridVersion.DEFAULT, + container_version=SliceContainerVersion.ROMBERG_DEFAULT, slice_grouping=SliceGrouping.UNIT, - slice_version=SliceVersion.ROMBERG_DEFAULT, - container_version=SliceContainerVersion.ROMBERG_DEFAULT): + slice_version=SliceVersion.ROMBERG_DEFAULT): self.boundary = boundary self.integrator = IntegratorArbitraryGridScalarProduct(self) self.a = a @@ -1204,9 +1206,10 @@ def __init__(self, a, b, boundary=True, modified_basis=False, assert not(modified_basis) assert boundary + self.grid_version = grid_version + self.container_version = container_version self.slice_grouping = slice_grouping self.slice_version = slice_version - self.container_version = container_version self.weight_cache = {} @@ -1223,11 +1226,16 @@ def compute_1D_quad_weights(self, grid_1D: Sequence[float], a: float, b: float, if self.do_cache and key in self.weight_cache: return self.weight_cache[key] - romberg_grid = ExtrapolationGrid(slice_grouping=self.slice_grouping, - slice_version=self.slice_version, - container_version=self.container_version) + romberg_grid = ExtrapolationGrid(grid_version=self.grid_version, + container_version=self.container_version, + slice_grouping=self.slice_grouping, + slice_version=self.slice_version) romberg_grid.set_grid(grid_1D, grid_levels_1D) + + # TODO Adapt function to one dimensional grid stripe (1D), because the function is multidimensional in the general case + # romberg_grid.set_function(self.function) + weights = romberg_grid.get_weights() if self.do_cache: @@ -1317,11 +1325,11 @@ def interpolate_grid(self, grid_points_for_dims: Sequence[Sequence[float]], comp evaluations = np.empty(self.dim, dtype=object) for d in range(self.dim): points_d = grid_points_for_dims[d] - evaluations1D = np.zeros((len(points_d), (len(self.splines[d])))) - for i, spline in enumerate(self.splines[d]): + evaluations1D = np.zeros((len(points_d), (len(self.basis[d])))) + for i, basis in enumerate(self.basis[d]): for j, p in enumerate(points_d): #print(p, spline(p)) - evaluations1D[j, i] = spline(p) + evaluations1D[j, i] = basis(p) evaluations[d] = evaluations1D #print(evaluations) @@ -1344,6 +1352,80 @@ def interpolate_grid(self, grid_points_for_dims: Sequence[Sequence[float]], comp #print(results) return results + def plot_basis_functions(self, support_points=None, interpolation_point=None, function=None, filename=None): + if self.basis is None: + print("No basis functions constructed") + return + + # Iterate over all dimensions and plot each basis + for d in range(self.dim): + plt.figure() + + if function is not None: + evaluation_points = np.linspace(self.a[d], self.b[d], len(support_points) * 100) + plt.plot(evaluation_points, [function([point]) for point in evaluation_points], + linestyle=':', color="black") + + # Plot basis functions + for i, basis in enumerate(self.basis[d]): + evaluation_points = np.linspace(self.a[d], self.b[d], len(basis.knots) * 100) + plt.plot(evaluation_points, [basis(point) for point in evaluation_points]) + + # Print interpolation grid + plt.plot(support_points, np.zeros(len(support_points)), 'bo', markersize=6, color="dimgray") + plt.title("Basis functions for the interpolation of point {}".format(interpolation_point)) + + plt.show() + + def plot_basis_functions_with_evaluations(self, support_points, evaluation_point, function=None): + # Iterate over all dimensions and plot each basis + # plt.figure() + + for d in range(self.dim): + # TODO plot dimension d + columns = 2 + fig, axes = plt.subplots(figsize=(10, 15), nrows=math.ceil(len(support_points) / columns), ncols=columns) + fig.subplots_adjust(hspace=0.5) + fig.suptitle("Basis functions for the interpolation of the point {}".format(evaluation_point)) + + axes = axes.ravel() + + a = support_points[0] + b = support_points[-1] + + for i, basis in enumerate(self.basis[d]): + if function is not None: + evaluation_points = np.linspace(a, b, len(support_points) * 100) + axes[i].plot(evaluation_points, [function([point]) for point in evaluation_points], + linestyle=':', color="black") + + # Plot basis functions + for j, basis_f in enumerate(self.basis[d]): + evaluation_points = np.linspace(self.a[d], self.b[d], len(basis.knots) * 100) + axes[i].plot(evaluation_points, [basis_f(point) for point in evaluation_points]) + + # Print interpolation grid + axes[i].plot(support_points, np.zeros(len(support_points)), 'bo', markersize=4, color="black") + + # Plot evaluation point + if evaluation_point is not None: + evaluation_point_value = basis(evaluation_point) + axes[i].plot(evaluation_point, evaluation_point_value, 'bo', markersize=4, color="red") + axes[i].plot([evaluation_point, evaluation_point], [0, evaluation_point_value], ':', + color="red") + + axes[i].set_title("Evaluation using basis function of point {}".format(support_points[i])) + + if len(support_points) < 15: + axes[i].set_xticks(support_points) + axes[i].set_xticklabels(support_points, fontsize=9) + + # Clear unused axes + for i in range(len(support_points), len(axes)): + fig.delaxes(axes[i]) + + fig.show() + class GlobalBSplineGrid(GlobalBasisGrid): def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, modified_basis: bool=False, p: int=3, chebyshev=False): @@ -1468,6 +1550,66 @@ def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, assert not(modified_basis) or not(boundary) self.surplus_values = {} + def compute_1D_quad_weights(self, grid_1D: Sequence[float], a: float, b: float, d: int, + grid_levels_1D: Sequence[int]=None) -> Sequence[float]: + max_level = max(grid_levels_1D) + grid_1D = list(grid_1D) + grid_levels_1D = np.asarray(grid_levels_1D) + level_coordinate_array = [[] for l in range(max_level+1)] + for i, p in enumerate(grid_1D): + level_coordinate_array[grid_levels_1D[i]].append(p) + + #print(level_coordinate_array) + self.start = a + self.end = b + # level = 0 + weights = np.zeros(len(grid_1D)) + starting_level = 0 if self.boundary else 1 + + for l in range(starting_level, max_level + 1): + for i in range(len(level_coordinate_array[l])): + x_basis = level_coordinate_array[l][i] + knots = grid_1D + + if self.modified_basis: + # spline = LagrangeBasisRestrictedModified(self.p, knots.index(x_basis), knots, a, b, l) + RuntimeError("Modified Basis is not supported") + + spline = LagrangeBasis(self.p, knots.index(x_basis), knots) + index = grid_1D.index(x_basis) + self.basis[d][index] = spline + weights[index] = self._get_spline_integral(spline, d) + + if not self.boundary: + self.basis[d] = self.basis[d][1:-1] + + for basis in self.basis[d]: + assert basis is not None + + return weights + + def _get_spline_integral(self, spline, d=None): + if self.coords_gauss is None: + self.coords_gauss, self.weights_gauss = legendre.leggauss(int(self.p/2) + 1) + + return spline.get_integral(self.start, self.end, self.coords_gauss, self.weights_gauss) + + +class GlobalHierarchicalLagrangeGrid(GlobalBasisGrid): + def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, modified_basis: bool=False, p: int=3): + self.boundary = boundary + self.integrator = IntegratorHierarchicalBasisFunctions(self) + self.a = a + self.b = b + self.dim = len(a) + self.length = np.array(b) - np.array(a) + self.modified_basis = modified_basis + assert p >= 1 + self.p = p + self.coords_gauss = None + assert not(modified_basis) or not(boundary) + self.surplus_values = {} + def compute_1D_quad_weights(self, grid_1D: Sequence[float], a: float, b: float, d: int, grid_levels_1D: Sequence[int]=None) -> Sequence[float]: max_level = max(grid_levels_1D) grid_1D = list(grid_1D) @@ -1559,7 +1701,7 @@ def get_parent(self, point_position: float, grid_1D: Sequence[float], grid_level assert False -class GlobalLagrangeGridWeighted(GlobalLagrangeGrid): +class GlobalLagrangeGridWeighted(GlobalHierarchicalLagrangeGrid): def __init__(self, a: Sequence[float], b: Sequence[float], uq_operation, boundary: bool=True, modified_basis: bool=False, p: int=3): super().__init__(a, b, boundary=boundary, modified_basis=modified_basis, p=p) self.distributions = uq_operation.get_distributions() diff --git a/sparseSpACE/spatiallyAdaptiveSingleDimension2.py b/sparseSpACE/spatiallyAdaptiveSingleDimension2.py index fbfa2cb..0da65b1 100644 --- a/sparseSpACE/spatiallyAdaptiveSingleDimension2.py +++ b/sparseSpACE/spatiallyAdaptiveSingleDimension2.py @@ -1374,7 +1374,7 @@ def sum_up_volumes_for_point_completely_vectorized(self, child_infos: NodeInfo, # on this slice. def calculate_surplusses(self, grid_points: Sequence[Sequence[float]], children_indices: Sequence[Sequence[int]], component_grid: ComponentGridInfo): tol = 10**-84 - if isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalLagrangeGrid): + if isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalHierarchicalLagrangeGrid): # grid_values = np.empty((self.f.output_length(), np.prod(self.grid.numPoints))) # points = self.grid.getPoints() # for i, point in enumerate(points): @@ -1383,7 +1383,7 @@ def calculate_surplusses(self, grid_points: Sequence[Sequence[float]], children_ for d in range(0, self.dim): k=0 refinement_dim = self.refinement.get_refinement_container_for_dim(d) - if isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalLagrangeGrid): + if isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalHierarchicalLagrangeGrid): hierarchization_operator = HierarchizationLSG(self.grid) surplusses_1d = hierarchization_operator.hierarchize_poles_for_dim(np.array(grid_values.T), self.grid.numPoints, d) surplus_pole = np.zeros((self.operation.point_output_length(), self.grid.numPoints[d])) @@ -1393,14 +1393,14 @@ def calculate_surplusses(self, grid_points: Sequence[Sequence[float]], children_ while i < np.prod(self.grid.numPoints): surplus_pole[:,j] += np.sum(abs(surplusses_1d[:,i:i+stride])) #* weights[i:i+stride])) i += stride * self.grid.numPoints[d] - if not (isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalLagrangeGrid)) and len(children_indices[d]) > 0: + if not (isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalHierarchicalLagrangeGrid)) and len(children_indices[d]) > 0: #print(children_indices) volumes, evaluations = self.sum_up_volumes_for_point_completely_vectorized(child_infos=children_indices[d], grid_points=grid_points, d=d, component_grid=component_grid) for i, child_info in enumerate(children_indices[d]): left_parent = child_info.left_parent right_parent = child_info.right_parent child = child_info.child - if isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalLagrangeGrid): + if isinstance(self.grid_surplusses, GlobalBSplineGrid) or isinstance(self.grid_surplusses, GlobalHierarchicalLagrangeGrid): index_child = grid_points[d].index(child) - int(not(self.grid.boundary)) volume = surplus_pole[:, index_child] / np.prod(self.grid.numPoints) * self.grid.numPoints[d] * self.grid.weights[d][index_child] evaluations = np.prod(self.grid.numPoints) / self.grid.numPoints[d] diff --git a/test/test_ExtrapolationInterpolatingGrid.py b/test/test_ExtrapolationInterpolatingGrid.py index 5f76860..1b1791b 100644 --- a/test/test_ExtrapolationInterpolatingGrid.py +++ b/test/test_ExtrapolationInterpolatingGrid.py @@ -2,11 +2,9 @@ import sparseSpACE from sparseSpACE.Extrapolation import SliceVersion, ExtrapolationGrid, SliceContainerVersion, \ - LagrangeRombergGridSliceContainer, SliceGrouping + LagrangeRombergGridSliceContainer, SliceGrouping, GridVersion from sparseSpACE.Function import Polynomial1d -# ----------------------------------------------------------------------------------------------------------------- -# --- Interpolating Grid class TestExtrapolationInterpolationGrid(unittest.TestCase): def test_interpolated_points(self): @@ -14,9 +12,10 @@ def test_interpolated_points(self): grid_levels = [0, 3, 4, 2, 3, 1, 2, 3, 0] # Containers: [0, 0.125], [0.125, .0.25], [0.25, 0.5], [0.5, 0.75], [0.75, 0.875] [0.875, 1] - romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, force_balanced_refinement_tree=False) romberg_grid.set_grid(grid, grid_levels) containers = romberg_grid.slice_containers @@ -38,9 +37,10 @@ def test_grid_without_interpolated_points(self): actual_grid = [] actual_grid_levels = [] - romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, force_balanced_refinement_tree=False) romberg_grid.set_grid(grid, grid_levels) containers = romberg_grid.slice_containers @@ -68,9 +68,10 @@ def test_interpolation_support_points_by_level(self): grid = [0.0, 0.125, 0.1875, 0.25, 0.5, 0.625, 0.75, 1] grid_levels = [0, 3, 4, 2, 1, 3, 2, 0] - romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, force_balanced_refinement_tree=False) romberg_grid.set_grid(grid, grid_levels) containers = romberg_grid.slice_containers @@ -83,9 +84,10 @@ def test_interpolation_support_points_geometrically(self): grid = [0.0, 0.125, 0.1875, 0.25, 0.5, 0.625, 0.75, 1] grid_levels = [0, 3, 4, 2, 1, 3, 2, 0] - romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, force_balanced_refinement_tree=False) romberg_grid.set_grid(grid, grid_levels) containers = romberg_grid.slice_containers @@ -97,23 +99,6 @@ def test_interpolation_support_points_geometrically(self): self.assertEqual([0.625, 0.75, 1], containers[2].get_support_points_for_interpolation_geometrically(0.875, 3)) - # Exactness tests - def test_exactness(self): - grid = [0, 0.5, 0.625, 0.75, 1] - grid_levels = [0, 1, 3, 2, 0] - - function = Polynomial1d([1, 0, 0, 2]) - - romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, - container_version=SliceContainerVersion.LAGRANGE_ROMBERG, - force_balanced_refinement_tree=True) - romberg_grid.set_grid(grid, grid_levels) - actual_result = romberg_grid.integrate(function) - - expected_result = 1.5 - self.assertAlmostEqual(expected_result, actual_result) - def test_weight_count(self): grid = [0.0, 0.0078125, 0.015625, 0.0234375, 0.03125, 0.0390625, 0.046875, 0.0625, 0.078125, 0.09375, 0.109375, 0.125, @@ -126,9 +111,10 @@ def test_weight_count(self): grid_levels = [0, 7, 6, 7, 5, 7, 6, 4, 6, 5, 6, 3, 6, 5, 6, 4, 6, 5, 6, 2, 6, 5, 6, 4, 6, 5, 6, 3, 6, 5, 4, 5, 1, 5, 4, 5, 3, 5, 4, 5, 2, 5, 4, 5, 3, 5, 4, 0] - romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, force_balanced_refinement_tree=False) romberg_grid.set_grid(grid, grid_levels) weights = romberg_grid.get_weights() @@ -143,9 +129,10 @@ def test_full_grid_interpolation(self): function = Polynomial1d([1, 0, 0, 2]) - romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_FULL_GRID, + container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, slice_version=SliceVersion.ROMBERG_DEFAULT, - container_version=SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG, force_balanced_refinement_tree=False) romberg_grid.set_grid(grid, grid_levels) actual_result = romberg_grid.integrate(function) @@ -153,6 +140,76 @@ def test_full_grid_interpolation(self): expected_result = 1.5 self.assertAlmostEqual(expected_result, actual_result) + def test_exactness_on_balanced_adaptive_interpolated_lagrange_grid(self): + grid = [0, 0.5, 0.625, 0.75, 1] + grid_levels = [0, 1, 3, 2, 0] + + function = Polynomial1d([1, 0, 0, 2]) + + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, + container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, + force_balanced_refinement_tree=True) + romberg_grid.set_grid(grid, grid_levels) + actual_result = romberg_grid.integrate(function) + + expected_result = 1.5 + self.assertAlmostEqual(expected_result, actual_result) + + def test_exactness_on_unbalanced_adaptive_interpolated_lagrange_grid(self): + grid = [0, 0.5, 0.625, 0.75, 1] + grid_levels = [0, 1, 3, 2, 0] + + function = Polynomial1d([1, 0, 0, 2]) + + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, + container_version=SliceContainerVersion.LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, + force_balanced_refinement_tree=False) + romberg_grid.set_grid(grid, grid_levels) + actual_result = romberg_grid.integrate(function) + + expected_result = 1.5 + + # TODO works only with force_balanced_refinement_tree = True + # self.assertAlmostEqual(expected_result, actual_result) + + def test_exactness_on_unbalanced_adaptive_interpolated_bspline_grid(self): + grid = [0, 0.5, 0.625, 0.75, 1] + grid_levels = [0, 1, 3, 2, 0] + + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, + container_version=SliceContainerVersion.BSPLINE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, + force_balanced_refinement_tree=False) + romberg_grid.set_grid(grid, grid_levels) + function = Polynomial1d([1, 0, 0, 2]) + + actual_result = romberg_grid.integrate(function) + expected_result = 1.5 + + # self.assertAlmostEqual(expected_result, actual_result) + + def test_exactness_on_unbalanced_adaptive_interpolated_hierarchical_lagrange_grid(self): + grid = [0, 0.5, 0.625, 0.75, 1] + grid_levels = [0, 1, 3, 2, 0] + + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, + container_version=SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG, + slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, + force_balanced_refinement_tree=True) + romberg_grid.set_grid(grid, grid_levels) + function = Polynomial1d([1, 0, 0, 2]) + + actual_result = romberg_grid.integrate(function) + expected_result = 1.5 + + # self.assertAlmostEqual(expected_result, actual_result) + # ----------------------------------------------------------------------------------------------------------------- # --- Unit Test diff --git a/test/test_Hierarchization.py b/test/test_Hierarchization.py index 3c83dad..d123964 100644 --- a/test/test_Hierarchization.py +++ b/test/test_Hierarchization.py @@ -14,7 +14,7 @@ def test_interpolation(self): a = -3 b = 6 for d in range(2, 5): - grid = GlobalLagrangeGrid(a*np.ones(d), b*np.ones(d), boundary= True, modified_basis = False, p = 1) + grid = GlobalHierarchicalLagrangeGrid(a*np.ones(d), b*np.ones(d), boundary= True, modified_basis = False, p = 1) for l in range(7 - d): f = FunctionLinear([10*(i+1) for i in range(d)]) grid_points = [np.linspace(a,b,2**(l+i)+ 1) for i in range(d)] diff --git a/test/test_Integrator.py b/test/test_Integrator.py index e055e6e..0f63aeb 100644 --- a/test/test_Integrator.py +++ b/test/test_Integrator.py @@ -36,7 +36,7 @@ def test_integrate_basis_functions(self): b = 6 for d in range(1, 5): for p in range(1,7): - grid = GlobalLagrangeGrid(a*np.ones(d), b*np.ones(d), boundary= True, modified_basis = False, p=p) + grid = GlobalHierarchicalLagrangeGrid(a*np.ones(d), b*np.ones(d), boundary= True, modified_basis = False, p=p) for l in range(p - 1, 8 - d): f = FunctionPolynomial([10*(i+1) for i in range(d)], degree=p) grid_points = [np.linspace(a,b,2**l+ 1) for _ in range(d)] diff --git a/test/test_spatiallyAdaptiveSingleDimension2.py b/test/test_spatiallyAdaptiveSingleDimension2.py index 71d4897..12430f3 100644 --- a/test/test_spatiallyAdaptiveSingleDimension2.py +++ b/test/test_spatiallyAdaptiveSingleDimension2.py @@ -31,7 +31,7 @@ def test_integrate(self): a = -3 b = 6 for d in range(2, 4): - grid = GlobalLagrangeGrid(a * np.ones(d), b * np.ones(d), boundary=True, modified_basis=False, p=2) + grid = GlobalHierarchicalLagrangeGrid(a * np.ones(d), b * np.ones(d), boundary=True, modified_basis=False, p=2) f = FunctionPolynomial([10 * (i + 1) for i in range(d)], degree=2) operation = Integration(f, grid=grid, dim=d, reference_solution=f.getAnalyticSolutionIntegral(a * np.ones(d), b * np.ones(d))) errorOperator = ErrorCalculatorSingleDimVolumeGuided() @@ -104,7 +104,7 @@ def test_interpolate(self): a = -1 b = 6 for d in range(2, 5): - grid = GlobalLagrangeGrid(a * np.ones(d), b * np.ones(d), boundary=True, modified_basis=False, p=2) + grid = GlobalHierarchicalLagrangeGrid(a * np.ones(d), b * np.ones(d), boundary=True, modified_basis=False, p=2) f = FunctionPolynomial([(i + 1) for i in range(d)], degree=2) operation = Integration(f, grid=grid, dim=d, reference_solution=f.getAnalyticSolutionIntegral(a * np.ones(d), b * np.ones(d))) errorOperator = ErrorCalculatorSingleDimVolumeGuided()