Add method to find analytically derivatives of SquarePolynomialVector

src/Classic/Fields/Interpolation/SquarePolynomialVector.cpp

@@ -114,7 +114,8 @@ void  SquarePolynomialVector::F(const double*   point,    double* value)
     for(unsigned int i=0; i<ValueDimension(); i++) value[i]  = valueV(i+1);
 }
 
-void  SquarePolynomialVector::F(const MVector<double>& point,   MVector<double>& value) const
+void  SquarePolynomialVector::F(const MVector<double>& point,
+                                MVector<double>& value) const
 {
     MVector<double> polyVector(_polyCoeffs.num_col(), 1);
     MakePolyVector(point, polyVector);
@@ -123,7 +124,9 @@ void  SquarePolynomialVector::F(const MVector<double>& point,   MVector<double>&
 }
 
 
-MVector<double>& SquarePolynomialVector::MakePolyVector(const MVector<double>& point, MVector<double>& polyVector) const
+MVector<double>& SquarePolynomialVector::MakePolyVector(
+                                const MVector<double>& point,
+                                MVector<double>& polyVector) const
 {
     for(unsigned int i=0; i<_polyCoeffs.num_col(); i++) {
         polyVector(i+1) = 1.;
@@ -161,7 +164,7 @@ void SquarePolynomialVector::IndexByPowerRecursive(std::vector<int> check, size_
 std::vector<int> SquarePolynomialVector::IndexByPower(int index, int point_dim) {
     if (point_dim < 1)
         throw(GeneralClassicException(
-            "PPSolveFactory::GetNearbyPointsSquares",
+            "SquarePolynomialVector::IndexByPower",
             "Point dimension must be > 0"
         ));
     while (int(_polyKeyByPower.size()) < point_dim)
@@ -253,5 +256,50 @@ unsigned int SquarePolynomialVector::PolynomialOrder() const {
     return maxPower;
 }
 
+SquarePolynomialVector SquarePolynomialVector::Deriv(const int* derivPower) const {
+    if (derivPower == NULL) {
+        throw(GeneralClassicException(
+            "SquarePolynomialVector::Deriv",
+            "Derivative points to NULL"
+        ));
+    }
+    MMatrix<double> newPolyCoeffs(_polyCoeffs.num_row(),  _polyCoeffs.num_col(), 0.);
+    std::vector< std::vector<int> > powerKey = _polyKeyByPower[_pointDim-1];
+    for (size_t j = 0; j < _polyCoeffs.num_col(); ++j) {
+        std::vector<int> thisPower = powerKey[j];
+        std::vector<int> newPower(_pointDim);
+        // f(x) = Product(x_i^t_i)
+        // d^n f(x)  / Product(d x_j^d_j) =
+        //       t_i >= d_j ... Product(x_i^(t_i - d_j) * t_i!/(t_i-d_j)!
+        //       t_i < d_j  ... 0
+        double newCoeff = 1.;
+        // calculate the coefficient
+        for (size_t k = 0; k < thisPower.size(); ++k) {
+            newPower[k] = thisPower[k] - derivPower[k];
+            if (newPower[k] < 0) {
+                newCoeff = 0.;
+                break;
+            }
+            newCoeff *= gsl_sf_fact(thisPower[k])/gsl_sf_fact(newPower[k]);
+        }
+        // put it in the matrix of coefficients
+        for (size_t k = 0; k < powerKey.size(); ++k) {
+            bool match = true;
+            for (size_t l = 0; l < powerKey[k].size(); ++l) {
+                if (powerKey[k][l] != newPower[l]) {
+                    match = false;
+                    break;
+                }
+            }
+            if (match) {
+                for (size_t i = 0; i < _polyCoeffs.num_row(); ++i) {
+                    newPolyCoeffs(i+1, k+1) = newCoeff*_polyCoeffs(i+1, j+1);
+                }
+            }
+        }
+    }
+    SquarePolynomialVector vec(_pointDim, newPolyCoeffs);
+    return vec;
+}
 
 }
\ No newline at end of file

src/Classic/Fields/Interpolation/SquarePolynomialVector.h

@@ -58,7 +58,7 @@ namespace interpolation {
  *  \f$ w_i = x_{i_1}x_{i_2} \ldots x_{i_n} \f$ \n
  *  \n
  *
- *  Nb: it is a SquarePolynomialVector because coefficients include all
+ *  Nb: it is a *Square*PolynomialVector because coefficients include all
  *  polynomial coefficients with \f$ i_ j <= n \f$; coefficients sit within an
  *  n-dimensional square. The distinction should be made with a PolynomialVector
  *  where coefficients include all polynomial coefficients with 
@@ -137,6 +137,17 @@ class SquarePolynomialVector {
      */
     void  F(const MVector<double>& point, MVector<double>& value) const;
 
+    /** Generate polynomial corresponding to the partial derivative of y
+     * 
+     *  @returns \f$ \frac{\partial^{j_1+j_2+\ldots} \vec{y}}
+     *                    {\partial x^j_1 \partial x^j_2 \ldots} \f$
+     *  at some set point \f$ \vec{x} \f$.
+     *
+     *  derivPower is an array with indices of the derivative in the 
+     *      "IndexByPower" style, of length PointDimension()
+     */
+    SquarePolynomialVector Deriv(const int* derivPower) const;
+
     /** Length of the input point (x) vector. */
     unsigned int PointDimension()         const {return _pointDim;}
 

tests/classic_src/Fields/Interpolation/SquarePolynomialVectorTest.cpp

@@ -108,4 +108,27 @@ TEST(SquarePolynomialVectorTest, TestF) {
     ref.F(&point3[0], &value[0]);
     EXPECT_EQ(value[0], refValue(1)); // sum (0, ..., 8)
     EXPECT_EQ(value[1], refValue(2)); // sum (9, ..., 17)
-}
+}
\ No newline at end of file
+
+TEST(SquarePolynomialVectorTest, TestDeriv) {
+    std::vector<double> data(32);
+    for (size_t i = 0; i < data.size(); ++i)
+        data[i] = i;
+    // Up to x^3
+    MMatrix<double> inputCoeffs(2, 16, &data[0]);
+    SquarePolynomialVector input(2, inputCoeffs);
+    std::vector<int> derivPower({1, 2});
+    SquarePolynomialVector test = input.Deriv(&derivPower[0]);
+    std::vector<double> refCoeffs({
+        28, 64, 156, 336, 132, 540, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
+        30, 68, 162, 348, 138, 558, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+    });
+    MMatrix<double> testCoeffs = test.GetCoefficientsAsMatrix();
+    size_t index = 0;
+    for (size_t i = 0; i < testCoeffs.num_row(); ++i) {
+        for (size_t j = 0; j < testCoeffs.num_col(); ++j) {
+            EXPECT_NEAR(testCoeffs(i+1, j+1), refCoeffs[index], 1e-12);
+            index++;
+        }
+    }
+}