LCOV - code coverage report
Current view: top level - chart2/source/tools - PolynomialRegressionCurveCalculator.cxx (source / functions) Hit Total Coverage
Test: commit e02a6cb2c3e2b23b203b422e4e0680877f232636 Lines: 0 113 0.0 %
Date: 2014-04-14 Functions: 0 7 0.0 %
Legend: Lines: hit not hit

          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           0 : PolynomialRegressionCurveCalculator::PolynomialRegressionCurveCalculator()
      34           0 : {}
      35             : 
      36           0 : PolynomialRegressionCurveCalculator::~PolynomialRegressionCurveCalculator()
      37           0 : {}
      38             : 
      39             : // ____ XRegressionCurveCalculator ____
      40           0 : 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           0 :     rtl::math::setNan(&m_fCorrelationCoeffitient);
      46             : 
      47             :     RegressionCalculationHelper::tDoubleVectorPair aValues(
      48           0 :         RegressionCalculationHelper::cleanup( aXValues, aYValues, RegressionCalculationHelper::isValid()));
      49             : 
      50           0 :     const sal_Int32 aNoValues = aValues.first.size();
      51             : 
      52           0 :     const sal_Int32 aNoPowers = mForceIntercept ? mDegree : mDegree + 1;
      53             : 
      54           0 :     mCoefficients.clear();
      55           0 :     mCoefficients.resize(aNoPowers, 0.0);
      56             : 
      57           0 :     double yAverage = 0.0;
      58             : 
      59           0 :     std::vector<double> aQRTransposed;
      60           0 :     aQRTransposed.resize(aNoValues * aNoPowers, 0.0);
      61             : 
      62           0 :     std::vector<double> yVector;
      63           0 :     yVector.resize(aNoValues, 0.0);
      64             : 
      65           0 :     for(sal_Int32 i = 0; i < aNoValues; i++)
      66             :     {
      67           0 :         double yValue = aValues.second[i];
      68           0 :         if (mForceIntercept)
      69           0 :             yValue -= mInterceptValue;
      70           0 :         yVector[i] = yValue;
      71           0 :         yAverage += yValue;
      72             :     }
      73           0 :     yAverage /= aNoValues;
      74             : 
      75           0 :     for(sal_Int32 j = 0; j < aNoPowers; j++)
      76             :     {
      77           0 :         sal_Int32 aPower = mForceIntercept ? j+1 : j;
      78           0 :         sal_Int32 aColumnIndex = j * aNoValues;
      79           0 :         for(sal_Int32 i = 0; i < aNoValues; i++)
      80             :         {
      81           0 :             double xValue = aValues.first[i];
      82           0 :             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           0 :     sal_Int32 aMinorSize = std::min(aNoValues, aNoPowers);
      88             : 
      89           0 :     std::vector<double> aDiagonal;
      90           0 :     aDiagonal.resize(aMinorSize, 0.0);
      91             : 
      92             :     // Calculate Householder reflectors
      93           0 :     for (sal_Int32 aMinor = 0; aMinor < aMinorSize; aMinor++)
      94             :     {
      95           0 :         double aNormSqr = 0.0;
      96           0 :         for (sal_Int32 x = aMinor; x < aNoValues; x++)
      97             :         {
      98           0 :             double c = aQRTransposed[x + aMinor * aNoValues];
      99           0 :             aNormSqr += c * c;
     100             :         }
     101             : 
     102             :         double a;
     103             : 
     104           0 :         if (aQRTransposed[aMinor + aMinor * aNoValues] > 0.0)
     105           0 :             a = -std::sqrt(aNormSqr);
     106             :         else
     107           0 :             a = std::sqrt(aNormSqr);
     108             : 
     109           0 :         aDiagonal[aMinor] = a;
     110             : 
     111           0 :         if (a != 0.0)
     112             :         {
     113           0 :             aQRTransposed[aMinor + aMinor * aNoValues] -= a;
     114             : 
     115           0 :             for (sal_Int32 aColumn = aMinor + 1; aColumn < aNoPowers; aColumn++)
     116             :             {
     117           0 :                 double alpha = 0.0;
     118           0 :                 for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
     119             :                 {
     120           0 :                     alpha -= aQRTransposed[aRow + aColumn * aNoValues] * aQRTransposed[aRow + aMinor * aNoValues];
     121             :                 }
     122           0 :                 alpha /= a * aQRTransposed[aMinor + aMinor * aNoValues];
     123             : 
     124           0 :                 for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
     125             :                 {
     126           0 :                     aQRTransposed[aRow + aColumn * aNoValues] -= alpha * aQRTransposed[aRow + aMinor * aNoValues];
     127             :                 }
     128             :             }
     129             :         }
     130             :     }
     131             : 
     132             :     // Solve the linear equation
     133           0 :     for (sal_Int32 aMinor = 0; aMinor < aMinorSize; aMinor++)
     134             :     {
     135           0 :         double aDotProduct = 0;
     136             : 
     137           0 :         for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
     138             :         {
     139           0 :             aDotProduct += yVector[aRow] * aQRTransposed[aRow + aMinor * aNoValues];
     140             :         }
     141           0 :         aDotProduct /= aDiagonal[aMinor] * aQRTransposed[aMinor + aMinor * aNoValues];
     142             : 
     143           0 :         for (sal_Int32 aRow = aMinor; aRow < aNoValues; aRow++)
     144             :         {
     145           0 :             yVector[aRow] += aDotProduct * aQRTransposed[aRow + aMinor * aNoValues];
     146             :         }
     147             : 
     148             :     }
     149             : 
     150           0 :     for (sal_Int32 aRow = aDiagonal.size() - 1; aRow >= 0; aRow--)
     151             :     {
     152           0 :         yVector[aRow] /= aDiagonal[aRow];
     153           0 :         double yRow = yVector[aRow];
     154           0 :         mCoefficients[aRow] = yRow;
     155             : 
     156           0 :         for (sal_Int32 i = 0; i < aRow; i++)
     157             :         {
     158           0 :             yVector[i] -= yRow * aQRTransposed[i + aRow * aNoValues];
     159             :         }
     160             :     }
     161             : 
     162           0 :     if(mForceIntercept)
     163             :     {
     164           0 :         mCoefficients.insert(mCoefficients.begin(), mInterceptValue);
     165             :     }
     166             : 
     167             :     // Calculate correlation coeffitient
     168           0 :     double aSumError = 0.0;
     169           0 :     double aSumTotal = 0.0;
     170           0 :     double aSumYpred2 = 0.0;
     171             : 
     172           0 :     for( sal_Int32 i = 0; i < aNoValues; i++ )
     173             :     {
     174           0 :         double xValue = aValues.first[i];
     175           0 :         double yActual = aValues.second[i];
     176           0 :         double yPredicted = getCurveValue( xValue );
     177           0 :         aSumTotal += (yActual - yAverage) * (yActual - yAverage);
     178           0 :         aSumError += (yActual - yPredicted) * (yActual - yPredicted);
     179           0 :         if(mForceIntercept)
     180           0 :             aSumYpred2 += (yPredicted - mInterceptValue) * (yPredicted - mInterceptValue);
     181             :     }
     182             : 
     183           0 :     double aRSquared = 0.0;
     184           0 :     if(mForceIntercept)
     185             :     {
     186           0 :         aRSquared = aSumYpred2 / (aSumError + aSumYpred2);
     187             :     }
     188             :     else
     189             :     {
     190           0 :         aRSquared = 1.0 - (aSumError / aSumTotal);
     191             :     }
     192             : 
     193           0 :     if (aRSquared > 0.0)
     194           0 :         m_fCorrelationCoeffitient = std::sqrt(aRSquared);
     195             :     else
     196           0 :         m_fCorrelationCoeffitient = 0.0;
     197           0 : }
     198             : 
     199           0 : double SAL_CALL PolynomialRegressionCurveCalculator::getCurveValue( double x )
     200             :     throw (lang::IllegalArgumentException,
     201             :            uno::RuntimeException, std::exception)
     202             : {
     203             :     double fResult;
     204           0 :     rtl::math::setNan(&fResult);
     205             : 
     206           0 :     if (mCoefficients.empty())
     207             :     {
     208           0 :         return fResult;
     209             :     }
     210             : 
     211           0 :     sal_Int32 aNoCoefficients = (sal_Int32) mCoefficients.size();
     212             : 
     213             :     // Horner's method
     214           0 :     fResult = 0.0;
     215           0 :     for (sal_Int32 i = aNoCoefficients - 1; i >= 0; i--)
     216             :     {
     217           0 :         fResult = mCoefficients[i] + (x * fResult);
     218             :     }
     219           0 :     return fResult;
     220             : }
     221             : 
     222           0 : 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           0 :     return RegressionCurveCalculator::getCurveValues( min, max, nPointCount, xScalingX, xScalingY, bMaySkipPointsInCalculation );
     232             : }
     233             : 
     234           0 : OUString PolynomialRegressionCurveCalculator::ImplGetRepresentation(
     235             :     const uno::Reference< util::XNumberFormatter >& xNumFormatter,
     236             :     sal_Int32 nNumberFormatKey ) const
     237             : {
     238           0 :     OUStringBuffer aBuf( "f(x) = ");
     239             : 
     240           0 :     sal_Int32 aLastIndex = mCoefficients.size() - 1;
     241           0 :     for (sal_Int32 i = aLastIndex; i >= 0; i--)
     242             :     {
     243           0 :         double aValue = mCoefficients[i];
     244           0 :         if (aValue == 0.0)
     245             :         {
     246           0 :             continue;
     247             :         }
     248           0 :         else if (aValue < 0.0)
     249             :         {
     250           0 :             aBuf.appendAscii( " - " );
     251             :         }
     252             :         else
     253             :         {
     254           0 :             if (i != aLastIndex)
     255           0 :                 aBuf.appendAscii( " + " );
     256             :         }
     257             : 
     258           0 :         aBuf.append( getFormattedString( xNumFormatter, nNumberFormatKey, std::abs( aValue ) ) );
     259             : 
     260           0 :         if(i > 0)
     261             :         {
     262           0 :             if (i == 1)
     263             :             {
     264           0 :                 aBuf.appendAscii( "x" );
     265             :             }
     266             :             else
     267             :             {
     268           0 :                 aBuf.appendAscii( "x^" );
     269           0 :                 aBuf.append(i);
     270             :             }
     271             :         }
     272             :     }
     273             : 
     274           0 :     return aBuf.makeStringAndClear();
     275             : }
     276             : 
     277             : } //  namespace chart
     278             : 
     279             : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */

Generated by: LCOV version 1.10