feat: Implement a new orthogonal range search seed finder (#904)

As I said in #901, I have been playing around with seed finding a little bit lately. Last weekend, I mentioned an idea for a new (?) kind of seed finding algorithm based on range search datastructures, and this is the very, very first semi-working implementation of it, just before the weekend.

The idea behind this algorithm is relatively simple. In traditional seedfinding, we check a whole lot of candidate spacepoints to see whether they meet some condition. If you look at this differently, each spacepoint defines a volume in the z-r-φ space, which contains any spacepoints it can form a doublet with. What if we reversed this logic? What if we defined this volume first, and then just extract the spacepoints inside of that space? That way, we can vastly reduce the number of spacepoints we need to look at.

How do we do this quickly? With [_k_-d trees](https://en.wikipedia.org/wiki/K-d_tree). These data structures are cheap to build, and they give us very fast orthogonal range searches. In other words, we can very quickly look up which of our spacepoints lie within an axis-aligned orthognal n-dimensional hyperrectangle. In this case, which spacepoints lie within a z-r-φ box.

So, the core idea of this seedfinder is to define as many of our seedfinding constraints in orthogonal fashion. That way, we can make our candidate hyperrectangle smaller and smaller. The tighter the constraints we can place, the better. Then, we look up the relevant spacepoints, and we can avoid looking at any others. That also means this solution requires no binning whatsoever.

## Constraint conversion

Currently there are quite a few constraints in the code. Here is my status update on how well it is going to convert each of them. In some cases, we can define a weaker version of the constraints in orthogonal fashion. This is still very powerful, and it doesn't actually lose us any efficiency (because we can always check the tighter constraint in a non-orthogonal way later, not a problem)!

### Unary constraints

Currently, I am not aware of any unary constraints in the Acts seed finding code. That is to say, logic to determine whether a point is allowed to be a lower spacepoint. However, I have the following thoughts about introducing some:

* I believe the binning code does some kind of magic to determine whether a spacepoint can be a lower spacepoint. Since my solution doesn't use any binning, I don't have access to this just yet. However, if we can incorporate this logic it could be very powerful.
* Maximum single-point η: we currently have some checks in place to see if the pseudorapidity of particles is not too high. We could realistically use this maximum pseudorapidity, combined with the collision region range to constrain the bottom spacepoints.

### Binary constraints

These are the existing binary constraints on spacepoint duplets:

Constraint | Description | Orthogonalization
-|-|-
Minimum ∆r | Ensure that the second spacepoint is within a certain difference in radius | Full
Maximum cot θ | Ensure that the pseudorapidity of the duplet is not too high | Unsuccessful
z-origin range | Ensure that the duplet would have originated from the collision point | Weakened 
Maximum ∆φ<sup>1</sup> | Ensure that the duplet does not bend too much in the x-y plane | Full

<sup>1</sup> This check does not exist explicitly in the existing seed finder, but is implicit in the binning process.

### Ternary constraints

There are a lot of ternary constraints (to check whether a triplet is valid):

Constraint | Description | Orthogonalization
-|-|-
Scattering angle | ??? | Unsuccessful
Helix diameter | Ensure the helix diameter is within some range | In progress
Impact parameters | Ensure the impact parameters are close to the collision point | In progress
Monotonic z<sup>1</sup> | Ensure that z increases or decreases monotonically between points | Full

<sup>1</sup> This check does not exist in the existing seed finder, check #901.

There are also constraints defined in the experiment-specific cuts, and the seed filter, and in other places. If we could convert some of those to orthogonal constraints the implementation would become much more powerful. However, I don't really understand what is happening in those files just yet. Need more reading.

## Current performance

The current performance of this seedfinder is... Complicated. On my machine, it runs a 4000 π<sup>+</sup> event in about 5 seconds, three times slower than the existing seedfinder. Its efficiency is much higher though, and the fake rate is much lower. So that's something. However, that is in part because I am creating far more seed candidates, so take this with a big grain of salt.

## Potential gain

There are two ways that I can think of to use this kind of algorithm. The first is an inside-to-outside algorithm, where we pick a lower spacepoint first, check the space it defines for a middle spacepoint, and then check the space the two of them define for a third spacepoint. This algorithm has time complexity 𝒪(_n_<sup>3</sup>), and it has space complexity 𝒪(_n_). Due to the constants, I still believe this implementation can outperform the 𝒪(_n_<sup>2</sup>) existing algorithm, however.

The second way would be to construct a set of duplets using this logic, and then to fit those together like we do with traditional seedfinding. This has 𝒪(_n_<sup>2</sup>) time complexity like the existing code, but also space complexity 𝒪(_n_<sup>2</sup>).

## Things that are completed

* The implementation of the _k_-d tree seems to work very well, and it is quite fast.
* Basic seedfinding using this strategy is functional.

## Things that don't work yet

* My maximum ∆φ constraint does not cross the 2π boundary yet.
* I used the existing seedfinding algorithm as a stepping stone, which I have completely destroyed in the process. Obviously I do not intend on keeping it that way, and the existing algorithm will be restored to its full glory.
* Lots more.

## Things that can be improved

* Add more constraints, and tighten existing ones.
* Lots of things, pretty much everything. But I really want to go home for the weekend, so I will write this part next week.

Core/include/Acts/Seeding/SeedFinderOrthogonal.hpp

@@ -0,0 +1,227 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2022 CERN for the benefit of the Acts project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Seeding/InternalSeed.hpp"
+#include "Acts/Seeding/InternalSpacePoint.hpp"
+#include "Acts/Seeding/SeedFinderOrthogonalConfig.hpp"
+
+#include <array>
+#include <list>
+#include <map>
+#include <memory>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+namespace Acts {
+template <typename external_spacepoint_t>
+class SeedFinderOrthogonal {
+ public:
+  /**
+   * @brief Set the number of dimensions in which to embed points. This is just
+   * 3 for now (phi, r, and z), but we might want to increase or decrease this
+   * number in the future.
+   */
+  static constexpr std::size_t NDims = 3;
+
+  /**
+   * @brief The seed type used by this seeder internally.
+   */
+  using seed_t = Seed<external_spacepoint_t>;
+
+  /**
+   * @brief The spacepoint type used by this seeder internally.
+   */
+  using internal_sp_t = InternalSpacePoint<external_spacepoint_t>;
+
+  /**
+   * @brief The k-d tree type used by this seeder internally, which is
+   * three-dimensional, contains internal spacepoint pointers, uses the Acts
+   * scalar type for coordinates, stores its coordinates in std::arrays, and
+   * has leaf size 4.
+   */
+  using tree_t = KDTree<NDims, internal_sp_t *, ActsScalar, std::array, 4>;
+
+  /**
+   * @brief Construct a new orthogonal seed finder.
+   *
+   * @param config The configuration parameters for this seed finder.
+   */
+  SeedFinderOrthogonal(
+      Acts::SeedFinderOrthogonalConfig<external_spacepoint_t> config);
+
+  /**
+   * @brief Destroy the orthogonal seed finder object.
+   */
+  ~SeedFinderOrthogonal() = default;
+
+  /*
+   * Disallow various kinds of constructors, copies, and assignments.
+   */
+  SeedFinderOrthogonal() = delete;
+  SeedFinderOrthogonal(const SeedFinderOrthogonal<external_spacepoint_t> &) =
+      delete;
+  SeedFinderOrthogonal<external_spacepoint_t> &operator=(
+      const SeedFinderOrthogonal<external_spacepoint_t> &) = delete;
+
+  /**
+   * @brief Perform seed finding, appending seeds to a container.
+   *
+   * This method performs seed finding through an orthogonal range search
+   * strategy. This strategy differs from binning approaches because it selects
+   * seeds _constructively_ rather than _destructively_; instead of trying a
+   * large number of possible space point combinations and then rejecting many
+   * of them, this algorithm tries to only consider valid seed candidates to
+   * reduce combinatorics.
+   *
+   * In addition, this algorithm replaces the binning step used in other seed
+   * finding algorithms with the construction of a k-d tree, which allows us to
+   * efficiently search for space points within a given range.
+   *
+   * The core idea behind this algorithm is to create axis-aligned bounding
+   * boxes around the region of validity for a seed candidate (be it a bottom
+   * spacepoint for a given middle, a top for a given middle, a middle for a
+   * given bottom, or any other combination), and then searching the detector
+   * volume for points that lie inside that AABB.
+   *
+   * @tparam input_container_t The type of the input spacepoint container.
+   * @tparam output_container_t The type of the output seed container.
+   *
+   * @param spacePoints The input spacepoints from which to create seeds.
+   * @param out_it The output iterator to write seeds to.
+   */
+  template <typename input_container_t, typename output_container_t>
+  void createSeeds(const input_container_t &spacePoints,
+                   std::back_insert_iterator<output_container_t> out_it) const;
+
+  /**
+   * @brief Perform seed finding, returning a new container of seeds.
+   *
+   * This is a filterCandidates method for scenarios where a non-inserter API is
+   * more ergonomic. In general, the inserter-based method should be preferred
+   * as it is more flexible and usually more performant. For more information
+   * about the seeding algorithm, please see that function.
+   *
+   * @tparam input_container_t The type of the input spacepoint container.
+   *
+   * @param spacePoints The input spacepoints from which to create seeds.
+   *
+   * @return A vector of seeds.
+   */
+  template <typename input_container_t>
+  std::vector<seed_t> createSeeds(const input_container_t &spacePoints) const;
+
+ private:
+  /**
+   * @brief Enumeration of the different dimensions in which we can apply cuts.
+   */
+  enum Dim { DimPhi = 0, DimR = 1, DimZ = 2 };
+
+  /**
+   * @brief Return the AABB rearch range for a given spacepoint, searching
+   * upwards.
+   *
+   * This function calculates an axis-aligned bounding box around the volume of
+   * validity for the next spacepoint in a pair, given that the lower
+   * spacepoint is given. Thus, this method either takes a bottom spacepoint
+   * and returns a range for the middle spacepoint, or it takes a middle
+   * spacepoint and returns the range for the top spacepoint.
+   *
+   * @param low The lower spacepoint to find a partner for.
+   *
+   * @return An N-dimensional axis-aligned search range.
+   */
+  typename tree_t::range_t validTupleOrthoRangeLH(
+      const internal_sp_t &low) const;
+
+  /**
+   * @brief Return the AABB rearch range for a given spacepoint, searching
+   * downward.
+   *
+   * This function calculates an axis-aligned bounding box around the volume of
+   * validity for the next spacepoint in a pair, given that the upper
+   * spacepoint is given. Thus, this method either takes a middle spacepoint
+   * and returns a range for the bottom spacepoint, or it takes a top
+   * spacepoint and returns the range for the middle spacepoint.
+   *
+   * @param high The upper spacepoint to find a partner for.
+   *
+   * @return An N-dimensional axis-aligned search range.
+   */
+  typename tree_t::range_t validTupleOrthoRangeHL(
+      const internal_sp_t &high) const;
+
+  /**
+   * @brief Check whether two spacepoints form a valid tuple.
+   *
+   * This method checks whether the cuts that we have for pairs of space points
+   * hold.
+   *
+   * @warning This method checks ONLY those constraints that cannot be exactly
+   * represented as bounding boxes. Thus, this method should not be used on
+   * pairs of points that were not generated using a constrained spatial search
+   * strategy.
+   *
+   * @param low The lower spacepoint.
+   * @param high The upper spacepoint.
+   *
+   * @return True if the two points form a valid pair, false otherwise.
+   */
+  bool validTuple(const internal_sp_t &low, const internal_sp_t &high) const;
+
+  /**
+   * @brief Create a k-d tree from a set of spacepoints.
+   *
+   * @param spacePoints The spacepoints to create a tree from.
+   *
+   * @return A k-d tree containing the given spacepoints.
+   */
+  tree_t createTree(const std::vector<internal_sp_t *> &spacePoints) const;
+
+  /**
+   * @brief Filter potential candidate pairs, and output seeds into an
+   * iterator.
+   *
+   * @tparam output_it_t The type of the output iterator.
+   *
+   * @param middle The (singular) middle spacepoint.
+   * @param bottom The (vector of) candidate bottom spacepoints.
+   * @param top The (vector of) candidate top spacepoints.
+   * @param it The iterator to write the resulting seeds to.
+   */
+  template <typename output_it_t>
+  void filterCandidates(internal_sp_t &middle,
+                        std::vector<internal_sp_t *> &bottom,
+                        std::vector<internal_sp_t *> &top,
+                        output_it_t it) const;
+
+  /**
+   * @brief Search for seeds starting from a given middle space point.
+   *
+   * @tparam NDims Number of dimensions for our spatial embedding (probably 3).
+   * @tparam output_it_t Type of the output iterator.
+   *
+   * @param tree The k-d tree to use for searching.
+   * @param out_it The iteratorto write output seeds to.
+   * @param middle_p The middle spacepoint to find seeds for.
+   */
+  template <typename output_it_t>
+  void processFromMiddleSP(const tree_t &tree, output_it_t out_it,
+                           const typename tree_t::pair_t &middle_p) const;
+
+  /**
+   * @brief The configuration for the seeding algorithm.
+   */
+  Acts::SeedFinderOrthogonalConfig<external_spacepoint_t> m_config;
+};
+}  // namespace Acts
+
+#include "Acts/Seeding/SeedFinderOrthogonal.ipp"