summaryrefslogtreecommitdiffstats
path: root/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
diff options
context:
space:
mode:
authorDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-07 09:06:44 +0000
committerDaniel Baumann <daniel.baumann@progress-linux.org>2024-04-07 09:06:44 +0000
commited5640d8b587fbcfed7dd7967f3de04b37a76f26 (patch)
tree7a5f7c6c9d02226d7471cb3cc8fbbf631b415303 /chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
parentInitial commit. (diff)
downloadlibreoffice-upstream.tar.xz
libreoffice-upstream.zip
Adding upstream version 4:7.4.7.upstream/4%7.4.7upstream
Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>
Diffstat (limited to 'chart2/source/tools/PolynomialRegressionCurveCalculator.cxx')
-rw-r--r--chart2/source/tools/PolynomialRegressionCurveCalculator.cxx392
1 files changed, 392 insertions, 0 deletions
diff --git a/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx b/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
new file mode 100644
index 000000000..3d7f1c076
--- /dev/null
+++ b/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
@@ -0,0 +1,392 @@
+/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+/*
+ * This file is part of the LibreOffice project.
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ *
+ * This file incorporates work covered by the following license notice:
+ *
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed
+ * with this work for additional information regarding copyright
+ * ownership. The ASF licenses this file to you under the Apache
+ * License, Version 2.0 (the "License"); you may not use this file
+ * except in compliance with the License. You may obtain a copy of
+ * the License at http://www.apache.org/licenses/LICENSE-2.0 .
+ */
+
+#include <PolynomialRegressionCurveCalculator.hxx>
+#include <RegressionCalculationHelper.hxx>
+
+#include <cmath>
+#include <limits>
+#include <rtl/math.hxx>
+#include <rtl/ustrbuf.hxx>
+
+#include <SpecialCharacters.hxx>
+
+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_fCorrelationCoefficient = std::sqrt(aRSquared);
+ else
+ m_fCorrelationCoefficient = 0.0;
+}
+
+// ____ XRegressionCurveCalculator ____
+void SAL_CALL PolynomialRegressionCurveCalculator::recalculateRegression(
+ const uno::Sequence< double >& aXValues,
+ const uno::Sequence< double >& aYValues )
+{
+ m_fCorrelationCoefficient = std::numeric_limits<double>::quiet_NaN();
+
+ RegressionCalculationHelper::tDoubleVectorPair aValues(
+ RegressionCalculationHelper::cleanup( aXValues, aYValues, RegressionCalculationHelper::isValid()));
+
+ const sal_Int32 aNoValues = aValues.first.size();
+
+ const sal_Int32 aNoPowers = mForceIntercept ? mDegree : mDegree + 1;
+
+ mCoefficients.clear();
+ mCoefficients.resize(aNoPowers, 0.0);
+
+ double yAverage = 0.0;
+
+ std::vector<double> yVector;
+ yVector.resize(aNoValues, 0.0);
+
+ for(sal_Int32 i = 0; i < aNoValues; i++)
+ {
+ double yValue = aValues.second[i];
+ if (mForceIntercept)
+ yValue -= mInterceptValue;
+ yVector[i] = yValue;
+ yAverage += yValue;
+ }
+ if (aNoValues != 0)
+ {
+ 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;
+ sal_Int32 aColumnIndex = j * aNoValues;
+ for(sal_Int32 i = 0; i < aNoValues; i++)
+ {
+ double xValue = aValues.first[i];
+ aQRTransposed[i + aColumnIndex] = std::pow(xValue, static_cast<int>(aPower));
+ }
+ }
+
+ // QR decomposition - based on org.apache.commons.math.linear.QRDecomposition from apache commons math (ASF)
+ sal_Int32 aMinorSize = std::min(aNoValues, aNoPowers);
+
+ std::vector<double> aDiagonal;
+ aDiagonal.resize(aMinorSize, 0.0);
+
+ // Calculate Householder reflectors
+ for (sal_Int32 aMinor = 0; aMinor < aMinorSize; aMinor++)
+ {
+ double aNormSqr = 0.0;
+ for (sal_Int32 x = aMinor; x < aNoValues; x++)
+ {
+ double c = aQRTransposed[x + aMinor * aNoValues];
+ aNormSqr += c * c;
+ }
+
+ double a;
+
+ if (aQRTransposed[aMinor + aMinor * aNoValues] > 0.0)
+ a = -std::sqrt(aNormSqr);
+ else
+ a = std::sqrt(aNormSqr);
+
+ aDiagonal[aMinor] = a;
+
+ if (a != 0.0)
+ {
+ aQRTransposed[aMinor + aMinor * aNoValues] -= a;
+
+ for (sal_Int32 aColumn = aMinor + 1; aColumn < aNoPowers; aColumn++)
+ {
+ double alpha = 0.0;
+ for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
+ {
+ alpha -= aQRTransposed[aRow + aColumn * aNoValues] * aQRTransposed[aRow + aMinor * aNoValues];
+ }
+ alpha /= a * aQRTransposed[aMinor + aMinor * aNoValues];
+
+ for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
+ {
+ aQRTransposed[aRow + aColumn * aNoValues] -= alpha * aQRTransposed[aRow + aMinor * aNoValues];
+ }
+ }
+ }
+ }
+
+ // Solve the linear equation
+ for (sal_Int32 aMinor = 0; aMinor < aMinorSize; aMinor++)
+ {
+ double aDotProduct = 0;
+
+ for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
+ {
+ aDotProduct += yVector[aRow] * aQRTransposed[aRow + aMinor * aNoValues];
+ }
+ aDotProduct /= aDiagonal[aMinor] * aQRTransposed[aMinor + aMinor * aNoValues];
+
+ for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
+ {
+ yVector[aRow] += aDotProduct * aQRTransposed[aRow + aMinor * aNoValues];
+ }
+
+ }
+
+ for (sal_Int32 aRow = aDiagonal.size() - 1; aRow >= 0; aRow--)
+ {
+ yVector[aRow] /= aDiagonal[aRow];
+ double yRow = yVector[aRow];
+ mCoefficients[aRow] = yRow;
+
+ for (sal_Int32 i = 0; i < aRow; i++)
+ {
+ yVector[i] -= yRow * aQRTransposed[i + aRow * aNoValues];
+ }
+ }
+
+ if(mForceIntercept)
+ {
+ mCoefficients.insert(mCoefficients.begin(), mInterceptValue);
+ }
+
+ // Calculate correlation coefficient
+ computeCorrelationCoefficient(aValues, aNoValues, yAverage);
+}
+
+double SAL_CALL PolynomialRegressionCurveCalculator::getCurveValue( double x )
+{
+ if (mCoefficients.empty())
+ return std::numeric_limits<double>::quiet_NaN();
+
+ sal_Int32 aNoCoefficients = static_cast<sal_Int32>(mCoefficients.size());
+
+ // Horner's method
+ double fResult = 0.0;
+ for (sal_Int32 i = aNoCoefficients - 1; i >= 0; i--)
+ {
+ fResult = mCoefficients[i] + (x * fResult);
+ }
+ return fResult;
+}
+
+OUString PolynomialRegressionCurveCalculator::ImplGetRepresentation(
+ const uno::Reference< util::XNumberFormatter >& xNumFormatter,
+ sal_Int32 nNumberFormatKey, sal_Int32* pFormulaMaxWidth /* = nullptr */ ) const
+{
+ OUStringBuffer aBuf( mYName + " = " );
+
+ sal_Int32 nValueLength=0;
+ sal_Int32 aLastIndex = mCoefficients.size() - 1;
+
+ if ( pFormulaMaxWidth && *pFormulaMaxWidth > 0 )
+ {
+ sal_Int32 nCharMin = aBuf.getLength(); // count characters different from coefficients
+ double nCoefficients = aLastIndex + 1.0; // number of coefficients
+ for (sal_Int32 i = aLastIndex; i >= 0; i--)
+ {
+ double aValue = mCoefficients[i];
+ if ( aValue == 0.0 )
+ { // do not count coefficient if it is 0
+ nCoefficients --;
+ continue;
+ }
+ if ( rtl::math::approxEqual( fabs( aValue ) , 1.0 ) )
+ { // do not count coefficient if it is 1
+ nCoefficients --;
+ if ( i == 0 ) // intercept = 1
+ nCharMin ++;
+ }
+ if ( i != aLastIndex )
+ nCharMin += 3; // " + "
+ if ( i > 0 )
+ {
+ nCharMin += mXName.getLength() + 1; // " x"
+ if ( i > 1 )
+ nCharMin +=1; // "^i"
+ if ( i >= 10 )
+ nCharMin ++; // 2 digits for i
+ }
+ }
+ nValueLength = ( *pFormulaMaxWidth - nCharMin ) / nCoefficients;
+ if ( nValueLength <= 0 )
+ nValueLength = 1;
+ }
+
+ bool bFindValue = false;
+ sal_Int32 nLineLength = aBuf.getLength();
+ for (sal_Int32 i = aLastIndex; i >= 0; i--)
+ {
+ double aValue = mCoefficients[i];
+ OUStringBuffer aTmpBuf(""); // temporary buffer
+ if (aValue == 0.0)
+ {
+ continue;
+ }
+ else if (aValue < 0.0)
+ {
+ if ( bFindValue ) // if it is not the first aValue
+ aTmpBuf.append( " " );
+ aTmpBuf.append( OUStringChar(aMinusSign) + " ");
+ aValue = - aValue;
+ }
+ else
+ {
+ if ( bFindValue ) // if it is not the first aValue
+ aTmpBuf.append( " + " );
+ }
+ bFindValue = true;
+
+ // if nValueLength not calculated then nullptr
+ sal_Int32* pValueLength = nValueLength ? &nValueLength : nullptr;
+ OUString aValueString = getFormattedString( xNumFormatter, nNumberFormatKey, aValue, pValueLength );
+ if ( i == 0 || aValueString != "1" ) // aValueString may be rounded to 1 if nValueLength is small
+ {
+ aTmpBuf.append( aValueString );
+ if ( i > 0 ) // insert blank between coefficient and x
+ aTmpBuf.append( " " );
+ }
+
+ if(i > 0)
+ {
+ aTmpBuf.append( mXName );
+ if (i > 1)
+ {
+ if (i < 10) // simple case if only one digit
+ aTmpBuf.append( aSuperscriptFigures[ i ] );
+ else
+ {
+ OUString aValueOfi = OUString::number( i );
+ for ( sal_Int32 n = 0; n < aValueOfi.getLength() ; n++ )
+ {
+ sal_Int32 nIndex = aValueOfi[n] - u'0';
+ aTmpBuf.append( aSuperscriptFigures[ nIndex ] );
+ }
+ }
+ }
+ }
+ addStringToEquation( aBuf, nLineLength, aTmpBuf, pFormulaMaxWidth );
+ }
+ if ( std::u16string_view(aBuf) == OUStringConcatenation( mYName + " = ") )
+ aBuf.append( "0" );
+
+ return aBuf.makeStringAndClear();
+}
+
+} // namespace chart
+
+/* vim:set shiftwidth=4 softtabstop=4 expandtab: */