unit2.cpp

程序采用可视化界面

//---------------------------------------------------------------------------

#include <vcl.h>
#pragma hdrstop

#include "Unit2.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm2 *Form2;
//---------------------------------------------------------------------------
__fastcall TForm2::TForm2(TComponent* Owner)
    : TForm(Owner)
{
}
//---------------------------------------------------------------------------



void __fastcall TForm2::FormCreate(TObject *Sender)
{
    StringGrid1->Cells[0][0]="序号";
    StringGrid1->Cells[1][0]="Y";
}
//---------------------------------------------------------------------------
node * TForm2::initiate()
{ node *h, *p, *q;
  h=(node *)malloc(sizeof(node));      
  (*h).data=(*h).flag=(*h).b=0;                          /* 此处表明头指针不是头结点（元素） */
  p=(node *)malloc(sizeof(node));
  h->next=p;                          // 一定要注意
  int i=0;
  while(1)
  { p->data=0;
    p->flag=0;
	p->b=0;
    if(i==M-1) break;
    q=(node *)malloc(sizeof(node));
    p->next=q;    
    p=q;
    i++;    
  }
  p->next=NULL;                       /* 这是良好的习惯 */
  return h;  
}
void TForm2::average(float **x, float y[], float ax[], float *ay)
{
  int i,j;
  float sum1, sum2;
  for(j=0;j<M;j++)
  { sum1=0;
	for(i=0;i<N;i++)  sum1+=x[i][j];
    ax[j]=sum1/N;
  }
  sum2=0;
  for(i=0;i<N;i++)  sum2+=y[i];
  *ay=sum2/N;
}
void TForm2::augmented_matrix(float **l, float **x, float y[], float ay, float ax[])
{
  int i,j,k;
  for(i=0;i<M;i++)
  { for(j=0;j<=i;j++)
	{ for(k=0;k<N;k++)
	  { l[i][j]+=(x[k][i]-ax[i])*(x[k][j]-ax[j]);            //    对称矩阵
        l[j][i]=l[i][j];            
	  }
	}
    for(k=0;k<N;k++)
        l[i][M]+=(x[k][i]-ax[i])*(y[k]-ay);
  }
  for(k=0;k<N;k++)
    l[M][M]+=(y[k]-ay)*(y[k]-ay);
}
void TForm2::extended_augmented_matrix(float **l, float **r)
{ 
  int i,j;
  for(i=0;i<M;i++)
  for(j=0;j<=i;j++)
  { r[i][j]=l[i][j]/(sqrt(l[i][i])*sqrt(l[j][j]));
    r[j][i]=r[i][j];
    r[i][M]=l[i][M]/(sqrt(l[i][i])*sqrt(l[M][M]));
  }
  r[M][M]=1;
}
void TForm2::symbol(node *h, int flag, int m)          // 引入或删除因子K
{ 
  int i;
  node *q=(*h).next;
  for(i=1;i<m;i++) q=(*q).next;
  (*q).flag=flag;
}
float TForm2::input_b(node *h, float **r, float **l, float ax[], float ay)
{
  int j=0;
  float sum=0;
  node *q;
  q=(*h).next;
  while(q)
  { (*q).b=r[j][M]*sqrt(l[M][M]/l[j][j]);
    if(q->flag)  sum+=(*q).b*ax[j];
    q=q->next;
    j++;
  }
  return (ay-sum);
}
void TForm2::input_p(node *h, float **r)
{
  int j=0;
  node *p;
  p=(*h).next;
  while(p)
  { (*p).data=r[j][M]*r[j][M]/r[j][j];
    p=p->next;
    j++;
  }
}
float TForm2::read(node *h, int flag, int *f, int *s0)
{ 
  int i,j;
  float t;
  node *q=(*h).next;
  i=j=0;
  while(q)
  { if((*q).flag!=flag)
    { if(!j) t=(*q).data;             // 找第一个符合要求的元素
      j++;                            // 记录集合{J}或{Jc}中元素的个数
	  q=(*q).next;
    }
    else  q=(*q).next;
  }
  if(flag) *f=M-j;
  else *f=j;
  q=(*h).next;
  while(q)
  { i++;    
    if(flag)
    { if(t<=(*q).data&(*q).flag!=flag)
      { t=(*q).data;
        *s0=i;
      }
    }
    else   
    { if(t>=(*q).data&(*q).flag!=flag)
      { t=(*q).data;
        *s0=i;
      }  
    }
    q=(*q).next;
  }  
  return t;             
}
void TForm2::elimination(float **r, int k)
{ 
  k=k-1;                         // 注意，此处 i，j，k 都是从0开始计数，故此处的i，j，k要比《随机数据处理》中的值小1。
  static float **s;                // 第(s+1)的r[][]必须由第(s)回的r[][]计算。
  int i,j;
  s=new float*[M+1];
  for(i=0;i<M+1;i++)
  {
    s[i]=new float[M+1];
  }
  for(i=0;i<M;i++)
  for(j=0;j<M;j++)
  {
	if(i!=k&j!=k) s[i][j]=r[i][j]-r[i][k]*r[k][j]/r[k][k];
    else if(i==k&j!=k) s[i][j]=r[k][j]/r[k][k];
         else if(i==k&j==k) s[i][j]=1/r[k][k];
              else s[i][j]=-r[i][k]/r[k][k];
    if(i==k) s[i][M]=r[k][M]/r[k][k];
    else s[i][M]=r[i][M]-r[i][k]*r[k][M]/r[k][k];
  }
  s[M][M]=r[M][M];                  // 在这里，s[M+1][M+1]充当r(s+1)，而r[M+1][M+1]充当r(s)。                   
  for(i=0;i<=M;i++)
  for(j=0;j<=M;j++)  r[i][j]=s[i][j];

   //输出矩阵r[][]，用以检验其正确性。

  for(i=0;i<M+1;i++)
  { for(j=0;j<M+1;j++)  printf("%f\t",r[i][j]);
    puts("");
  }
  puts("");
  for(i=0;i<M+1;i++)
  {
    delete[] s[i];
  }
  delete[] s;
}
float TForm2::in(node *h, float **r)//二维数组用二维指针代替
{
  int fc, f, k0;
  float f1, p1, residual;
  input_p(h,r);
  p1=read(h, 1, &f, &k0);
  fc=M-f;
  if(!fc)  return 0;
  residual=r[M][M]-p1;             // 重点：注意此处算法在（引入）与（剔除）操作中的不同
  r[M][M]=residual;
  f1=p1*(N-f-2)/residual;
  if(f1<FIN)  return 0;        // 若不能引入,则返回0
  else
  { elimination(r, k0);
    symbol(h, 1, k0);
    return 1;
  }
}
int TForm2::out(node *h, float **r)
{
  int j0, f, sign;
  float f0, p0, residual;
  sign=0;                 // 若不能剔除,则sign为0
  int count;
  float pp=0;               // 用于记录对应于上一步所剔除因子的 P0 值
  count=1;
  while(1)
  {
	input_p(h,r);
	p0=read(h, 0, &f, &j0);
    if(!f)  break;
    if(count==1)  residual=r[M][M];
    else                               // 重点：注意此处算法在（引入）与（剔除）操作中的不同
	{ residual=r[M][M]+pp;
	  r[M][M]=residual;
	}
	pp=p0;                            // 用于记录对应于“上一步”所剔除因子的 P0 值，第一次删除因子除外
	f0=p0*(N-f-1)/residual;
    if(f0>FOUT)  break;
    sign++;
    elimination(r, j0);
    symbol(h, 0, j0);
	count++;
  }
  return sign;
}
void __fastcall TForm2::Button1Click(TObject *Sender)
{
    float b0, ay, *ax,**x,*y;
    ax=new float[M];
  int i, s1, s0, s=0,j;
  static float **l,**r;
  l=new float*[M+1]; //动态申请指针数组 ,先申请行再申请列
  r=new float*[M+1];
  for(i=0;i<M+1;i++)
  {
    l[i]=new float[M+1];
    r[i]=new float[M+1];
  }
  x=new float*[N];
  y=new float[N];
  for(i=0;i<N;i++)
  {
    x[i]=new float[M];
  }
  for(i=0;i<N;i++)
  {
    for(j=0;j<M;j++)
    {
        x[i][j]=(float)StrToFloat(StringGrid1->Cells[j+2][i+1]);
    }
    y[i]=(float)StrToFloat(StringGrid1->Cells[1][i+1]);
  }
  average(x, y, ax, &ay);
  augmented_matrix(l, x, y, ay, ax);
  extended_augmented_matrix(l, r);
  node *h;
  h=initiate();
  while(s<3)         // 经证明，第(s+1)和(s+1)步引入的变量，不可能在第(s+3)步中被剔除，所以此处先连续引入三个因子。
  { s1=in(h,r);
    if(!s1)  break;
    s++;
  }
  while(1)
  {
    s0=out(h,r);
    s1=in(h,r);
    if(!s0&!s1)  break;
  }
  b0=input_b(h, r, l, ax, ay);
  node *p=h->next;
  /*
  printf("\tb0=%f\n",b0);
  for(i=0;i<M;i++)
  {  printf("\tflag=%d\t\tb(%d)=%f\n",p->flag,i+1,p->b);
     p=p->next;
  }
  */
  AnsiString Text;
  Edit3->Text=Text.sprintf("Y = %.5f",b0);
  for(i=0;i<M;i++)
  {
    if(p->flag)
    {
        if(p->b<0)
        {
            Edit3->Text=Edit3->Text+Text.sprintf(" %.5f X(%d)",p->b,i+1);
        }
        else
        {
            Edit3->Text=Edit3->Text+Text.sprintf(" +%.5f X(%d)",p->b,i+1);
        }
    }
    p=p->next;
  }
  float sc, sh, sz;
  sc=l[M][M]*r[M][M];
  sh=l[M][M]*(1-sc);
  sz=l[M][M];
//  printf("lyy=%f\tSc=%f\tSh=%f\tSz=%f\n",l[M][M],sc,sh,sz);
    Edit4->Text=FloatToStr(l[M][M]);
    Edit5->Text=FloatToStr(sc);
    Edit6->Text=FloatToStr(sh);
    Edit7->Text=FloatToStr(sz);
    delete[] ax;
    for(i=0;i<M+1;i++)
    {
        delete[] l[i];
        delete[] r[i];
    }
    delete[] l;
    delete[] r;
    for(i=0;i<M;i++)
    {
        delete[] x[i];
    }
    delete[] x;
    delete[] y;
}
//---------------------------------------------------------------------------
void __fastcall TForm2::Edit1Exit(TObject *Sender)
{
    int m;
    m=StrToInt(Edit1->Text);
    StringGrid1->ColCount=m+2;
    int i;
    for(i=0;i<m;i++)
    {
        StringGrid1->Cells[i+2][0]="X"+IntToStr(i+1);
    }
    M=m;
}
//---------------------------------------------------------------------------

void __fastcall TForm2::Edit2Exit(TObject *Sender)
{
    int n;
    n=StrToInt(Edit2->Text);
    StringGrid1->RowCount=n+1;
    int i;
    for(i=0;i<n;i++)
    {
        StringGrid1->Cells[0][i+1]=IntToStr(i+1);   //Cells是列优先排列
    }
    N=n;

}
//---------------------------------------------------------------------------
#include <fstream.h>
void __fastcall TForm2::Button3Click(TObject *Sender)
{
    ofstream out;
    AnsiString FileName;
    float TempFloat;
    if(SaveDialog1->Execute())
    {
        FileName=SaveDialog1->FileName;
        out.open(FileName.c_str());
        if(!out)
        {
            return;
        }
        for(int i=0;i<N;i++)
        {
            for(int j=0;j<M+1;j++)
            {
                TempFloat=(float)StrToFloat(StringGrid1->Cells[j+1][i+1]);
                out<<TempFloat<<'\t';
            }
            out<<endl;
        }
        out.close();
    }

}
//---------------------------------------------------------------------------

void __fastcall TForm2::Button4Click(TObject *Sender)
{
    ifstream in;
    AnsiString FileName;
    float TempFloat;
    if(OpenDialog1->Execute())
    {
        FileName=OpenDialog1->FileName;
        in.open(FileName.c_str());
        if(!in)
        {
            return;
        }
        for(int i=0;i<N;i++)
        {
            for(int j=0;j<M+1;j++)
            {
                in>>TempFloat;
                StringGrid1->Cells[j+1][i+1]=FloatToStr(TempFloat);
            }
        }
        in.close();
    }
}
//---------------------------------------------------------------------------