Add Moving Least Squares example #919

mrlag31 · 2023-08-01T16:11:46Z

This pull request adds a new example for ArborX, showing how it can be used to compute and approximate a function from a point cloud to another, using the Moving Least Squares method. This example is taken from DataTransferKit's MeshFree package and removes all usage of DTK. This is a draft PR as some extra work is still needed.

Some important differences from DTK's code:

Views are multidimensional instead of data vectors. I feel like it helps better understand the layout of each object and the operations done on them.
The inverse is computed using Gaussian elimination instead of SVD as a start.
Data is baked within the example.

aprokop

Really good start. Quick partial review to address some of the common issues.

examples/moving_least_squares/moving_least_squares.cpp

aprokop · 2023-08-01T17:48:02Z

examples/moving_least_squares/moving_least_squares.cpp

+  template <typename Double3D>
+  KOKKOS_INLINE_FUNCTION Kokkos::Array<double, size>
+  operator()(Double3D const &p) const


I think it should be templated on a Point. We can assert that the point is 3D (can we extend the basis to other dimensions)? Also, would want to switch away from using double.

This was taken directly from DTK's code of different polynomial bases. I will strip out the templates and switch double to float on the example. As I am using ArborX::Point, it is 3D.

I think it is possible to generate that code for any number of dimensions and polynomial degree, but it will most likely require a lot of templated code and classes.

examples/moving_least_squares/moving_least_squares.cpp

aprokop · 2023-08-01T17:56:28Z

examples/moving_least_squares/moving_least_squares.cpp

+  auto radii = Kokkos::View<double*, MemorySpace>(
+    "radii", target_points_num);
+  constexpr double epsilon = std::numeric_limits<double>::epsilon();


Probably want float here, it's up for a discussion.

It depends a lot on how the values should be computed. Should the coefficients be handled as doubles or floats?

It may make sense to do the calculations based on the floating type of the data. We have Point<DIM, FloatingPoint>, so should probably use the same type as the second template argument of the Point. Note that we need #918 with GeometricTraits::coordinate_type<Point>::type to figure out that type.

aprokop · 2023-08-01T17:57:24Z

examples/moving_least_squares/moving_least_squares.cpp

+  auto phi = Kokkos::View<double**, MemorySpace>(
+    "phi", target_points_num, num_neighbors);
+  Kokkos::parallel_for(
+    "phi_computation",


We also try to prefix the kernels to indicate the origin, like Example::phi_computation.

aprokop · 2023-08-01T17:59:14Z

examples/moving_least_squares/moving_least_squares.cpp

+  for (int i = 0; i < target_points_num; i++) {
+    std::cout << "====\n"
+              << target_values_host(i) << '\n'
+              << target_values_exact_host(i) << "\n====\n";
+  }


Probably want to compute and print an error.

aprokop · 2023-08-01T19:44:35Z

@mrlag31 You'll need to use clang-format 14 and fix the styling (e.g., using scripts/check_format_cpp.sh) to be able to test the PR.

mrlag31 · 2023-08-01T20:08:39Z

I've fixed some of the issues you found.

No templates (replaced with Point)
Now using only floats
Proper source point generation (but should also be done for target points)

examples/moving_least_squares/CMakeLists.txt

examples/moving_least_squares/moving_least_squares.cpp

aprokop · 2023-08-02T00:41:24Z

examples/moving_least_squares/moving_least_squares.cpp

+  constexpr std::size_t source_points_num = cube_side * cube_side * cube_side;
+  constexpr std::size_t target_points_num = 4;
+
+  auto source_points = Kokkos::View<ArborX::Point *, MemorySpace>(


It would be better to not use the default execution space.
Instead, it is preferable to create ExecutionSpace space;, and pass it to all view allocations, deep copies and kernels. This way, we can better keep track of the kernels, and also allow running MLS on multiple CUDA streams at the same time.

Also not sure why you would prefer to write it as auto blah = Type() instead of Type blah(), if you spell out the type anyway.

I'll add the execution space on those calls.

Regarding why I write auto blah = Type(), it's some habits that I have when declaring a variable. The other way is too similar to a function call in my opinion, and I think it helps readability. But I can change the declarations to use Type blah() instead.

examples/moving_least_squares/moving_least_squares.cpp

aprokop · 2023-08-02T00:45:41Z

examples/moving_least_squares/moving_least_squares.cpp

+  auto queries = Kokkos::View<ArborX::Nearest<ArborX::Point> *, MemorySpace>(
+      "MLS_EX::queries", target_points_num);
+  Kokkos::parallel_for(
+      "MLS_EX::make_queries",
+      Kokkos::RangePolicy<ExecutionSpace>(0, target_points_num),
+      KOKKOS_LAMBDA(int const i) {
+        queries(i) = ArborX::nearest(target_points(i), num_neighbors);
+      });


We should probably steer user to a better approach in the example. Instead of creating a view of queries, requiring an extra memory, we could provide a specialization of AccessTraits for a structure that contains target_points, and create queries on the fly.
For an example, see struct ArborX::AccessTraits<NearestToOrigin, ArborX::PredicatesTag> in examples/callback/example_callback.cpp.

Not a blocker, though.

That's how it was done in DTK, so I copied. But I agree that using AccessTraits would be a more "correct" way to do it.

That's how it was done in DTK

We've come a long way since DTK. DTK is pretty outdated in the way it uses ArborX. It is a good starting point, but we want and should move more towards a better and more efficient use of ArborX.

examples/moving_least_squares/moving_least_squares.cpp

aprokop · 2023-08-02T14:29:34Z

examples/moving_least_squares/moving_least_squares.cpp

+      KOKKOS_LAMBDA(int const i) {
+        float radius = 10. * epsilon;
+
+        for (int j = 0; j < num_neighbors; j++)


This seems wrong. What if a query returned less than num_neighbors? For example, if the number of source points is less? This loop would then use uninitialized data.

If the number of source points is less than the number of requested neighbors, I would most likely consider that an input error from the side of the user. Is there other reasons on why would a query not return the requested number of neighbors?

Yes, it is the only case where it would be not num_neighbors.

mrlag31 · 2023-08-07T18:39:01Z

I have switched from using Gaussian method to the SVD pseudo-inverse. I have also made adjustments from DTK's version by using the fact that the moment matrix is symmetric.

aprokop

I think DTK has a pretty decent interface for MLS, and we should follow it more closely here, with some differences:

Switch DeviceType to MemorySpace (do we need it?)
DeviceType is really outdated in the usage, and we prefer to use MemorySpace and ExecutionSpace separately
Introduce ExecutionSpace to the constructor and apply
They don't have to be the same (e.g., using different CUDA streams should be ok)
Duplicate communicator (similar to how it's done in the DistributedTree)
We want to make sure our communication does not overlap with users. Duplicating communicator is the safe way to do that.
(maybe) store Distributor
DTK stores indices, ranks, etc. Maybe we can simplify that.

Other things:

Anything regarding data setup should probably belong to the example .cpp (e.g., TargetPoints)
I think MLSComputation is really a collection of free functions, which should be in ArborX_DetailsMovingLeastSquares.hpp in the Details namespace (similar to DTK)
Likely don't want to hardcode user data to be Views (or, at least, point data), and switch to using AccessTraits to access source and target points
This would allow the interface to be a lot more flexible, and user to avoid converting their data to Kokkos Views.
We prefer function to follow camel case (and other styling things, see here)
num_neighbors should not be a parameter, it is determined from the polynomial basis size

aprokop

Another round of comments I had. Still haven't looked through all the code, though.

aprokop · 2023-08-16T14:52:43Z

examples/moving_least_squares/common.hpp

+template <typename T>
+using inner_value_t = std::decay_t<std::invoke_result_t<
+    decltype(ArborX::AccessTraits<T, ArborX::PrimitivesTag>::get), T const &,
+    int>>;


We have a utility to do that, AccessTraitsHelper. You can use it as AccessTraitsHelper<Access>::type. So you don't need this.

aprokop · 2023-08-16T14:57:04Z

examples/moving_least_squares/mls.hpp

+template <typename ValueType, typename PolynomialBasis, typename RBF,
+          typename MemorySpace>


Imho, this should be

Suggested change

template <typename ValueType, typename PolynomialBasis, typename RBF,

typename MemorySpace>

template <typename MemorySpace, typename PolynomialBasis = <some default>, typename RBF = <some default>, typename FloatingPoint = float>

I think we can find some reasonable default to arguments. We can debate on what arguments a user may more likely to change, which would affect the order.

Not sure what to call that FloatingPoint. That's the type that the calculations are done in, correct?

ValueType (and FloatingPoint in your example) is the type of the calculations and the stored matrices, but it might be different from the coordinate type and the evaluation type. FloatingCalculationType would be a better name?

FloatingCalculationType would be a better name?

May be a bit better than a generic ValueType, but it still does not quite express the intent. There must be a better name.

aprokop · 2023-08-16T14:58:14Z

examples/moving_least_squares/mls.hpp

+
+template <typename ValueType, typename PolynomialBasis, typename RBF,
+          typename MemorySpace>
+class MLS


Maybe rename to MovingLeastSquares, and then do an MLS alias?

aprokop · 2023-08-16T14:59:29Z

examples/moving_least_squares/mls.hpp

+  MLS(ExecutionSpace const &space, MPI_Comm comm, Points const &source_points,
+      Points const &target_points,


Do we want to assume that source and target points are in a similar container? Or should they be allowed to have different access?

I can't recall the order of execution space and mpi comm in the distributed tree. We should probably match that here, as we wouldn't want to have different order in different places.

Both the source_points and target_points are coming from user space and represent array of points. But I guess they can still have different types, so I will probably add more types to the templat.

aprokop · 2023-08-16T15:02:47Z

examples/moving_least_squares/mls.hpp

+        "Example::MLS::local_indices", _tgt_size * _num_neighbors);
+    Kokkos::View<int *, MemorySpace> local_ranks("Example::MLS::local_ranks",
+                                                 _tgt_size * _num_neighbors);


These views don't need initializations, as they will be set in the parallel loop below. So, please use WithoutInitializing.

aprokop · 2023-08-16T15:10:02Z

examples/moving_least_squares/mls.hpp

+                                                     source_points);
+
+    // Perform the query
+    Kokkos::View<Kokkos::pair<int, int> *, MemorySpace> index_ranks(


We have a type for that, I think, PairIndexRank from src/details/ArborX_PairIndexRank.hpp.

aprokop · 2023-08-16T15:13:33Z

examples/moving_least_squares/symmetric_pseudoinverse_svd.hpp

+// a sign matrix (only 1 or -1 in the diagonal, 0 elsewhere).
+// Thus A = U.E.S.U^T
+template <class ValueType, typename MemorySpace>
+class SymmPseudoInverseSVD


Maybe we should move this in a separate PR, and then add some tests for it?

I'm also a bit confused. Why is this a class, and why do we need to store data? It's called exactly once from MLS, and I don't think we use any of the stored _es, _inv, etc, anywhere in MLS. Why not make it a free function?

Maybe we should move this in a separate PR, and then add some tests for it?

Why not, as it is one of the more complex set of functions and could actually be used for other things.

I'm also a bit confused. Why is this a class, and why do we need to store data? It's called exactly once from MLS, and I don't think we use any of the stored _es, _inv, etc, anywhere in MLS. Why not make it a free function?

This is usually how I create complex algorithms that need temporary variables. The user cannot instantiate it so these variables are dropped as soon as the static method returns. However, I can transform it so it uses only free functions.

I have modified this part so it resembles what one could find in ArborX::Details. The inverse is now computed using only free functions and with less restrictions on the input.

aprokop · 2023-08-16T15:18:11Z

examples/moving_least_squares/mpi_comms.hpp

+#include <mpi.h>
+
+template <typename MemorySpace>
+class MPIComms


I haven't closely looked at this yet.

masterleinad · 2023-08-17T19:57:42Z

I'm concerned that this is getting too large. I would rather have a node-level (non-MPI) version and the MPI examples #732 separately and possibly an extension using concepts in common with #737.

mrlag31 · 2023-08-23T13:09:32Z

I have attempted to add support for hypergeometry (i.e. points that are not 3D) but DistributedTree does not support geometries that are not 3D. This seems to be dealt in PR #907 so I won't be able to test them with something else other than 3D points.

masterleinad · 2023-08-23T13:20:25Z

I have attempted to add support for hypergeometry (i.e. points that are not 3D) but DistributedTree does not support geometries that are not 3D. This seems to be dealt in PR #907 so I won't be able to test them with something else other than 3D points.

I would just not use ExperimentalHyperGeometry::Point then. Let's not make this pull request even more complex by suggesting to use functionality that is not available.

aprokop · 2023-08-23T13:24:28Z

I would just not use ExperimentalHyperGeometry::Point then. Let's not make this pull request even more complex by suggesting to use functionality that is not available.

I think I suggested it before MPI was added.

…cation

mrlag31 · 2023-08-24T12:24:49Z

I'm concerned that this is getting too large. I would rather have a node-level (non-MPI) version and the MPI examples #732 separately and possibly an extension using concepts in common with #737.

It is true that this starts to be quite a large example (although it can be used as an ArborX add-on). I think it would be interesting to discuss how can this PR be split.

aprokop · 2023-08-24T17:20:09Z

I think it would be interesting to discuss how can this PR be split.

One of the ways we could proceed (requires general consensus) is like this:

Create arborx/interpolation directory
Move SVD to arborx/interpolation/src/details/ArborX_DetailsSVD.hpp, and add testing to arborx/interpolation/test
Add radial basis functions and tests
Add non-distributed MLS and tests
Add non-distributed MLS example to arborx/interpolation/examples
Add distributed MLS and tests
Add distributed MLS example

This is just an idea. It really depends on whether we want to grow the set of the interpolation algorithms in the future and make them easily available to a user.

Rombur · 2023-08-24T17:26:34Z

It really depends on whether we want to grow the set of the interpolation algorithms in the future and make them easily available to a user.

Next year I have a project where I will need to couple two codes together (I don't remember the name of the codes), so I may use this.

aprokop · 2024-01-05T20:24:10Z

Can this be closed?

aprokop added the examples Anything to do with examples label Aug 1, 2023

aprokop reviewed Aug 1, 2023

View reviewed changes

aprokop changed the title ~~Moving Least Squares example (from DataTransferKit)~~ Add Moving Least Squares example Aug 2, 2023

aprokop reviewed Aug 2, 2023

View reviewed changes

mrlag31 force-pushed the dtk-moving-least-squares-example branch 2 times, most recently from 05272fb to efd1abe Compare August 14, 2023 18:54

aprokop reviewed Aug 14, 2023

View reviewed changes

aprokop reviewed Aug 16, 2023

View reviewed changes

mrlag31 added 15 commits August 23, 2023 16:24

New example program

12105e8

Source and target point creation

9289fc6

Completion of MLS (not tested)

8c97beb

Small fixes and values computation

ea3a2b4

Small fixup (memory and kernel names, removing templates, ...)

357bda3

source points generation (20x20x20 cube)

4c22a22

double to float

2586400

clang format

da4671b

Fixup (memory and kernel names, floats, error, exec space)

4f44e1c

Modifies predicates array into AccessTraits

15b150b

Correct declaration convention

d69f348

Typo fix, execution spaces in deep copies and range policies simplifi…

f6cd686

…cation

Switching from gaussian inverse to SVD

793a5f2

Specifying ExecutionSpace in RangePolicies

15ca7a1

Fixing wrong SVD calculation

62de5ed

mrlag31 added 21 commits August 23, 2023 16:24

Switching from std::size_t to int const and removing missing floats

4edbe19

Templation of the proper MLS computation

3a9afcf

CMake MPI check

68e199a

Templated Moving Least Squares

894e80f

Better RBF

1e6d4a5

Removed DeviceType, duplicated communicators

52fe5e1

Style fixes and num_neighbors as an optional arg

db39716

Moving ExecutionSpace templates

03b600c

Swicthing to AccessTraits for user inputs (attempt)

dec46be

Simplification of traits access

210243d

Assertions for public interfaces

b1267dd

Replacing code using ArborX's code and convention

9315a4c

Better symmetric pseudo inverse (free function and better template)

daf9822

Better MLS computations, more permissive templates

0c6af81

Proper usage of AccessTraits for MPI comms

a1bd291

Proper MLS public interface

a7650db

Extra RBFs

8fc8a75

NVCC and CUDA compliance (compilation errors)

f3556de

Automatic polynomial basis generation and better rbf interface

db133c8

Hypergeometry (only works in 3D)

bc388f8

Using point clouds creation and gathering of the maximum error

220407c

mrlag31 force-pushed the dtk-moving-least-squares-example branch from 0d2efb8 to 2065168 Compare August 24, 2023 12:17

Back and forth MLS

dcdedd3

mrlag31 force-pushed the dtk-moving-least-squares-example branch from 2065168 to dcdedd3 Compare August 24, 2023 12:20

mrlag31 mentioned this pull request Sep 7, 2023

Moving least squares interpolation #946

Merged

mrlag31 closed this Jan 11, 2024

		template <typename ValueType, typename PolynomialBasis, typename RBF,
		typename MemorySpace>

	template <typename ValueType, typename PolynomialBasis, typename RBF,
	typename MemorySpace>
	template <typename MemorySpace, typename PolynomialBasis = <some default>, typename RBF = <some default>, typename FloatingPoint = float>

		MLS(ExecutionSpace const &space, MPI_Comm comm, Points const &source_points,
		Points const &target_points,

Add Moving Least Squares example #919

Add Moving Least Squares example #919

Conversation

mrlag31 commented Aug 1, 2023 • edited Loading

aprokop 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

aprokop commented Aug 1, 2023

mrlag31 commented Aug 1, 2023

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 Aug 2, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mrlag31 commented Aug 7, 2023

aprokop left a comment • edited Loading

Choose a reason for hiding this comment

aprokop 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

mrlag31 Aug 17, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

masterleinad commented Aug 17, 2023

mrlag31 commented Aug 23, 2023 • edited Loading

masterleinad commented Aug 23, 2023

aprokop commented Aug 23, 2023

mrlag31 commented Aug 24, 2023

aprokop commented Aug 24, 2023

Rombur commented Aug 24, 2023

aprokop commented Jan 5, 2024

mrlag31 commented Aug 1, 2023 •

edited

Loading

mrlag31 Aug 2, 2023 •

edited

Loading

aprokop left a comment •

edited

Loading

mrlag31 Aug 17, 2023 •

edited

Loading

mrlag31 commented Aug 23, 2023 •

edited

Loading