From 2bbb2a83d4ba22d1a11aedeb6e1526269988715f Mon Sep 17 00:00:00 2001 From: xuri Date: Sun, 17 Apr 2022 23:49:37 +0800 Subject: [PATCH] ref #65, initial formula functions: GROWTH and TREND --- calc.go | 825 +++++++++++++++++++++++++++++++++++++++++++++++---- calc_test.go | 85 +++++- 2 files changed, 848 insertions(+), 62 deletions(-) diff --git a/calc.go b/calc.go index 03d467b851..f4ab9c4a72 100644 --- a/calc.go +++ b/calc.go @@ -462,6 +462,7 @@ type formulaFuncs struct { // GCD // GEOMEAN // GESTEP +// GROWTH // HARMEAN // HEX2BIN // HEX2DEC @@ -682,6 +683,7 @@ type formulaFuncs struct { // TINV // TODAY // TRANSPOSE +// TREND // TRIM // TRIMMEAN // TRUE @@ -960,7 +962,11 @@ func (f *File) evalInfixExpFunc(sheet, cell string, token, nextToken efp.Token, argsStack.Peek().(*list.List).PushBack(arg) } } else { - opdStack.Push(efp.Token{TValue: arg.Value(), TType: efp.TokenTypeOperand, TSubType: efp.TokenSubTypeNumber}) + val := arg.Value() + if arg.Type == ArgMatrix && len(arg.Matrix) > 0 && len(arg.Matrix[0]) > 0 { + val = arg.Matrix[0][0].Value() + } + opdStack.Push(efp.Token{TValue: val, TType: efp.TokenTypeOperand, TSubType: efp.TokenSubTypeNumber}) } return nil } @@ -2015,7 +2021,6 @@ func (fn *formulaFuncs) COMPLEX(argsList *list.List) formulaArg { // cmplx2str replace complex number string characters. func cmplx2str(num complex128, suffix string) string { - c := fmt.Sprint(num) realPart, imagPart := fmt.Sprint(real(num)), fmt.Sprint(imag(num)) isNum, i := isNumeric(realPart) if isNum && i > 15 { @@ -2025,7 +2030,7 @@ func cmplx2str(num complex128, suffix string) string { if isNum && i > 15 { imagPart = roundPrecision(imagPart, -1) } - c = realPart + c := realPart if imag(num) > 0 { c += "+" } @@ -4073,7 +4078,8 @@ func newFormulaArgMatrix(numMtx [][]float64) (arg [][]formulaArg) { for r, row := range numMtx { arg = append(arg, make([]formulaArg, len(row))) for c, cell := range row { - arg[r][c] = newNumberFormulaArg(cell) + decimal, _ := big.NewFloat(cell).SetPrec(15).Float64() + arg[r][c] = newNumberFormulaArg(decimal) } } return @@ -5819,12 +5825,12 @@ func (fn *formulaFuncs) BETAdotDIST(argsList *list.List) formulaArg { if a.Number == b.Number { return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) } - fScale := b.Number - a.Number - x.Number = (x.Number - a.Number) / fScale + scale := b.Number - a.Number + x.Number = (x.Number - a.Number) / scale if cumulative.Number == 1 { return newNumberFormulaArg(getBetaDist(x.Number, alpha.Number, beta.Number)) } - return newNumberFormulaArg(getBetaDistPDF(x.Number, alpha.Number, beta.Number) / fScale) + return newNumberFormulaArg(getBetaDistPDF(x.Number, alpha.Number, beta.Number) / scale) } // BETADIST function calculates the cumulative beta probability density @@ -6665,12 +6671,12 @@ func getLogGamma(fZ float64) float64 { // getLowRegIGamma returns lower regularized incomplete gamma function. func getLowRegIGamma(fA, fX float64) float64 { - fLnFactor := fA*math.Log(fX) - fX - getLogGamma(fA) - fFactor := math.Exp(fLnFactor) + lnFactor := fA*math.Log(fX) - fX - getLogGamma(fA) + factor := math.Exp(lnFactor) if fX > fA+1 { - return 1 - fFactor*getGammaContFraction(fA, fX) + return 1 - factor*getGammaContFraction(fA, fX) } - return fFactor * getGammaSeries(fA, fX) + return factor * getGammaSeries(fA, fX) } // getChiSqDistCDF returns left tail for the Chi-Square distribution. @@ -7610,6 +7616,703 @@ func (fn *formulaFuncs) GEOMEAN(argsList *list.List) formulaArg { return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) } +// getNewMatrix create matrix by given columns and rows. +func getNewMatrix(c, r int) (matrix [][]float64) { + for i := 0; i < c; i++ { + for j := 0; j < r; j++ { + for x := len(matrix); x <= i; x++ { + matrix = append(matrix, []float64{}) + } + for y := len(matrix[i]); y <= j; y++ { + matrix[i] = append(matrix[i], 0) + } + matrix[i][j] = 0 + } + } + return +} + +// approxSub subtract two values, if signs are identical and the values are +// equal, will be returns 0 instead of calculating the subtraction. +func approxSub(a, b float64) float64 { + if ((a < 0 && b < 0) || (a > 0 && b > 0)) && math.Abs(a-b) < 2.22045e-016 { + return 0 + } + return a - b +} + +// matrixClone return a copy of all elements of the original matrix. +func matrixClone(matrix [][]float64) (cloneMatrix [][]float64) { + for i := 0; i < len(matrix); i++ { + for j := 0; j < len(matrix[i]); j++ { + for x := len(cloneMatrix); x <= i; x++ { + cloneMatrix = append(cloneMatrix, []float64{}) + } + for k := len(cloneMatrix[i]); k <= j; k++ { + cloneMatrix[i] = append(cloneMatrix[i], 0) + } + cloneMatrix[i][j] = matrix[i][j] + } + } + return +} + +// trendGrowthMatrixInfo defined matrix checking result. +type trendGrowthMatrixInfo struct { + trendType, nCX, nCY, nRX, nRY, M, N int + mtxX, mtxY [][]float64 +} + +// prepareTrendGrowthMtxX is a part of implementation of the trend growth prepare. +func prepareTrendGrowthMtxX(mtxX [][]float64) [][]float64 { + var mtx [][]float64 + for i := 0; i < len(mtxX); i++ { + for j := 0; j < len(mtxX[i]); j++ { + if mtxX[i][j] == 0 { + return nil + } + for x := len(mtx); x <= j; x++ { + mtx = append(mtx, []float64{}) + } + for y := len(mtx[j]); y <= i; y++ { + mtx[j] = append(mtx[j], 0) + } + mtx[j][i] = mtxX[i][j] + } + } + return mtx +} + +// prepareTrendGrowthMtxY is a part of implementation of the trend growth prepare. +func prepareTrendGrowthMtxY(bLOG bool, mtxY [][]float64) [][]float64 { + var mtx [][]float64 + for i := 0; i < len(mtxY); i++ { + for j := 0; j < len(mtxY[i]); j++ { + if mtxY[i][j] == 0 { + return nil + } + for x := len(mtx); x <= j; x++ { + mtx = append(mtx, []float64{}) + } + for y := len(mtx[j]); y <= i; y++ { + mtx[j] = append(mtx[j], 0) + } + mtx[j][i] = mtxY[i][j] + } + } + if bLOG { + var pNewY [][]float64 + for i := 0; i < len(mtxY); i++ { + for j := 0; j < len(mtxY[i]); j++ { + fVal := mtxY[i][j] + if fVal <= 0 { + return nil + } + for x := len(pNewY); x <= j; x++ { + pNewY = append(pNewY, []float64{}) + } + for y := len(pNewY[j]); y <= i; y++ { + pNewY[j] = append(pNewY[j], 0) + } + pNewY[j][i] = math.Log(fVal) + } + } + mtx = pNewY + } + return mtx +} + +// prepareTrendGrowth check and return the result. +func prepareTrendGrowth(bLOG bool, mtxX, mtxY [][]float64) (*trendGrowthMatrixInfo, formulaArg) { + var nCX, nRX, M, N, trendType int + nRY, nCY := len(mtxY), len(mtxY[0]) + cntY := nCY * nRY + newY := prepareTrendGrowthMtxY(bLOG, mtxY) + if newY == nil { + return nil, newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) + } + var newX [][]float64 + if len(mtxX) != 0 { + nRX, nCX = len(mtxX), len(mtxX[0]) + if newX = prepareTrendGrowthMtxX(mtxX); newX == nil { + return nil, newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) + } + if nCX == nCY && nRX == nRY { + trendType, M, N = 1, 1, cntY // simple regression + } else if nCY != 1 && nRY != 1 { + return nil, newErrorFormulaArg(formulaErrorREF, formulaErrorREF) + } else if nCY == 1 { + if nRX != nRY { + return nil, newErrorFormulaArg(formulaErrorREF, formulaErrorREF) + } + trendType, M, N = 2, nCX, nRY + } else if nCX != nCY { + return nil, newErrorFormulaArg(formulaErrorREF, formulaErrorREF) + } else { + trendType, M, N = 3, nRX, nCY + } + } else { + newX = getNewMatrix(nCY, nRY) + nCX, nRX = nCY, nRY + num := 1.0 + for i := 0; i < nRY; i++ { + for j := 0; j < nCY; j++ { + newX[j][i] = num + num++ + } + } + trendType, M, N = 1, 1, cntY + } + return &trendGrowthMatrixInfo{ + trendType: trendType, + nCX: nCX, + nCY: nCY, + nRX: nRX, + nRY: nRY, + M: M, + N: N, + mtxX: newX, + mtxY: newY, + }, newEmptyFormulaArg() +} + +// calcPosition calculate position for matrix by given index. +func calcPosition(mtx [][]float64, idx int) (row, col int) { + rowSize := len(mtx[0]) + col = idx + if rowSize > 1 { + col = idx / rowSize + } + row = idx - col*rowSize + return +} + +// getDouble returns float64 data type value in the matrix by given index. +func getDouble(mtx [][]float64, idx int) float64 { + row, col := calcPosition(mtx, idx) + return mtx[col][row] +} + +// putDouble set a float64 data type value in the matrix by given index. +func putDouble(mtx [][]float64, idx int, val float64) { + row, col := calcPosition(mtx, idx) + mtx[col][row] = val +} + +// calcMeanOverAll returns mean of the given matrix by over all element. +func calcMeanOverAll(mtx [][]float64, n int) float64 { + var sum float64 + for i := 0; i < len(mtx); i++ { + for j := 0; j < len(mtx[i]); j++ { + sum += mtx[i][j] + } + } + return sum / float64(n) +} + +// calcSumProduct returns uses the matrices as vectors of length M over all +// element. +func calcSumProduct(mtxA, mtxB [][]float64, m int) float64 { + sum := 0.0 + for i := 0; i < m; i++ { + sum += getDouble(mtxA, i) * getDouble(mtxB, i) + } + return sum +} + +// calcColumnMeans calculates means of the columns of matrix. +func calcColumnMeans(mtxX, mtxRes [][]float64, c, r int) { + for i := 0; i < c; i++ { + var sum float64 + for k := 0; k < r; k++ { + sum += mtxX[i][k] + } + putDouble(mtxRes, i, sum/float64(r)) + } + return +} + +// calcColumnsDelta calculates subtract of the columns of matrix. +func calcColumnsDelta(mtx, columnMeans [][]float64, c, r int) { + for i := 0; i < c; i++ { + for k := 0; k < r; k++ { + mtx[i][k] = approxSub(mtx[i][k], getDouble(columnMeans, i)) + } + } +} + +// calcSign returns sign by given value, no mathematical signum, but used to +// switch between adding and subtracting. +func calcSign(val float64) float64 { + if val > 0 { + return 1 + } + return -1 +} + +// calcColsMaximumNorm is a special version for use within QR +// decomposition. Maximum norm of column index c starting in row index r; +// matrix A has count n rows. +func calcColsMaximumNorm(mtxA [][]float64, c, r, n int) float64 { + var norm float64 + for row := r; row < n; row++ { + if norm < math.Abs(mtxA[c][row]) { + norm = math.Abs(mtxA[c][row]) + } + } + return norm +} + +// calcFastMult returns multiply n x m matrix A with m x l matrix B to n x l matrix R. +func calcFastMult(mtxA, mtxB, mtxR [][]float64, n, m, l int) { + var sum float64 + for row := 0; row < n; row++ { + for col := 0; col < l; col++ { + sum = 0.0 + for k := 0; k < m; k++ { + sum += mtxA[k][row] * mtxB[col][k] + } + mtxR[col][row] = sum + } + } +} + +// calcRowsEuclideanNorm is a special version for use within QR +// decomposition. Euclidean norm of column index c starting in row index r; +// matrix a has count n rows. +func calcRowsEuclideanNorm(mtxA [][]float64, c, r, n int) float64 { + var norm float64 + for row := r; row < n; row++ { + norm += mtxA[c][row] * mtxA[c][row] + } + return math.Sqrt(norm) +} + +// calcRowsSumProduct is a special version for use within QR decomposition. +// starting in row index r; +// a and b are indices of columns, matrices A and B have count n rows. +func calcRowsSumProduct(mtxA [][]float64, a int, mtxB [][]float64, b, r, n int) float64 { + var result float64 + for row := r; row < n; row++ { + result += mtxA[a][row] * mtxB[b][row] + } + return result +} + +// calcSolveWithUpperRightTriangle solve for X in R*X=S using back substitution. +func calcSolveWithUpperRightTriangle(mtxA [][]float64, vecR []float64, mtxS [][]float64, k int, bIsTransposed bool) { + var row int + for rowp1 := k; rowp1 > 0; rowp1-- { + row = rowp1 - 1 + sum := getDouble(mtxS, row) + for col := rowp1; col < k; col++ { + if bIsTransposed { + sum -= mtxA[row][col] * getDouble(mtxS, col) + } else { + sum -= mtxA[col][row] * getDouble(mtxS, col) + } + } + putDouble(mtxS, row, sum/vecR[row]) + } +} + +// calcRowQRDecomposition calculates a QR decomposition with Householder +// reflection. +func calcRowQRDecomposition(mtxA [][]float64, vecR []float64, k, n int) bool { + for col := 0; col < k; col++ { + scale := calcColsMaximumNorm(mtxA, col, col, n) + if scale == 0 { + return false + } + for row := col; row < n; row++ { + mtxA[col][row] = mtxA[col][row] / scale + } + euclid := calcRowsEuclideanNorm(mtxA, col, col, n) + factor := 1.0 / euclid / (euclid + math.Abs(mtxA[col][col])) + signum := calcSign(mtxA[col][col]) + mtxA[col][col] = mtxA[col][col] + signum*euclid + vecR[col] = -signum * scale * euclid + // apply Householder transformation to A + for c := col + 1; c < k; c++ { + sum := calcRowsSumProduct(mtxA, col, mtxA, c, col, n) + for row := col; row < n; row++ { + mtxA[c][row] = mtxA[c][row] - sum*factor*mtxA[col][row] + } + } + } + return true +} + +// calcApplyColsHouseholderTransformation transposed matrices A and Y. +func calcApplyColsHouseholderTransformation(mtxA [][]float64, r int, mtxY [][]float64, n int) { + denominator := calcColsSumProduct(mtxA, r, mtxA, r, r, n) + numerator := calcColsSumProduct(mtxA, r, mtxY, 0, r, n) + factor := 2 * (numerator / denominator) + for col := r; col < n; col++ { + putDouble(mtxY, col, getDouble(mtxY, col)-factor*mtxA[col][r]) + } +} + +// calcRowMeans calculates means of the rows of matrix. +func calcRowMeans(mtxX, mtxRes [][]float64, c, r int) { + for k := 0; k < r; k++ { + var fSum float64 + for i := 0; i < c; i++ { + fSum += mtxX[i][k] + } + mtxRes[k][0] = fSum / float64(c) + } +} + +// calcRowsDelta calculates subtract of the rows of matrix. +func calcRowsDelta(mtx, rowMeans [][]float64, c, r int) { + for k := 0; k < r; k++ { + for i := 0; i < c; i++ { + mtx[i][k] = approxSub(mtx[i][k], rowMeans[k][0]) + } + } +} + +// calcColumnMaximumNorm returns maximum norm of row index R starting in col +// index C; matrix A has count N columns. +func calcColumnMaximumNorm(mtxA [][]float64, r, c, n int) float64 { + var norm float64 + for col := c; col < n; col++ { + if norm < math.Abs(mtxA[col][r]) { + norm = math.Abs(mtxA[col][r]) + } + } + return norm +} + +// calcColsEuclideanNorm returns euclidean norm of row index R starting in +// column index C; matrix A has count N columns. +func calcColsEuclideanNorm(mtxA [][]float64, r, c, n int) float64 { + var norm float64 + for col := c; col < n; col++ { + norm += (mtxA[col][r]) * (mtxA[col][r]) + } + return math.Sqrt(norm) +} + +// calcColsSumProduct returns sum product for given matrix. +func calcColsSumProduct(mtxA [][]float64, a int, mtxB [][]float64, b, c, n int) float64 { + var result float64 + for col := c; col < n; col++ { + result += mtxA[col][a] * mtxB[col][b] + } + return result +} + +// calcColQRDecomposition same with transposed matrix A, N is count of +// columns, k count of rows. +func calcColQRDecomposition(mtxA [][]float64, vecR []float64, k, n int) bool { + var sum float64 + for row := 0; row < k; row++ { + // calculate vector u of the householder transformation + scale := calcColumnMaximumNorm(mtxA, row, row, n) + if scale == 0 { + return false + } + for col := row; col < n; col++ { + mtxA[col][row] = mtxA[col][row] / scale + } + euclid := calcColsEuclideanNorm(mtxA, row, row, n) + factor := 1 / euclid / (euclid + math.Abs(mtxA[row][row])) + signum := calcSign(mtxA[row][row]) + mtxA[row][row] = mtxA[row][row] + signum*euclid + vecR[row] = -signum * scale * euclid + // apply Householder transformation to A + for r := row + 1; r < k; r++ { + sum = calcColsSumProduct(mtxA, row, mtxA, r, row, n) + for col := row; col < n; col++ { + mtxA[col][r] = mtxA[col][r] - sum*factor*mtxA[col][row] + } + } + } + return true +} + +// calcApplyRowsHouseholderTransformation applies a Householder transformation to a +// column vector Y with is given as Nx1 Matrix. The vector u, from which the +// Householder transformation is built, 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 calcRowQRDecomposition. +func calcApplyRowsHouseholderTransformation(mtxA [][]float64, c int, mtxY [][]float64, n int) { + denominator := calcRowsSumProduct(mtxA, c, mtxA, c, c, n) + numerator := calcRowsSumProduct(mtxA, c, mtxY, 0, c, n) + factor := 2 * (numerator / denominator) + for row := c; row < n; row++ { + putDouble(mtxY, row, getDouble(mtxY, row)-factor*mtxA[c][row]) + } +} + +// calcTrendGrowthSimpleRegression calculate simple regression for the calcTrendGrowth. +func calcTrendGrowthSimpleRegression(bConstant, bGrowth bool, mtxY, mtxX, newX, mtxRes [][]float64, meanY float64, N int) { + var meanX float64 + if bConstant { + meanX = calcMeanOverAll(mtxX, N) + for i := 0; i < len(mtxX); i++ { + for j := 0; j < len(mtxX[i]); j++ { + mtxX[i][j] = approxSub(mtxX[i][j], meanX) + } + } + } + sumXY := calcSumProduct(mtxX, mtxY, N) + sumX2 := calcSumProduct(mtxX, mtxX, N) + slope := sumXY / sumX2 + var help float64 + var intercept float64 + if bConstant { + intercept = meanY - slope*meanX + for i := 0; i < len(mtxRes); i++ { + for j := 0; j < len(mtxRes[i]); j++ { + help = newX[i][j]*slope + intercept + if bGrowth { + mtxRes[i][j] = math.Exp(help) + } else { + mtxRes[i][j] = help + } + } + } + } else { + for i := 0; i < len(mtxRes); i++ { + for j := 0; j < len(mtxRes[i]); j++ { + help = newX[i][j] * slope + if bGrowth { + mtxRes[i][j] = math.Exp(help) + } else { + mtxRes[i][j] = help + } + } + } + } +} + +// calcTrendGrowthMultipleRegressionPart1 calculate multiple regression for the +// calcTrendGrowth. +func calcTrendGrowthMultipleRegressionPart1(bConstant, bGrowth bool, mtxY, mtxX, newX, mtxRes [][]float64, meanY float64, RXN, K, N int) { + vecR := make([]float64, N) // for QR decomposition + means := getNewMatrix(K, 1) // mean of each column + slopes := getNewMatrix(1, K) // from b1 to bK + if len(means) == 0 || len(slopes) == 0 { + return + } + if bConstant { + calcColumnMeans(mtxX, means, K, N) + calcColumnsDelta(mtxX, means, K, N) + } + if !calcRowQRDecomposition(mtxX, vecR, K, N) { + return + } + // Later on we will divide by elements of vecR, so make sure that they aren't zero. + bIsSingular := false + for row := 0; row < K && !bIsSingular; row++ { + bIsSingular = bIsSingular || vecR[row] == 0 + } + if bIsSingular { + return + } + for col := 0; col < K; col++ { + calcApplyRowsHouseholderTransformation(mtxX, col, mtxY, N) + } + for col := 0; col < K; col++ { + putDouble(slopes, col, getDouble(mtxY, col)) + } + calcSolveWithUpperRightTriangle(mtxX, vecR, slopes, K, false) + // Fill result matrix + calcFastMult(newX, slopes, mtxRes, RXN, K, 1) + if bConstant { + intercept := meanY - calcSumProduct(means, slopes, K) + for row := 0; row < RXN; row++ { + mtxRes[0][row] = mtxRes[0][row] + intercept + } + } + if bGrowth { + for i := 0; i < RXN; i++ { + putDouble(mtxRes, i, math.Exp(getDouble(mtxRes, i))) + } + } +} + +// calcTrendGrowthMultipleRegressionPart2 calculate multiple regression for the +// calcTrendGrowth. +func calcTrendGrowthMultipleRegressionPart2(bConstant, bGrowth bool, mtxY, mtxX, newX, mtxRes [][]float64, meanY float64, nCXN, K, N int) { + vecR := make([]float64, N) // for QR decomposition + means := getNewMatrix(K, 1) // mean of each row + slopes := getNewMatrix(K, 1) // row from b1 to bK + if len(means) == 0 || len(slopes) == 0 { + return + } + if bConstant { + calcRowMeans(mtxX, means, N, K) + calcRowsDelta(mtxX, means, N, K) + } + if !calcColQRDecomposition(mtxX, vecR, K, N) { + return + } + // later on we will divide by elements of vecR, so make sure that they aren't zero + bIsSingular := false + for row := 0; row < K && !bIsSingular; row++ { + bIsSingular = bIsSingular || vecR[row] == 0 + } + if bIsSingular { + return + } + for row := 0; row < K; row++ { + calcApplyColsHouseholderTransformation(mtxX, row, mtxY, N) + } + for col := 0; col < K; col++ { + putDouble(slopes, col, getDouble(mtxY, col)) + } + calcSolveWithUpperRightTriangle(mtxX, vecR, slopes, K, true) + // fill result matrix + calcFastMult(slopes, newX, mtxRes, 1, K, nCXN) + if bConstant { + fIntercept := meanY - calcSumProduct(means, slopes, K) + for col := 0; col < nCXN; col++ { + mtxRes[col][0] = mtxRes[col][0] + fIntercept + } + } + if bGrowth { + for i := 0; i < nCXN; i++ { + putDouble(mtxRes, i, math.Exp(getDouble(mtxRes, i))) + } + } +} + +// calcTrendGrowthRegression is a part of implementation of the calcTrendGrowth. +func calcTrendGrowthRegression(bConstant, bGrowth bool, trendType, nCXN, nRXN, K, N int, mtxY, mtxX, newX, mtxRes [][]float64) { + if len(mtxRes) == 0 { + return + } + var meanY float64 + if bConstant { + copyX, copyY := matrixClone(mtxX), matrixClone(mtxY) + mtxX, mtxY = copyX, copyY + meanY = calcMeanOverAll(mtxY, N) + for i := 0; i < len(mtxY); i++ { + for j := 0; j < len(mtxY[i]); j++ { + mtxY[i][j] = approxSub(mtxY[i][j], meanY) + } + } + } + switch trendType { + case 1: + calcTrendGrowthSimpleRegression(bConstant, bGrowth, mtxY, mtxX, newX, mtxRes, meanY, N) + break + case 2: + calcTrendGrowthMultipleRegressionPart1(bConstant, bGrowth, mtxY, mtxX, newX, mtxRes, meanY, nRXN, K, N) + break + default: + calcTrendGrowthMultipleRegressionPart2(bConstant, bGrowth, mtxY, mtxX, newX, mtxRes, meanY, nCXN, K, N) + } +} + +// calcTrendGrowth returns values along a predicted exponential trend. +func calcTrendGrowth(mtxY, mtxX, newX [][]float64, bConstant, bGrowth bool) ([][]float64, formulaArg) { + getMatrixParams, errArg := prepareTrendGrowth(bGrowth, mtxX, mtxY) + if errArg.Type != ArgEmpty { + return nil, errArg + } + trendType := getMatrixParams.trendType + nCX := getMatrixParams.nCX + nRX := getMatrixParams.nRX + K := getMatrixParams.M + N := getMatrixParams.N + mtxX = getMatrixParams.mtxX + mtxY = getMatrixParams.mtxY + // checking if data samples are enough + if (bConstant && (N < K+1)) || (!bConstant && (N < K)) || (N < 1) || (K < 1) { + return nil, errArg + } + // set the default newX if necessary + nCXN, nRXN := nCX, nRX + if len(newX) == 0 { + newX = matrixClone(mtxX) // mtxX will be changed to X-meanX + } else { + nRXN, nCXN = len(newX[0]), len(newX) + if (trendType == 2 && K != nCXN) || (trendType == 3 && K != nRXN) { + return nil, errArg + } + } + var mtxRes [][]float64 + switch trendType { + case 1: + mtxRes = getNewMatrix(nCXN, nRXN) + break + case 2: + mtxRes = getNewMatrix(1, nRXN) + break + default: + mtxRes = getNewMatrix(nCXN, 1) + } + calcTrendGrowthRegression(bConstant, bGrowth, trendType, nCXN, nRXN, K, N, mtxY, mtxX, newX, mtxRes) + return mtxRes, errArg +} + +// trendGrowth is an implementation of the formula functions GROWTH and TREND. +func (fn *formulaFuncs) trendGrowth(name string, argsList *list.List) formulaArg { + if argsList.Len() < 1 { + return newErrorFormulaArg(formulaErrorVALUE, fmt.Sprintf("%s requires at least 1 argument", name)) + } + if argsList.Len() > 4 { + return newErrorFormulaArg(formulaErrorVALUE, fmt.Sprintf("%s allows at most 4 arguments", name)) + } + var knowY, knowX, newX [][]float64 + var errArg formulaArg + constArg := newBoolFormulaArg(true) + knowY, errArg = newNumberMatrix(argsList.Front().Value.(formulaArg), false) + if errArg.Type == ArgError { + return errArg + } + if argsList.Len() > 1 { + knowX, errArg = newNumberMatrix(argsList.Front().Next().Value.(formulaArg), false) + if errArg.Type == ArgError { + return errArg + } + } + if argsList.Len() > 2 { + newX, errArg = newNumberMatrix(argsList.Front().Next().Next().Value.(formulaArg), false) + if errArg.Type == ArgError { + return errArg + } + } + if argsList.Len() > 3 { + if constArg = argsList.Back().Value.(formulaArg).ToBool(); constArg.Type != ArgNumber { + return constArg + } + } + var mtxNewX [][]float64 + for i := 0; i < len(newX); i++ { + for j := 0; j < len(newX[i]); j++ { + for x := len(mtxNewX); x <= j; x++ { + mtxNewX = append(mtxNewX, []float64{}) + } + for k := len(mtxNewX[j]); k <= i; k++ { + mtxNewX[j] = append(mtxNewX[j], 0) + } + mtxNewX[j][i] = newX[i][j] + } + } + mtx, errArg := calcTrendGrowth(knowY, knowX, mtxNewX, constArg.Number == 1, name == "GROWTH") + if errArg.Type != ArgEmpty { + return errArg + } + return newMatrixFormulaArg(newFormulaArgMatrix(mtx)) +} + +// GROWTH function calculates the exponential growth curve through a given set +// of y-values and (optionally), one or more sets of x-values. The function +// then extends the curve to calculate additional y-values for a further +// supplied set of new x-values. The syntax of the function is: +// +// GROWTH(known_y's,[known_x's],[new_x's],[const]) +// +func (fn *formulaFuncs) GROWTH(argsList *list.List) formulaArg { + return fn.trendGrowth("GROWTH", argsList) +} + // HARMEAN function calculates the harmonic mean of a supplied set of values. // The syntax of the function is: // @@ -9652,92 +10355,98 @@ func (fn *formulaFuncs) TINV(argsList *list.List) formulaArg { return fn.TdotINVdot2T(argsList) } +// TREND function calculates the linear trend line through a given set of +// y-values and (optionally), a given set of x-values. The function then +// extends the linear trendline to calculate additional y-values for a further +// supplied set of new x-values. The syntax of the function is: +// +// TREND(known_y's,[known_x's],[new_x's],[const]) +// +func (fn *formulaFuncs) TREND(argsList *list.List) formulaArg { + return fn.trendGrowth("TREND", argsList) +} + // tTest calculates the probability associated with the Student's T Test. -func tTest(bTemplin bool, pMat1, pMat2 [][]formulaArg, nC1, nC2, nR1, nR2 int, fT, fF float64) (float64, float64, bool) { - var fCount1, fCount2, fSum1, fSumSqr1, fSum2, fSumSqr2 float64 +func tTest(bTemplin bool, mtx1, mtx2 [][]formulaArg, c1, c2, r1, r2 int, fT, fF float64) (float64, float64, bool) { + var cnt1, cnt2, sum1, sumSqr1, sum2, sumSqr2 float64 var fVal formulaArg - for i := 0; i < nC1; i++ { - for j := 0; j < nR1; j++ { - fVal = pMat1[i][j].ToNumber() + for i := 0; i < c1; i++ { + for j := 0; j < r1; j++ { + fVal = mtx1[i][j].ToNumber() if fVal.Type == ArgNumber { - fSum1 += fVal.Number - fSumSqr1 += fVal.Number * fVal.Number - fCount1++ + sum1 += fVal.Number + sumSqr1 += fVal.Number * fVal.Number + cnt1++ } } } - for i := 0; i < nC2; i++ { - for j := 0; j < nR2; j++ { - fVal = pMat2[i][j].ToNumber() + for i := 0; i < c2; i++ { + for j := 0; j < r2; j++ { + fVal = mtx2[i][j].ToNumber() if fVal.Type == ArgNumber { - fSum2 += fVal.Number - fSumSqr2 += fVal.Number * fVal.Number - fCount2++ + sum2 += fVal.Number + sumSqr2 += fVal.Number * fVal.Number + cnt2++ } } } - if fCount1 < 2.0 || fCount2 < 2.0 { + if cnt1 < 2.0 || cnt2 < 2.0 { return 0, 0, false } if bTemplin { - fS1 := (fSumSqr1 - fSum1*fSum1/fCount1) / (fCount1 - 1) / fCount1 - fS2 := (fSumSqr2 - fSum2*fSum2/fCount2) / (fCount2 - 1) / fCount2 + fS1 := (sumSqr1 - sum1*sum1/cnt1) / (cnt1 - 1) / cnt1 + fS2 := (sumSqr2 - sum2*sum2/cnt2) / (cnt2 - 1) / cnt2 if fS1+fS2 == 0 { return 0, 0, false } c := fS1 / (fS1 + fS2) - fT = math.Abs(fSum1/fCount1-fSum2/fCount2) / math.Sqrt(fS1+fS2) - fF = 1 / (c*c/(fCount1-1) + (1-c)*(1-c)/(fCount2-1)) + fT = math.Abs(sum1/cnt1-sum2/cnt2) / math.Sqrt(fS1+fS2) + fF = 1 / (c*c/(cnt1-1) + (1-c)*(1-c)/(cnt2-1)) return fT, fF, true } - fS1 := (fSumSqr1 - fSum1*fSum1/fCount1) / (fCount1 - 1) - fS2 := (fSumSqr2 - fSum2*fSum2/fCount2) / (fCount2 - 1) - fT = math.Abs(fSum1/fCount1-fSum2/fCount2) / math.Sqrt((fCount1-1)*fS1+(fCount2-1)*fS2) * math.Sqrt(fCount1*fCount2*(fCount1+fCount2-2)/(fCount1+fCount2)) - fF = fCount1 + fCount2 - 2 + fS1 := (sumSqr1 - sum1*sum1/cnt1) / (cnt1 - 1) + fS2 := (sumSqr2 - sum2*sum2/cnt2) / (cnt2 - 1) + fT = math.Abs(sum1/cnt1-sum2/cnt2) / math.Sqrt((cnt1-1)*fS1+(cnt2-1)*fS2) * math.Sqrt(cnt1*cnt2*(cnt1+cnt2-2)/(cnt1+cnt2)) + fF = cnt1 + cnt2 - 2 return fT, fF, true } // tTest is an implementation of the formula function TTEST. -func (fn *formulaFuncs) tTest(pMat1, pMat2 [][]formulaArg, fTails, fTyp float64) formulaArg { +func (fn *formulaFuncs) tTest(mtx1, mtx2 [][]formulaArg, fTails, fTyp float64) formulaArg { var fT, fF float64 - nC1 := len(pMat1) - nC2 := len(pMat2) - nR1 := len(pMat1[0]) - nR2 := len(pMat2[0]) - ok := true + c1, c2, r1, r2, ok := len(mtx1), len(mtx2), len(mtx1[0]), len(mtx2[0]), true if fTyp == 1 { - if nC1 != nC2 || nR1 != nR2 { + if c1 != c2 || r1 != r2 { return newErrorFormulaArg(formulaErrorNA, formulaErrorNA) } - var fCount, fSum1, fSum2, fSumSqrD float64 + var cnt, sum1, sum2, sumSqrD float64 var fVal1, fVal2 formulaArg - for i := 0; i < nC1; i++ { - for j := 0; j < nR1; j++ { - fVal1 = pMat1[i][j].ToNumber() - fVal2 = pMat2[i][j].ToNumber() + for i := 0; i < c1; i++ { + for j := 0; j < r1; j++ { + fVal1, fVal2 = mtx1[i][j].ToNumber(), mtx2[i][j].ToNumber() if fVal1.Type != ArgNumber || fVal2.Type != ArgNumber { continue } - fSum1 += fVal1.Number - fSum2 += fVal2.Number - fSumSqrD += (fVal1.Number - fVal2.Number) * (fVal1.Number - fVal2.Number) - fCount++ + sum1 += fVal1.Number + sum2 += fVal2.Number + sumSqrD += (fVal1.Number - fVal2.Number) * (fVal1.Number - fVal2.Number) + cnt++ } } - if fCount < 1 { + if cnt < 1 { return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) } - fSumD := fSum1 - fSum2 - fDivider := fCount*fSumSqrD - fSumD*fSumD - if fDivider == 0 { + sumD := sum1 - sum2 + divider := cnt*sumSqrD - sumD*sumD + if divider == 0 { return newErrorFormulaArg(formulaErrorDIV, formulaErrorDIV) } - fT = math.Abs(fSumD) * math.Sqrt((fCount-1)/fDivider) - fF = fCount - 1 + fT = math.Abs(sumD) * math.Sqrt((cnt-1)/divider) + fF = cnt - 1 } else if fTyp == 2 { - fT, fF, ok = tTest(false, pMat1, pMat2, nC1, nC2, nR1, nR2, fT, fF) + fT, fF, ok = tTest(false, mtx1, mtx2, c1, c2, r1, r2, fT, fF) } else { - fT, fF, ok = tTest(true, pMat1, pMat2, nC1, nC2, nR1, nR2, fT, fF) + fT, fF, ok = tTest(true, mtx1, mtx2, c1, c2, r1, r2, fT, fF) } if !ok { return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) diff --git a/calc_test.go b/calc_test.go index 6cae4a3ff7..c09174734d 100644 --- a/calc_test.go +++ b/calc_test.go @@ -572,9 +572,9 @@ func TestCalcCellValue(t *testing.T) { "=IMPRODUCT(COMPLEX(5,2),COMPLEX(0,1))": "-2+5i", "=IMPRODUCT(A1:C1)": "4", // MINVERSE - "=MINVERSE(A1:B2)": "", + "=MINVERSE(A1:B2)": "-0", // MMULT - "=MMULT(A4:A4,A4:A4)": "", + "=MMULT(A4:A4,A4:A4)": "0", // MOD "=MOD(6,4)": "2", "=MOD(6,3)": "0", @@ -597,7 +597,7 @@ func TestCalcCellValue(t *testing.T) { `=MULTINOMIAL("",3,1,2,5)`: "27720", "=MULTINOMIAL(MULTINOMIAL(1))": "1", // _xlfn.MUNIT - "=_xlfn.MUNIT(4)": "", + "=_xlfn.MUNIT(4)": "1", // ODD "=ODD(22)": "23", "=ODD(1.22)": "3", @@ -4444,6 +4444,83 @@ func TestCalcFORMULATEXT(t *testing.T) { } } +func TestCalcGROWTHandTREND(t *testing.T) { + cellData := [][]interface{}{ + {"known_x's", "known_y's", 0, -1}, + {1, 10, 1}, + {2, 20, 1}, + {3, 40}, + {4, 80}, + {}, + {"new_x's", "new_y's"}, + {5}, + {6}, + {7}, + } + f := prepareCalcData(cellData) + formulaList := map[string]string{ + "=GROWTH(A2:B2)": "1", + "=GROWTH(B2:B5,A2:A5,A8:A10)": "160", + "=GROWTH(B2:B5,A2:A5,A8:A10,FALSE)": "467.84375", + "=GROWTH(A4:A5,A2:B3,A8:A10,FALSE)": "", + "=GROWTH(A3:A5,A2:B4,A2:B3)": "2", + "=GROWTH(A4:A5,A2:B3)": "", + "=GROWTH(A2:B2,A2:B3)": "", + "=GROWTH(A2:B2,A2:B3,A2:B3,FALSE)": "1.28399658203125", + "=GROWTH(A2:B2,A4:B5,A4:B5,FALSE)": "1", + "=GROWTH(A3:C3,A2:C3,A2:B3)": "2", + "=TREND(A2:B2)": "1", + "=TREND(B2:B5,A2:A5,A8:A10)": "95", + "=TREND(B2:B5,A2:A5,A8:A10,FALSE)": "81.66796875", + "=TREND(A4:A5,A2:B3,A8:A10,FALSE)": "", + "=TREND(A4:A5,A2:B3,A2:B3,FALSE)": "1.5", + "=TREND(A3:A5,A2:B4,A2:B3)": "2", + "=TREND(A4:A5,A2:B3)": "", + "=TREND(A2:B2,A2:B3)": "", + "=TREND(A2:B2,A2:B3,A2:B3,FALSE)": "1", + "=TREND(A2:B2,A4:B5,A4:B5,FALSE)": "1", + "=TREND(A3:C3,A2:C3,A2:B3)": "2", + } + for formula, expected := range formulaList { + assert.NoError(t, f.SetCellFormula("Sheet1", "C1", formula)) + result, err := f.CalcCellValue("Sheet1", "C1") + assert.NoError(t, err, formula) + assert.Equal(t, expected, result, formula) + } + calcError := map[string]string{ + "=GROWTH()": "GROWTH requires at least 1 argument", + "=GROWTH(B2:B5,A2:A5,A8:A10,TRUE,0)": "GROWTH allows at most 4 arguments", + "=GROWTH(A1:B1,A2:A5,A8:A10,TRUE)": "strconv.ParseFloat: parsing \"known_x's\": invalid syntax", + "=GROWTH(B2:B5,A1:B1,A8:A10,TRUE)": "strconv.ParseFloat: parsing \"known_x's\": invalid syntax", + "=GROWTH(B2:B5,A2:A5,A1:B1,TRUE)": "strconv.ParseFloat: parsing \"known_x's\": invalid syntax", + "=GROWTH(B2:B5,A2:A5,A8:A10,\"\")": "strconv.ParseBool: parsing \"\": invalid syntax", + "=GROWTH(A2:B3,A4:B4)": "#REF!", + "=GROWTH(A4:B4,A2:A2)": "#REF!", + "=GROWTH(A2:A2,A4:A5)": "#REF!", + "=GROWTH(C1:C1,A2:A3)": "#NUM!", + "=GROWTH(D1:D1,A2:A3)": "#NUM!", + "=GROWTH(A2:A3,C1:C1)": "#NUM!", + "=TREND()": "TREND requires at least 1 argument", + "=TREND(B2:B5,A2:A5,A8:A10,TRUE,0)": "TREND allows at most 4 arguments", + "=TREND(A1:B1,A2:A5,A8:A10,TRUE)": "strconv.ParseFloat: parsing \"known_x's\": invalid syntax", + "=TREND(B2:B5,A1:B1,A8:A10,TRUE)": "strconv.ParseFloat: parsing \"known_x's\": invalid syntax", + "=TREND(B2:B5,A2:A5,A1:B1,TRUE)": "strconv.ParseFloat: parsing \"known_x's\": invalid syntax", + "=TREND(B2:B5,A2:A5,A8:A10,\"\")": "strconv.ParseBool: parsing \"\": invalid syntax", + "=TREND(A2:B3,A4:B4)": "#REF!", + "=TREND(A4:B4,A2:A2)": "#REF!", + "=TREND(A2:A2,A4:A5)": "#REF!", + "=TREND(C1:C1,A2:A3)": "#NUM!", + "=TREND(D1:D1,A2:A3)": "#REF!", + "=TREND(A2:A3,C1:C1)": "#NUM!", + } + for formula, expected := range calcError { + assert.NoError(t, f.SetCellFormula("Sheet1", "C1", formula)) + result, err := f.CalcCellValue("Sheet1", "C1") + assert.EqualError(t, err, expected, formula) + assert.Equal(t, "", result, formula) + } +} + func TestCalcHLOOKUP(t *testing.T) { cellData := [][]interface{}{ {"Example Result Table"}, @@ -4953,7 +5030,7 @@ func TestCalcMODE(t *testing.T) { formulaList := map[string]string{ "=MODE(A1:A10)": "3", "=MODE(B1:B6)": "2", - "=MODE.MULT(A1:A10)": "", + "=MODE.MULT(A1:A10)": "3", "=MODE.SNGL(A1:A10)": "3", "=MODE.SNGL(B1:B6)": "2", }