From 28fdf00b32010b48d01bf15ebc4c4b10516139af Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Fri, 16 Oct 2020 14:31:27 +0200 Subject: [PATCH 01/14] Add BSpline Extrapolation. --- ipynb/Extrapolation/Extrapolation_Test.ipynb | 8 +- src/Extrapolation.py | 209 ++++++++++++++----- src/Grid.py | 6 +- 3 files changed, 162 insertions(+), 61 deletions(-) diff --git a/ipynb/Extrapolation/Extrapolation_Test.ipynb b/ipynb/Extrapolation/Extrapolation_Test.ipynb index 32c64ad..85115b9 100644 --- a/ipynb/Extrapolation/Extrapolation_Test.ipynb +++ b/ipynb/Extrapolation/Extrapolation_Test.ipynb @@ -271,7 +271,7 @@ " 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", @@ -306,13 +306,13 @@ "pycharm": { "stem_cell": { "cell_type": "raw", + "source": [], "metadata": { "collapsed": false - }, - "source": [] + } } } }, "nbformat": 4, "nbformat_minor": 1 -} +} \ No newline at end of file diff --git a/src/Extrapolation.py b/src/Extrapolation.py index c3d7445..1d019a5 100644 --- a/src/Extrapolation.py +++ b/src/Extrapolation.py @@ -543,6 +543,7 @@ class SliceContainerVersion(Enum): 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 + BSPLINE_ROMBERG = 5 class ExtrapolationGrid: @@ -593,7 +594,8 @@ def __init__(self, def interpolation_is_enabled(self) -> bool: return self.container_version == SliceContainerVersion.LAGRANGE_ROMBERG or \ - self.container_version == SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG + self.container_version == SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG or \ + self.container_version == SliceContainerVersion.BSPLINE_ROMBERG # Integration def integrate(self, function: Function = None) -> float: @@ -2254,7 +2256,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 +2278,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 +2414,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 @@ -2473,10 +2473,109 @@ def populate_interpolated_point_weights(self, factory: RombergWeights, normalized_max_level) for support_point in support_points: - interpolated_extrapolated_weight = weight * self.get_langrange_basis(support_points, support_point, - interp_point) + interpolated_extrapolated_weight = weight * self.get_interpolation_weight(support_points, support_point, + interp_point) weight_dictionary[support_point].append(interpolated_extrapolated_weight) + 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_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 + + +class BSplineRombergGridSliceContainer(InterpolationGridSliceContainer): + + """ + This class interpolates (BSplines) 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(BSplineRombergGridSliceContainer, self).__init__(function) + self.max_interpolation_support_points = 7 + + def get_interpolation_weight(self, support_points, basis_point, evaluation_point) -> float: + # Prevent circular dependencies + from Grid import GlobalBSplineGrid + from ComponentGridInfo import ComponentGridInfo + + # Interpolation grid is computed beforehand + bspline_grid = GlobalBSplineGrid([support_points[0]], [support_points[-1]], boundary=True) + + support_point_grid_levels = InterpolationGridSliceContainer.get_levels_for_support_point_grid(support_points) + bspline_grid.set_grid([support_points], [support_point_grid_levels]) + + level_vec = [1] + # component_grid = ComponentGridInfo(level_vec, 1) + # print("Check surplusses") + + # interpolation_indicator = self.get_interpolated_grid_points_indicator() + # interp_points = [self.get_grid()[i] for i, indicator in enumerate(interpolation_indicator) if indicator] + + # bspline_grid.integrate(self.function, level_vec, self.left_point * np.ones(1), self.right_point * np.ones(1)) + bspline_grid.integrate(self.function, level_vec, support_points[0] * np.ones(1), support_points[-1] * np.ones(1)) + + # Compute weight + assert len(bspline_grid.basis) == 1 + basis = bspline_grid.basis[0] + assert len(basis) == len(support_points) + + i = support_points.index(basis_point) # Element must be found + + weight = basis[i](evaluation_point) + + # interp = bspline_grid.interpolate([(point, 0) for point in interp_points], component_grid) + + return weight + class TrapezoidalGridSlice(ExtrapolationGridSlice): """ @@ -2597,6 +2696,8 @@ def get_grid_slice_container(self, function: Function = None) -> ExtrapolationGr 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) else: raise RuntimeError("Wrong ContainerVersion provided.") diff --git a/src/Grid.py b/src/Grid.py index 4aa1c77..70b4aa0 100644 --- a/src/Grid.py +++ b/src/Grid.py @@ -1271,11 +1271,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) From 8214a270573a07d7201e07d3087aa50565eae117 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Fri, 16 Oct 2020 15:35:15 +0200 Subject: [PATCH 02/14] Add hierarchical lagrange interpolated extrapolation. Generalized adapter for interpolated extrapolation using GlobalBasisGrid. --- src/Extrapolation.py | 106 ++++++++++++++---- .../test_ExtrapolationInterpolatingGrid.py | 32 ++++++ 2 files changed, 114 insertions(+), 24 deletions(-) diff --git a/src/Extrapolation.py b/src/Extrapolation.py index 1d019a5..9606193 100644 --- a/src/Extrapolation.py +++ b/src/Extrapolation.py @@ -544,6 +544,7 @@ class SliceContainerVersion(Enum): 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 BSPLINE_ROMBERG = 5 + HIERARCHICAL_LAGRANGE_ROMBERG = 6 class ExtrapolationGrid: @@ -595,7 +596,8 @@ def __init__(self, def interpolation_is_enabled(self) -> bool: return self.container_version == SliceContainerVersion.LAGRANGE_ROMBERG or \ self.container_version == SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG or \ - self.container_version == SliceContainerVersion.BSPLINE_ROMBERG + self.container_version == SliceContainerVersion.BSPLINE_ROMBERG or \ + self.container_version == SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG # Integration def integrate(self, function: Function = None) -> float: @@ -2524,12 +2526,13 @@ def get_langrange_basis(support_points: Sequence[float], basis_point: float, eva return evaluated_basis -class BSplineRombergGridSliceContainer(InterpolationGridSliceContainer): +class GlobalBasisGridRombergGridSliceContainerAdapter(InterpolationGridSliceContainer): """ - This class interpolates (BSplines) as many missing points as possible inside the container and executes a default - Romberg method on this interpolated partial full grid. + 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. """ @@ -2539,44 +2542,97 @@ def __init__(self, function: Function = None): :param function: for error computation. """ - super(BSplineRombergGridSliceContainer, self).__init__(function) + super(GlobalBasisGridRombergGridSliceContainerAdapter, self).__init__(function) self.max_interpolation_support_points = 7 - def get_interpolation_weight(self, support_points, basis_point, evaluation_point) -> float: - # Prevent circular dependencies - from Grid import GlobalBSplineGrid - from ComponentGridInfo import ComponentGridInfo + @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 - bspline_grid = GlobalBSplineGrid([support_points[0]], [support_points[-1]], boundary=True) + interpolation_grid = self.construct_interpolation_grid(support_points) support_point_grid_levels = InterpolationGridSliceContainer.get_levels_for_support_point_grid(support_points) - bspline_grid.set_grid([support_points], [support_point_grid_levels]) + interpolation_grid.set_grid([support_points], [support_point_grid_levels]) level_vec = [1] - # component_grid = ComponentGridInfo(level_vec, 1) - # print("Check surplusses") - - # interpolation_indicator = self.get_interpolated_grid_points_indicator() - # interp_points = [self.get_grid()[i] for i, indicator in enumerate(interpolation_indicator) if indicator] - - # bspline_grid.integrate(self.function, level_vec, self.left_point * np.ones(1), self.right_point * np.ones(1)) - bspline_grid.integrate(self.function, level_vec, support_points[0] * np.ones(1), support_points[-1] * np.ones(1)) + interpolation_grid.integrate(self.function, level_vec, + support_points[0] * np.ones(1), support_points[-1] * np.ones(1)) # Compute weight - assert len(bspline_grid.basis) == 1 - basis = bspline_grid.basis[0] + 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) - # interp = bspline_grid.interpolate([(point, 0) for point in interp_points], component_grid) - return weight +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 Grid import GlobalBSplineGrid + + # Interpolation grid is computed beforehand + return GlobalBSplineGrid([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 Grid import GlobalLagrangeGrid + + # Interpolation grid is computed beforehand + return GlobalLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) + + class TrapezoidalGridSlice(ExtrapolationGridSlice): """ This type of GridSlice executes no extrapolation. It produces default trapezoidal weights. @@ -2698,6 +2754,8 @@ def get_grid_slice_container(self, function: Function = None) -> ExtrapolationGr 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/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py b/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py index af1b0b4..80bd25b 100644 --- a/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py +++ b/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py @@ -154,6 +154,38 @@ def test_full_grid_interpolation(self): expected_result = 1.5 self.assertAlmostEqual(expected_result, actual_result) + def test_bspline_interpolation(self): + grid = [0, 0.5, 0.625, 0.75, 1] + grid_levels = [0, 1, 3, 2, 0] + + romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, + container_version=SliceContainerVersion.BSPLINE_ROMBERG, + 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_hierarchical_lagrange_interpolation(self): + grid = [0, 0.5, 0.625, 0.75, 1] + grid_levels = [0, 1, 3, 2, 0] + + romberg_grid = ExtrapolationGrid(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + slice_version=SliceVersion.ROMBERG_DEFAULT, + container_version=SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG, + 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) + # ----------------------------------------------------------------------------------------------------------------- # --- Unit Test From 9708bf34223369adc8b17a8137e2b34ea44df751 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Tue, 27 Oct 2020 13:18:09 +0100 Subject: [PATCH 03/14] Fix interpolation adapter for extrapolation --- ipynb/Extrapolation/Extrapolation_Test.ipynb | 32 +++++++++++++++---- src/Extrapolation.py | 10 ++++-- src/Grid.py | 4 +++ .../test_ExtrapolationInterpolatingGrid.py | 16 ++++++++++ 4 files changed, 53 insertions(+), 9 deletions(-) diff --git a/ipynb/Extrapolation/Extrapolation_Test.ipynb b/ipynb/Extrapolation/Extrapolation_Test.ipynb index 85115b9..de43bdb 100644 --- a/ipynb/Extrapolation/Extrapolation_Test.ipynb +++ b/ipynb/Extrapolation/Extrapolation_Test.ipynb @@ -42,11 +42,11 @@ " force_balanced_refinement_tree=force_balanced_refinement_tree)\n", "\n", "# Settings\n", - "dim = 5\n", + "dim = 2\n", "a = np.zeros(dim)\n", "b = np.ones(dim)\n", "max_tol = 20 # set to max, for better comparison\n", - "max_evaluations = 10 ** 5\n", + "max_evaluations = 10 ** 3\n", "evaluation_points = evaluation_points(a, b, dim)\n", "functions = []\n", "sasd_version = 6\n", @@ -93,9 +93,9 @@ " 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", + " f_genz_disc\n", "]\n", "\n", "# Iterate over all functions\n", @@ -202,6 +202,22 @@ " container_version=SliceContainerVersion.LAGRANGE_FULL_GRID_ROMBERG,\n", " force_balanced_refinement_tree=False)\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", + " \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", + " \n", " # Romberg Constant Subtraction\n", " r_constant_subtraction = get_single_dim_romberg(a, b, f, reference_solution,\n", " sasd_version=sasd_version,\n", @@ -255,16 +271,20 @@ " # (r_grouped_optimized_trapezoidal_balanced, 1, 2, errorOperator, 'Extrapolation Grid (Optimized, Trapezoid, Romberg, Balanced)'),\n", " # \n", " # (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", + " # (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_full_grid, 1, 2, errorOperator, 'Extrapolation Grid (Grouped Optimized, Romberg, Lagrange Full Romberg)'),\n", + " \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", diff --git a/src/Extrapolation.py b/src/Extrapolation.py index 9606193..cb200ac 100644 --- a/src/Extrapolation.py +++ b/src/Extrapolation.py @@ -2566,9 +2566,13 @@ def get_interpolation_weight(self, support_points, basis_point, evaluation_point support_point_grid_levels = InterpolationGridSliceContainer.get_levels_for_support_point_grid(support_points) interpolation_grid.set_grid([support_points], [support_point_grid_levels]) - level_vec = [1] - interpolation_grid.integrate(self.function, level_vec, - support_points[0] * np.ones(1), support_points[-1] * np.ones(1)) + # level_vec = [max(support_point_grid_levels)] + # left_boundary = support_points[0] + # right_boundary = support_points[-1] + + # TODO Adapt function to one dimensional stripe (1D) + # interpolation_grid.integrate(self.function, level_vec, + # left_boundary * np.ones(1), right_boundary * np.ones(1)) # Compute weight assert len(interpolation_grid.basis) == 1 diff --git a/src/Grid.py b/src/Grid.py index 6f99b29..7892eb3 100644 --- a/src/Grid.py +++ b/src/Grid.py @@ -1234,6 +1234,10 @@ def compute_1D_quad_weights(self, grid_1D: Sequence[float], a: float, b: float, container_version=self.container_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: diff --git a/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py b/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py index 80bd25b..8d9e3b4 100644 --- a/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py +++ b/test/Extrapolation/test_ExtrapolationInterpolatingGrid.py @@ -154,6 +154,22 @@ def test_full_grid_interpolation(self): expected_result = 1.5 self.assertAlmostEqual(expected_result, actual_result) + def test_langrange_interpolation(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=False) + 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_bspline_interpolation(self): grid = [0, 0.5, 0.625, 0.75, 1] grid_levels = [0, 1, 3, 2, 0] From a47abfacd87684b485d3b238fb97320b9860af42 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Tue, 27 Oct 2020 13:18:37 +0100 Subject: [PATCH 04/14] Remove comments --- src/Extrapolation.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/Extrapolation.py b/src/Extrapolation.py index cb200ac..63fefa5 100644 --- a/src/Extrapolation.py +++ b/src/Extrapolation.py @@ -2566,14 +2566,6 @@ def get_interpolation_weight(self, support_points, basis_point, evaluation_point support_point_grid_levels = InterpolationGridSliceContainer.get_levels_for_support_point_grid(support_points) interpolation_grid.set_grid([support_points], [support_point_grid_levels]) - # level_vec = [max(support_point_grid_levels)] - # left_boundary = support_points[0] - # right_boundary = support_points[-1] - - # TODO Adapt function to one dimensional stripe (1D) - # interpolation_grid.integrate(self.function, level_vec, - # left_boundary * np.ones(1), right_boundary * np.ones(1)) - # Compute weight assert len(interpolation_grid.basis) == 1 basis = interpolation_grid.basis[0] From 16fc5e173afbfb3cb6eef11456aa0b4fd7e79143 Mon Sep 17 00:00:00 2001 From: Pascal Date: Sun, 27 Dec 2020 14:18:41 +0100 Subject: [PATCH 05/14] Fix import bug --- sparseSpACE/Extrapolation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index 5ee70ba..88df763 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -2597,7 +2597,7 @@ def __init__(self, function: Function = None): def construct_interpolation_grid(self, support_points): # Prevent circular dependencies - from Grid import GlobalBSplineGrid + from sparseSpACE.Grid import GlobalBSplineGrid # Interpolation grid is computed beforehand return GlobalBSplineGrid([support_points[0]], [support_points[-1]], boundary=True) @@ -2623,7 +2623,7 @@ def __init__(self, function: Function = None): def construct_interpolation_grid(self, support_points): # Prevent circular dependencies - from Grid import GlobalLagrangeGrid + from sparseSpACE.Grid import GlobalLagrangeGrid # Interpolation grid is computed beforehand return GlobalLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) From d9b9aa937bfd4dbac8bc4010668cc5a8f223f0fc Mon Sep 17 00:00:00 2001 From: Pascal Date: Sun, 3 Jan 2021 19:17:56 +0100 Subject: [PATCH 06/14] Update Extrapolation grid: Differentiate between no/subgrid/full interpolation --- ipynb/Extrapolation/Extrapolation_Test.ipynb | 34 ++++---- sparseSpACE/Extrapolation.py | 55 +++++++----- sparseSpACE/Grid.py | 18 ++-- sparseSpACE/TestInterpolationAdapter.py | 1 + test/test_ExtrapolationInterpolatingGrid.py | 91 +++++++++++--------- 5 files changed, 110 insertions(+), 89 deletions(-) create mode 100644 sparseSpACE/TestInterpolationAdapter.py diff --git a/ipynb/Extrapolation/Extrapolation_Test.ipynb b/ipynb/Extrapolation/Extrapolation_Test.ipynb index 7f7d5c0..6c25afc 100644 --- a/ipynb/Extrapolation/Extrapolation_Test.ipynb +++ b/ipynb/Extrapolation/Extrapolation_Test.ipynb @@ -5,8 +5,8 @@ "execution_count": null, "metadata": { "pycharm": { - "is_executing": false, - "name": "#%%\n" + "name": "#%%\n", + "is_executing": true } }, "outputs": [], @@ -31,11 +31,12 @@ "# 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", + " grid_version=grid_version,\n", + " container_version=container_version,\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", @@ -194,14 +195,16 @@ " 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", @@ -209,7 +212,8 @@ " 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", + " 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", @@ -217,7 +221,8 @@ " 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", + " 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", @@ -323,15 +328,6 @@ "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.9" - }, - "pycharm": { - "stem_cell": { - "cell_type": "raw", - "source": [], - "metadata": { - "collapsed": false - } - } } }, "nbformat": 4, diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index 88df763..55e4b23 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,25 +532,27 @@ 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 = 3 # Interpolated grid points are interpolated using Lagrange + HIERARCHICAL_LAGRANGE_ROMBERG = 4 # Interpolated grid points are interpolated using hierarchical Lagrange + BSPLINE_ROMBERG = 5 -class SliceContainerVersion(Enum): + +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 - BSPLINE_ROMBERG = 5 - HIERARCHICAL_LAGRANGE_ROMBERG = 6 + 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: @@ -561,9 +569,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 @@ -578,26 +587,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 or \ - self.container_version == SliceContainerVersion.BSPLINE_ROMBERG or \ - self.container_version == SliceContainerVersion.HIERARCHICAL_LAGRANGE_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: @@ -2743,8 +2755,7 @@ 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: return LagrangeRombergGridSliceContainer(function) elif self.slice_container_version == SliceContainerVersion.SIMPSON_ROMBERG: return SimpsonRombergGridSliceContainer(function) diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index 91fc6a0..71ce0d3 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): @@ -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,9 +1226,10 @@ 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) diff --git a/sparseSpACE/TestInterpolationAdapter.py b/sparseSpACE/TestInterpolationAdapter.py new file mode 100644 index 0000000..7e98437 --- /dev/null +++ b/sparseSpACE/TestInterpolationAdapter.py @@ -0,0 +1 @@ +# TODO Implement \ No newline at end of file diff --git a/test/test_ExtrapolationInterpolatingGrid.py b/test/test_ExtrapolationInterpolatingGrid.py index 1f04d23..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,29 +140,50 @@ def test_full_grid_interpolation(self): expected_result = 1.5 self.assertAlmostEqual(expected_result, actual_result) - def test_langrange_interpolation(self): + 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(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, + 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_bspline_interpolation(self): + 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(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, + 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]) @@ -185,14 +193,15 @@ def test_bspline_interpolation(self): # self.assertAlmostEqual(expected_result, actual_result) - def test_hierarchical_lagrange_interpolation(self): + 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(slice_grouping=SliceGrouping.GROUPED_OPTIMIZED, - slice_version=SliceVersion.ROMBERG_DEFAULT, + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, container_version=SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG, - force_balanced_refinement_tree=False) + 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]) From 8f8c1aa77abaff44638e099ff910c59b933f0c35 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Wed, 17 Mar 2021 17:35:46 +0100 Subject: [PATCH 07/14] Add possibility to plot basis functions for GlobalBasisGrids --- ipynb/Diss/Figures_Diss.ipynb | 2 +- ipynb/Extrapolation/Extrapolation_Test.ipynb | 4 +- sparseSpACE/BasisFunctions.py | 10 +- sparseSpACE/Extrapolation.py | 18 +++- sparseSpACE/Grid.py | 36 ++++++- sparseSpACE/TestInterpolationAdapter.py | 103 ++++++++++++++++++- 6 files changed, 163 insertions(+), 10 deletions(-) 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 6c25afc..57005d3 100644 --- a/ipynb/Extrapolation/Extrapolation_Test.ipynb +++ b/ipynb/Extrapolation/Extrapolation_Test.ipynb @@ -313,9 +313,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "name": "pycharm-fcc6cfaa", "language": "python", - "name": "python3" + "display_name": "PyCharm (bachelorarbeit)" }, "language_info": { "codemirror_mode": { 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 55e4b23..41c64ca 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -2487,12 +2487,17 @@ def populate_interpolated_point_weights(self, factory: RombergWeights, normalized_max_level) for support_point in support_points: - interpolated_extrapolated_weight = weight * self.get_interpolation_weight(support_points, support_point, - interp_point) + interpolation_weight = self.get_interpolation_weight(support_points, support_point, interp_point) + interpolated_extrapolated_weight = weight * interpolation_weight weight_dictionary[support_point].append(interpolated_extrapolated_weight) - def get_interpolation_weight(self, support_points: Sequence[float], basis_point: float, evaluation_point: float)\ - -> float: + self.plot_interpolation_basis(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(self, support_points, evaluation_point): pass @@ -2556,6 +2561,7 @@ def __init__(self, function: Function = None): """ super(GlobalBasisGridRombergGridSliceContainerAdapter, self).__init__(function) self.max_interpolation_support_points = 7 + self.interpolation_grid = None @abstractmethod def construct_interpolation_grid(self, support_points): @@ -2574,6 +2580,7 @@ def construct_interpolation_grid(self, support_points): 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]) @@ -2588,6 +2595,9 @@ def get_interpolation_weight(self, support_points, basis_point, evaluation_point return weight + def plot_interpolation_basis(self, support_points, evaluation_point): + self.interpolation_grid.plot_basis(support_points, evaluation_point, self.function) + class BSplineRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter): diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index 71ce0d3..dc14710 100644 --- a/sparseSpACE/Grid.py +++ b/sparseSpACE/Grid.py @@ -353,7 +353,7 @@ def __init__(self, a: float, b: float, boundary: bool=True, p: int=3, modified_b self.boundary = boundary self.a = a self.b = b - assert len(a) == len(b) + assert len(a) == len(b) # TODO a, b are not floats self.dim = len(a) self.integrator = IntegratorHierarchicalBasisFunctions(self) self.grids = [LagrangeGrid1D(a=a[d], b=b[d], boundary=self.boundary, p=p, modified_basis=modified_basis) for d in range(self.dim)] @@ -1352,6 +1352,37 @@ def interpolate_grid(self, grid_points_for_dims: Sequence[Sequence[float]], comp #print(results) return results + def plot_basis(self, support_points=None, evaluation_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") + + # Plit 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") + + # Plot evaluation point + if evaluation_point is not None: + evaluation_point_value = basis(evaluation_point) + plt.plot(evaluation_point, evaluation_point_value, 'bo', markersize=6, color="red") + plt.plot([evaluation_point, evaluation_point], [0, evaluation_point_value], ':', + color="red") + + plt.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): @@ -1461,6 +1492,9 @@ def get_full_level_hierarchy(self, grid_1D, grid_levels_1D, a, b): return level_coordinate_array_complete +# TODO Default Global Lagrange grid with set grid?? + +# Hierarchical Lagrange Grid class GlobalLagrangeGrid(GlobalBasisGrid): def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, modified_basis: bool=False, p: int=3): self.boundary = boundary diff --git a/sparseSpACE/TestInterpolationAdapter.py b/sparseSpACE/TestInterpolationAdapter.py index 7e98437..40c4fc2 100644 --- a/sparseSpACE/TestInterpolationAdapter.py +++ b/sparseSpACE/TestInterpolationAdapter.py @@ -1 +1,102 @@ -# TODO Implement \ No newline at end of file +import sparseSpACE +from sparseSpACE.Extrapolation import ExtrapolationGrid, GridVersion, SliceContainerVersion, SliceGrouping, SliceVersion +from sparseSpACE.Function import Polynomial1d + +# TODO refactor in a test case + + +def run_extrapolation(grid, grid_levels, function, container_version: SliceContainerVersion): + print("with {}".format(container_version)) + romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, + container_version=container_version, + 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 + + # romberg_grid.plot_grid_with_containers() + print("- expected result: {}".format(expected_result)) + print("- actual result: {}".format(actual_result)) + print("- delta: {}".format(actual_result - expected_result)) + print() + + +# -------------------------------------------------------- +# Balanced adaptive grid +grid = [0, 0.5, 0.625, 0.75, 1] +grid_levels = [0, 1, 3, 2, 0] +function = Polynomial1d([1, 0, 0, 2]) +print("Balanced adaptive grid") +print(grid) +print(grid_levels) +print() + +# -- with Lagrange interpolation +# run_extrapolation(grid, grid_levels, function, SliceContainerVersion.LAGRANGE_ROMBERG) + +# -- with hierarchical lagrange grid +run_extrapolation(grid, grid_levels, function, SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG) + +# -- with BSpline interpolation +run_extrapolation(grid, grid_levels, function, SliceContainerVersion.BSPLINE_ROMBERG) + + +# -------------------------------------------------------- +# Unbalanced grid + +# 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) \ No newline at end of file From 6647702f3a2792dd6e01810c42bac910efacbb19 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Wed, 17 Mar 2021 19:10:14 +0100 Subject: [PATCH 08/14] Fix plots of basis functions --- sparseSpACE/Extrapolation.py | 44 ++++++++++++++++++++++--- sparseSpACE/Grid.py | 9 +++-- sparseSpACE/TestInterpolationAdapter.py | 6 ++-- 3 files changed, 47 insertions(+), 12 deletions(-) diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index 41c64ca..fd57130 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -2486,18 +2486,22 @@ 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: 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) - self.plot_interpolation_basis(support_points, interp_point) + # TODO add a flag for plotting this + # Plot basis functions + self.plot_interpolation_basis(support_points, support_point, interp_point) def get_interpolation_weight(self, support_points: Sequence[float], basis_point: float, evaluation_point: float) -> float: pass - def plot_interpolation_basis(self, support_points, evaluation_point): + def plot_interpolation_basis(self, support_points: Sequence[float], basis_point: float, evaluation_point: float): pass @@ -2542,6 +2546,37 @@ def get_langrange_basis(support_points: Sequence[float], basis_point: float, eva return evaluated_basis + # TODO use plots from GlobalBasisGrid + def plot_interpolation_basis(self, support_points, basis_point, evaluation_point): + # Iterate over all dimensions and plot each basis + plt.figure() + + a = support_points[0] + b = support_points[-1] + + if self.function is not None: + evaluation_points = np.linspace(a, b, len(support_points) * 100) + plt.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) + plt.plot(evaluation_points, [self.get_langrange_basis(support_points, support_point, point) + for point in evaluation_points]) + + # Print interpolation grid + plt.plot(support_points, np.zeros(len(support_points)), 'bo', markersize=6, color="dimgray") + + # Plot evaluation point + if evaluation_point is not None: + evaluation_point_value = self.get_langrange_basis(support_points, basis_point, evaluation_point) + plt.plot(evaluation_point, evaluation_point_value, 'bo', markersize=6, color="red") + plt.plot([evaluation_point, evaluation_point], [0, evaluation_point_value], ':', + color="red") + + plt.show() + class GlobalBasisGridRombergGridSliceContainerAdapter(InterpolationGridSliceContainer): @@ -2592,11 +2627,12 @@ def get_interpolation_weight(self, support_points, basis_point, evaluation_point 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(self, support_points, evaluation_point): - self.interpolation_grid.plot_basis(support_points, evaluation_point, self.function) + def plot_interpolation_basis(self, support_points, basis_point, evaluation_point): + self.interpolation_grid.plot_basis(support_points, basis_point, evaluation_point, self.function) class BSplineRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter): diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index dc14710..6e27281 100644 --- a/sparseSpACE/Grid.py +++ b/sparseSpACE/Grid.py @@ -1352,7 +1352,7 @@ def interpolate_grid(self, grid_points_for_dims: Sequence[Sequence[float]], comp #print(results) return results - def plot_basis(self, support_points=None, evaluation_point=None, function=None, filename=None): + def plot_basis(self, support_points=None, basis_point=None, evaluation_point=None, function=None, filename=None): if self.basis is None: print("No basis functions constructed") return @@ -1366,7 +1366,7 @@ def plot_basis(self, support_points=None, evaluation_point=None, function=None, plt.plot(evaluation_points, [function([point]) for point in evaluation_points], linestyle=':', color="black") - # Plit basis functions + # 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]) @@ -1376,6 +1376,7 @@ def plot_basis(self, support_points=None, evaluation_point=None, function=None, # Plot evaluation point if evaluation_point is not None: + basis = self.basis[d][support_points.index(basis_point)] evaluation_point_value = basis(evaluation_point) plt.plot(evaluation_point, evaluation_point_value, 'bo', markersize=6, color="red") plt.plot([evaluation_point, evaluation_point], [0, evaluation_point_value], ':', @@ -1492,9 +1493,7 @@ def get_full_level_hierarchy(self, grid_1D, grid_levels_1D, a, b): return level_coordinate_array_complete -# TODO Default Global Lagrange grid with set grid?? - -# Hierarchical Lagrange Grid +# TODO Rename into GlobalHierarchicalLagrangeGrid class GlobalLagrangeGrid(GlobalBasisGrid): def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, modified_basis: bool=False, p: int=3): self.boundary = boundary diff --git a/sparseSpACE/TestInterpolationAdapter.py b/sparseSpACE/TestInterpolationAdapter.py index 40c4fc2..cbd6303 100644 --- a/sparseSpACE/TestInterpolationAdapter.py +++ b/sparseSpACE/TestInterpolationAdapter.py @@ -34,13 +34,13 @@ def run_extrapolation(grid, grid_levels, function, container_version: SliceConta print() # -- with Lagrange interpolation -# run_extrapolation(grid, grid_levels, function, SliceContainerVersion.LAGRANGE_ROMBERG) +run_extrapolation(grid, grid_levels, function, SliceContainerVersion.LAGRANGE_ROMBERG) # -- with hierarchical lagrange grid -run_extrapolation(grid, grid_levels, function, SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG) +# run_extrapolation(grid, grid_levels, function, SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG) # -- with BSpline interpolation -run_extrapolation(grid, grid_levels, function, SliceContainerVersion.BSPLINE_ROMBERG) +# run_extrapolation(grid, grid_levels, function, SliceContainerVersion.BSPLINE_ROMBERG) # -------------------------------------------------------- From a858b8fa9826d5e7b440eec36347e8d06c214f3e Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Wed, 17 Mar 2021 20:29:36 +0100 Subject: [PATCH 09/14] Optimize plots of basis functions --- sparseSpACE/Extrapolation.py | 80 ++++++++++++++++++++++++------------ sparseSpACE/Grid.py | 60 +++++++++++++++++++++++---- 2 files changed, 105 insertions(+), 35 deletions(-) diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index fd57130..c3e1e64 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -2493,15 +2493,20 @@ def populate_interpolated_point_weights(self, factory: RombergWeights, 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(support_points, support_point, interp_point) + # 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(self, support_points: Sequence[float], basis_point: float, evaluation_point: float): + 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 @@ -2547,35 +2552,55 @@ def get_langrange_basis(support_points: Sequence[float], basis_point: float, eva return evaluated_basis # TODO use plots from GlobalBasisGrid - def plot_interpolation_basis(self, support_points, basis_point, evaluation_point): + # TODO add function to plot one basis point + # TODO add function to plot only basis functions + def plot_interpolation_basis_functions_with_evaluations(self, support_points, evaluation_point): # Iterate over all dimensions and plot each basis - plt.figure() + # 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] - if self.function is not None: - evaluation_points = np.linspace(a, b, len(support_points) * 100) - plt.plot(evaluation_points, [self.function([point]) for point in evaluation_points], - linestyle=':', color="black") + 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]) - # Plot basis functions - for support_point in support_points: - evaluation_points = np.linspace(a, b, len(support_points) * 100) - plt.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") - # Print interpolation grid - plt.plot(support_points, np.zeros(len(support_points)), 'bo', markersize=6, color="dimgray") + # 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") - # Plot evaluation point - if evaluation_point is not None: - evaluation_point_value = self.get_langrange_basis(support_points, basis_point, evaluation_point) - plt.plot(evaluation_point, evaluation_point_value, 'bo', markersize=6, color="red") - plt.plot([evaluation_point, evaluation_point], [0, evaluation_point_value], ':', - color="red") + axes[i].set_title("Evaluation using basis function of point {}".format(basis_point)) - plt.show() + 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): @@ -2631,8 +2656,11 @@ def get_interpolation_weight(self, support_points, basis_point, evaluation_point return weight - def plot_interpolation_basis(self, support_points, basis_point, evaluation_point): - self.interpolation_grid.plot_basis(support_points, basis_point, evaluation_point, self.function) + 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): diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index 6e27281..1e16b4e 100644 --- a/sparseSpACE/Grid.py +++ b/sparseSpACE/Grid.py @@ -1352,7 +1352,7 @@ def interpolate_grid(self, grid_points_for_dims: Sequence[Sequence[float]], comp #print(results) return results - def plot_basis(self, support_points=None, basis_point=None, evaluation_point=None, function=None, filename=None): + 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 @@ -1373,17 +1373,59 @@ def plot_basis(self, support_points=None, basis_point=None, evaluation_point=Non # Print interpolation grid plt.plot(support_points, np.zeros(len(support_points)), 'bo', markersize=6, color="dimgray") - - # Plot evaluation point - if evaluation_point is not None: - basis = self.basis[d][support_points.index(basis_point)] - evaluation_point_value = basis(evaluation_point) - plt.plot(evaluation_point, evaluation_point_value, 'bo', markersize=6, color="red") - plt.plot([evaluation_point, evaluation_point], [0, evaluation_point_value], ':', - color="red") + 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): From 591de88ee1c101da125d4b7f3d3cd2159d9f6a30 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Thu, 25 Mar 2021 13:42:21 +0100 Subject: [PATCH 10/14] Add GlobalLagrangeGrid without hierarchical construction --- sparseSpACE/Extrapolation.py | 211 +++++++++++++++++++---------------- sparseSpACE/Grid.py | 61 ++++++++++ 2 files changed, 178 insertions(+), 94 deletions(-) diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index c3e1e64..04af4d8 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -2495,8 +2495,8 @@ def populate_interpolated_point_weights(self, factory: RombergWeights, # 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) + # 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: @@ -2509,98 +2509,95 @@ def plot_interpolation_basis_functions_with_evaluations(self, support_points: Se evaluation_point: 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_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 - - # TODO use plots from GlobalBasisGrid - # TODO add function to plot one basis point - # TODO add function to plot only basis functions - 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() +# This class was replaced with LagrangeRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter) +# 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_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): @@ -2689,6 +2686,32 @@ def construct_interpolation_grid(self, support_points): 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 GlobalDefaultLagrangeGrid + + # Interpolation grid is computed beforehand + return GlobalDefaultLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) + + class HierarchicalLagrangeRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter): """ diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index 1e16b4e..c40f3e4 100644 --- a/sparseSpACE/Grid.py +++ b/sparseSpACE/Grid.py @@ -1535,6 +1535,67 @@ def get_full_level_hierarchy(self, grid_1D, grid_levels_1D, a, b): return level_coordinate_array_complete +# TODO Rename into GlobalLagrangeGrid +class GlobalDefaultLagrangeGrid(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) + 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) + + # TODO Rename into GlobalHierarchicalLagrangeGrid class GlobalLagrangeGrid(GlobalBasisGrid): def __init__(self, a: Sequence[float], b: Sequence[float], boundary: bool=True, modified_basis: bool=False, p: int=3): From 5ee818c57f849bb98082eb68ee16b939c6b23d55 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Thu, 25 Mar 2021 13:52:45 +0100 Subject: [PATCH 11/14] Update class names for global Lagrange grids --- sparseSpACE/Extrapolation.py | 8 ++++---- sparseSpACE/Grid.py | 8 +++----- sparseSpACE/spatiallyAdaptiveSingleDimension2.py | 8 ++++---- test/test_Hierarchization.py | 2 +- test/test_Integrator.py | 2 +- test/test_spatiallyAdaptiveSingleDimension2.py | 4 ++-- 6 files changed, 15 insertions(+), 17 deletions(-) diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index 04af4d8..144ca86 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -2706,10 +2706,10 @@ def __init__(self, function: Function = None): def construct_interpolation_grid(self, support_points): # Prevent circular dependencies - from sparseSpACE.Grid import GlobalDefaultLagrangeGrid + from sparseSpACE.Grid import GlobalLagrangeGrid # Interpolation grid is computed beforehand - return GlobalDefaultLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) + return GlobalLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) class HierarchicalLagrangeRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter): @@ -2732,10 +2732,10 @@ def __init__(self, function: Function = None): def construct_interpolation_grid(self, support_points): # Prevent circular dependencies - from sparseSpACE.Grid import GlobalLagrangeGrid + from sparseSpACE.Grid import GlobalHierarchicalLagrangeGrid # Interpolation grid is computed beforehand - return GlobalLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) + return GlobalHierarchicalLagrangeGrid([support_points[0]], [support_points[-1]], boundary=True) class TrapezoidalGridSlice(ExtrapolationGridSlice): diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index c40f3e4..1cdc18f 100644 --- a/sparseSpACE/Grid.py +++ b/sparseSpACE/Grid.py @@ -1535,8 +1535,7 @@ def get_full_level_hierarchy(self, grid_1D, grid_levels_1D, a, b): return level_coordinate_array_complete -# TODO Rename into GlobalLagrangeGrid -class GlobalDefaultLagrangeGrid(GlobalBasisGrid): +class GlobalLagrangeGrid(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) @@ -1596,8 +1595,7 @@ def _get_spline_integral(self, spline, d=None): return spline.get_integral(self.start, self.end, self.coords_gauss, self.weights_gauss) -# TODO Rename into GlobalHierarchicalLagrangeGrid -class GlobalLagrangeGrid(GlobalBasisGrid): +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) @@ -1703,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 5e9f7bf..0f19d71 100644 --- a/sparseSpACE/spatiallyAdaptiveSingleDimension2.py +++ b/sparseSpACE/spatiallyAdaptiveSingleDimension2.py @@ -1365,7 +1365,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): @@ -1374,7 +1374,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])) @@ -1384,14 +1384,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_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() From 42f749b3c6ba66f4a8cf78a93aa63f9cd376e7d2 Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Thu, 25 Mar 2021 19:18:45 +0100 Subject: [PATCH 12/14] Setup error plots for interpolation adapter comparison --- ipynb/Extrapolation/Extrapolation_Test.ipynb | 32 ++-- sparseSpACE/Extrapolation.py | 188 ++++++++++--------- 2 files changed, 117 insertions(+), 103 deletions(-) diff --git a/ipynb/Extrapolation/Extrapolation_Test.ipynb b/ipynb/Extrapolation/Extrapolation_Test.ipynb index 57005d3..c204b8b 100644 --- a/ipynb/Extrapolation/Extrapolation_Test.ipynb +++ b/ipynb/Extrapolation/Extrapolation_Test.ipynb @@ -5,8 +5,7 @@ "execution_count": null, "metadata": { "pycharm": { - "name": "#%%\n", - "is_executing": true + "name": "#%%\n" } }, "outputs": [], @@ -44,11 +43,11 @@ " force_balanced_refinement_tree=force_balanced_refinement_tree)\n", "\n", "# Settings\n", - "dim = 2\n", + "dim = 5\n", "a = np.zeros(dim)\n", "b = np.ones(dim)\n", "max_tol = 20 # set to max, for better comparison\n", - "max_evaluations = 10 ** 3\n", + "max_evaluations = 10 ** 5\n", "evaluation_points = evaluation_points(a, b, dim)\n", "functions = []\n", "sasd_version = 6\n", @@ -92,12 +91,12 @@ "# Manually override functions array:\n", "functions = [\n", " f_genz_corner,\n", - " f_genz_product,\n", - " f_genz_cont,\n", - " f_genz_gaussian,\n", - " f_exp_var,\n", - " f_genz_osz,\n", - " f_genz_disc\n", + " # f_genz_product,\n", + " # f_genz_cont,\n", + " # f_genz_gaussian,\n", + " # f_exp_var,\n", + " # f_genz_osz,\n", + " # f_genz_disc\n", "]\n", "\n", "# Iterate over all functions\n", @@ -190,6 +189,14 @@ " 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", @@ -279,9 +286,10 @@ " # (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", + "\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", diff --git a/sparseSpACE/Extrapolation.py b/sparseSpACE/Extrapolation.py index 144ca86..0e4b7f8 100644 --- a/sparseSpACE/Extrapolation.py +++ b/sparseSpACE/Extrapolation.py @@ -541,9 +541,12 @@ class SliceContainerVersion(Enum): SIMPSON_ROMBERG = 2 # Romberg extrapolation with simpson rule as base quadrature rule # Interpolation before extrapolation - LAGRANGE_ROMBERG = 3 # Interpolated grid points are interpolated using Lagrange - HIERARCHICAL_LAGRANGE_ROMBERG = 4 # Interpolated grid points are interpolated using hierarchical Lagrange - BSPLINE_ROMBERG = 5 + LAGRANGE_ROMBERG_OLD = 3 # Interpolated grid points are interpolated using Lagrange + + # 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): @@ -2509,95 +2512,96 @@ def plot_interpolation_basis_functions_with_evaluations(self, support_points: Se evaluation_point: float): pass + # This class was replaced with LagrangeRombergGridSliceContainer(GlobalBasisGridRombergGridSliceContainerAdapter) -# 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_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 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): @@ -2852,6 +2856,8 @@ 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_OLD: + return LagrangeRombergGridSliceContainerOld(function) elif self.slice_container_version == SliceContainerVersion.LAGRANGE_ROMBERG: return LagrangeRombergGridSliceContainer(function) elif self.slice_container_version == SliceContainerVersion.SIMPSON_ROMBERG: From c0ea48c8715a1676798e16df0f63a7eb2bbf7496 Mon Sep 17 00:00:00 2001 From: Pascal Date: Tue, 6 Apr 2021 15:04:40 +0200 Subject: [PATCH 13/14] Remove test file --- sparseSpACE/TestInterpolationAdapter.py | 102 ------------------------ 1 file changed, 102 deletions(-) delete mode 100644 sparseSpACE/TestInterpolationAdapter.py diff --git a/sparseSpACE/TestInterpolationAdapter.py b/sparseSpACE/TestInterpolationAdapter.py deleted file mode 100644 index cbd6303..0000000 --- a/sparseSpACE/TestInterpolationAdapter.py +++ /dev/null @@ -1,102 +0,0 @@ -import sparseSpACE -from sparseSpACE.Extrapolation import ExtrapolationGrid, GridVersion, SliceContainerVersion, SliceGrouping, SliceVersion -from sparseSpACE.Function import Polynomial1d - -# TODO refactor in a test case - - -def run_extrapolation(grid, grid_levels, function, container_version: SliceContainerVersion): - print("with {}".format(container_version)) - romberg_grid = ExtrapolationGrid(grid_version=GridVersion.INTERPOLATE_SUB_GRIDS, - container_version=container_version, - 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 - - # romberg_grid.plot_grid_with_containers() - print("- expected result: {}".format(expected_result)) - print("- actual result: {}".format(actual_result)) - print("- delta: {}".format(actual_result - expected_result)) - print() - - -# -------------------------------------------------------- -# Balanced adaptive grid -grid = [0, 0.5, 0.625, 0.75, 1] -grid_levels = [0, 1, 3, 2, 0] -function = Polynomial1d([1, 0, 0, 2]) -print("Balanced adaptive grid") -print(grid) -print(grid_levels) -print() - -# -- with Lagrange interpolation -run_extrapolation(grid, grid_levels, function, SliceContainerVersion.LAGRANGE_ROMBERG) - -# -- with hierarchical lagrange grid -# run_extrapolation(grid, grid_levels, function, SliceContainerVersion.HIERARCHICAL_LAGRANGE_ROMBERG) - -# -- with BSpline interpolation -# run_extrapolation(grid, grid_levels, function, SliceContainerVersion.BSPLINE_ROMBERG) - - -# -------------------------------------------------------- -# Unbalanced grid - -# 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) \ No newline at end of file From 5182586ca58bad210f21314bd871f114679fafbe Mon Sep 17 00:00:00 2001 From: Pascal Resch Date: Sun, 25 Apr 2021 18:06:01 +0200 Subject: [PATCH 14/14] Fix some Grid type hints --- sparseSpACE/Grid.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/sparseSpACE/Grid.py b/sparseSpACE/Grid.py index 84f3329..d5836a6 100644 --- a/sparseSpACE/Grid.py +++ b/sparseSpACE/Grid.py @@ -349,11 +349,11 @@ 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 - assert len(a) == len(b) # TODO a, b are not floats + assert len(a) == len(b) self.dim = len(a) self.integrator = IntegratorHierarchicalBasisFunctions(self) self.grids = [LagrangeGrid1D(a=a[d], b=b[d], boundary=self.boundary, p=p, modified_basis=modified_basis) for d in range(self.dim)] @@ -367,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) @@ -483,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 @@ -500,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