summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--chart2/source/inc/PolynomialRegressionCurveCalculator.hxx3
-rw-r--r--chart2/source/tools/PolynomialRegressionCurveCalculator.cxx149
-rw-r--r--chart2/source/tools/gauss.hxx166
-rw-r--r--chart2/source/view/charttypes/VSeriesPlotter.cxx5
4 files changed, 109 insertions, 214 deletions
diff --git a/chart2/source/inc/PolynomialRegressionCurveCalculator.hxx b/chart2/source/inc/PolynomialRegressionCurveCalculator.hxx
index e3b43a26095f..6c1d990f9d58 100644
--- a/chart2/source/inc/PolynomialRegressionCurveCalculator.hxx
+++ b/chart2/source/inc/PolynomialRegressionCurveCalculator.hxx
@@ -25,6 +25,7 @@
namespace chart
{
+
class PolynomialRegressionCurveCalculator : public RegressionCurveCalculator
{
public:
@@ -57,7 +58,7 @@ private:
throw (com::sun::star::lang::IllegalArgumentException,
com::sun::star::uno::RuntimeException);
- std::vector<double> mResult;
+ std::vector<double> mCoefficients;
};
} // namespace chart
diff --git a/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx b/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
index 19d224d83cfb..1c84731141e4 100644
--- a/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
+++ b/chart2/source/tools/PolynomialRegressionCurveCalculator.cxx
@@ -24,11 +24,9 @@
#include <cmath>
#include <rtl/math.hxx>
#include <rtl/ustrbuf.hxx>
-#include "gauss.hxx"
using namespace com::sun::star;
-
namespace chart
{
@@ -49,80 +47,138 @@ void SAL_CALL PolynomialRegressionCurveCalculator::recalculateRegression(
RegressionCalculationHelper::tDoubleVectorPair aValues(
RegressionCalculationHelper::cleanup( aXValues, aYValues, RegressionCalculationHelper::isValid()));
- sal_Int32 aNoElements = mForceIntercept ? mDegree : mDegree + 1;
- sal_Int32 aNumberOfPowers = 2 * aNoElements - 1;
+ const sal_Int32 aNoValues = aValues.first.size();
- std::vector<double> aPowers;
- aPowers.resize(aNumberOfPowers, 0.0);
+ const sal_Int32 aNoPowers = mForceIntercept ? mDegree : mDegree + 1;
- sal_Int32 aNoColumns = aNoElements;
- sal_Int32 aNoRows = aNoElements + 1;
+ mCoefficients.clear();
+ mCoefficients.resize(aNoPowers, 0.0);
- std::vector<double> aMatrix;
- aMatrix.resize(aNoColumns * aNoRows, 0.0);
+ double yAverage = 0.0;
- const size_t aNoValues = aValues.first.size();
+ std::vector<double> aQRTransposed;
+ aQRTransposed.resize(aNoValues * aNoPowers, 0.0);
- double yAverage = 0.0;
+ std::vector<double> yVector;
+ yVector.resize(aNoValues, 0.0);
- for( size_t i = 0; i < aNoValues; ++i )
+ for(sal_Int32 i = 0; i < aNoValues; i++)
{
- double x = aValues.first[i];
- double y = aValues.second[i];
+ double yValue = aValues.second[i];
+ if (mForceIntercept)
+ yValue -= mInterceptValue;
+ yVector[i] = yValue;
+ yAverage += yValue;
+ }
+ yAverage /= aNoValues;
- for (sal_Int32 j = 0; j < aNumberOfPowers; j++)
+ 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++)
{
- if (mForceIntercept)
- aPowers[j] += std::pow(x, (int) j + 2);
- else
- aPowers[j] += std::pow(x, (int) j);
+ double xValue = aValues.first[i];
+ aQRTransposed[i + aColumnIndex] = std::pow(xValue, 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);
- for (sal_Int32 j = 0; j < aNoElements; j++)
+ // Calculate Householder reflectors
+ for (sal_Int32 aMinor = 0; aMinor < aMinorSize; aMinor++)
+ {
+ double aNormSqr = 0.0;
+ for (sal_Int32 x = aMinor; x < aNoValues; x++)
{
- if (mForceIntercept)
- aMatrix[j * aNoRows + aNoElements] += std::pow(x, (int) j + 1) * ( y - mInterceptValue );
- else
- aMatrix[j * aNoRows + aNoElements] += std::pow(x, (int) j) * y;
+ double c = aQRTransposed[x + aMinor * aNoValues];
+ aNormSqr += c * c;
}
- yAverage += y;
- }
+ double a;
+
+ if (aQRTransposed[aMinor + aMinor * aNoValues] > 0.0)
+ a = -std::sqrt(aNormSqr);
+ else
+ a = std::sqrt(aNormSqr);
+
+ aDiagonal[aMinor] = a;
- yAverage = yAverage / aNoValues;
+ if (a != 0.0)
+ {
+ aQRTransposed[aMinor + aMinor * aNoValues] -= a;
- for (sal_Int32 y = 0; y < aNoElements; y++)
+ 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++)
{
- for (sal_Int32 x = 0; x < aNoElements; x++)
+ 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++)
{
- aMatrix[y * aNoRows + x] = aPowers[y + x];
+ yVector[aRow] += aDotProduct * aQRTransposed[aRow + aMinor * aNoValues];
}
+
}
- mResult.clear();
- mResult.resize(aNoElements, 0.0);
+ for (sal_Int32 aRow = aDiagonal.size() - 1; aRow >= 0; aRow--)
+ {
+ yVector[aRow] /= aDiagonal[aRow];
+ double yRow = yVector[aRow];
+ mCoefficients[aRow] = yRow;
- solve(aMatrix, aNoColumns, aNoRows, mResult, 1.0e-20);
+ for (sal_Int32 i = 0; i < aRow; i++)
+ {
+ yVector[i] -= yRow * aQRTransposed[i + aRow * aNoValues];
+ }
+ }
- // Set intercept value if force intercept is enabled
- if (mForceIntercept) {
- mResult.insert( mResult.begin(), mInterceptValue );
+ if(mForceIntercept)
+ {
+ mCoefficients.insert(mCoefficients.begin(), mInterceptValue);
}
// Calculate correlation coeffitient
double aSumError = 0.0;
double aSumTotal = 0.0;
- for( size_t i = 0; i < aNoValues; ++i )
+ for( sal_Int32 i = 0; i < aNoValues; i++ )
{
- double x = aValues.first[i];
+ double xValue = aValues.first[i];
double yActual = aValues.second[i];
- double yPredicted = getCurveValue( x );
+ double yPredicted = getCurveValue( xValue );
aSumTotal += (yActual - yAverage) * (yActual - yAverage);
aSumError += (yActual - yPredicted) * (yActual - yPredicted);
}
double aRSquared = 1.0 - (aSumError / aSumTotal);
+
if (aRSquared > 0.0)
m_fCorrelationCoeffitient = std::sqrt(aRSquared);
else
@@ -136,15 +192,18 @@ double SAL_CALL PolynomialRegressionCurveCalculator::getCurveValue( double x )
double fResult;
rtl::math::setNan(&fResult);
- if (mResult.empty())
+ if (mCoefficients.empty())
{
return fResult;
}
+ sal_Int32 aNoCoefficients = (sal_Int32) mCoefficients.size();
+
+ // Horner's method
fResult = 0.0;
- for (size_t i = 0; i<mResult.size(); i++)
+ for (sal_Int32 i = aNoCoefficients - 1; i >= 0; i--)
{
- fResult += mResult[i] * std::pow(x, (int) i);
+ fResult = mCoefficients[i] + (x * fResult);
}
return fResult;
}
@@ -167,10 +226,10 @@ OUString PolynomialRegressionCurveCalculator::ImplGetRepresentation(
{
OUStringBuffer aBuf( "f(x) = ");
- sal_Int32 aLastIndex = mResult.size() - 1;
+ sal_Int32 aLastIndex = mCoefficients.size() - 1;
for (sal_Int32 i = aLastIndex; i >= 0; i--)
{
- double aValue = mResult[i];
+ double aValue = mCoefficients[i];
if (aValue == 0.0)
{
continue;
diff --git a/chart2/source/tools/gauss.hxx b/chart2/source/tools/gauss.hxx
deleted file mode 100644
index 6cd63ce30691..000000000000
--- a/chart2/source/tools/gauss.hxx
+++ /dev/null
@@ -1,166 +0,0 @@
-/* -*- 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 .
- */
-
-/** This method eliminates elements below main diagonal in the given
- matrix by gaussian elimination.
-
- @param matrix
- The matrix to operate on. Last column is the result vector (right
- hand side of the linear equation). After successful termination,
- the matrix is upper triangular. The matrix is expected to be in
- row major order.
-
- @param rows
- Number of rows in matrix
-
- @param cols
- Number of columns in matrix
-
- @param minPivot
- If the pivot element gets lesser than minPivot, this method fails,
- otherwise, elimination succeeds and true is returned.
-
- @return true, if elimination succeeded.
- */
-template <class Matrix, typename BaseType>
-bool eliminate( Matrix& matrix,
- int rows,
- int cols,
- const BaseType& minPivot )
-{
- BaseType temp;
- int max, i, j, k; /* *must* be signed, when looping like: j>=0 ! */
-
- /* eliminate below main diagonal */
- for(i=0; i<cols-1; ++i)
- {
- /* find best pivot */
- max = i;
- for(j=i+1; j<rows; ++j)
- if( fabs(matrix[ j*cols + i ]) > fabs(matrix[ max*cols + i ]) )
- max = j;
-
- /* check pivot value */
- if( fabs(matrix[ max*cols + i ]) < minPivot )
- return false; /* pivot too small! */
-
- /* interchange rows 'max' and 'i' */
- for(k=0; k<cols; ++k)
- {
- temp = matrix[ i*cols + k ];
- matrix[ i*cols + k ] = matrix[ max*cols + k ];
- matrix[ max*cols + k ] = temp;
- }
-
- /* eliminate column */
- for(j=i+1; j<rows; ++j)
- for(k=cols-1; k>=i; --k)
- matrix[ j*cols + k ] -= matrix[ i*cols + k ] *
- matrix[ j*cols + i ] / matrix[ i*cols + i ];
- }
-
- /* everything went well */
- return true;
-}
-
-
-/** Retrieve solution vector of linear system by substituting backwards.
-
- This operation _relies_ on the previous successful
- application of eliminate()!
-
- @param matrix
- Matrix in upper diagonal form, as e.g. generated by eliminate()
-
- @param rows
- Number of rows in matrix
-
- @param cols
- Number of columns in matrix
-
- @param result
- Result vector. Given matrix must have space for one column (rows entries).
-
- @return true, if back substitution was possible (i.e. no division
- by zero occurred).
- */
-template <class Matrix, class Vector, typename BaseType>
-bool substitute( const Matrix& matrix,
- int rows,
- int cols,
- Vector& result )
-{
- BaseType temp;
- int j,k; /* *must* be signed, when looping like: j>=0 ! */
-
- /* substitute backwards */
- for(j=rows-1; j>=0; --j)
- {
- temp = 0.0;
- for(k=j+1; k<cols-1; ++k)
- temp += matrix[ j*cols + k ] * result[k];
-
- if( matrix[ j*cols + j ] == 0.0 )
- return false; /* imminent division by zero! */
-
- result[j] = (matrix[ j*cols + cols-1 ] - temp) / matrix[ j*cols + j ];
- }
-
- /* everything went well */
- return true;
-}
-
-
-/** This method determines solution of given linear system, if any
-
- This is a wrapper for eliminate and substitute, given matrix must
- contain right side of equation as the last column.
-
- @param matrix
- The matrix to operate on. Last column is the result vector (right
- hand side of the linear equation). After successful termination,
- the matrix is upper triangular. The matrix is expected to be in
- row major order.
-
- @param rows
- Number of rows in matrix
-
- @param cols
- Number of columns in matrix
-
- @param minPivot
- If the pivot element gets lesser than minPivot, this method fails,
- otherwise, elimination succeeds and true is returned.
-
- @return true, if elimination succeeded.
- */
-template <class Matrix, class Vector, typename BaseType>
-bool solve( Matrix& matrix,
- int rows,
- int cols,
- Vector& result,
- BaseType minPivot )
-{
- if( eliminate<Matrix,BaseType>(matrix, rows, cols, minPivot) )
- return substitute<Matrix,Vector,BaseType>(matrix, rows, cols, result);
-
- return false;
-}
-
-/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
diff --git a/chart2/source/view/charttypes/VSeriesPlotter.cxx b/chart2/source/view/charttypes/VSeriesPlotter.cxx
index 6a9a62a33704..40bf94cab050 100644
--- a/chart2/source/view/charttypes/VSeriesPlotter.cxx
+++ b/chart2/source/view/charttypes/VSeriesPlotter.cxx
@@ -1025,12 +1025,13 @@ void VSeriesPlotter::createRegressionCurvesShapes( VDataSeries& rVDataSeries,
fPointScale = (fMaxX - fMinX) / (fChartMaxX - fChartMinX);
}
-
xCalculator->setRegressionProperties(aDegree, aForceIntercept, aInterceptValue, aPeriod);
xCalculator->recalculateRegression( rVDataSeries.getAllX(), rVDataSeries.getAllY() );
-
sal_Int32 nPointCount = 100 * fPointScale;
+ if ( nPointCount < 2 )
+ nPointCount = 2;
+
drawing::PolyPolygonShape3D aRegressionPoly;
aRegressionPoly.SequenceX.realloc(1);
aRegressionPoly.SequenceY.realloc(1);
generated by cgit (git 2.34.1) at 2024-08-12 22:21:20 +0000