// 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-