最小二乘法C++封装类

// FitFormula.cpp: implementation of the CFitFormula class. // ////////////////////////////////////////////////////////////////////// #include "stdafx.h" #include "MCTool.h" #include "FitFormula.h" #ifdef _DEBUG #undef THIS_FILE static char THIS_FILE[]=__FILE__; #define new DEBUG_NEW #endif ////////////////////////////////////////////////////////////////////// // Construction/Destruction ////////////////////////////////////////////////////////////////////// CFitFormula::CFitFormula() { m_dbMaxVal = 0.0f; m_dbMinVal = 0.0f; } CFitFormula::~CFitFormula() { } int CFitFormula::GetFormula(CString strPath, int nPower, CString &strFormula) { int nRetCode = 0; CStdioFile file; int nPtCount = 0; CString strBuffer; CStringArray arrPt; CString strTemp; int i = 0, j = 0; CDoubleArray X, Y, A; if (file.Open(strPath, CFile::modeRead | CFile::typeText)) { while (file.ReadString(strBuffer) != FALSE) { if (strBuffer != "") { arrPt.Add(strBuffer); nPtCount ++; } } file.Close(); } else { return -1; } int nSize = arrPt.GetSize(); if (nSize > 0) { j = arrPt[0].FindOneOf(","); if (j == -1) { return -1; } else { m_dbMaxVal = int(atof(arrPt[0].Mid(j + 1)) * 10000.0f + 0.5f) / 10000.0f; m_dbMinVal = m_dbMaxVal; } } else { return -1; } for (i = 0; i < nSize; i ++) { arrPt[i].Remove(' '); j = arrPt[i].Find(",", 0); double d1; __int64 n1 = 0; d1 = 0.0f; n1 = __int64(atof(arrPt[i].Left(j)) * 10000.0f + 0.5f); d1 = n1 / 10000.0f; Y.SetAtGrow(i, d1); d1 = 0.0f; n1 = 0; n1 = __int64(atof(arrPt[i].Mid(j + 1)) * 10000.0f + 0.5f); d1 = n1 / 10000.0f; X.SetAtGrow(i, d1); if (X[i] > m_dbMaxVal) { m_dbMaxVal = X[i]; } if (X[i] < m_dbMinVal) { m_dbMinVal = X[i]; } } //------------------------------ A.SetSize(nPower); GetMINNHvalue(X.GetData(), Y.GetData(), nPtCount, A.GetData(), nPower); //CalculateCurveParameter(&X, &Y, nPower, nPtCount, &A); strFormula = ""; for (i = nPower - 1; i > -1; i --) { strTemp.Format("(%e)", A[i]); if (strTemp.Replace("e", "*10^(") != 0) { strTemp += ")"; } strFormula += strTemp; strTemp.Format("*X^%d", i); strFormula += strTemp; if (i != 0) { strFormula += "+"; } } return nRetCode; } double CFitFormula::GetHLValue(double **n, long wide) { double nt[256]; long z1, z2; double z3, z4, z5, z6; long double d1 = 0.0f, d2 = 0.0f, d3 = 0.0f; if (wide > 2) { z5 = 0.0f; for (z1 = 0; z1 < wide; z1 ++) { z4 = n[z1][wide]; if (z4 == 0.0f) { continue; } for (z2 = 0; z2 <= wide; z2 ++) { nt[z2] = n[z1][z2]; n[z1][z2] = n[wide][z2]; } z3 = GetHLValue(n, wide - 1); for (z2 = 0; z2 <= wide; z2 ++) { n[z1][z2] = nt[z2]; } if (wide % 2 != 0) { z5 += z3 * z4; } else { z5 -= z3 * z4; } z3 = GetHLValue(n, wide - 1); if (wide % 2 == 0) { z6 = z5 + n[wide][wide] * z3; } else { z6 = z5 - n[wide][wide] * z3; } } return z6; } else if (wide == 2) { z6 = 0.0f; for (z2 = 0; z2 <= wide; z2 ++) { z3 = 1.0f; z4 = 1.0f; for (z1 = 0; z1 <= wide; z1 ++) { z3 = z3 * n[(z2 + z1) % (wide + 1)][z1]; z4 = z4 * n[(z2 - z1 + wide + 1) % (wide + 1)][z1]; } z6 += z3 - z4; } return z6; } else if (wide == 1) { d1 = n[0][0] * n[1][1]; d2 = n[0][1] * n[1][0]; d3 = d1 - d2; return double(d3); } else if (wide == 0) { return n[0][0]; } return -2; } long CFitFormula::GetZZXCvalue(long N1, long N2) { int i; N1 = abs(N1); N2 = abs(N2); i = min(N1, N2); if (i == 0 || N1 + N2 == 0) { return 1; } while (TRUE) { if (N1 < N2) { N2 %= N1; if (N2 == 0) { return N1; } } else { N1 %= N2; if (N1 == 0) { return N2; } } } } BOOL CFitFormula::ExplainEquation(double **fc, double *proot, long power) { long z1, z2; proot[power + 1] = GetHLValue(fc, power); for (z1 = 0; z1 <= power; z1 ++) { for (z2 = 0; z2 <= power; z2 ++) { fc[power + 2][z2] = fc[z1][z2]; fc[z1][z2] = fc[power + 1][z2]; } proot[z1] = GetHLValue(fc, power); for (z2 = 0; z2 <= power; z2 ++) { fc[z1][z2] = fc[power + 2][z2]; } } for (z1 = 0; z1 <= power + 1; z1 ++) { if (proot[z1] != int(proot[z1])) { return FALSE; } } z2 = GetZZXCvalue(long(proot[0]), long(proot[1])); for (z1 = 0; z1 <= power; z1 ++) { z2 = GetZZXCvalue(z2, long(proot[z1 + 1])); } if (proot[power + 1] < 0) { z2 = -z2; } for (z1 = 0; z1 <= power + 1; z1 ++) { proot[z1] = proot[z1] / z2; } return TRUE; } long CFitFormula::GetMINNHvalue(double *x, double *y, long data_len, double *hroot, long max_power/*3 == x^2+x^1+x^0*/) { double d1, d2, d3, d4, d5, d6; __int64 n4, n5, n6; long z1, z2, z3; long zfx = max(max_power - 1, 0), zfy = data_len; double *pdata = new double[data_len]; double *proot = new double[max_power]; double **sum = NULL, **fc = NULL; CMCToolApp::MKCreatArray(sizeof(double), max_power, data_len, reinterpret_cast<void ***>(&sum)); CMCToolApp::MKCreatArray(sizeof(double), max_power + 1, max_power, reinterpret_cast<void ***>(&fc)); memcpy(pdata, x, sizeof(double) * zfy); memcpy(sum[zfx], y, sizeof(double) * data_len); z3 = 1; for (z2 = 0; z2 < zfx; z2 ++) { for (z1 = 0; z1 < data_len; z1 ++) { d1 = pdata[z1]; sum[z2][z1] = pow(d1, z3); } z3 ++; } for (z2 = 0; z2 < max_power; z2 ++) { for (z1 = 0; z1 < max_power; z1 ++) { d1 = 0.0f; d2 = 0.0f; d3 = 0.0f; for (z3 = 0; z3 < data_len; z3 ++) { d1 += sum[z1][z3] * sum[z2][z3]; d2 += sum[z1][z3]; d3 += sum[z2][z3]; } d4 = 0.0f; d5 = 0.0f; d6 = 0.0f; d4 = d2 * d3; d5 = data_len / 1.0f; d6 = d4 / d5; n4 = 0; n5 = 0; n6 = 0; n4 = __int64(d1 * 10000000.0f + 0.5f); n5 = __int64(d6 * 10000000.0f + 0.5f); n6 = n4 - n5; fc[z1][z2] = n6 / 10000000.0f; } } ExplainEquation(fc, proot, zfx - 1); if (proot[zfx] == 0) { return -1; } d1 = 0.0f; for (z1 = 0; z1 < data_len; z1 ++) { d1 += sum[zfx][z1]; } for (z2 = 0; z2 < zfx; z2 ++) { d2 = 0.0f; for (z1 = 0; z1 < data_len; z1 ++) { d2 += sum[z2][z1]; } d1 -= d2 * proot[z2] / proot[zfx]; } hroot[0] = d1 / zfy; z1 = 0; for (z2 = 1; z2 < max_power; z2 ++) { hroot[z2] = proot[z1] / proot[zfx]; z1 ++; } d1 = 0.0f; for (z1 = 0; z1 < max_power; z1 ++) { d1 += proot[z1] * fc[zfx][z1]; } z1 = long(sqrt(d1 / hroot[zfx] / fc[zfx][zfx])); CMCToolApp::MKFreeArray(reinterpret_cast<void ***>(&sum)); CMCToolApp::MKFreeArray(reinterpret_cast<void ***>(&fc)); delete [] proot; delete [] pdata; return z1; } BOOL CFitFormula::CalculateCurveParameter(CDoubleArray *X,CDoubleArray *Y,long M,long N,CDoubleArray *A) { //X,Y -- X,Y两轴的坐标 //M -- 结果变量组数 //N -- 采样数目 //A -- 结果参数 register long i,j,k; double Z,D1,D2,C,P,G,Q; CDoubleArray B,T,S; B.SetSize(N); T.SetSize(N); S.SetSize(N); if(M>N)M=N; for(i=0;i<M;i++) (*A)[i]=0; Z=0; B[0]=1; D1=N; P=0; C=0; for(i=0;i<N;i++) { P=P+(*X)[i]-Z; C=C+(*Y)[i]; } C=C/D1; P=P/D1; (*A)[0]=C*B[0]; if(M>1) { T[1]=1; T[0]=-P; D2=0; C=0; G=0; for(i=0;i<N;i++) { Q=(*X)[i]-Z-P; D2=D2+Q*Q; C=(*Y)[i]*Q+C; G=((*X)[i]-Z)*Q*Q+G; } C=C/D2; P=G/D2; Q=D2/D1; D1=D2; (*A)[1]=C*T[1]; (*A)[0]=C*T[0]+(*A)[0]; } for(j=2;j<M;j++) { S[j]=T[j-1]; S[j-1]=-P*T[j-1]+T[j-2]; if(j>=3) { for(k=j-2;k>=1;k--) S[k]=-P*T[k]+T[k-