Hi Regina,
On Sun, 2010-10-24 at 18:13 +0200, Regina Henschel wrote:
I'm currently working on LINEST and have attached a draft to issue
Cool ! this is an awesome patch :-)
There is no mathematical problem, but I'm uncertain about coding style.
And your coding style looks beautiful, a great improvement over what
was there before.
The algorithms work on matrices. They have a lot of parts which are
nearly identical but the matrices are transposed. How to handle that?
Good question; and one to which I don't know the answer (sadly) -
whatever uses the least code, and yet performs well I suppose.
And I do not know how much comments are needed for those, who have to
maintain the code later, and if a separate documentation is needed.
Your level of commenting is great.
What are your opinions?
First - I'm sorry it took so long to get back to you; sad. Secondly
I've turned your changes into a patch (which I append). I split your
changes out into a new module - so that this (nice new) piece can be
licensed under the LGPLv3+/MPL combination - if you're happy with that.
Finally - you updated the CheckMatrix signature, but I didn't see the
header change for that; it'd be great to expand on the diff to include
those bits & return it.
It'd be wonderful to have your changes included - though we have a
feature freeze in 2 days now ;-)
Many thanks & looking forward to working with you [ do you hang out on
IRC ? it'd be great to introduce yourself there #libreoffice on
irc.freenode.net ].
Thanks !
Michael.
PS. do you think a self contained regression test in sc/qa/unit/ would
be good for this ? or does that belong better in a spreadsheet we can
load and compare results in ? [ it would be good to get some of them
into qa/unit so we can load / calculate and check the answers during
build I suspect ].
--
michael.meeks@novell.com <><, Pseudo Engineer, itinerant idiot
diff --git a/sc/source/core/tool/interpr5.cxx b/sc/source/core/tool/interpr5.cxx
index ff72817..6fe6546 100644
--- a/sc/source/core/tool/interpr5.cxx
+++ b/sc/source/core/tool/interpr5.cxx
@@ -34,7 +34,6 @@
#include <rtl/math.hxx>
#include <rtl/logfile.hxx>
#include <string.h>
-#include <math.h>
#include <stdio.h>
#if OSL_DEBUG_LEVEL > 1
@@ -2060,302 +2059,6 @@ void ScInterpreter::ScRGP()
RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::ScRGP" );
CalulateRGPRKP(FALSE);
}
-bool ScInterpreter::CheckMatrix(BOOL _bLOG,BOOL _bTrendGrowth,BYTE& nCase,SCSIZE& nCX,SCSIZE&
nCY,SCSIZE& nRX,SCSIZE& nRY,SCSIZE& M,SCSIZE& N,ScMatrixRef& pMatX,ScMatrixRef& pMatY)
-{
- nCX = 0;
- nCY = 0;
- nRX = 0;
- nRY = 0;
- M = 0;
- N = 0;
- pMatY->GetDimensions(nCY, nRY);
- const SCSIZE nCountY = nCY * nRY;
- for ( SCSIZE i = 0; i < nCountY; i++ )
- {
- if (!pMatY->IsValue(i))
- {
- PushIllegalArgument();
- return false;
- }
- }
-
- if ( _bLOG )
- {
- ScMatrixRef pNewY = pMatY->CloneIfConst();
- for (SCSIZE nElem = 0; nElem < nCountY; nElem++)
- {
- const double fVal = pNewY->GetDouble(nElem);
- if (fVal <= 0.0)
- {
- PushIllegalArgument();
- return false;
- }
- else
- pNewY->PutDouble(log(fVal), nElem);
- }
- pMatY = pNewY;
- }
-
- if (pMatX)
- {
- pMatX->GetDimensions(nCX, nRX);
- const SCSIZE nCountX = nCX * nRX;
- for ( SCSIZE i = 0; i < nCountX; i++ )
- if (!pMatX->IsValue(i))
- {
- PushIllegalArgument();
- return false;
- }
- if (nCX == nCY && nRX == nRY)
- nCase = 1; // einfache Regression
- else if (nCY != 1 && nRY != 1)
- {
- PushIllegalArgument();
- return false;
- }
- else if (nCY == 1)
- {
- if (nRX != nRY)
- {
- PushIllegalArgument();
- return false;
- }
- else
- {
- nCase = 2; // zeilenweise
- N = nRY;
- M = nCX;
- }
- }
- else if (nCX != nCY)
- {
- PushIllegalArgument();
- return false;
- }
- else
- {
- nCase = 3; // spaltenweise
- N = nCY;
- M = nRX;
- }
- }
- else
- {
- pMatX = GetNewMat(nCY, nRY);
- if ( _bTrendGrowth )
- {
- nCX = nCY;
- nRX = nRY;
- }
- if (!pMatX)
- {
- PushIllegalArgument();
- return false;
- }
- for ( SCSIZE i = 1; i <= nCountY; i++ )
- pMatX->PutDouble((double)i, i-1);
- nCase = 1;
- }
- return true;
-}
-void ScInterpreter::CalulateRGPRKP(BOOL _bRKP)
-{
- RTL_LOGFILE_CONTEXT_AUTHOR( aLogger, "sc", "er", "ScInterpreter::CheckMatrix" );
- BYTE nParamCount = GetByte();
- if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
- return;
- BOOL bConstant, bStats;
- if (nParamCount == 4)
- bStats = GetBool();
- else
- bStats = FALSE;
- if (nParamCount >= 3)
- bConstant = GetBool();
- else
- bConstant = TRUE;
- ScMatrixRef pMatX;
- ScMatrixRef pMatY;
- if (nParamCount >= 2)
- pMatX = GetMatrix();
- else
- pMatX = NULL;
- pMatY = GetMatrix();
- if (!pMatY)
- {
- PushIllegalParameter();
- return;
- } // if (!pMatY)
- BYTE nCase; // 1 = normal, 2,3 = mehrfach
- SCSIZE nCX, nCY;
- SCSIZE nRX, nRY;
- SCSIZE M = 0, N = 0;
- if ( !CheckMatrix(_bRKP,FALSE,nCase,nCX,nCY,nRX,nRY,M,N,pMatX,pMatY) )
- return;
-
- ScMatrixRef pResMat;
- if (nCase == 1)
- {
- if (!bStats)
- pResMat = GetNewMat(2,1);
- else
- pResMat = GetNewMat(2,5);
- if (!pResMat)
- {
- PushIllegalArgument();
- return;
- }
- double fCount = 0.0;
- double fSumX = 0.0;
- double fSumSqrX = 0.0;
- double fSumY = 0.0;
- double fSumSqrY = 0.0;
- double fSumXY = 0.0;
- double fValX, fValY;
- for (SCSIZE i = 0; i < nCY; i++)
- for (SCSIZE j = 0; j < nRY; j++)
- {
- fValX = pMatX->GetDouble(i,j);
- fValY = pMatY->GetDouble(i,j);
- fSumX += fValX;
- fSumSqrX += fValX * fValX;
- fSumY += fValY;
- fSumSqrY += fValY * fValY;
- fSumXY += fValX*fValY;
- fCount++;
- }
- if (fCount < 1.0)
- PushNoValue();
- else
- {
- double f1 = fCount*fSumXY-fSumX*fSumY;
- double fX = fCount*fSumSqrX-fSumX*fSumX;
- double b, m;
- if (bConstant)
- {
- b = fSumY/fCount - f1/fX*fSumX/fCount;
- m = f1/fX;
- }
- else
- {
- b = 0.0;
- m = fSumXY/fSumSqrX;
- }
- pResMat->PutDouble(_bRKP ? exp(m) : m, 0, 0);
- pResMat->PutDouble(_bRKP ? exp(b) : b, 1, 0);
- if (bStats)
- {
- double fY = fCount*fSumSqrY-fSumY*fSumY;
- double fSyx = fSumSqrY-b*fSumY-m*fSumXY;
- double fR2 = f1*f1/(fX*fY);
- pResMat->PutDouble (fR2, 0, 2);
- if (fCount < 3.0)
- {
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 1 );
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 1, 1 );
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 1, 2 );
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 3 );
- }
- else
- {
- pResMat->PutDouble(sqrt(fSyx*fCount/(fX*(fCount-2.0))), 0, 1);
- pResMat->PutDouble(sqrt(fSyx*fSumSqrX/fX/(fCount-2.0)), 1, 1);
- pResMat->PutDouble(
- sqrt((fCount*fSumSqrY - fSumY*fSumY - f1*f1/fX)/
- (fCount*(fCount-2.0))), 1, 2);
- if (fR2 == 1.0)
- pResMat->PutString(ScGlobal::GetRscString(STR_NO_VALUE), 0, 3 );
- else
- pResMat->PutDouble(fR2*(fCount-2.0)/(1.0-fR2), 0, 3);
- }
- pResMat->PutDouble(((double)(nCY*nRY))-2.0, 1, 3);
- pResMat->PutDouble(fY/fCount-fSyx, 0, 4);
- pResMat->PutDouble(fSyx, 1, 4);
- }
- }
- } // if (nCase == 1)
- if ( nCase != 1 )
- {
- SCSIZE i, j, k;
- if (!bStats)
- pResMat = GetNewMat(M+1,1);
- else
- pResMat = GetNewMat(M+1,5);
- if (!pResMat)
- {
- PushIllegalArgument();
- return;
- }
- ScMatrixRef pQ = GetNewMat(M+1, M+2);
- ScMatrixRef pE = GetNewMat(M+2, 1);
- ScMatrixRef pV = GetNewMat(M+1, 1);
- pE->PutDouble(0.0, M+1);
- pQ->FillDouble(0.0, 0, 0, M, M+1);
- if (nCase == 2)
- {
- for (k = 0; k < N; k++)
- {
- double Yk = pMatY->GetDouble(k);
- pE->PutDouble( pE->GetDouble(M+1)+Yk*Yk, M+1 );
- double sumYk = pQ->GetDouble(0, M+1) + Yk;
- pQ->PutDouble( sumYk, 0, M+1 );
- pE->PutDouble( sumYk, 0 );
- for (i = 0; i < M; i++)
- {
- double Xik = pMatX->GetDouble(i,k);
- double sumXik = pQ->GetDouble(0, i+1) + Xik;
- pQ->PutDouble( sumXik, 0, i+1);
- pQ->PutDouble( sumXik, i+1, 0);
- double sumXikYk = pQ->GetDouble(i+1, M+1) + Xik * Yk;
- pQ->PutDouble( sumXikYk, i+1, M+1);
- pE->PutDouble( sumXikYk, i+1);
- for (j = i; j < M; j++)
- {
- const double fVal = pMatX->GetDouble(j,k);
- double sumXikXjk = pQ->GetDouble(j+1, i+1) +
- Xik * fVal;
- pQ->PutDouble( sumXikXjk, j+1, i+1);
- pQ->PutDouble( sumXikXjk, i+1, j+1);
- }
- }
- }
- }
- else
- {
- for (k = 0; k < N; k++)
- {
- double Yk = pMatY->GetDouble(k);
- pE->PutDouble( pE->GetDouble(M+1)+Yk*Yk, M+1 );
- double sumYk = pQ->GetDouble(0, M+1) + Yk;
- pQ->PutDouble( sumYk, 0, M+1 );
- pE->PutDouble( sumYk, 0 );
- for (i = 0; i < M; i++)
- {
- double Xki = pMatX->GetDouble(k,i);
- double sumXki = pQ->GetDouble(0, i+1) + Xki;
- pQ->PutDouble( sumXki, 0, i+1);
- pQ->PutDouble( sumXki, i+1, 0);
- double sumXkiYk = pQ->GetDouble(i+1, M+1) + Xki * Yk;
- pQ->PutDouble( sumXkiYk, i+1, M+1);
- pE->PutDouble( sumXkiYk, i+1);
- for (j = i; j < M; j++)
- {
- const double fVal = pMatX->GetDouble(k,j);
- double sumXkiXkj = pQ->GetDouble(j+1, i+1) +
- Xki * fVal;
- pQ->PutDouble( sumXkiXkj, j+1, i+1);
- pQ->PutDouble( sumXkiXkj, i+1, j+1);
- }
- }
- }
- }
- if ( !Calculate4(_bRKP,pResMat,pQ,bConstant,N,M) )
- return;
-
- if (bStats)
- Calculate(pResMat,pE,pQ,pV,pMatX,bConstant,N,M,nCase);
- }
- PushMatrix(pResMat);
-}
void ScInterpreter::ScRKP()
{
diff --git a/sc/source/core/tool/interpr7.cxx b/sc/source/core/tool/interpr7.cxx
new file mode 100644
index 0000000..4b2c8b0
--- /dev/null
+++ b/sc/source/core/tool/interpr7.cxx
@@ -0,0 +1,1013 @@
+// MARKER(update_precomp.py): autogen include statement, do not remove
+#include "precompiled_sc.hxx"
+
+#include <rtl/math.hxx>
+#include "interpre.hxx"
+#include "globstr.hrc"
+
+// -----------------------------------------------------------------------------
+// Helper methods for LINEST/LOGEST and TREND/GROWTH
+// All matrices must already exist and have the needed size, no control tests
+// done. Those methodes, which names start with lcl_T, are adapted to case 3,
+// where Y (=observed values) is given as row.
+// Remember, ScMatrix matrices are zero based, index access (column,row).
+// -----------------------------------------------------------------------------
+
+namespace {
+// Multiply n x m Mat A with m x l Mat B to n x l Mat R
+void lcl_MFastMult(ScMatrixRef pA, ScMatrixRef pB, ScMatrixRef pR,
+ SCSIZE n, SCSIZE m, SCSIZE l)
+
+{
+ double sum;
+ for (SCSIZE row = 0; row < n; row++)
+ {
+ for (SCSIZE col = 0; col < l; col++)
+ { // result element(col, row) =sum[ (row of A) * (column of B)]
+ sum = 0.0;
+ for (SCSIZE k = 0; k < m; k++)
+ sum += pA->GetDouble(k,row) * pB->GetDouble(col,k);
+ pR->PutDouble(sum, col, row);
+ }
+ }
+}
+
+// <A;B> over all elements; uses the matrices as vectors of length M
+double lcl_GetSumProduct(ScMatrixRef pMatA, ScMatrixRef pMatB, SCSIZE nM)
+{
+ double fSum = 0.0;
+ for (SCSIZE i=0; i<nM; i++)
+ fSum += pMatA->GetDouble(i) * pMatB->GetDouble(i);
+ return fSum;
+}
+
+// Special version for use within QR decomposition.
+// Euclidean norm of column index C starting in row index R;
+// matrix A has count N rows.
+double lcl_GetColumnEuclideanNorm(ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE row=nR; row<nN; row++)
+ fNorm += (pMatA->GetDouble(nC,row)) * (pMatA->GetDouble(nC,row));
+ return sqrt(fNorm);
+}
+
+// Euclidean norm of row index R starting in column index C;
+// matrix A has count N columns.
+double lcl_TGetColumnEuclideanNorm(ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE col=nC; col<nN; col++)
+ fNorm += (pMatA->GetDouble(col,nR)) * (pMatA->GetDouble(col,nR));
+ return sqrt(fNorm);
+}
+
+// Special version for use within QR decomposition.
+// Maximum norm of column index C starting in row index R;
+// matrix A has count N rows.
+double lcl_GetColumnMaximumNorm(ScMatrixRef pMatA, SCSIZE nC, SCSIZE nR, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE row=nR; row<nN; row++)
+ if (fNorm < fabs(pMatA->GetDouble(nC,row)))
+ fNorm = fabs(pMatA->GetDouble(nC,row));
+ return fNorm;
+}
+
+// Maximum norm of row index R starting in col index C;
+// matrix A has count N columns.
+double lcl_TGetColumnMaximumNorm(ScMatrixRef pMatA, SCSIZE nR, SCSIZE nC, SCSIZE nN)
+{
+ double fNorm = 0.0;
+ for (SCSIZE col=nC; col<nN; col++)
+ if (fNorm < fabs(pMatA->GetDouble(col,nR)))
+ fNorm = fabs(pMatA->GetDouble(col,nR));
+ return fNorm;
+}
+
+// Special version for use within QR decomposition.
+// <A(Ca);B(Cb)> starting in row index R;
+// Ca and Cb are indices of columns, matrices A and B have count N rows.
+double lcl_GetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nCa,
+ ScMatrixRef pMatB, SCSIZE nCb, SCSIZE nR, SCSIZE nN)
+{
+ double fResult = 0.0;
+ for (SCSIZE row=nR; row<nN; row++)
+ fResult += pMatA->GetDouble(nCa,row) * pMatB->GetDouble(nCb,row);
+ return fResult;
+}
+
+// <A(Ra);B(Rb)> starting in column index C;
+// Ra and Rb are indices of rows, matrices A and B have count N columns.
+double lcl_TGetColumnSumProduct(ScMatrixRef pMatA, SCSIZE nRa,
+ ScMatrixRef pMatB, SCSIZE nRb, SCSIZE nC, SCSIZE nN)
+{
+ double fResult = 0.0;
+ for (SCSIZE col=nC; col<nN; col++)
+ fResult += pMatA->GetDouble(col,nRa) * pMatB->GetDouble(col,nRb);
+ return fResult;
+}
+
+double lcl_GetSign(double fValue)
+{
+ if (fValue < 0.0)
+ return -1.0;
+ else if (fValue > 0.0)
+ return 1.0;
+ else
+ return 0.0;
+}
+
+/* Calculates a QR decomposition with Householder reflection.
+ * For each NxK matrix A exists a decomposition A=Q*R with an orthogonal
+ * NxN matrix Q and a NxK matrix R.
+ * Q=H1*H2*...*Hk with Householder matrices H. Such a householder matrix can
+ * be build from a vector u by H=I-(2/u'u)*(u u'). This vectors u are returned
+ * in the columns of matrix A, overwriting the old content.
+ * The matrix R has a quadric upper part KxK with values in the upper right
+ * triangle and zeros in all other elements. Here the diagonal elements of R
+ * are stored in the vector R and the other upper right elements in the upper
+ * right of the matrix A.
+ * The function returns false, if calculation breaks. But because of round-off
+ * errors singularity is often not detected.
+ */
+bool lcl_CalculateQRdecomposition(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
+{
+ double fScale ;
+ double fEuclid ;
+ double fFactor ;
+ double fSignum ;
+ double fSum ;
+ // ScMatrix matrices are zero based, index access (column,row)
+ for (SCSIZE col = 0; col <nK; col++)
+ {
+ // calculate vector u of the householder transformation
+ fScale = lcl_GetColumnMaximumNorm(pMatA, col, col, nN);
+ if (fScale == 0.0)
+ {
+ // A is singular
+ return false;
+ }
+ for (SCSIZE row = col; row <nN; row++)
+ pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
+
+ fEuclid = lcl_GetColumnEuclideanNorm(pMatA, col, col, nN);
+ fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(col,col)));
+ fSignum = lcl_GetSign(pMatA->GetDouble(col,col));
+ pMatA->PutDouble( pMatA->GetDouble(col,col) + fSignum*fEuclid, col,col);
+ pVecR[col] = -fSignum * fScale * fEuclid;
+
+ // apply Householder transformation to A
+ for (SCSIZE c=col+1; c<nK; c++)
+ {
+ fSum =lcl_GetColumnSumProduct(pMatA, col, pMatA, c, col, nN);
+ for (SCSIZE row = col; row <nN; row++)
+ pMatA->PutDouble( pMatA->GetDouble(c,row)
+ - fSum * fFactor * pMatA->GetDouble(col,row), c, row);
+ }
+ }
+ return true;
+}
+
+// same with transposed matrix A, N is count of columns, K count of rows
+bool lcl_TCalculateQRdecomposition(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, SCSIZE nK, SCSIZE nN)
+{
+ double fScale ;
+ double fEuclid ;
+ double fFactor ;
+ double fSignum ;
+ double fSum ;
+ // ScMatrix matrices are zero based, index access (column,row)
+ for (SCSIZE row = 0; row <nK; row++)
+ {
+ // calculate vector u of the householder transformation
+ fScale = lcl_TGetColumnMaximumNorm(pMatA, row, row, nN);
+ if (fScale == 0.0)
+ {
+ // A is singular
+ return false;
+ }
+ for (SCSIZE col = row; col <nN; col++)
+ pMatA->PutDouble( pMatA->GetDouble(col,row)/fScale, col, row);
+
+ fEuclid = lcl_TGetColumnEuclideanNorm(pMatA, row, row, nN);
+ fFactor = 1.0/fEuclid/(fEuclid + fabs(pMatA->GetDouble(row,row)));
+ fSignum = lcl_GetSign(pMatA->GetDouble(row,row));
+ pMatA->PutDouble( pMatA->GetDouble(row,row) + fSignum*fEuclid, row,row);
+ pVecR[row] = -fSignum * fScale * fEuclid;
+
+ // apply Householder transformation to A
+ for (SCSIZE r=row+1; r<nK; r++)
+ {
+ fSum =lcl_TGetColumnSumProduct(pMatA, row, pMatA, r, row, nN);
+ for (SCSIZE col = row; col <nN; col++)
+ pMatA->PutDouble( pMatA->GetDouble(col,r)
+ - fSum * fFactor * pMatA->GetDouble(col,row), col, r);
+ }
+ }
+ return true;
+}
+
+
+/* Applies a Householder transformation to a column vector Y with is given as
+ * Nx1 Matrix. The Vektor u, from which the Householder transformation is build,
+ * is the column part in matrix A, with column index C, starting with row
+ * index C. A is the result of the QR decomposition as obtained from
+ * lcl_CaluclateQRdecomposition.
+ */
+void lcl_ApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nC,
+ ScMatrixRef pMatY, SCSIZE nN)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fDenominator = lcl_GetColumnSumProduct(pMatA, nC, pMatA, nC, nC, nN);
+ double fNumerator = lcl_GetColumnSumProduct(pMatA, nC, pMatY, 0, nC, nN);
+ double fFactor = 2.0 * (fNumerator/fDenominator);
+ for (SCSIZE row = nC; row < nN; row++)
+ pMatY->PutDouble(
+ pMatY->GetDouble(row) - fFactor * pMatA->GetDouble(nC,row), row);
+}
+
+// Same with transposed matrices A and Y.
+void lcl_TApplyHouseholderTransformation(ScMatrixRef pMatA, SCSIZE nR,
+ ScMatrixRef pMatY, SCSIZE nN)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fDenominator = lcl_TGetColumnSumProduct(pMatA, nR, pMatA, nR, nR, nN);
+ double fNumerator = lcl_TGetColumnSumProduct(pMatA, nR, pMatY, 0, nR, nN);
+ double fFactor = 2.0 * (fNumerator/fDenominator);
+ for (SCSIZE col = nR; col < nN; col++)
+ pMatY->PutDouble(
+ pMatY->GetDouble(col) - fFactor * pMatA->GetDouble(col,nR), col);
+}
+
+/* Solve for X in R*X=S using back substitution. The solution X overwrites S.
+ * Uses R from the result of the QR decomposition of a NxK matrix A.
+ * S is a column vector given as matrix, with at least elements on index
+ * 0 to K-1; elements on index>=K are ignored. Vector R must not have zero
+ * elements, no check is done.
+ */
+void lcl_SolveWithUpperRightTriangle(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, ScMatrixRef pMatS,
+ SCSIZE nK, bool bIsTransposed)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fSum;
+ SCSIZE row;
+ // SCSIZE is never negative, therefore test with rowp1=row+1
+ for (SCSIZE rowp1 = nK; rowp1>0; rowp1--)
+ {
+ row = rowp1-1;
+ fSum = pMatS->GetDouble(row);
+ for (SCSIZE col = rowp1; col<nK ; col++)
+ if (bIsTransposed)
+ fSum -= pMatA->GetDouble(row,col) * pMatS->GetDouble(col);
+ else
+ fSum -= pMatA->GetDouble(col,row) * pMatS->GetDouble(col);
+ pMatS->PutDouble( fSum / pVecR[row] , row);
+ }
+}
+
+/* Solve for X in R' * X= T using forward substitution. The solution X
+ * overwrites T. Uses R from the result of the QR decomposition of a NxK
+ * matrix A. T is a column vectors given as matrix, with at least elements on
+ * index 0 to K-1; elements on index>=K are ignored. Vector R must not have
+ * zero elements, no check is done.
+ */
+void lcl_SolveWithLowerLeftTriangle(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, ScMatrixRef pMatT,
+ SCSIZE nK, bool bIsTransposed)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fSum;
+ for (SCSIZE row = 0; row < nK; row++)
+ {
+ fSum = pMatT -> GetDouble(row);
+ for (SCSIZE col=0; col < row; col++)
+ {
+ if (bIsTransposed)
+ fSum -= pMatA->GetDouble(col,row) * pMatT->GetDouble(col);
+ else
+ fSum -= pMatA->GetDouble(row,col) * pMatT->GetDouble(col);
+ }
+ pMatT->PutDouble( fSum / pVecR[row] , row);
+ }
+}
+
+/* Calculates Z = R * B
+ * R is given in matrix A and vector VecR as obtained from the QR
+ * decompostion in lcl_CalculateQRdecomposition. B and Z are column vectors
+ * given as matrix with at least index 0 to K-1; elements on index>=K are
+ * not used.
+ */
+void lcl_ApplyUpperRightTriangle(ScMatrixRef pMatA,
+ ::std::vector< double>& pVecR, ScMatrixRef pMatB,
+ ScMatrixRef pMatZ, SCSIZE nK, bool bIsTransposed)
+{
+ // ScMatrix matrices are zero based, index access (column,row)
+ double fSum;
+ for (SCSIZE row = 0; row < nK; row++)
+ {
+ fSum = pVecR[row] * pMatB->GetDouble(row);
+ for (SCSIZE col = row+1; col < nK; col++)
+ if (bIsTransposed)
+ fSum += pMatA->GetDouble(row,col) * pMatB->GetDouble(col);
+ else
+ fSum += pMatA->GetDouble(col,row) * pMatB->GetDouble(col);
+ pMatZ->PutDouble( fSum, row);
+ }
+}
+
+double lcl_GetMeanOverAll(ScMatrixRef pMat, SCSIZE nN)
+{
+ double fSum = 0.0;
+ for (SCSIZE i=0 ; i<nN; i++)
+ fSum += pMat->GetDouble(i);
+ return fSum/static_cast<double>(nN);
+}
+
+// Calculates means of the columns of matrix X. X is a RxC matrix;
+// ResMat is a 1xC matrix (=row).
+void lcl_CalculateColumnMeans(ScMatrixRef pX, ScMatrixRef pResMat,
+ SCSIZE nC, SCSIZE nR)
+{
+ double fSum = 0.0;
+ for (SCSIZE i=0; i < nC; i++)
+ {
+ fSum =0.0;
+ for (SCSIZE k=0; k < nR; k++)
+ fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
+ pResMat ->PutDouble( fSum/static_cast<double>(nR),i);
+ }
+}
+
+// Calculates means of the rows of matrix X. X is a RxC matrix;
+// ResMat is a Rx1 matrix (=column).
+void lcl_CalculateRowMeans(ScMatrixRef pX, ScMatrixRef pResMat,
+ SCSIZE nC, SCSIZE nR)
+{
+ double fSum = 0.0;
+ for (SCSIZE k=0; k < nR; k++)
+ {
+ fSum =0.0;
+ for (SCSIZE i=0; i < nC; i++)
+ fSum += pX->GetDouble(i,k); // GetDouble(Column,Row)
+ pResMat ->PutDouble( fSum/static_cast<double>(nC),k);
+ }
+}
+
+void lcl_CalculateColumnsDelta(ScMatrixRef pMat, ScMatrixRef pColumnMeans,
+ SCSIZE nC, SCSIZE nR)
+{
+ for (SCSIZE i = 0; i < nC; i++)
+ for (SCSIZE k = 0; k < nR; k++)
+ pMat->PutDouble( ::rtl::math::approxSub
+ (pMat->GetDouble(i,k) , pColumnMeans->GetDouble(i) ) , i, k);
+}
+
+void lcl_CalculateRowsDelta(ScMatrixRef pMat, ScMatrixRef pRowMeans,
+ SCSIZE nC, SCSIZE nR)
+{
+ for (SCSIZE k = 0; k < nR; k++)
+ for (SCSIZE i = 0; i < nC; i++)
+ pMat->PutDouble( ::rtl::math::approxSub
+ ( pMat->GetDouble(i,k) , pRowMeans->GetDouble(k) ) , i, k);
+}
+
+// Case1 = simple regression
+// MatX = X - MeanX, MatY = Y - MeanY, y - haty = (y - MeanY) - (haty - MeanY)
+// = (y-MeanY)-((slope*x+a)-(slope*MeanX+a)) = (y-MeanY)-slope*(x-MeanX)
+double lcl_GetSSresid(ScMatrixRef pMatX, ScMatrixRef pMatY, double fSlope,
+ SCSIZE nN)
+{
+ double fSum = 0.0;
+ double fTemp = 0.0;
+ for (SCSIZE i=0; i<nN; i++)
+ {
+ fTemp = pMatY->GetDouble(i) - fSlope * pMatX->GetDouble(i);
+ fSum += fTemp * fTemp;
+ }
+ return fSum;
+}
+} // private namespace
+
+void ScInterpreter::CalulateRGPRKP(BOOL _bRKP)
+{
+ BYTE nParamCount = GetByte();
+ if ( !MustHaveParamCount( nParamCount, 1, 4 ) )
+ return;
+ bool bConstant, bStats;
+
+// optional forth parameter
+ if (nParamCount == 4)
+ bStats = GetBool();
+ else
+ bStats = false;
+
+// The third parameter may not be missing in ODF, if the forth parameter
+// is present.
+ if (nParamCount >= 3)
+ {
+ if (IsMissing())
+ {
+ PushIllegalParameter();
+ return;
+ }
+ else
+ bConstant = GetBool();
+ }
+ else
+ bConstant = TRUE;
+
+ ScMatrixRef pMatX;
+ if (nParamCount >= 2)
+ {
+ if (IsMissing())
+ { //In ODF1.2 empty second parameter (which is two ;; ) is allowed
+ Pop();
+ pMatX = NULL;
+ }
+ else
+ {
+ pMatX = GetMatrix();
+ }
+ }
+ else
+ pMatX = NULL;
+
+ ScMatrixRef pMatY;
+ pMatY = GetMatrix();
+ if (!pMatY)
+ {
+ PushIllegalParameter();
+ return;
+ }
+
+ // 1 = simple; 2 = multiple with Y as column; 3 = multiple with Y as row
+ BYTE nCase;
+
+ SCSIZE nCX, nCY; // number of columns
+ SCSIZE nRX, nRY; //number of rows
+ SCSIZE K = 0, N = 0; // K=number of variables X, N=number of data samples
+
+ // Fill default values in matrix X, transform Y to log(Y) in case _bLOGEST,
+ // determine sizes of matrices X and Y, determine kind of regression, clone
+ // Y in case bLOGEST, if constant.
+ if ( !CheckMatrix(_bRKP,nCase,nCX,nCY,nRX,nRY,K,N,pMatX,pMatY) )
+ {
+ PushIllegalParameter();
+ return;
+ }
+
+ // Enough data samples?
+ if ( (bConstant && (N<K+1)) || (!bConstant && (N<K)) || (N<1) || (K<1) )
+ {
+ PushIllegalParameter();
+ return;
+ }
+
+ ScMatrixRef pResMat;
+ if (bStats)
+ pResMat = GetNewMat(K+1,5);
+ else
+ pResMat = GetNewMat(K+1,1);
+ if (!pResMat)
+ {
+ PushError(errStackOverflow);
+ return;
+ }
+ // Fill unused cells in pResMat; order (column,row)
+ if (bStats)
+ {
+ for (SCSIZE i=2; i<K+1; i++)
+ {
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 2 );
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 3 );
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), i, 4 );
+ }
+ }
+
+ // Uses sum(x-MeanX)^2 and not [sum x^2]-N * MeanX^2 in case bConstant.
+ // Clone constant matrices, so that Mat = Mat - Mean is possible.
+ double fMeanY = 0.0;
+ if (bConstant)
+ {
+ ScMatrixRef pNewX = pMatX->CloneIfConst();
+ ScMatrixRef pNewY = pMatY->CloneIfConst();
+ if (!pNewX || !pNewY)
+ {
+ PushError(errStackOverflow);
+ return;
+ }
+ pMatX = pNewX;
+ pMatY = pNewY;
+ // DeltaY is possible here; DeltaX depends on nCase, so later
+ fMeanY = lcl_GetMeanOverAll(pMatY, N);
+ for (SCSIZE i=0; i<N; i++)
+ {
+ pMatY->PutDouble( ::rtl::math::approxSub(pMatY->GetDouble(i),fMeanY), i );
+ }
+ }
+
+ if (nCase==1)
+ {
+ // calculate simple regression
+ double fMeanX = 0.0;
+ if (bConstant)
+ { // Mat = Mat - Mean
+ fMeanX = lcl_GetMeanOverAll(pMatX, N);
+ for (SCSIZE i=0; i<N; i++)
+ {
+ pMatX->PutDouble( ::rtl::math::approxSub(pMatX->GetDouble(i),fMeanX), i );
+ }
+ }
+ double fSumXY = lcl_GetSumProduct(pMatX,pMatY,N);
+ double fSumX2 = lcl_GetSumProduct(pMatX,pMatX,N);
+ if (fSumX2==0.0)
+ {
+ PushNoValue(); // all x-values are identical
+ return;
+ }
+ double fSlope = fSumXY / fSumX2;
+ double fIntercept = 0.0;
+ if (bConstant)
+ fIntercept = fMeanY - fSlope * fMeanX;
+ pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, 1, 0); //order (column,row)
+ pResMat->PutDouble(_bRKP ? exp(fSlope) : fSlope, 0, 0);
+
+ if (bStats)
+ {
+ double fSSreg = fSlope * fSlope * fSumX2;
+ pResMat->PutDouble(fSSreg, 0, 4);
+
+ double fDegreesFreedom =static_cast<double>( (bConstant) ? N-2 : N-1 );
+ pResMat->PutDouble(fDegreesFreedom, 1, 3);
+
+ double fSSresid = lcl_GetSSresid(pMatX,pMatY,fSlope,N);
+ pResMat->PutDouble(fSSresid, 1, 4);
+
+ if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
+ { // exact fit; test SSreg too, because SSresid might be
+ // unequal zero due to round of errors
+ pResMat->PutDouble(0.0, 1, 4); // SSresid
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
+ pResMat->PutDouble(0.0, 1, 2); // RMSE
+ pResMat->PutDouble(0.0, 0, 1); // SigmaSlope
+ if (bConstant)
+ pResMat->PutDouble(0.0, 1, 1); //SigmaIntercept
+ else
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
+ pResMat->PutDouble(1.0, 0, 2); // R^2
+ }
+ else
+ {
+ double fFstatistic = (fSSreg / static_cast<double>(K))
+ / (fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fFstatistic, 0, 3);
+
+ // standard error of estimate
+ double fRMSE = sqrt(fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fRMSE, 1, 2);
+
+ double fSigmaSlope = fRMSE / sqrt(fSumX2);
+ pResMat->PutDouble(fSigmaSlope, 0, 1);
+
+ if (bConstant)
+ {
+ double fSigmaIntercept = fRMSE
+ * sqrt(fMeanX*fMeanX/fSumX2 + 1.0/static_cast<double>(N));
+ pResMat->PutDouble(fSigmaIntercept, 1, 1);
+ }
+ else
+ {
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 1, 1);
+ }
+
+ double fR2 = fSSreg / (fSSreg + fSSresid);
+ pResMat->PutDouble(fR2, 0, 2);
+ }
+ }
+ PushMatrix(pResMat);
+ }
+ else // calculate multiple regression;
+ {
+ // Uses a QR decomposition X = QR. The solution B = (X'X)^(-1) * X' * Y
+ // becomes B = R^(-1) * Q' * Y
+ if (nCase ==2) // Y is column
+ {
+ ::std::vector< double> aVecR(N); // for QR decomposition
+ // Enough memory for needed matrices?
+ ScMatrixRef pMeans = GetNewMat(K, 1); // mean of each column
+ ScMatrixRef pMatZ; // for Q' * Y , inter alia
+ if (bStats)
+ pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
+ else
+ pMatZ = pMatY; // Y can be overwritten
+ ScMatrixRef pSlopes = GetNewMat(1,K); // from b1 to bK
+ if (!pMeans || !pMatZ || !pSlopes)
+ {
+ PushError(errStackOverflow);
+ return;
+ }
+ if (bConstant)
+ {
+ lcl_CalculateColumnMeans(pMatX, pMeans, K, N);
+ lcl_CalculateColumnsDelta(pMatX, pMeans, K, N);
+ }
+ if (!lcl_CalculateQRdecomposition(pMatX, aVecR, K, N))
+ {
+ PushNoValue();
+ return;
+ }
+ // Later on we will divide by elements of aVecR, so make sure
+ // that they aren't zero.
+ bool bIsSingular=false;
+ for (SCSIZE row=0; row < K && !bIsSingular; row++)
+ bIsSingular = bIsSingular || aVecR[row]==0.0;
+ if (bIsSingular)
+ {
+ PushNoValue();
+ return;
+ }
+ // Z = Q' Y;
+ for (SCSIZE col = 0; col < K; col++)
+ {
+ lcl_ApplyHouseholderTransformation(pMatX, col, pMatZ, N);
+ }
+ // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
+ // result Z should have zeros for index>=K; if not, ignore values
+ for (SCSIZE col = 0; col < K ; col++)
+ {
+ pSlopes->PutDouble( pMatZ->GetDouble(col), col);
+ }
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, false);
+ double fIntercept = 0.0;
+ if (bConstant)
+ fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
+ // Fill first line in result matrix
+ pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
+ for (SCSIZE i = 0; i < K; i++)
+ pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
+ : pSlopes->GetDouble(i) , K-1-i, 0);
+
+
+ if (bStats)
+ {
+ double fSSreg = 0.0;
+ double fSSresid = 0.0;
+ // re-use memory of Z;
+ pMatZ->FillDouble(0.0, 0, 0, 0, N-1);
+ // Z = R * Slopes
+ lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, false);
+ // Z = Q * Z, that is Q * R * Slopes = X * Slopes
+ for (SCSIZE colp1 = K; colp1 > 0; colp1--)
+ {
+ lcl_ApplyHouseholderTransformation(pMatX, colp1-1, pMatZ,N);
+ }
+ fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
+ // re-use Y for residuals, Y = Y-Z
+ for (SCSIZE row = 0; row < N; row++)
+ pMatY->PutDouble(pMatY->GetDouble(row) - pMatZ->GetDouble(row), row);
+ fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
+ pResMat->PutDouble(fSSreg, 0, 4);
+ pResMat->PutDouble(fSSresid, 1, 4);
+
+ double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K );
+ pResMat->PutDouble(fDegreesFreedom, 1, 3);
+
+ if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
+ { // exact fit; incl. observed values Y are identical
+ pResMat->PutDouble(0.0, 1, 4); // SSresid
+ // F = (SSreg/K) / (SSresid/df) = #DIV/0!
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
+ // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
+ pResMat->PutDouble(0.0, 1, 2); // RMSE
+ // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
+ for (SCSIZE i=0; i<K; i++)
+ pResMat->PutDouble(0.0, K-1-i, 1);
+
+ // SigmaIntercept = RMSE * sqrt(...) = 0
+ if (bConstant)
+ pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
+ else
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
+
+ // R^2 = SSreg / (SSreg + SSresid) = 1.0
+ pResMat->PutDouble(1.0, 0, 2); // R^2
+ }
+ else
+ {
+ double fFstatistic = (fSSreg / static_cast<double>(K))
+ / (fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fFstatistic, 0, 3);
+
+ // standard error of estimate = root mean SSE
+ double fRMSE = sqrt(fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fRMSE, 1, 2);
+
+ // standard error of slopes
+ // = RMSE * sqrt(diagonal element of (R' R)^(-1) )
+ // standard error of intercept
+ // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
+ // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
+ // a whole matrix, but iterate over unit vectors.
+ double fSigmaSlope = 0.0;
+ double fSigmaIntercept = 0.0;
+ for (SCSIZE col = 0; col < K; col++)
+ {
+ //re-use memory of MatZ
+ pMatZ->FillDouble(0.0,0,0,0,K-1); // Z = unit vector e
+ pMatZ->PutDouble(1.0, col);
+ //Solve R' * Z = e
+ lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, false);
+ // Solve R * Znew = Zold
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, false);
+ // now Z is column col in (R' R)^(-1)
+ fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(col));
+ pResMat->PutDouble(fSigmaSlope, K-1-col, 1);
+ // (R' R) ^(-1) is symmetric
+ if (bConstant)
+ fSigmaIntercept += lcl_GetSumProduct(pMeans, pMatZ, K);
+ }
+ if (bConstant)
+ {
+ fSigmaIntercept = fRMSE
+ * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+ pResMat->PutDouble(fSigmaIntercept, K, 1);
+ }
+ else
+ {
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
+ }
+
+ double fR2 = fSSreg / (fSSreg + fSSresid);
+ pResMat->PutDouble(fR2, 0, 2);
+ } // end not exact fit
+ } //end bStats
+ PushMatrix(pResMat);
+ } //end nCase == 2
+ else // nCase == 3, Y is row, all matrices are transposed
+ {
+ ::std::vector< double> aVecR(N); // for QR decomposition
+ // Enough memory for needed matrices?
+ ScMatrixRef pMeans = GetNewMat(1, K); // mean of each row
+ ScMatrixRef pMatZ; // for Q' * Y , inter alia
+ if (bStats)
+ pMatZ = pMatY->Clone(); // Y is used in statistic, keep it
+ else
+ pMatZ = pMatY; // Y can be overwritten
+ ScMatrixRef pSlopes = GetNewMat(K,1); // from b1 to bK
+ if (!pMeans || !pMatZ || !pSlopes)
+ {
+ PushError(errStackOverflow);
+ return;
+ }
+ if (bConstant)
+ {
+ lcl_CalculateRowMeans(pMatX, pMeans, N, K);
+ lcl_CalculateRowsDelta(pMatX, pMeans, N, K);
+ }
+
+ if (!lcl_TCalculateQRdecomposition(pMatX, aVecR, K, N))
+ {
+ PushNoValue();
+ return;
+ }
+
+ // Later on we will divide by elements of aVecR, so make sure
+ // that they aren't zero.
+ bool bIsSingular=false;
+ for (SCSIZE row=0; row < K && !bIsSingular; row++)
+ bIsSingular = bIsSingular || aVecR[row]==0.0;
+ if (bIsSingular)
+ {
+ PushNoValue();
+ return;
+ }
+ // Z = Q' Y
+ for (SCSIZE row = 0; row < K; row++)
+ {
+ lcl_TApplyHouseholderTransformation(pMatX, row, pMatZ, N);
+ }
+ // B = R^(-1) * Q' * Y <=> B = R^(-1) * Z <=> R * B = Z
+ // result Z should have zeros for index>=K; if not, ignore values
+ for (SCSIZE col = 0; col < K ; col++)
+ {
+ pSlopes->PutDouble( pMatZ->GetDouble(col), col);
+ }
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pSlopes, K, true);
+ double fIntercept = 0.0;
+ if (bConstant)
+ fIntercept = fMeanY - lcl_GetSumProduct(pMeans,pSlopes,K);
+ // Fill first line in result matrix
+ pResMat->PutDouble(_bRKP ? exp(fIntercept) : fIntercept, K, 0 );
+ for (SCSIZE i = 0; i < K; i++)
+ pResMat->PutDouble(_bRKP ? exp(pSlopes->GetDouble(i))
+ : pSlopes->GetDouble(i) , K-1-i, 0);
+
+
+ if (bStats)
+ {
+ double fSSreg = 0.0;
+ double fSSresid = 0.0;
+ // re-use memory of Z;
+ pMatZ->FillDouble(0.0, 0, 0, N-1, 0);
+ // Z = R * Slopes
+ lcl_ApplyUpperRightTriangle(pMatX, aVecR, pSlopes, pMatZ, K, true);
+ // Z = Q * Z, that is Q * R * Slopes = X * Slopes
+ for (SCSIZE rowp1 = K; rowp1 > 0; rowp1--)
+ {
+ lcl_TApplyHouseholderTransformation(pMatX, rowp1-1, pMatZ,N);
+ }
+ fSSreg =lcl_GetSumProduct(pMatZ, pMatZ, N);
+ // re-use Y for residuals, Y = Y-Z
+ for (SCSIZE col = 0; col < N; col++)
+ pMatY->PutDouble(pMatY->GetDouble(col) - pMatZ->GetDouble(col), col);
+ fSSresid = lcl_GetSumProduct(pMatY, pMatY, N);
+ pResMat->PutDouble(fSSreg, 0, 4);
+ pResMat->PutDouble(fSSresid, 1, 4);
+
+ double fDegreesFreedom =static_cast<double>( (bConstant) ? N-K-1 : N-K );
+ pResMat->PutDouble(fDegreesFreedom, 1, 3);
+
+ if (fDegreesFreedom == 0.0 || fSSresid == 0.0 || fSSreg == 0.0)
+ { // exact fit; incl. case observed values Y are identical
+ pResMat->PutDouble(0.0, 1, 4); // SSresid
+ // F = (SSreg/K) / (SSresid/df) = #DIV/0!
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), 0, 3); // F
+ // RMSE = sqrt(SSresid / df) = sqrt(0 / df) = 0
+ pResMat->PutDouble(0.0, 1, 2); // RMSE
+ // SigmaSlope[i] = RMSE * sqrt(matrix[i,i]) = 0 * sqrt(...) = 0
+ for (SCSIZE i=0; i<K; i++)
+ pResMat->PutDouble(0.0, K-1-i, 1);
+
+ // SigmaIntercept = RMSE * sqrt(...) = 0
+ if (bConstant)
+ pResMat->PutDouble(0.0, K, 1); //SigmaIntercept
+ else
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
+
+ // R^2 = SSreg / (SSreg + SSresid) = 1.0
+ pResMat->PutDouble(1.0, 0, 2); // R^2
+ }
+ else
+ {
+ double fFstatistic = (fSSreg / static_cast<double>(K))
+ / (fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fFstatistic, 0, 3);
+
+ // standard error of estimate = root mean SSE
+ double fRMSE = sqrt(fSSresid / fDegreesFreedom);
+ pResMat->PutDouble(fRMSE, 1, 2);
+
+ // standard error of slopes
+ // = RMSE * sqrt(diagonal element of (R' R)^(-1) )
+ // standard error of intercept
+ // = RMSE * sqrt( Xmean * (R' R)^(-1) * Xmean' + 1/N)
+ // (R' R)^(-1) = R^(-1) * (R')^(-1). Do not calculate it as
+ // a whole matrix, but iterate over unit vectors.
+ double fSigmaSlope = 0.0;
+ double fSigmaIntercept = 0.0;
+ for (SCSIZE row = 0; row < K; row++)
+ {
+ //re-use memory of MatZ
+ pMatZ->FillDouble(0.0,0,0,K-1,0); // Z = unit vector e
+ pMatZ->PutDouble(1.0, row);
+ //Solve R' * Z = e
+ lcl_SolveWithLowerLeftTriangle(pMatX, aVecR, pMatZ, K, true);
+ // Solve R * Znew = Zold
+ lcl_SolveWithUpperRightTriangle(pMatX, aVecR, pMatZ, K, true);
+ // now Z is column col in (R' R)^(-1)
+ fSigmaSlope = fRMSE * sqrt(pMatZ->GetDouble(row));
+ pResMat->PutDouble(fSigmaSlope, K-1-row, 1);
+ // (R' R) ^(-1) is symmetric
+ if (bConstant)
+ fSigmaIntercept += lcl_GetSumProduct(pMeans, pMatZ, K);
+ }
+ if (bConstant)
+ {
+ fSigmaIntercept = fRMSE
+ * sqrt(fSigmaIntercept + 1.0 / static_cast<double>(N));
+ pResMat->PutDouble(fSigmaIntercept, K, 1);
+ }
+ else
+ {
+ pResMat->PutString(ScGlobal::GetRscString(STR_NV_STR), K, 1);
+ }
+
+ double fR2 = fSSreg / (fSSreg + fSSresid);
+ pResMat->PutDouble(fR2, 0, 2);
+ } // end not exact fit
+ } //end bStats
+ PushMatrix(pResMat);
+ } //end nCase == 3
+ } // end multiple regression
+ return;
+}
+
+
+// Changed CheckMatrix, which allways returns correct values for M and N
+bool ScInterpreter::CheckMatrix(BOOL _bLOG,BYTE& nCase,SCSIZE& nCX,SCSIZE& nCY,SCSIZE& nRX,SCSIZE&
nRY,SCSIZE& M,SCSIZE& N,ScMatrixRef& pMatX,ScMatrixRef& pMatY)
+{
+ nCX = 0;
+ nCY = 0;
+ nRX = 0;
+ nRY = 0;
+ M = 0;
+ N = 0;
+ pMatY->GetDimensions(nCY, nRY);
+ const SCSIZE nCountY = nCY * nRY;
+ for ( SCSIZE i = 0; i < nCountY; i++ )
+ {
+ if (!pMatY->IsValue(i))
+ {
+ PushIllegalArgument();
+ return false;
+ }
+ }
+
+
+ if ( _bLOG )
+ {
+ ScMatrixRef pNewY = pMatY->CloneIfConst();
+ for (SCSIZE nElem = 0; nElem < nCountY; nElem++)
+ {
+ const double fVal = pNewY->GetDouble(nElem);
+ if (fVal <= 0.0)
+ {
+ PushIllegalArgument();
+ return false;
+ }
+ else
+ pNewY->PutDouble(log(fVal), nElem);
+ }
+ pMatY = pNewY;
+ }
+
+
+ if (pMatX)
+ {
+ pMatX->GetDimensions(nCX, nRX);
+ const SCSIZE nCountX = nCX * nRX;
+ for ( SCSIZE i = 0; i < nCountX; i++ )
+ if (!pMatX->IsValue(i))
+ {
+ PushIllegalArgument();
+ return false;
+ }
+ if (nCX == nCY && nRX == nRY)
+ {
+ nCase = 1; // simple regression
+ M = 1;
+ N = nCountY;
+ }
+ else if (nCY != 1 && nRY != 1)
+ {
+ PushIllegalArgument();
+ return false;
+ }
+ else if (nCY == 1)
+ {
+ if (nRX != nRY)
+ {
+ PushIllegalArgument();
+ return false;
+ }
+ else
+ {
+ nCase = 2; // Y is column
+ N = nRY;
+ M = nCX;
+ }
+ }
+ else if (nCX != nCY)
+ {
+ PushIllegalArgument();
+ return false;
+ }
+ else
+ {
+ nCase = 3; // Y is row
+ N = nCY;
+ M = nRX;
+ }
+ }
+ else
+ {
+ pMatX = GetNewMat(nCY, nRY);
+ nCX = nCY;
+ nRX = nRY;
+ if (!pMatX)
+ {
+ PushIllegalArgument();
+ return false;
+ }
+ for ( SCSIZE i = 1; i <= nCountY; i++ )
+ pMatX->PutDouble((double)i, i-1);
+ nCase = 1;
+ N = nCountY;
+ M = 1;
+ }
+ return true;
+}
diff --git a/sc/source/core/tool/makefile.mk b/sc/source/core/tool/makefile.mk
index b1cad58..bd241b3 100644
--- a/sc/source/core/tool/makefile.mk
+++ b/sc/source/core/tool/makefile.mk
@@ -82,6 +82,7 @@ SLOFILES = \
$(SLO)$/interpr4.obj \
$(SLO)$/interpr5.obj \
$(SLO)$/interpr6.obj \
+ $(SLO)$/interpr7.obj \
$(SLO)$/lookupcache.obj \
$(SLO)$/navicfg.obj \
$(SLO)$/odffmap.obj \
Context
Privacy Policy |
Impressum (Legal Info) |
Copyright information: Unless otherwise specified, all text and images
on this website are licensed under the
Creative Commons Attribution-Share Alike 3.0 License.
This does not include the source code of LibreOffice, which is
licensed under the Mozilla Public License (
MPLv2).
"LibreOffice" and "The Document Foundation" are
registered trademarks of their corresponding registered owners or are
in actual use as trademarks in one or more countries. Their respective
logos and icons are also subject to international copyright laws. Use
thereof is explained in our
trademark policy.