欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

最小二乘法实现

程序员文章站 2022-06-04 10:52:32
...
#include <windows.h>
#include <stdio.h>
#include <math.h>
#include "COMMCTRL.H"
#define true TRUE
#define false FALSE
double Sigma[6];
static HFONT font;
static HWND mainHwnd;
static HDC mainHDC;
const char *szAppName="最小二乘法" ;
double x_value[]={36.9,46.7,63.7,77.8,84.0,87.5};
double y_value[]={181.0,197.0,235.0,270.0,283.0,292.0};

double det(double *A,int n)
{//|A| 矩阵A所有元素构成的行列式   
    double f=0.0f;
    double *B=(double*)malloc(sizeof(double)*(n-1)*(n-1));
    int i,j,k;
    if(n>2){
       for(j=0;j<n;j++){ 
          for(i=0;i<n;i++){
             if(i<j){
               for(k=0;k<n-1;k++){
                  B[k*(n-1)+i]=A[(k+1)*n+i];
               }
             }else if(i>j){
               for(k=0;k<n-1;k++){
                  B[k*(n-1)+i-1]=A[(k+1)*n+i];
               }
             }        
          }    
          if(j%2==1)f+=A[j]*det(B,n-1)*(-1);
          else      f+=A[j]*det(B,n-1); 
               
       }
          
    }
    else  f=A[0]*A[3]-A[1]*A[2];
  
    free(B);
    return f;
}

void one_variable_cubic(double *x,double *y,int n)
{
   int i,t;
   double X_pow[8]={0.0f};
   double *M=(double *)malloc(sizeof(double)*n);
   for(i=0;i<n;i++) {
   
   
   
   }
   free(M);
}
BOOL ShowSigma(HWND hwnd,double *value,int row,int col)
{  //有bug 只能拆开???  
   int i;
   char *s;
   for(i=0;i<col;i++){  
       s=(char*)calloc(12*sizeof(char));
       sprintf(s,"%.2f",value[i]);
       ListView_SetItemText(hwnd,row,i+1,s) ;
       free(s);
   }
}
BOOL ShowValue(HWND hwnd,double *value,int row,int col)
{ 
   int i;
   char *s;
   for(i=0;i<col;i++){  
       s=(char*)calloc(12*sizeof(char));
       sprintf(s,"%.2f",value[i]);
       ListView_SetItemText(hwnd,i,row,s) ;
       free(s);
   }
}
BOOL AddListViewItem(HWND hwnd,int row)
{
  LVITEM lvi;
  DWORD index;
  ZeroMemory(&lvi,sizeof(lvi));
  lvi.mask=LVIF_TEXT|LVIF_IMAGE|LVIF_PARAM|LVIF_STATE;
  lvi.state=0;
  lvi.stateMask=0;
  char *s;
  int length=sizeof(x_value)/sizeof(double);
  if(sizeof(y_value)/sizeof(double)!=length) {
      MessageBox(0,"x&y length is not equal","error",0);
      return false;
  }
  for(index=1;index<length+1;index++){
     s=(char*)calloc(10*sizeof(char));
     sprintf(s,"%d",index);
     lvi.iItem=index;
     lvi.iImage=0;
     lvi.iSubItem=row;
     lvi.pszText=s;//"undefine";
     lvi.cchTextMax=lstrlen(lvi.pszText) ;
     if(ListView_InsertItem(hwnd,&lvi)==-1)
        return false;
     free(s);
  }
  lvi.iItem=index;
  lvi.iImage=0;
  lvi.iSubItem=row;
  lvi.pszText="∑";
  lvi.cchTextMax=lstrlen(lvi.pszText) ;
  if(ListView_InsertItem(hwnd,&lvi)==-1)
       return false;
  return true;
}
void InitListViewColumns(HWND hwnd)
{
  char szText[256];
  LVCOLUMN lvc;
  DWORD i=0;
  lvc.mask=LVCF_FMT|LVCF_WIDTH|LVCF_TEXT|LVCF_SUBITEM;
  lvc.pszText=szText;
  lvc.cx=90;
  lvc.iImage=1;
  lvc.fmt=LVCF_FMT;
  char *s[]={" ","x","y","x^2","xy"}; 
  for(i=0;i<5;i++){
     lvc.pszText= s[i];//ColNames[i];
     lvc.iSubItem=i;
     ListView_InsertColumn(hwnd,i,&lvc);  
     free(s);
  }      
  if(!AddListViewItem(hwnd,0))
         MessageBox(0,"create listview failed","error",0); 
   ShowValue(hwnd,x_value,1,sizeof(x_value)/sizeof(double));
   ShowValue(hwnd,y_value,2,sizeof(x_value)/sizeof(double));
   mainHwnd=hwnd;
}
void double_mul(double *a,double *b,double *c,int len)
{
  int i;
  for(i=0;i<len;i++)c[i]=a[i]*b[i];
}
double Matrix(double *mat,int len,int a)
{ //y=ax+b,a=1计算a,a=0计算b;
    double ret,div;
    div=mat[2]*len-mat[0]*mat[0]+0.0001 ;
    if(a==1) ret=mat[3]*len-mat[0]*mat[1];
    else     ret=mat[2]*mat[1]-mat[0]*mat[3];
    ret=ret/div;
    return ret;
}
void fitting(double *sigma)
{     
      int i;
      int len=sizeof(x_value)/sizeof(double);      
      double *tmp=(double *)malloc(sizeof(double)*(len*2)); 
      double_mul(x_value,x_value,tmp,len);
      for(i=0;i<len;i++){
        sigma[0]+= x_value[i];
        sigma[1]+= y_value[i];
        sigma[2]+= tmp[i];
      }
      ShowValue(mainHwnd,tmp,3,len);
      double_mul(x_value,y_value,tmp,len);
      ShowValue(mainHwnd,tmp,4,len);
      for(i=0;i<len;i++)
            sigma[3]+= tmp[i];
      ShowSigma(mainHwnd,sigma,len,4);
      free(tmp);     
      sigma[4]=Matrix(sigma,len,1);
      sigma[5]=Matrix(sigma,len,0);      
}
void draw_piont(HDC hdc,int x,int y)
{
   int xLeft, yTop, xRight, yBottom;
   xLeft=x-3;
   xRight=x+3;
   yTop=y-3;
   yBottom=y+3;
   Ellipse (hdc, xLeft, yTop, xRight, yBottom) ;
}
void show_fit_line(HDC hdc,double *sigma)
{
   char s[32]={0};
   int x,y,i,len;
   double bx=0,by=0;
   HPEN pen=CreatePen(PS_SOLID,0,0X101010);
   HPEN bluepen=CreatePen(PS_SOLID,0,0Xff0000);
   HPEN redpen=CreatePen(PS_SOLID,0,0X0000ff);
   SetTextColor (hdc, 0x0000ff) ;
   SetBkColor (hdc, 0xf2f2f2) ;
   sprintf(s,"y=%.2fx+%.2f",sigma[4],sigma[5]);
   TextOut(hdc,640,360,s,strlen(s)) ;
     
   SelectObject(hdc,pen) ; 
   MoveToEx(hdc,620,300,NULL);
   LineTo(hdc,960,300) ;
   MoveToEx(hdc,620,300,NULL);
   LineTo(hdc,620,10) ; 
   
   
    
   SelectObject(hdc,bluepen) ;  
   len=sizeof(x_value)/sizeof(double);  
   for(i=0;i<len;i++){
     if(bx<x_value[i])bx=x_value[i]+0.001;
     if(by<y_value[i])by=y_value[i]+0.001;   
   }
   for(i=0;i<len;i++){
     x=(int)(x_value[i]*340/bx)+620;
     y=300-(int)(y_value[i]*290/by);     
     draw_piont(hdc,x,y);
   }   
   SelectObject(hdc,redpen) ; 
   POINT *porabi=(POINT*)malloc(1024*sizeof(POINT));   
   for(i=0;i<1024;i++){
     porabi[i].x=(int)(i*340/bx)+620;
     porabi[i].y=300-(int)((sigma[4]*i+ sigma[5])*290/by);   
   }   
   Polyline (hdc,porabi,1024);
   free(porabi);    
}       
void CreateButton(HWND hwnd,LPARAM lParam)
{
    static HWND hwndButton[3] ;
    hwndButton[0]=CreateWindow(TEXT("button"),"拟合",
              WS_CHILD|WS_VISIBLE|BS_FLAT,640,480,70,43,
              hwnd,(HMENU)0,((LPCREATESTRUCT)lParam)->hInstance,NULL);        
    hwndButton[1] =CreateWindow(WC_LISTVIEW,"dd",
             WS_CHILDWINDOW|WS_VISIBLE|LVS_REPORT|LVS_EDITLABELS,
             0,0,585,525,
            hwnd,(HMENU)1,((LPCREATESTRUCT)lParam)->hInstance,NULL);
    hwndButton[2]=CreateWindow(TEXT("button"),"一元三次",
              WS_CHILD|WS_VISIBLE|BS_FLAT,720,480,90,43,
              hwnd,(HMENU)2,((LPCREATESTRUCT)lParam)->hInstance,NULL); 
   InitListViewColumns(hwndButton[1]);
}
LRESULT CALLBACK SetChildFont(HWND hwnd, LPARAM lParam)
{
	SendMessage(hwnd,WM_SETFONT,(WPARAM)lParam,true);
	return true;
}
LRESULT CALLBACK WndProc(HWND hwnd,UINT message,WPARAM wParam,LPARAM lParam)
{
 int i;
 HDC hdc;
 PAINTSTRUCT ps;
 RECT  rect ,dataRect;
 static HANDLE hThread;
 static int flag=0;
 switch(message){
   case WM_CREATE:    
       CreateButton(hwnd,lParam);
       font=CreateFont(24,9,0,0,FW_NORMAL,false,false,
           0,ANSI_CHARSET,OUT_DEFAULT_PRECIS,CLIP_DEFAULT_PRECIS,
           DEFAULT_QUALITY,DEFAULT_PITCH|FF_SWISS,"微软雅黑");
       EnumChildWindows(hwnd,SetChildFont,(LPARAM)font);  break;    
   case WM_PAINT :
        hdc = BeginPaint (hwnd, &ps) ;
        mainHDC=hdc;
        GetClientRect (hwnd, &rect) ;
        SelectObject(hdc,font) ;
        if(flag)show_fit_line(hdc,Sigma);
        ReleaseDC(hwnd, hdc);
        DeleteDC(hdc);
        EndPaint (hwnd, &ps) ;break; 
   case WM_DESTROY:
        exit(0);
   case   WM_COMMAND :   
        switch(LOWORD (wParam)){ //对应按钮执行的任务
          case 0:fitting(Sigma); 
                 flag=1;
                InvalidateRect(hwnd,NULL,true);break;
          default:break;
        }
 }
 return DefWindowProc(hwnd,message,wParam,lParam);
}
int WINAPI WinMain(HINSTANCE hInstance,HINSTANCE hPrevInstance,
                                   PSTR szCmdLine,int nCmdShow)
{
   HBRUSH brush=CreateSolidBrush(0xf2f2f2);
   HWND hwnd=NULL;
   HICON hicon;
   MSG msg;
   WNDCLASS wc;
   hicon=LoadIcon(hInstance,MAKEINTRESOURCE(101));
   ZeroMemory(&wc,sizeof wc);
   wc.hInstance     = hInstance;
   wc.lpszClassName = szAppName;
   wc.lpfnWndProc   = (WNDPROC)WndProc;
   wc.style         = CS_DBLCLKS|CS_VREDRAW|CS_HREDRAW;
   wc.hbrBackground = brush;//(HBRUSH)GetStockObject(WHITE_BRUSH);
   wc.hIcon         = hicon;
   wc.hCursor       = LoadCursor(NULL,IDC_ARROW);
   if (false==RegisterClass(&wc))return 0;
   hwnd = CreateWindow (szAppName,szAppName,
           WS_OVERLAPPEDWINDOW,500,100,1000,570,
           HWND_DESKTOP,NULL,hInstance,NULL);
   ShowWindow(hwnd,nCmdShow);UpdateWindow(hwnd); 
   while(GetMessage(&msg,NULL,0,0)>0){
        TranslateMessage(&msg);
        DispatchMessage(&msg);
    }
   return msg.wParam;
}

最小二乘法实现

相关标签: math