Line data    Source code

       1             : #pragma once
       2             : 
       3             : #include "vec.hpp"
       4             : 
       5             : namespace pa {
       6             : 
       7         127 : struct Box {
       8             :     vec upperbound;
       9             :     vec lowerbound;
      10             : };
      11             : 
      12             : /// Project a vector onto a box.
      13             : /// @f[ \Pi_C(v) @f]
      14             : template <class Vec>
      15        7717 : inline auto project(const Vec &v,  ///< [in] The vector to project
      16             :                     const Box &box ///< [in] The box to project onto
      17             : ) {
      18             :     using binary_real_f = real_t (*)(real_t, real_t);
      19       15434 :     return v.binaryExpr(box.lowerbound, binary_real_f(std::fmax))
      20        7717 :         .binaryExpr(box.upperbound, binary_real_f(std::fmin));
      21             : }
      22             : 
      23             : /// Get the difference between the given vector and its projection.
      24             : /// @f[ v - \Pi_C(v) @f]
      25             : /// @warning    Beware catastrophic cancellation!
      26             : template <class Vec>
      27             : inline auto
      28        1826 : projecting_difference(const Vec &v,  ///< [in] The vector to project
      29             :                       const Box &box ///< [in] The box to project onto
      30             : ) {
      31        1826 :     return v - project(v, box);
      32             : }
      33             : 
      34             : /// Get the distance squared between the given vector and its projection.
      35             : /// @f[ \left\| v - \Pi_C(v) \right\|_2^2 @f]
      36             : /// @warning    Beware catastrophic cancellation!
      37             : inline real_t dist_squared(crvec v,       ///< [in] The vector to project
      38             :                            const Box &box ///< [in] The box to project onto
      39             : ) {
      40             :     return projecting_difference(v, box).squaredNorm();
      41             : }
      42             : 
      43             : /// Get the distance squared between the given vector and its projection in the
      44             : /// Σ norm.
      45             : /// @f[ \left\| v - \Pi_C(v) \right\|_\Sigma^2
      46             : /// = \left(v - \Pi_C(v)\right)^\top \Sigma \left(v - \Pi_C(v)\right) @f]
      47             : /// @warning    Beware catastrophic cancellation!
      48        1888 : inline real_t dist_squared(crvec v,        ///< [in] The vector to project
      49             :                            const Box &box, ///< [in] The box to project onto
      50             :                            crvec Σ ///< [in] Diagonal matrix defining norm
      51             : ) {
      52             :     // TODO: Does this allocate?
      53             :     //       Does it have dangling references to temporaries?
      54        1888 :     auto d = v - project(v, box);
      55        1888 :     return d.dot(Σ.asDiagonal() * d);
      56        1888 : }
      57             : 
      58             : } // namespace pa