fitpro.cpp

专门处理曲线拟合数值计算的C++类

// FitPro.cpp : implementation file
//

#include "stdafx.h"
#include "dbfetch.h"
#include "FitPro.h"
#include "math.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// CFitPro dialog


CFitPro::CFitPro(CWnd* pParent /*=NULL*/)
	: CDialog(CFitPro::IDD, pParent)
{
	//{{AFX_DATA_INIT(CFitPro)
	m_stepno = 0;
	m_pvalue = 1.0;
	m_errcheck = FALSE;
	m_diffcheck = FALSE;
	m_diffe = 10.0;
	m_erre = 0.0;
	m_loopcheck = FALSE;
	m_loope = 1000;
	m_aveerr = 1.0;
	//}}AFX_DATA_INIT
	LoopFlag=0;
}


void CFitPro::DoDataExchange(CDataExchange* pDX)
{
	CDialog::DoDataExchange(pDX);
	//{{AFX_DATA_MAP(CFitPro)
	DDX_Control(pDX, IDC_PARA, m_paralist);
	DDX_Text(pDX, IDC_STEPNO, m_stepno);
	DDX_Text(pDX, IDC_PVALUE, m_pvalue);
	DDX_Check(pDX, IDC_ERROR, m_errcheck);
	DDX_Check(pDX, IDC_DIFF, m_diffcheck);
	DDX_Text(pDX, IDC_DIFFE, m_diffe);
	DDX_Text(pDX, IDC_ERRORE, m_erre);
	DDX_Check(pDX, IDC_NLOOP, m_loopcheck);
	DDX_Text(pDX, IDC_NLOOPE, m_loope);
	DDV_MinMaxUInt(pDX, m_loope, 0, 4294967295);
	DDX_Text(pDX, IDC_AVEERR, m_aveerr);
	//}}AFX_DATA_MAP
}


BEGIN_MESSAGE_MAP(CFitPro, CDialog)
	//{{AFX_MSG_MAP(CFitPro)
	ON_BN_CLICKED(IDC_START, OnStart)
	ON_BN_CLICKED(IDC_STOP, OnStop)
	ON_BN_CLICKED(IDC_PAUSE, OnPause)
	ON_BN_CLICKED(IDC_SAVE, OnSave)
	ON_LBN_SELCHANGE(IDC_PARA, OnSelchangePara)
	ON_EN_KILLFOCUS(IDC_PVALUE, OnKillfocusPvalue)
	ON_BN_CLICKED(IDC_ERROR, OnError)
	ON_BN_CLICKED(IDC_APPLY, OnApply)
	ON_BN_CLICKED(IDC_DIFF, OnDiff)
	ON_BN_CLICKED(IDC_CONTINUE, OnContinue)
	ON_BN_CLICKED(IDC_NLOOP, OnNloop)
	//}}AFX_MSG_MAP
END_MESSAGE_MAP()

/////////////////////////////////////////////////////////////////////////////
// CFitPro message handlers

void CFitPro::OnStart() 
{
	GetDlgItem(IDC_START)->EnableWindow(FALSE);
	GetDlgItem(IDC_SAVE)->EnableWindow(FALSE);
	GetDlgItem(IDC_CONTINUE)->EnableWindow(FALSE);
	GetDlgItem(IDOK)->EnableWindow(FALSE);
	//
	GetDlgItem(IDC_STOP)->EnableWindow(TRUE);
	GetDlgItem(IDC_PAUSE)->EnableWindow(TRUE);
	m_stepno=0;
	LoopFlag=0;
	OldPar=new double[parnum];
	dPar=new double[parnum];
	Par=new double[parnum];
	DPar=new double[parnum];
	DPai=new double[parnum];
	NegPai=new double[parnum];
	PosPai=new double[parnum];
	DDPai=new double[parnum*parnum];
	Pai=0;
	for(int i=0;i<parnum;i++)
	{
		DPar[i]=1.0;
	}
	UpdateData(FALSE);
	Iteration();
}

void CFitPro::OnStop() 
{
	GetDlgItem(IDC_START)->EnableWindow(TRUE);
	GetDlgItem(IDC_SAVE)->EnableWindow(TRUE);
	GetDlgItem(IDC_PAUSE)->EnableWindow(FALSE);
	GetDlgItem(IDC_STOP)->EnableWindow(FALSE);
	GetDlgItem(IDC_CONTINUE)->EnableWindow(FALSE);
	GetDlgItem(IDOK)->EnableWindow(TRUE);
	delete DDPai;
	delete PosPai;
	delete NegPai;
	delete OldPar;
	delete DPar;
	delete dPar;
	delete Par;
	delete DPai;
	LoopFlag=2;
}

void CFitPro::OnPause() 
{
	GetDlgItem(IDC_SAVE)->EnableWindow(TRUE);
	GetDlgItem(IDC_CONTINUE)->EnableWindow(TRUE);
	GetDlgItem(IDC_PAUSE)->EnableWindow(FALSE);
	LoopFlag=2;
}
void CFitPro::OnContinue() 
{
	GetDlgItem(IDC_SAVE)->EnableWindow(FALSE);
	GetDlgItem(IDC_CONTINUE)->EnableWindow(FALSE);
	GetDlgItem(IDC_PAUSE)->EnableWindow(TRUE);
	LoopFlag=1;
	Iteration();
}
void CFitPro::OnSave() 
{
	CFileDialog fdlg( FALSE,"txt", NULL, OFN_HIDEREADONLY, "Text Files(*.*)|*.txt" );//OFN_OVERWRITEPROMPT
	if(fdlg.DoModal()==IDOK)
	{
		CString fname, str,str1, strt;
		fname=fdlg.GetPathName();
		if(fname.IsEmpty())
		{
			AfxMessageBox("No file name is saved");
		}
		CStdioFile stdf;
		if(stdf.Open(fname,CFile::modeReadWrite|CFile::modeCreate|CFile::modeNoTruncate))
		{
			while(stdf.ReadString(str))
			{;}
			CTime ctm=CTime::GetCurrentTime();
			str=ctm.Format("%H:%M:%S, %A, %B %d, %Y");
			strt.Format("\r\n");
			str+=strt;
			stdf.WriteString(str);
			for(int i=0;i<parnum;i++)
			{
				str.Format("P(%d)=\t%#10.4G\r\n",i,ParValues[i]);
				stdf.WriteString(str);
			}
			stdf.Close();
		}
		else
		{
			AfxMessageBox("Open file failure");
			return;
		}
	}
}

BOOL CFitPro::OnInitDialog() 
{
	CDialog::OnInitDialog();
	
	for(int i=0;i<parnum;i++)
	{
		m_paralist.AddString(ParNames[i]);
		m_paralist.SetCurSel(0);
	}
	if(m_loopcheck)
	{
		GetDlgItem(IDC_NLOOPE)->EnableWindow(TRUE);
	}else
	{
		GetDlgItem(IDC_NLOOPE)->EnableWindow(FALSE);
	}
	if(m_errcheck)
	{
		GetDlgItem(IDC_ERRORE)->EnableWindow(TRUE);
	}else
	{
		GetDlgItem(IDC_ERRORE)->EnableWindow(FALSE);
	}
	if(m_diffcheck)
	{
		GetDlgItem(IDC_DIFFE)->EnableWindow(TRUE);
	}else
	{
		GetDlgItem(IDC_DIFFE)->EnableWindow(FALSE);
	}
	if(m_errcheck||m_loopcheck||m_diffcheck)
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(TRUE);
	}
	else
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(FALSE);
	}
	UpdateData(FALSE);
	GetDlgItem(IDC_STOP)->EnableWindow(FALSE);
	GetDlgItem(IDC_PAUSE)->EnableWindow(FALSE);
	GetDlgItem(IDC_SAVE)->EnableWindow(FALSE);
	GetDlgItem(IDC_CONTINUE)->EnableWindow(FALSE);
	
	return TRUE;  // return TRUE unless you set the focus to a control
	              // EXCEPTION: OCX Property Pages should return FALSE
}

void CFitPro::OnSelchangePara() 
{
	int index;
	index=m_paralist.GetCurSel();
	m_pvalue=ParValues[index];
	UpdateData(FALSE);
}

void CFitPro::OnKillfocusPvalue() 
{
	int index;
	index=m_paralist.GetCurSel();
	UpdateData();
	if(index>=0)
	ParValues[index]=m_pvalue;
}

void CFitPro::OnNloop() 
{
	UpdateData();
	if(m_loopcheck)
	{
		GetDlgItem(IDC_NLOOPE)->EnableWindow(TRUE);
	}else
	{
		GetDlgItem(IDC_NLOOPE)->EnableWindow(FALSE);
	}
	if(m_errcheck||m_loopcheck||m_diffcheck)
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(TRUE);
	}
	else
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(FALSE);
	}
}
void CFitPro::OnError() 
{
	UpdateData();
	if(m_errcheck)
	{
		GetDlgItem(IDC_ERRORE)->EnableWindow(TRUE);
	}else
	{
		GetDlgItem(IDC_ERRORE)->EnableWindow(FALSE);
	}
	if(m_errcheck||m_loopcheck||m_diffcheck)
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(TRUE);
	}
	else
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(FALSE);
	}
}
void CFitPro::OnDiff() 
{
	UpdateData();
	if(m_diffcheck)
	{
		GetDlgItem(IDC_DIFFE)->EnableWindow(TRUE);
	}else
	{
		GetDlgItem(IDC_DIFFE)->EnableWindow(FALSE);
	}
	if(m_errcheck||m_loopcheck||m_diffcheck)
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(TRUE);
	}
	else
	{
		GetDlgItem(IDC_APPLY)->EnableWindow(FALSE);
	}
}


void CFitPro::OnApply() 
{
	UpdateData();
}
double CFitPro::FuncValue(CExpression *pexp,CString *names,double *values, int count)
{
	double value;
	pexp->NumOfVar=count;
	pexp->VarNames=names;
	pexp->VarValues=values;
	pexp->Value(value);
	return value;
}
double CFitPro::PaiValue(double * pars)
{
	int i,j;
	double tmp,sum=0;
	for(i=0;i<datanum;i++)
	{
		for(j=0;j<varnum;j++)
		{
			ExpVValue[j]=VarValues[i][j+1];
		}
		for(j=0;j<parnum;j++)
		{
			ExpVValue[j+varnum]=pars[j];
		}
		tmp=VarValues[i][0]-FuncValue(m_pexp,ExpVName,ExpVValue,varnum+parnum);
		sum+=tmp*tmp;
	}
	return sum;
}
void CFitPro::Iteration()
{
	int i,j;
	double paitmp,dpmax,fc;
	do
	{
		m_stepno++;
		OldPai=Pai;
		for(i=0;i<parnum;i++)
		{
			OldPar[i]=ParValues[i];
			if((m_stepno/2)==((m_stepno+1)/2))
			{
				dPar[i]=fabs(ParValues[i])*1e-5;
			}else
			{
				dPar[i]=fabs(ParValues[i])*2e-1;
			}
			if(dPar[i]<1e-3)
			{
				dPar[i]=1e-3;
			}
		}
		Pai=PaiValue(ParValues);
		for(j=0;j<parnum;j++)
		{
			Par[j]=ParValues[j];
		}
		for(i=0;i<parnum;i++)
		{
			Par[i]=ParValues[i]+dPar[i];
			PosPai[i]=PaiValue(Par);
			Par[i]=ParValues[i]-dPar[i];
			NegPai[i]=PaiValue(Par);
			Par[i]=ParValues[i];
		}
		dpmax=0;
		for(i=0;i<parnum;i++)
		{
			DPai[i]=(PosPai[i]-NegPai[i])/(dPar[i]*2.0);
			DPar[i]=0;
			dpmax+=DPai[i]*DPai[i];
		}
		dpmax=sqrt(dpmax);
		for(i=0;i<parnum;i++)
		{
			for(j=0;j<i;j++)
			{
				Par[i]=ParValues[i]+dPar[i];
				Par[j]=ParValues[j]+dPar[j];
				paitmp=PaiValue(Par);
				DDPai[i*parnum+j]=((paitmp-PosPai[j])-(PosPai[i]-Pai))/dPar[i]/dPar[j];
				DDPai[j*parnum+i]=DDPai[i*parnum+j];
				Par[i]=ParValues[i];
				Par[j]=ParValues[j];
			}
			DDPai[i*parnum+i]=(PosPai[i]+NegPai[i]-Pai*2.0)/dPar[i]/dPar[i];
		}
		if((Pai/dpmax)>dpmax)
		{
			for(i=0;i<parnum;i++)
			{
				DPar[i]=DPai[i];
			}
		}else
		{
			for(i=0;i<parnum;i++)
			{
				DPar[i]=DPai[i]*Pai/dpmax/dpmax;
			}
		}
		fc=Factor(ParValues,DPar,parnum,1.0);
		if(fc<1e-3)
		{
			double tt=fabs((Range(DDPai,parnum)));
			if((tt/parnum)>1e-5)
			{
				if(parnum==1)
				{
					if(DDPai!=0)
					{
						DPar[0]=DPai[0]/fabs(DDPai[0]);
					}else
					{
						DPar[0]=0;
					}
				}else
				{
					CGause(DDPai, DPai, parnum, DPar);
				}
			}
			fc=Factor(ParValues,DPar,parnum,1.0);
		}
		for(i=0;i<parnum;i++)
		{
			ParValues[i]-=DPar[i]*fc;
		}
		m_aveerr=fabs(Pai)/datanum;
		if((m_aveerr>1e10)||(m_aveerr<-1e10))
		{
			m_aveerr=1e10;
		}
//		if((m_stepno/10)!=((m_stepno-1)/10))
		{
			int index=m_paralist.GetCurSel();
			if(index>=0)
			m_pvalue=ParValues[index];
			UpdateData(FALSE);
		}
		if(m_loopcheck&&(m_stepno>=m_loope))
		{
			OnStop();
		}
		if(m_errcheck)
		{
			if(paitmp<=m_erre)
			{
				OnStop();
			}
		}
		if(m_diffcheck)
		{
			double tmp;
			paitmp=0;
			for(i=0;i<parnum;i++)
			{
				if(fabs(ParValues[i])>1e-5)
				{
					tmp=fabs((ParValues[i]-OldPar[i])/ParValues[i])*1e2;
				}else
				{
					tmp=fabs((ParValues[i]-OldPar[i]))*1e2;
				}
				paitmp=tmp>paitmp?tmp:paitmp;
			}
			if(paitmp<=m_diffe)
			{
				OnStop();
			}
		}
		MSG msg;
		if(GetMessage(&msg, NULL, 0, 0)) 
		{
			TranslateMessage(&msg);
		    DispatchMessage(&msg);
		}
	}
	while(LoopFlag==0||LoopFlag==1);
}
int CFitPro::CGause(double *a, double *b, int n, double *x)
{
	int *js,l,k,i,j,is,p,q;
	double d,t;
	js=(int*)malloc(n*sizeof(int));
	l=1;
    for (k=0;k<=n-2;k++)
	{
		d=0.0;
        for (i=k;i<=n-1;i++)
		{
			for (j=k;j<=n-1;j++)
            {
				t=fabs(a[i*n+j]);
				if (t>d)
				{
					d=t; 
					js[k]=j;
					is=i;
				}
            }
			if (d+1.0==1.0)
			{
				l=0;
			}
			else
			{
				if (js[k]!=k)
					for (i=0;i<=n-1;i++)
					{
						p=i*n+k;
						q=i*n+js[k];
						t=a[p];
						a[p]=a[q];
						a[q]=t;
					}
				if (is!=k)
				{
					for (j=k;j<=n-1;j++)
					{
						p=k*n+j;
						q=is*n+j;
						t=a[p];
						a[p]=a[q];
						a[q]=t;
					}
					t=b[k];
					b[k]=b[is];
					b[is]=t;
				}
			}
			if (l==0)
			{
				free(js); printf("fail\n");
				return(0);
			}
			d=a[k*n+k];
			for (j=k+1;j<=n-1;j++)
			{
				p=k*n+j; a[p]=a[p]/d;
			}
			b[k]=b[k]/d;
			for (i=k+1;i<=n-1;i++)
			{
				for (j=k+1;j<=n-1;j++)
				{
					p=i*n+j;
					a[p]=a[p]-a[i*n+k]*a[k*n+j];
				}
				b[i]=b[i]-a[i*n+k]*b[k];
			}
		}
		d=a[(n-1)*n+n-1];
		if (fabs(d)+1.0==1.0)
		{
			free(js); printf("fail\n");
			return(0);
		}
		x[n-1]=b[n-1]/d;
		for (i=n-2;i>=0;i--)
		{
			t=0.0;
			for (j=i+1;j<=n-1;j++)
			{
				t=t+a[i*n+j]*x[j];
			}
			x[i]=b[i]-t;
		}
		js[n-1]=n-1;
	}
	for (k=n-1;k>=0;k--)
	{
		if (js[k]!=k)
		{
			t=x[k];
			x[k]=x[js[k]];
			x[js[k]]=t;
		}
	}
	free(js);
	return(1);
}
double CFitPro::Range(double *A, int n)
{
	if(n==1) return A[0];
	int i,j,k,pnt;
	double p=1,v=0;
	double* B=new double[(n-1)*(n-1)];
	for(i=0;i<n;i++)
	{
		pnt=0;
		for(j=0;j<n;j++)
		{
			if(j!=i)
			{
				for(k=1;k<n;k++)
				{
					B[pnt]=A[j*n+k];
					pnt++;
				}
			}
		}
		v=v+p*A[i*n]*Range(B,n-1);
		p=p*(-1.0);
	}
	delete [] B;
	return v;
}

double CFitPro::Factor(double *pars, double *dpars, int num, double fact)
{
	int i;
	double pai0,pai1,pai2;
	double * part=new double [num];
	for(i=0;i<num;i++)
	{
		part[i]=pars[i]-fact*dpars[i];
	}
	pai1=PaiValue(part);
	for(i=0;i<num;i++)
	{
		part[i]=pars[i]-0.5*fact*dpars[i];
	}
	pai2=PaiValue(part);
	pai0=PaiValue(pars);
	delete [] part;
	if(pai1<pai2&&pai1<pai0&&pai1>0)
	{
		return fact;
	}else
	{
		if(fact<1e-5)
		{
			if(pai2<PaiValue(pars))
			{
				return fact;
			}else
			{
				return 0;
			}
		}
		return Factor(pars,dpars,num,fact*0.5);
	}
}