diff --git a/src/Classic/Fields/Interpolation/SquarePolynomialVector.cpp b/src/Classic/Fields/Interpolation/SquarePolynomialVector.cpp
index b375fe636ca700accf7d41bfe16ed7f5687514ff..bc352bb28b1e662b9a3f71124c4ee999cad5d2bf 100644
--- a/src/Classic/Fields/Interpolation/SquarePolynomialVector.cpp
+++ b/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
diff --git a/src/Classic/Fields/Interpolation/SquarePolynomialVector.h b/src/Classic/Fields/Interpolation/SquarePolynomialVector.h
index 0ff53fca824b5d908c3156c0e73a64a70eb30963..46ffdfd38abf77737c24949ae442eef7441a5a26 100644
--- a/src/Classic/Fields/Interpolation/SquarePolynomialVector.h
+++ b/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;}
diff --git a/tests/classic_src/Fields/Interpolation/SquarePolynomialVectorTest.cpp b/tests/classic_src/Fields/Interpolation/SquarePolynomialVectorTest.cpp
index bad61c86ede53d21a68219086bb769f3ca571fa4..dcda97b99294a2a89dc938f4049c3668046a9415 100644
--- a/tests/classic_src/Fields/Interpolation/SquarePolynomialVectorTest.cpp
+++ b/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++;
+ }
+ }
+}