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.me...@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 \
_______________________________________________ LibreOffice mailing list LibreOffice@lists.freedesktop.org http://lists.freedesktop.org/mailman/listinfo/libreoffice