Line data Source code
1 : /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 : /*
3 : * This file is part of the LibreOffice project.
4 : *
5 : * This Source Code Form is subject to the terms of the Mozilla Public
6 : * License, v. 2.0. If a copy of the MPL was not distributed with this
7 : * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 : *
9 : * This file incorporates work covered by the following license notice:
10 : *
11 : * Licensed to the Apache Software Foundation (ASF) under one or more
12 : * contributor license agreements. See the NOTICE file distributed
13 : * with this work for additional information regarding copyright
14 : * ownership. The ASF licenses this file to you under the Apache
15 : * License, Version 2.0 (the "License"); you may not use this file
16 : * except in compliance with the License. You may obtain a copy of
17 : * the License at http://www.apache.org/licenses/LICENSE-2.0 .
18 : */
19 :
20 : #include "PolynomialRegressionCurveCalculator.hxx"
21 : #include "macros.hxx"
22 : #include "RegressionCalculationHelper.hxx"
23 :
24 : #include <cmath>
25 : #include <rtl/math.hxx>
26 : #include <rtl/ustrbuf.hxx>
27 :
28 : using namespace com::sun::star;
29 :
30 : namespace chart
31 : {
32 :
33 28 : PolynomialRegressionCurveCalculator::PolynomialRegressionCurveCalculator()
34 28 : {}
35 :
36 42 : PolynomialRegressionCurveCalculator::~PolynomialRegressionCurveCalculator()
37 42 : {}
38 :
39 : // ____ XRegressionCurveCalculator ____
40 28 : void SAL_CALL PolynomialRegressionCurveCalculator::recalculateRegression(
41 : const uno::Sequence< double >& aXValues,
42 : const uno::Sequence< double >& aYValues )
43 : throw (uno::RuntimeException, std::exception)
44 : {
45 28 : rtl::math::setNan(&m_fCorrelationCoeffitient);
46 :
47 : RegressionCalculationHelper::tDoubleVectorPair aValues(
48 28 : RegressionCalculationHelper::cleanup( aXValues, aYValues, RegressionCalculationHelper::isValid()));
49 :
50 28 : const sal_Int32 aNoValues = aValues.first.size();
51 :
52 28 : const sal_Int32 aNoPowers = mForceIntercept ? mDegree : mDegree + 1;
53 :
54 28 : mCoefficients.clear();
55 28 : mCoefficients.resize(aNoPowers, 0.0);
56 :
57 28 : double yAverage = 0.0;
58 :
59 56 : std::vector<double> aQRTransposed;
60 28 : aQRTransposed.resize(aNoValues * aNoPowers, 0.0);
61 :
62 56 : std::vector<double> yVector;
63 28 : yVector.resize(aNoValues, 0.0);
64 :
65 224 : for(sal_Int32 i = 0; i < aNoValues; i++)
66 : {
67 196 : double yValue = aValues.second[i];
68 196 : if (mForceIntercept)
69 98 : yValue -= mInterceptValue;
70 196 : yVector[i] = yValue;
71 196 : yAverage += yValue;
72 : }
73 28 : yAverage /= aNoValues;
74 :
75 98 : for(sal_Int32 j = 0; j < aNoPowers; j++)
76 : {
77 70 : sal_Int32 aPower = mForceIntercept ? j+1 : j;
78 70 : sal_Int32 aColumnIndex = j * aNoValues;
79 560 : for(sal_Int32 i = 0; i < aNoValues; i++)
80 : {
81 490 : double xValue = aValues.first[i];
82 490 : aQRTransposed[i + aColumnIndex] = std::pow(xValue, (int) aPower);
83 : }
84 : }
85 :
86 : // QR decomposition - based on org.apache.commons.math.linear.QRDecomposition from apache commons math (ASF)
87 28 : sal_Int32 aMinorSize = std::min(aNoValues, aNoPowers);
88 :
89 56 : std::vector<double> aDiagonal;
90 28 : aDiagonal.resize(aMinorSize, 0.0);
91 :
92 : // Calculate Householder reflectors
93 98 : for (sal_Int32 aMinor = 0; aMinor < aMinorSize; aMinor++)
94 : {
95 70 : double aNormSqr = 0.0;
96 504 : for (sal_Int32 x = aMinor; x < aNoValues; x++)
97 : {
98 434 : double c = aQRTransposed[x + aMinor * aNoValues];
99 434 : aNormSqr += c * c;
100 : }
101 :
102 : double a;
103 :
104 70 : if (aQRTransposed[aMinor + aMinor * aNoValues] > 0.0)
105 28 : a = -std::sqrt(aNormSqr);
106 : else
107 42 : a = std::sqrt(aNormSqr);
108 :
109 70 : aDiagonal[aMinor] = a;
110 :
111 70 : if (a != 0.0)
112 : {
113 70 : aQRTransposed[aMinor + aMinor * aNoValues] -= a;
114 :
115 126 : for (sal_Int32 aColumn = aMinor + 1; aColumn < aNoPowers; aColumn++)
116 : {
117 56 : double alpha = 0.0;
118 434 : for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
119 : {
120 378 : alpha -= aQRTransposed[aRow + aColumn * aNoValues] * aQRTransposed[aRow + aMinor * aNoValues];
121 : }
122 56 : alpha /= a * aQRTransposed[aMinor + aMinor * aNoValues];
123 :
124 434 : for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
125 : {
126 378 : aQRTransposed[aRow + aColumn * aNoValues] -= alpha * aQRTransposed[aRow + aMinor * aNoValues];
127 : }
128 : }
129 : }
130 : }
131 :
132 : // Solve the linear equation
133 98 : for (sal_Int32 aMinor = 0; aMinor < aMinorSize; aMinor++)
134 : {
135 70 : double aDotProduct = 0;
136 :
137 504 : for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
138 : {
139 434 : aDotProduct += yVector[aRow] * aQRTransposed[aRow + aMinor * aNoValues];
140 : }
141 70 : aDotProduct /= aDiagonal[aMinor] * aQRTransposed[aMinor + aMinor * aNoValues];
142 :
143 504 : for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
144 : {
145 434 : yVector[aRow] += aDotProduct * aQRTransposed[aRow + aMinor * aNoValues];
146 : }
147 :
148 : }
149 :
150 98 : for (sal_Int32 aRow = aDiagonal.size() - 1; aRow >= 0; aRow--)
151 : {
152 70 : yVector[aRow] /= aDiagonal[aRow];
153 70 : double yRow = yVector[aRow];
154 70 : mCoefficients[aRow] = yRow;
155 :
156 126 : for (sal_Int32 i = 0; i < aRow; i++)
157 : {
158 56 : yVector[i] -= yRow * aQRTransposed[i + aRow * aNoValues];
159 : }
160 : }
161 :
162 28 : if(mForceIntercept)
163 : {
164 14 : mCoefficients.insert(mCoefficients.begin(), mInterceptValue);
165 : }
166 :
167 : // Calculate correlation coeffitient
168 28 : double aSumError = 0.0;
169 28 : double aSumTotal = 0.0;
170 28 : double aSumYpred2 = 0.0;
171 :
172 224 : for( sal_Int32 i = 0; i < aNoValues; i++ )
173 : {
174 196 : double xValue = aValues.first[i];
175 196 : double yActual = aValues.second[i];
176 196 : double yPredicted = getCurveValue( xValue );
177 196 : aSumTotal += (yActual - yAverage) * (yActual - yAverage);
178 196 : aSumError += (yActual - yPredicted) * (yActual - yPredicted);
179 196 : if(mForceIntercept)
180 98 : aSumYpred2 += (yPredicted - mInterceptValue) * (yPredicted - mInterceptValue);
181 : }
182 :
183 28 : double aRSquared = 0.0;
184 28 : if(mForceIntercept)
185 : {
186 14 : aRSquared = aSumYpred2 / (aSumError + aSumYpred2);
187 : }
188 : else
189 : {
190 14 : aRSquared = 1.0 - (aSumError / aSumTotal);
191 : }
192 :
193 28 : if (aRSquared > 0.0)
194 28 : m_fCorrelationCoeffitient = std::sqrt(aRSquared);
195 : else
196 28 : m_fCorrelationCoeffitient = 0.0;
197 28 : }
198 :
199 1274 : double SAL_CALL PolynomialRegressionCurveCalculator::getCurveValue( double x )
200 : throw (lang::IllegalArgumentException,
201 : uno::RuntimeException, std::exception)
202 : {
203 : double fResult;
204 1274 : rtl::math::setNan(&fResult);
205 :
206 1274 : if (mCoefficients.empty())
207 : {
208 0 : return fResult;
209 : }
210 :
211 1274 : sal_Int32 aNoCoefficients = (sal_Int32) mCoefficients.size();
212 :
213 : // Horner's method
214 1274 : fResult = 0.0;
215 6118 : for (sal_Int32 i = aNoCoefficients - 1; i >= 0; i--)
216 : {
217 4844 : fResult = mCoefficients[i] + (x * fResult);
218 : }
219 1274 : return fResult;
220 : }
221 :
222 14 : uno::Sequence< geometry::RealPoint2D > SAL_CALL PolynomialRegressionCurveCalculator::getCurveValues(
223 : double min, double max, sal_Int32 nPointCount,
224 : const uno::Reference< chart2::XScaling >& xScalingX,
225 : const uno::Reference< chart2::XScaling >& xScalingY,
226 : sal_Bool bMaySkipPointsInCalculation )
227 : throw (lang::IllegalArgumentException,
228 : uno::RuntimeException, std::exception)
229 : {
230 :
231 14 : return RegressionCurveCalculator::getCurveValues( min, max, nPointCount, xScalingX, xScalingY, bMaySkipPointsInCalculation );
232 : }
233 :
234 28 : OUString PolynomialRegressionCurveCalculator::ImplGetRepresentation(
235 : const uno::Reference< util::XNumberFormatter >& xNumFormatter,
236 : sal_Int32 nNumberFormatKey ) const
237 : {
238 28 : OUStringBuffer aBuf( "f(x) = ");
239 :
240 28 : sal_Int32 aLastIndex = mCoefficients.size() - 1;
241 112 : for (sal_Int32 i = aLastIndex; i >= 0; i--)
242 : {
243 84 : double aValue = mCoefficients[i];
244 84 : if (aValue == 0.0)
245 : {
246 0 : continue;
247 : }
248 84 : else if (aValue < 0.0)
249 : {
250 56 : aBuf.appendAscii( " - " );
251 : }
252 : else
253 : {
254 28 : if (i != aLastIndex)
255 28 : aBuf.appendAscii( " + " );
256 : }
257 :
258 84 : aBuf.append( getFormattedString( xNumFormatter, nNumberFormatKey, std::abs( aValue ) ) );
259 :
260 84 : if(i > 0)
261 : {
262 56 : if (i == 1)
263 : {
264 28 : aBuf.appendAscii( "x" );
265 : }
266 : else
267 : {
268 28 : aBuf.appendAscii( "x^" );
269 28 : aBuf.append(i);
270 : }
271 : }
272 : }
273 :
274 28 : return aBuf.makeStringAndClear();
275 : }
276 :
277 : } // namespace chart
278 :
279 : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|
