/***************************************************
*  Parabolic Least Squares Demonstration Program   *
* ------------------------------------------------ *
* Reference: BASIC Scientific Subroutines, Vol. II *
* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
*                                                  *
*              C++ version by J-P Moreau, Paris    *
*              with possibility of a graph.        *
*    (to be used with lstsqr2.mak and graph2d.cpp) *
*                     (www.jpmoreau.fr)            *
* ------------------------------------------------ *
* This program calculates a parabolic least squares*
* fit to a given data set.                         *
*                                                  *
* INSTRUCTIONS                                     *
* ------------                                     *
*                                                  *
* The number of data coordinates provided must be  *
* greater than three. Otherwise, a divide by zero  *
* error may result.                                *
*                                                  *
* SAMPLE DATA:                                     *
*                                                  *
* Number of data points : 4                        *
*                                                  *
*   X[0]=1  Y[0]=1                                 *
*   X[1]=2  Y[1]=4                                 *
*   X[2]=3  Y[2]=9                                 *
*   X[3]=5  Y[3]=24.95                             *
*                                                  *
* Fitted equation is:                              *
*   Y = -0.017727 + 0.022045 X + 0.994318 X^2      *
*                                                  *
* Standard deviation of fit: 0.004767              *
*                                                  *
* (A graph is displayed).                          *
***************************************************/
#include <windows.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>


HDC  hdc;
RECT rect;
HPEN hpen, hpenOld;

//Functions used here of Graph2D.cpp
void Conversion(double,double,int *,int *);
void MinMax(int,float *,double *,double *);
void Initwindow(int,double,double,double,double);
void MoveXY(double,double);
void LineXY(double,double);
void Legendes (int,char *,char *,char *,int,int,char *,int,int,char *);


//"home made" graphic commands for hdc environment used by above functions
void DrawPixel(int ix,int iy) {
  //sorry, no other available command found
  Rectangle(hdc,rect.left+ix,rect.top+iy,
	            rect.left+ix+2,rect.top+iy+1);
}

void Swap(int *i1,int *i2) {
  int it;
  it=*i1; *i1=*i2; *i2=it;
}

void DrawLine(int ix1,int iy1,int ix2,int iy2) {
  int i,il,ix,iy;
  double dx,dy;
  if (ix2<ix1) {
	Swap(&ix1,&ix2); Swap(&iy1,&iy2);
  }
  il=(int) sqrt((ix2-ix1)*(ix2-ix1)+(iy2-iy1)*(iy2-iy1));
  dx=(double) (ix2-ix1)/il; 
  dy=(double) (iy2-iy1)/il;
  ix=ix1; iy=iy1;
  for (i=1;i<il+1; i++) {
	DrawPixel(ix,iy);
	ix=(int) (ix1+i*dx); 
	iy=(int) (iy1+i*dy);
  }
}

void OutText(int ix,int iy,char *s) {
  RECT rct;
  rct.left=ix; rct.right=ix+10*strlen(s);
  rct.top=iy-5; rct.bottom=iy+15;
  DrawText(hdc, s, -1, &rct,
                DT_SINGLELINE | DT_CENTER | DT_VCENTER);
}


#define SIZE  100

double a,b,c,d;
int    n;
float  X[SIZE], Y[SIZE];      // Real vectors


//Values may be changed by user
void Data() {
  n=4;
  X[0]=1; Y[0]=(float) 1.00;
  X[1]=2; Y[1]=(float) 4.00;
  X[2]=3; Y[2]=(float) 9.00;
  X[3]=5; Y[3]=(float) 24.95;
}

// Draw a cross at physical point (x,y)
void CrossXY(double x,double y) {
  int jx,jy;
  Conversion(x,y,&jx,&jy);
  DrawLine(jx-5,jy,jx+5,jy);
  DrawLine(jx,jy-5,jx,jy+5);
}

/****************************************************************
*           Parabolic least squares subroutine                  *
* ------------------------------------------------------------- *
* In:  integer n := number of points                            *
*      n values x[i), y(i) shared with main                     *
* Out: coefficients a,b,c of fit (a+b*x+c*x^2) shared with main *
*      standard deviation d shared with main                    *
****************************************************************/
void Least_Square(int n) {
  double a0,a1,a2,a3,a4,b0,b1,b2,d1;
  int m;
  a0 = 1;a1 = 0;a2 = 0;a3 = 0;a4 = 0;
  b0 = 0;b1 = 0;b2 = 0;
  for (m = 0; m<n; m++) {
    a1 = a1 + X[m];
    a2 = a2 + X[m] * X[m];
    a3 = a3 + X[m] * X[m] * X[m];
    a4 = a4 + X[m] * X[m] * X[m] * X[m];
    b0 = b0 + Y[m];
    b1 = b1 + Y[m] * X[m];
    b2 = b2 + Y[m] * X[m] * X[m];
  }
  a1 = a1 / n;a2 = a2 / n;a3 = a3 / n;a4 = a4 / n;
  b0 = b0 / n;b1 = b1 / n;b2 = b2 / n;
  d = a0 * (a2 * a4 - a3 * a3) - a1 * (a1 * a4 - a2 * a3) + a2 * (a1 * a3 - a2 * a2);
  a = b0 * (a2 * a4 - a3 * a3) + b1 * (a2 * a3 - a1 * a4) + b2 * (a1 * a3 - a2 * a2);
  a = a / d;
  b = b0 * (a2 * a3 - a1 * a4) + b1 * (a0 * a4 - a2 * a2) + b2 * (a1 * a2 - a0 * a3);
  b = b / d;
  c = b0 * (a1 * a3 - a2 * a2) + b1 * (a2 * a1 - a0 * a3) + b2 * (a0 * a2 - a1 * a1);
  c = c / d;
  // Evaluation of standard deviation d
  d = 0;
  for (m = 0; m<n; m++) {
    d1 = Y[m] - a - b * X[m] - c * X[m] * X[m];
    d += d1 * d1;
  }
  d = sqrt(d / (n - 3));
}

void Draw_graph() {
  double dx,x,xmn,xmx,ymn,ymx;
  extern MaxX, MaxY;
  char   buf[40],titre[80];
  int    m,npoints=25;
  //client drawing zone in pixels
  MaxX=(rect.right-rect.left);
  MaxY=(rect.bottom-rect.top);
  //seek minimum, maximum of table Y
  MinMax(n,X,&xmn,&xmx);
  MinMax(n,Y,&ymn,&ymx);
  Initwindow(10,xmn,xmx,ymn,ymx);
  //Crosses represent input data points
  for (m = 0; m<n; m++)
    CrossXY(X[m],Y[m]);
  //Draw parabolic line (25 points) 
  dx=(X[n-1]-X[0])/(npoints-1);
  x=X[0]-dx;
  MoveXY(X[0],c*X[0]*X[0]+b*X[0]+a);
  for (m=0; m<npoints; m++) {
    x+=dx;
    LineXY(x,c*x*x+b*x+a);
  }
  //Graph caption preparation (skip null coefficients)
  strcpy(titre,"  Y = ");
  if (a != 0) {sprintf(buf,"%f ",a); strcat(titre,buf);}
  if (b != 0) {
    if (b > 0 && a != 0) strcat(titre,"+ ");
    if (fabs(b) != 1)  {
      sprintf(buf,"%f X ",b); 
      strcat(titre,buf);
    }
    else if (b > 0) strcat(titre," X ");
    else strcat(titre," -X ");
  }
  if (c != 0) {
    if (c > 0 && (a != 0 || b != 0)) strcat(titre,"+ ");
    if (fabs(c) != 1) {
      sprintf(buf,"%f X^2 ",c); 
      strcat(titre,buf); 
    }
    else if (c > 0) strcat(buf," X^2 ");
    else strcat(titre," -X^2 ");
  }
  //Display graph caption, names of axes and 
  //standard deviation, d
  sprintf(buf," Standard deviation = %f",d);
  Legendes(10,titre,"X","Y",50,50,buf,0,0,"");
}

void Lstsqr2() {
  Data();
  Least_Square(n);
  Draw_graph();
}


//Handle window Paint and Destroy messages
long FAR PASCAL /*_export*/ WndProc(HWND hwnd,
    UINT message, UINT wParam, LONG lParam) 
{
   
    PAINTSTRUCT ps;

    switch (message) {
        case WM_PAINT:
	    hdc = BeginPaint(hwnd, &ps);
	    GetClientRect(hwnd, &rect);
            hpen = CreatePen(PS_SOLID, 1, RGB(0, 0, 255));
            hpenOld = SelectObject(hdc,hpen);

            Lstsqr2();
               
            SelectObject(hdc,hpenOld);
            DeleteObject(hpen);
            EndPaint(hwnd, &ps);
            return 0;
        case WM_DESTROY:
            PostQuitMessage(0);
            return 0;
    }
    return DefWindowProc(hwnd, message, wParam, lParam);
}

// main program
int PASCAL WinMain(HANDLE hInstance, HANDLE hPrevInstance,
    LPSTR lpszCmdParam, int nCmdShow)
{
    static char szAppName[] = "Lstsqr2";
    HWND        hwnd;
    MSG         msg;
    WNDCLASS    wndclass;

    if (!hPrevInstance) {
        wndclass.style         = CS_HREDRAW | CS_VREDRAW;
        wndclass.lpfnWndProc   = WndProc;
        wndclass.cbClsExtra    = 0;
        wndclass.cbWndExtra    = 0;
        wndclass.hInstance     = hInstance;
        wndclass.hIcon         = LoadIcon(NULL, IDI_APPLICATION);
        wndclass.hCursor       = LoadCursor(NULL, IDC_ARROW);
        wndclass.hbrBackground = GetStockObject(WHITE_BRUSH);
        wndclass.lpszMenuName  = NULL;
        wndclass.lpszClassName = szAppName;

        RegisterClass(&wndclass);
    }

    hwnd = CreateWindow(szAppName,  // window class name
        " PARABOLIC LEAST SQUARES", // window caption
        WS_OVERLAPPEDWINDOW,        // window style                       
        CW_USEDEFAULT,              // initial x position
        CW_USEDEFAULT,              // initial y position
        CW_USEDEFAULT,              // initial x size
        CW_USEDEFAULT,              // initial y size
        NULL,                       // parent window handle
        NULL,                       // window menu handle
        hInstance,                  // program instance handle
        NULL);                      // creation parameters

   ShowWindow(hwnd, nCmdShow);
   UpdateWindow(hwnd);

   while (GetMessage(&msg, NULL, 0, 0)) {
       TranslateMessage(&msg);
       DispatchMessage(&msg);
   }
   return msg.wParam;
}

//end of file lstsqr2.cpp