polyval operator

Author: cliffburdick

docs_input/api/polynomials/polyval.rst

@@ -0,0 +1,18 @@
+.. _polyval_func:
+
+polyval
+=======
+
+Evaluate a polynomial given an input sequence and coefficients
+
+.. doxygenfunction:: polyval(const Op &op, const Coeffs &coeffs)
+
+Examples
+~~~~~~~~
+
+.. literalinclude:: ../../../test/00_operators/OperatorTests.cu
+   :language: cpp
+   :start-after: example-begin polyval-test-1
+   :end-before: example-end polyval-test-1
+   :dedent:
+

include/matx/operators/operators.h

@@ -78,6 +78,7 @@
 #include "matx/operators/percentile.h"
 #include "matx/operators/permute.h"
 #include "matx/operators/planar.h"
+#include "matx/operators/polyval.h"
 #include "matx/operators/pwelch.h"
 #include "matx/operators/qr.h"
 #include "matx/operators/r2c.h"

include/matx/operators/polyval.h

@@ -0,0 +1,101 @@
+////////////////////////////////////////////////////////////////////////////////
+// BSD 3-Clause License
+//
+// Copyright (c) 2021, NVIDIA Corporation
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// 1. Redistributions of source code must retain the above copyright notice, this
+//    list of conditions and the following disclaimer.
+//
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+//
+// 3. Neither the name of the copyright holder nor the names of its
+//    contributors may be used to endorse or promote products derived from
+//    this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+/////////////////////////////////////////////////////////////////////////////////
+
+#pragma once
+
+
+#include "matx/core/type_utils.h"
+#include "matx/operators/base_operator.h"
+
+namespace matx
+{
+
+  /**
+   * Returns the polynomial evaluated at each point
+   */
+  namespace detail {
+    template <typename Op, typename Coeffs>
+    class PolyvalOp : public BaseOp<PolyvalOp<Op, Coeffs>>
+    {
+      private:
+        Op op_;
+        Coeffs coeffs_;
+
+      public:
+        using matxop = bool;
+        using scalar_type = typename Op::scalar_type;
+
+        __MATX_INLINE__ std::string str() const { return "polyval()"; }
+        __MATX_INLINE__ PolyvalOp(const Op &op, const Coeffs &coeffs) : op_(op), coeffs_(coeffs) {
+          MATX_STATIC_ASSERT_STR(coeffs.Rank() == 1, matxInvalidDim, "Coefficient must be rank 1");
+          MATX_STATIC_ASSERT_STR(op.Rank() == 1, matxInvalidDim, "Input operator must be rank 1");
+        };
+
+        __MATX_INLINE__ __MATX_DEVICE__ __MATX_HOST__ auto operator()(index_t idx) const
+        {
+          // Horner's method for computing polynomial
+          scalar_type ttl{coeffs_(0)};
+          for(int i = 1; i < coeffs_.Size(0); i++) {
+              ttl = ttl * op_(idx) + coeffs_(i);
+          }
+
+          return ttl;
+        }
+
+        static __MATX_INLINE__ constexpr __MATX_HOST__ __MATX_DEVICE__ int32_t Rank()
+        {
+          return 1;
+        }
+        constexpr __MATX_INLINE__ __MATX_HOST__ __MATX_DEVICE__ index_t Size([[maybe_unused]] int dim) const
+        {
+          return op_.Size(0);
+        }
+    };
+  }
+
+
+  /**
+   * @brief Evaluate a polynomial
+   * 
+   * Currently only allows 1D input and coefficients
+   * 
+   * @tparam Op Type of input values to evaluate
+   * @tparam Coeffs Type of coefficients
+   * @param op Input values to evaluate
+   * @param coeffs Coefficient values
+   * @return polyval operator 
+   */
+  template <typename Op, typename Coeffs>
+  __MATX_INLINE__ auto polyval(const Op &op, const Coeffs &coeffs) {
+    return detail::PolyvalOp(op, coeffs);
+  }
+} // end namespace matx

test/00_operators/OperatorTests.cu

@@ -3274,6 +3274,33 @@ TYPED_TEST(OperatorTestsFloatNonComplexAllExecs, Sphere2Cart)
   MATX_EXIT_HANDLER();
 }
 
+TYPED_TEST(OperatorTestsFloatNonComplexNonHalfAllExecs, PolyVal)
+{
+  MATX_ENTER_HANDLER();
+  using TestType = std::tuple_element_t<0, TypeParam>;
+  using ExecType = std::tuple_element_t<1, TypeParam>;
+
+  auto pb = std::make_unique<detail::MatXPybind>();
+  pb->InitAndRunTVGenerator<TestType>("00_operators", "polyval_operator", "run", {4, 100});
+
+  ExecType exec{};
+  auto x = make_tensor<TestType>({100});
+  auto c = make_tensor<TestType>({4});
+  auto out = make_tensor<TestType>({100});
+
+  pb->NumpyToTensorView(x, "x");
+  pb->NumpyToTensorView(c, "c");
+
+  // example-begin polyval-test-1
+  (out = polyval(x, c)).run();
+  // example-end polyval-test-1
+  cudaStreamSynchronize(0);
+  MATX_TEST_ASSERT_COMPARE(pb, out, "out", 0.01);
+
+  MATX_EXIT_HANDLER();
+}
+
+
 TYPED_TEST(OperatorTestsNumericAllExecs, Upsample)
 {
   MATX_ENTER_HANDLER();

test/test_vectors/generators/00_operators.py

@@ -4,6 +4,21 @@
 import scipy.signal as ss
 from typing import Dict, List
 
+class polyval_operator:
+    def __init__(self, dtype: str, size: List[int]):
+        self.size = size
+        pass
+
+    def run(self) -> Dict[str, np.array]:
+        c = np.random.rand(self.size[0])
+        x = np.random.rand(self.size[1])
+
+        return {
+            'c': c,
+            'x': x,
+            'out': np.polyval(c, x),
+        }
+
 
 class kron_operator:
     def __init__(self, dtype: str, size: List[int]):