tdf#128995: Special case for single variable regression...

like in LINEST implementation in Calc. Use a straightforward
regression solver in this case, so that it is easier to
handle the numerical error in the intercept term using
::rtl::math::approxSub().

Change-Id: I627c0c48e377cac5385a85050c4f472ee963f3d6
Reviewed-on: https://gerrit.libreoffice.org/85222
Tested-by: Jenkins
Reviewed-by: Tomaž Vajngerl <quikee@gmail.com>

1 files changed, 101 insertions, 36 deletions

diff --git a/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx b/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
index 8658f6004c1e..c1e17594a316 100644
--- a/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
+++ b/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
@@ -31,12 +31,60 @@ using namespace com::sun::star;
 namespace chart
 {
 
+static double lcl_GetDotProduct(std::vector<double>& aVec1, std::vector<double>& aVec2)
+{
+    double fResult = 0.0;
+    assert(aVec1.size() == aVec2.size());
+    for (size_t i = 0; i < aVec1.size(); ++i)
+        fResult += aVec1[i] * aVec2[i];
+    return fResult;
+}
+
 PolynomialRegressionCurveCalculator::PolynomialRegressionCurveCalculator()
 {}
 
 PolynomialRegressionCurveCalculator::~PolynomialRegressionCurveCalculator()
 {}
 
+void PolynomialRegressionCurveCalculator::computeCorrelationCoefficient(
+    RegressionCalculationHelper::tDoubleVectorPair& rValues,
+    const sal_Int32 aNoValues,
+    double yAverage )
+{
+    double aSumError = 0.0;
+    double aSumTotal = 0.0;
+    double aSumYpred2 = 0.0;
+
+    for( sal_Int32 i = 0; i < aNoValues; i++ )
+    {
+        double xValue = rValues.first[i];
+        double yActual = rValues.second[i];
+        double yPredicted = getCurveValue( xValue );
+        aSumTotal += (yActual - yAverage) * (yActual - yAverage);
+        aSumError += (yActual - yPredicted) * (yActual - yPredicted);
+        if(mForceIntercept)
+            aSumYpred2 += (yPredicted - mInterceptValue) * (yPredicted - mInterceptValue);
+    }
+
+    double aRSquared = 0.0;
+    if(mForceIntercept)
+    {
+        if (auto const div = aSumError + aSumYpred2)
+        {
+            aRSquared = aSumYpred2 / div;
+        }
+    }
+    else if (aSumTotal != 0.0)
+    {
+        aRSquared = 1.0 - (aSumError / aSumTotal);
+    }
+
+    if (aRSquared > 0.0)
+        m_fCorrelationCoeffitient = std::sqrt(aRSquared);
+    else
+        m_fCorrelationCoeffitient = 0.0;
+}
+
 // ____ XRegressionCurveCalculator ____
 void SAL_CALL PolynomialRegressionCurveCalculator::recalculateRegression(
     const uno::Sequence< double >& aXValues,
@@ -56,9 +104,6 @@ void SAL_CALL PolynomialRegressionCurveCalculator::recalculateRegression(
 
     double yAverage = 0.0;
 
-    std::vector<double> aQRTransposed;
-    aQRTransposed.resize(aNoValues * aNoPowers, 0.0);
-
     std::vector<double> yVector;
     yVector.resize(aNoValues, 0.0);
 
@@ -75,6 +120,57 @@ void SAL_CALL PolynomialRegressionCurveCalculator::recalculateRegression(
         yAverage /= aNoValues;
     }
 
+    // Special case for single variable regression like in LINEST
+    // implementation in Calc.
+    if (mDegree == 1)
+    {
+        std::vector<double> xVector;
+        xVector.resize(aNoValues, 0.0);
+        double xAverage = 0.0;
+
+        for(sal_Int32 i = 0; i < aNoValues; ++i)
+        {
+            double xValue = aValues.first[i];
+            xVector[i] = xValue;
+            xAverage += xValue;
+        }
+        if (aNoValues != 0)
+        {
+            xAverage /= aNoValues;
+        }
+
+        if (!mForceIntercept)
+        {
+            for (sal_Int32 i = 0; i < aNoValues; ++i)
+            {
+                xVector[i] -= xAverage;
+                yVector[i] -= yAverage;
+            }
+        }
+        double fSumXY = lcl_GetDotProduct(xVector, yVector);
+        double fSumX2 = lcl_GetDotProduct(xVector, xVector);
+
+        double fSlope = fSumXY / fSumX2;
+
+        if (!mForceIntercept)
+        {
+            mInterceptValue = ::rtl::math::approxSub(yAverage, fSlope * xAverage);
+            mCoefficients[0] = mInterceptValue;
+            mCoefficients[1] = fSlope;
+        }
+        else
+        {
+            mCoefficients[0] = fSlope;
+            mCoefficients.insert(mCoefficients.begin(), mInterceptValue);
+        }
+
+        computeCorrelationCoefficient(aValues, aNoValues, yAverage);
+        return;
+    }
+
+    std::vector<double> aQRTransposed;
+    aQRTransposed.resize(aNoValues * aNoPowers, 0.0);
+
     for(sal_Int32 j = 0; j < aNoPowers; j++)
     {
         sal_Int32 aPower = mForceIntercept ? j+1 : j;
@@ -167,39 +263,8 @@ void SAL_CALL PolynomialRegressionCurveCalculator::recalculateRegression(
         mCoefficients.insert(mCoefficients.begin(), mInterceptValue);
     }
 
-    // Calculate correlation coeffitient
-    double aSumError = 0.0;
-    double aSumTotal = 0.0;
-    double aSumYpred2 = 0.0;
-
-    for( sal_Int32 i = 0; i < aNoValues; i++ )
-    {
-        double xValue = aValues.first[i];
-        double yActual = aValues.second[i];
-        double yPredicted = getCurveValue( xValue );
-        aSumTotal += (yActual - yAverage) * (yActual - yAverage);
-        aSumError += (yActual - yPredicted) * (yActual - yPredicted);
-        if(mForceIntercept)
-            aSumYpred2 += (yPredicted - mInterceptValue) * (yPredicted - mInterceptValue);
-    }
-
-    double aRSquared = 0.0;
-    if(mForceIntercept)
-    {
-        if (auto const div = aSumError + aSumYpred2)
-        {
-            aRSquared = aSumYpred2 / div;
-        }
-    }
-    else if (aSumTotal != 0.0)
-    {
-        aRSquared = 1.0 - (aSumError / aSumTotal);
-    }
-
-    if (aRSquared > 0.0)
-        m_fCorrelationCoeffitient = std::sqrt(aRSquared);
-    else
-        m_fCorrelationCoeffitient = 0.0;
+    // Calculate correlation coefficient
+    computeCorrelationCoefficient(aValues, aNoValues, yAverage);
 }
 
 double SAL_CALL PolynomialRegressionCurveCalculator::getCurveValue( double x )