Moving least squares interpolation #946

mrlag31 · 2023-09-07T13:29:47Z

This is a follow-up from #919. It properly adds the moving least square interpolation method as well as some tests.

dalg24

Some feedback on the example
(Haven't looked at the rest yet)

dalg24 · 2023-10-17T04:19:47Z

examples/CMakeLists.txt

 find_package(Boost COMPONENTS program_options)
 if(Boost_FOUND)
  add_subdirectory(viz)
  add_subdirectory(raytracing)
  add_subdirectory(brute_force)
-endif()
+endif()


dalg24 · 2023-10-17T13:49:30Z

examples/moving_least_squares/moving_least_squares.cpp

+template <typename T>
+KOKKOS_INLINE_FUNCTION double step(T const &p)


Why is the argument templated but not the return type?

Old code when I tested different point dimensions. Simply forgot to remove the template.

dalg24 · 2023-10-17T13:52:49Z

examples/moving_least_squares/moving_least_squares.cpp

+template <typename T>
+KOKKOS_INLINE_FUNCTION double step(T const &p)
+{
+  return Kokkos::signbit(p[0]) ? 0 : 1;


You could also return 1 - signbit(p[0]) but I think

auto const x = p[0]; return x >= 0 ? 1 : 0;

might be more readable.

dalg24 · 2023-10-17T13:56:12Z

examples/moving_least_squares/moving_least_squares.cpp

+  Kokkos::ScopeGuard guard(argc, argv);
+  ExecutionSpace space{};
+
+  static constexpr std::size_t num_points = 1000;


Did you consider making it the default value and optionally taking it from the command line arguments?
This would be handy when studying convergence as the number of points increases.

I made that change. You can input the number of points using the command line arguments with --points.

examples/moving_least_squares/moving_least_squares.cpp

dalg24 · 2023-10-17T14:05:10Z

examples/moving_least_squares/moving_least_squares.cpp

+
+  auto approx_values = mls.apply(space, source_values);
+
+  double max_error = 0.;


We usually suggest not to initialize the value you reduce into because it may give the wrong impression to users,
that the reduced value will be contributed to an initial value or whatnot, when Kokkos will actually ignore that value and overwrite it with the result of the reduction.

Suggested change

double max_error = 0.;

double max_error;

dalg24 · 2023-10-17T14:07:58Z

examples/moving_least_squares/moving_least_squares.cpp

+        loc_error = Kokkos::max(
+            loc_error, Kokkos::abs(target_values(i) - approx_values(i)));
+      },
+      Kokkos::Max<double>(max_error));


Why did you go for max absolute difference rather than L1 or L2 error?

I used the first norm that came to my mind. I have changed it to use an L2 norm.

dalg24 · 2023-10-17T14:08:55Z

examples/moving_least_squares/moving_least_squares.cpp

+      },
+      Kokkos::Max<double>(max_error));
+
+  std::cout << "Error: " << max_error << '\n';


Did you have a chance to study the convergence of the error as the number of points increases?

A bit using the new L2 norm. It seems the more points, the better the error.

masterleinad · 2023-10-17T18:03:28Z

examples/moving_least_squares/moving_least_squares.cpp

+  ArborX::Interpolation::MovingLeastSquares<MemorySpace, double> mls(
+      space, source_points, target_points);
+
+  auto approx_values = mls.apply(space, source_values);


What does mls.apply mean? I would prefer something more expressive like mls.interpolate.

'apply' was the function name in DTK and I kept it, I will switch it

src/interpolation/ArborX_Interp.hpp

masterleinad · 2023-10-17T19:11:17Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+// This is done to avoid a clash with another predicates access trait
+template <typename Points>
+struct MLSTargetPointsPredicateWrapper
+{
+  Points target_points;
+  int num_neighbors;
+};


Which would that be?
What do the possible types for Points look like? I assume these are Views or would the implementation also allow for something different?

The implementation allows any predicates access traits as long as "get" resolves to a point. However, I could not really test if Points is valid there as ArborX::Details::check_valid_access_traits and ArborX::GeometryTraits::check_valid_geometry_traits does not return a boolean.

masterleinad · 2023-10-17T19:14:20Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+class MovingLeastSquares
+{
+public:
+  // If num_neighbors is 0 or negative, it will instead be a default value


What is the default value/behavior? Of course, the behavior should be documented in the Wiki.
I would move that comment to line 116. Also, did you consider using std::optional?

The default value is the size of the polynomial basis. And using optional is a better alternative than using zero or a negative.

src/geometry/ArborX_GeometryTraits.hpp

masterleinad · 2023-10-18T17:33:12Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+    int const num_targets = tgt_acc::size(target_points);
+    _source_size = source_points.extent(0);
+    // There must be enough source points
+    ARBORX_ASSERT(num_neighbors <= _source_size);


I'm not sure a ArborX::SearchException makes sense here.

Suggested change

ARBORX_ASSERT(num_neighbors <= _source_size);

KOKKOS_ASSERT(num_neighbors <= _source_size);

As this code is part of ArborX, it makes more sense to use ArborX's assert than Kokkos one. It is also used elsewhere for non-search errors as well.

masterleinad · 2023-10-18T17:36:36Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+                     SourcePoints const &source_points,
+                     TargetPoints const &target_points)
+      : MovingLeastSquares(space, source_points, target_points,
+                           CRBF::Wendland<2>{}, PolynomialDegree<2>{})


We need to document default values. Why did you choose these?

They were the ones that yielded the best results in my firsts small examples. We could discuss on which default would actually be best.

masterleinad · 2023-10-18T17:37:59Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+        Kokkos::view_alloc(Kokkos::WithoutInitializing,
+                           "ArborX::MovingLeastSquares::source_view"),
+        num_targets, num_neighbors);
+    auto const values_indices = _values_indices;


KOKKOS_CLASS_LAMBDA?

Same as previously, it was to get around nvcc's lambda restrictions.

src/interpolation/details/ArborX_InterpDetailsPolynomialBasis.hpp

masterleinad · 2023-10-18T17:51:32Z

test/tstInterpDetailsMLSC.cpp

+  Kokkos::View<double *, MemorySpace> tgtv0("Testing::tgtv0", 3);
+  Kokkos::View<double **, MemorySpace> coeffs0("Testing::coeffs0", 0, 0);
+  Kokkos::parallel_for(
+      "for", Kokkos::RangePolicy<ExecutionSpace>(space, 0, 3),


Suggested change

"for", Kokkos::RangePolicy<ExecutionSpace>(space, 0, 3),

"Testing::mls_coefficients_case_1", Kokkos::RangePolicy<ExecutionSpace>(space, 0, 3),

etc.

masterleinad · 2023-10-23T15:32:08Z

examples/moving_least_squares/moving_least_squares.cpp

+      Kokkos::view_alloc(Kokkos::WithoutInitializing, "Example::target_values"),
+      num_points);
+  filledBoxRandom(0.5, source_points);
+  filledBoxRandom(0.5, target_points);


What kind of convergence behavior would you expect based on the number (and positions) of the source points?
IMHO, the interpretation of the error would be easier if at least the target_points were evenly distributed.

Currently, both source and target points are uniformely distributed. However, I may try to implement a way to use any regular mesh.

I would expect, depending on the number of source points, to converge as an inverse exponential. I do not really know regarding position of the source point.

masterleinad · 2023-10-23T18:09:56Z

examples/moving_least_squares/moving_least_squares.cpp

+                     << '\n';
+  }
+
+  dump_file_stream.close();


Not necessary. std::fstream will close in the destructor.

Suggested change

dump_file_stream.close();

masterleinad · 2023-10-23T18:13:50Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+    int const num_targets = tgt_acc::size(target_points);
+    _source_size = source_points.extent(0);
+    // There must be enough source points
+    ARBORX_ASSERT(0 < num_neighbors_val && num_neighbors_val <= _source_size);


Suggested change

ARBORX_ASSERT(0 < num_neighbors_val && num_neighbors_val <= _source_size);

KOKKOS_ASSERT(0 < num_neighbors_val && num_neighbors_val <= _source_size);

masterleinad · 2023-10-23T18:14:29Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+
+    // Source values must be a valuation on the points so must be as big as the
+    // original input
+    ARBORX_ASSERT(_source_size == source_values.extent_int(0));


Suggested change

ARBORX_ASSERT(_source_size == source_values.extent_int(0));

KOKKOS_ASSERT(_source_size == source_values.extent_int(0));

mrlag31 · 2023-11-21T16:06:03Z

nvcc 11.0.3 used to segfault on the previous CI (and probably still segfaults). Except that, it should be ready to be merged.

masterleinad

What does a typical profile look like?

masterleinad · 2023-11-21T19:41:29Z

examples/moving_least_squares/moving_least_squares.cpp

+// Randomly fills a 2 * half_edge box centered at 0
+void filledBoxRandom(double half_edge,
+                     Kokkos::View<Point *, MemorySpace> points)
+{
+  int n = points.extent(0);
+  auto points_host = Kokkos::create_mirror_view(points);
+
+  std::uniform_real_distribution<double> dist(-half_edge, half_edge);
+  std::random_device rd;
+  std::default_random_engine gen(rd());
+  auto random = [&dist, &gen]() { return dist(gen); };
+  for (int i = 0; i < n; ++i)
+    for (int d = 0; d < DIM; ++d)
+      points_host(i)[d] = random();
+
+  Kokkos::deep_copy(points, points_host);
+}


Do we really need multiple variants of creating points in an example?

It is meant to show the user that they can use any set of points, either random or a regular mesh. It can be discussed if it should be kept or not.

masterleinad · 2023-11-21T19:43:00Z

examples/moving_least_squares/moving_least_squares.cpp

+  Kokkos::deep_copy(points, points_host);
+}
+
+// Switches the set of points to the new base


Why is this necessary?

By default, everything is generated in a box. This lets skew and rotate the box however the user wants to. Like the previous comment, it can be removed, as well as the creation of the basis matrix.

masterleinad · 2023-11-21T19:45:12Z

examples/moving_least_squares/moving_least_squares.cpp

+  // Basis change
+  Kokkos::View<Point[DIM], MemorySpace> source_basis("Example::source_basis");
+  Kokkos::View<Point[DIM], MemorySpace> target_basis("Example::target_basis");
+  Kokkos::parallel_for(
+      "Example::fill_basis", Kokkos::RangePolicy<ExecutionSpace>(space, 0, 1),
+      KOKKOS_LAMBDA(int const) {
+        for (int i = 0; i < DIM; i++)
+        {
+          source_basis(i) = Point{};
+          target_basis(i) = Point{};
+          for (int j = 0; j < DIM; j++)
+          {
+            source_basis(i)[j] = double(i == j);
+            target_basis(i)[j] = double(i == j);
+          }
+        }
+      });


If this is just the identity why do we need to change the basis?

Mostly to show that the user can do it. I also used it to build a skeded/triangular mesh in my report. It can be removed with the basis change function.

masterleinad · 2023-11-21T19:55:55Z

examples/moving_least_squares/moving_least_squares.cpp

+      num_points);
+
+  // Sets the meshes
+  filledBoxEven(0.5, {side_len, side_len}, source_points);


Can't we just have this function return the source_points (and similar for the other Views/functions)?

That can be done. But I don't know which of "returning the view" or "giving an uninitialized view" is the better design choice.

masterleinad · 2023-11-21T20:04:59Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+        (!num_neighbors)
+            ? Details::polynomialBasisSize<dimension, PolynomialDegree::value>()
+            : *num_neighbors;


Suggested change

(!num_neighbors)

? Details::polynomialBasisSize<dimension, PolynomialDegree::value>()

: *num_neighbors;

num_neighbors

? *num_neighbors

: Details::polynomialBasisSize<dimension, PolynomialDegree::value>();

to improve readability some.

masterleinad · 2023-11-21T20:08:50Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+    auto const source_view = fillValuesIndicesAndGetSourceView(
+        space, indices, offsets, num_targets, num_neighbors_val, source_points);
+
+    // Compute the Moving Least Squares


Suggested change

// Compute the Moving Least Squares

// Compute the Moving Least Squares coefficients

masterleinad · 2023-11-21T20:15:13Z

src/interpolation/ArborX_InterpMovingLeastSquares.hpp

+    Details::movingLeastSquaresCoefficients<CRBF, PolynomialDegree>(
+        space, source_view, target_points, _coeffs);


Is there a good reason not to return the coefficients?

Like previously, I do not know if it is better to return the coefficients or to give an uninitialized array to it. Querying a tree requires to give out empty views for example.

mrlag31 · 2023-11-22T12:57:06Z

What does a typical profile look like?

What do you mean by "profile"?

masterleinad · 2023-11-22T14:29:36Z

What do you mean by "profile"?

What does the output for space-time-stack output look like for sufficiently big problems?

mrlag31 · 2023-11-22T14:51:50Z

Using Cuda and my laptop's Nvidia A1000, here are the results for both 4096 points and 1000000 points (source and target).

sts4k.txt
sts1m.txt

In short, to compute the coefficients, it lasts 3ms for 4096 points and 200ms for 1000000 points. And, as I had to do the same measures in my report, the duration scales linearly with the number of points.

masterleinad · 2023-11-27T20:22:16Z

In short, to compute the coefficients, it lasts 3ms for 4096 points and 200ms for 1000000 points. And, as I had to do the same measures in my report, the duration scales linearly with the number of points.

So ArborX::SymmetricPseudoInverseSVD::computations or ArborX::BVH::query::nearest would be where to look first if we wanted to improve performance.

Rombur · 2023-11-27T20:23:37Z

nvcc 11.0.3 used to segfault on the previous CI (and probably still segfaults). Except that, it should be ready to be merged.

This PR doesn't compile with CUDA 11.0.3 https://cloud.cees.ornl.gov/jenkins-ci/job/arborx/job/PR-946/4/pipeline-console/?selected-node=44

Rombur · 2023-11-27T20:26:50Z

Retest this please.

dalg24 · 2023-11-27T20:37:13Z

retest this please

masterleinad · 2023-11-27T21:56:24Z

I can compile and run the example successfully on Summit using

$ module list 

Currently Loaded Modules:
  1) xl/16.1.1-10                     3) lsf-tools/2.0   5) darshan-runtime/3.4.0-lite   7) DefApps                   9) nsight-systems/2021.3.1.54  11) cmake/3.23.2
  2) spectrum-mpi/10.4.0.3-20210112   4) hsi/5.0.2.p5    6) xalt/1.2.1                   8) nsight-compute/2021.2.1  10) cuda/11.0.3                 12) boost/1.77.0

masterleinad · 2023-11-27T21:57:18Z

@mrlag31 Can you post the issues you are seeing here if you can find them?

mrlag31 · 2023-11-28T13:14:56Z

For the moment, the crash appeared after updating my code to the new callback interface. I will investigate it.

examples/moving_least_squares/moving_least_squares.cpp

masterleinad

I'm fine with cleaning up in a follow-up.

aprokop · 2023-12-21T16:44:58Z

Couple builds did not run (HIP, SYCL, Cuda-Clang), but that's fine.

mrlag31 added a commit to mrlag31/ArborX that referenced this pull request Sep 11, 2023

SVD from arborx#946

a9b9eea

mrlag31 mentioned this pull request Sep 11, 2023

Pseudo-inverse of symmetric matrices using SVD / Utility for moving least squares #950

Merged

mrlag31 mentioned this pull request Sep 27, 2023

Compact radial basis functions and generic polynomial basis / Utility for moving least squares #954

Merged

mrlag31 force-pushed the moving-least-squares branch 4 times, most recently from f637b0f to 7c3f0e1 Compare October 5, 2023 13:33

mrlag31 force-pushed the moving-least-squares branch from 7c3f0e1 to 354efe9 Compare October 12, 2023 16:52

mrlag31 marked this pull request as ready for review October 13, 2023 12:49

dalg24 reviewed Oct 17, 2023

View reviewed changes

aprokop mentioned this pull request Oct 17, 2023

Do not initialize the value we reduce to #960

Merged

mrlag31 force-pushed the moving-least-squares branch from 0d9c491 to fdc710a Compare October 18, 2023 13:15

masterleinad reviewed Oct 18, 2023

View reviewed changes

masterleinad reviewed Oct 23, 2023

View reviewed changes

mrlag31 force-pushed the moving-least-squares branch from 75b03ba to 76cf2fe Compare November 9, 2023 15:28

mrlag31 force-pushed the moving-least-squares branch from 76cf2fe to 8b2e464 Compare November 21, 2023 14:35

masterleinad reviewed Nov 21, 2023

View reviewed changes

aprokop mentioned this pull request Nov 27, 2023

Consider increasing the minimum nvcc version from 11.0.3 to 11.4+ #974

Closed

mrlag31 force-pushed the moving-least-squares branch from b585725 to f6a7b1e Compare November 28, 2023 15:12

mrlag31 added 18 commits December 20, 2023 13:05

KOKKOS_ASSERT and auto file closing

01271f6

Regular meshes and basis changes

4e8c483

Adding neighbors to example and changing defaults

03afffc

New callback api fix

c0bd551

Readability and moving coefficients creation space

07ba1e0

Naming typo and better cl args

a743130

Update some tests (avoid repeated evaluation of 0s)

d2653ba

nvcc 11.0.3 / empty example

b07d0f0

Example rework

e2854ab

Interface update

afc0072

std::optional nvcc fix

a2e5a74

Using Kokkos scoped regions

d13afe5

Reworked example

ef2bebb

MLS Interface rework

da7c6ee

Renaming test

b945686

Polynomial degree alias

728458a

Usage of new BVH API

b1aca0c

Moving spaces inside main

ff7e6c7

mrlag31 force-pushed the moving-least-squares branch from 88b439f to ff7e6c7 Compare December 20, 2023 18:50

aprokop reviewed Dec 20, 2023

View reviewed changes

examples/moving_least_squares/moving_least_squares.cpp Outdated Show resolved Hide resolved

Removing hard-coded values in example

46ae5fd

dalg24 reviewed Dec 21, 2023

View reviewed changes

examples/moving_least_squares/moving_least_squares.cpp Outdated Show resolved Hide resolved

Comment update and better naming

37161b9

aprokop approved these changes Dec 21, 2023

View reviewed changes

Get source size through AccessTraits interface

6529ca0

masterleinad reviewed Dec 21, 2023

View reviewed changes

aprokop merged commit a522958 into arborx:master Dec 21, 2023
1 of 2 checks passed

mrlag31 mentioned this pull request Dec 22, 2023

Moving least squares (extra fixes) #992

Merged

aprokop added the enhancement New feature or request label Dec 23, 2023

		template <typename T>
		KOKKOS_INLINE_FUNCTION double step(T const &p)


		auto approx_values = mls.apply(space, source_values);

		double max_error = 0.;

	ARBORX_ASSERT(num_neighbors <= _source_size);
	KOKKOS_ASSERT(num_neighbors <= _source_size);

	"for", Kokkos::RangePolicy<ExecutionSpace>(space, 0, 3),
	"Testing::mls_coefficients_case_1", Kokkos::RangePolicy<ExecutionSpace>(space, 0, 3),

	ARBORX_ASSERT(0 < num_neighbors_val && num_neighbors_val <= _source_size);
	KOKKOS_ASSERT(0 < num_neighbors_val && num_neighbors_val <= _source_size);

	ARBORX_ASSERT(_source_size == source_values.extent_int(0));
	KOKKOS_ASSERT(_source_size == source_values.extent_int(0));

	// Compute the Moving Least Squares
	// Compute the Moving Least Squares coefficients

		Details::movingLeastSquaresCoefficients<CRBF, PolynomialDegree>(
		space, source_view, target_points, _coeffs);

Moving least squares interpolation #946

Moving least squares interpolation #946

Conversation

mrlag31 commented Sep 7, 2023 • edited Loading

dalg24 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 Oct 18, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 Oct 24, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 commented Nov 21, 2023

masterleinad left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 commented Nov 22, 2023

masterleinad commented Nov 22, 2023

mrlag31 commented Nov 22, 2023 • edited Loading

masterleinad commented Nov 27, 2023

Rombur commented Nov 27, 2023

Rombur commented Nov 27, 2023

dalg24 commented Nov 27, 2023

masterleinad commented Nov 27, 2023

masterleinad commented Nov 27, 2023 • edited Loading

mrlag31 commented Nov 28, 2023

masterleinad left a comment

Choose a reason for hiding this comment

aprokop commented Dec 21, 2023

mrlag31 commented Sep 7, 2023 •

edited

Loading

mrlag31 Oct 18, 2023 •

edited

Loading

mrlag31 Oct 24, 2023 •

edited

Loading

mrlag31 commented Nov 22, 2023 •

edited

Loading

masterleinad commented Nov 27, 2023 •

edited

Loading