/******************************************************
*  PROGRAM TO DEMONSTRATE ONE DIMENSIONAL OPERATION   *
*    OF THE MULTI-NONLINEAR REGRESSION SUBROUTINE     *
* --------------------------------------------------- *
*  Ref.: BASIC Scientific Subroutines Vol. II,        *
*  By F.R. Ruckdeschel, Byte/McGRAW-HILL,1981 [1].    *
*                                                     *
*                  C++ version by J-P Moreau, Paris   *
*                         (www.jpmoreau.fr)           *
* --------------------------------------------------- *
* SAMPLE RUN:                                         *
*                                                     *
* MULTI-NON LINEAR REGRESSION                         *
*                                                     *
* Number of data points (mini 3, maxi 25)........: 11 *
* Degree of the polynomial to be fitted (maxi 10): 4  *
*                                                     *
* Input the data in (x,y) pairs as prompted:          *
*  1   X  Y = 1,1                                     *
*  2   X  Y = 2,2                                     *
*  3   X  Y = 3,3                                     *
*  4   X  Y = 4,4                                     *
*  5   X  Y = 5,5                                     *
*  6   X  Y = 6,6                                     *
*  7   X  Y = 7,7                                     *
*  8   X  Y = 8,8                                     *
*  9   X  Y = 9,9                                     *
*  10   X  Y = 0,0                                    *
*  11   X  Y = 1,1                                    *
*                                                     *
* The calculated coefficients are:                    *
*  0    -0.000000                                     *
*  1    1.000000                                      *
*  2    -0.000000                                     *
*  3    0.000000                                      *
*  4    -0.000000                                     *
*                                                     *
* Standard deviation:  0.000000                       *
*                                                     *
******************************************************/
#include <stdio.h>
#include <math.h>

#define MMAX  25
#define NMAX  10

double  D[NMAX+2], X[MMAX+1], Y[MMAX+1]; 
int     i,m,n;
double  dd,dx,xmn,xmx,xx,ymn,ymx,yy;

double  Z[MMAX+1][NMAX+2];      // maxi m=25 data points
double  A[MMAX+1][MMAX+1];      // maxi n=10 order of regression
double  B[MMAX+1][2*MMAX+1];
double  C[MMAX+1][MMAX+1];


// Matrix transpose B = Tr(A)
// A and B are n by m matrices
void Transpose()  {
  int i,j;
  for (i=1; i<n+1; i++) 
    for (j=1; j<m+1; j++) 
      B[i][j] = A[j][i];
}

// Matrix save (A in C) 
void A_IN_C(int n1,int n2) {
  int i1,i2;
  if (n1 * n2 != 0)  
    for (i1 = 1; i1<n1+1; i1++) 
      for (i2 = 1; i2<n2+1; i2++) 
        C[i1][i2] = A[i1][i2];
}

// Matrix save (B in A) 
void B_IN_A(int n1,int n2) {
  int i1,i2;
  if (n1 * n2 != 0)  
    for (i1 = 1; i1<n1+1; i1++) 
      for (i2 = 1; i2<n2+1; i2++) 
        A[i1][i2] = B[i1][i2];
}

// Matrix save (C in B) 
void C_IN_B(int n1,int n2) {
  int i1,i2;
  if (n1 * n2 != 0)  
    for (i1 = 1; i1<n1+1; i1++) 
      for (i2 = 1; i2<n2+1; i2++) 
        B[i1][i2] = C[i1][i2];
}

// Matrix save (C in A) 
void C_IN_A(int n1,int n2) {
  int i1,i2;
  if (n1 * n2 != 0)  
    for (i1 = 1; i1<n1+1; i1++) 
      for (i2 = 1; i2<n2+1; i2++) 
        A[i1][i2] = C[i1][i2];
}

// Matrix multiplication C = A * B
// A is m1 by n1, B is m2 by n2 ( m2 = n1 )
// Result C is m1 by n2.
void Matmult(int m1,int n1,int n2) {
  int i,j,k;
  for (i = 1; i < m1+1; i++) 
    for (j = 1; j < n2+1; j++) { 
      C[i][j] = 0;
      for (k = 1; k < n1+1; k++)
        C[i][j] += A[i][k] * B[k][j];
    }
}

// Matrix inversion: B = Inv(A) by Gauss-Jordan method
// A and B are n by n matrices
void Matinv()  {
  // Labels: E10, E20, E30;
  int i,j,k; double bb;
  for (i = 1; i < n+1; i++) {
    for (j = 1; j < n+1; j++) {
      B[i][j + n] = 0;
      B[i][j] = A[i][j];
    }
    B[i][i + n] = 1;
  }

  for (k = 1; k < n+1; k++) {
    if (k == n) goto E10;
    m = k;
    for (i = k + 1; i < n+1; i++)
      if (fabs(B[i][k]) > fabs(B[m][k]))  m = i;
    if (m == k) goto E10;
    for (j = k; j < 2*n+1; j++) {
      bb = B[k][j];
      B[k][j] = B[m][j];
      B[m][j] = bb;
    }
E10:for (j = k + 1; j < 2*n+1; j++)
      B[k][j] /= B[k][k];
    if (k == 1) goto E20;
    for (i = 1; i < k; i++)
      for (j = k + 1; j < 2*n+1; j++)
        B[i][j] -= B[i][k] * B[k][j];
    if (k == n) goto E30;
E20:for (i = k + 1; i < n+1; i++)
      for (j = k + 1; j < 2*n+1; j++)
        B[i][j] -= B[i][k] * B[k][j];
  } // k loop
E30:for (i = 1; i < n+1; i++)
      for (j = 1; j < n+1; j++)
        B[i][j] = B[i][j + n];
} // Matinv()

// Coefficient matrix generation
void Gene_Coeffs() {
  int i,j; double b;
  for (i = 1; i < m+1; i++) {
    b = 1;
    for (j = 1; j < n + 2; j++) {
      Z[i][j] = b;
      b = b * X[i];
    }
  }
} // Matinv()

// Least squares multidimensional fitting subroutine
void LstSqrN()  {
  int i,j,m4,n4;

  m4 = m;
  n4 = n;
  for (i = 1; i < m+1; i++)
    for (j = 1; j < n+1; j++)
      A[i][j] = Z[i][j];
  Transpose();                           // B=Transpose(A) 
  A_IN_C(m,n);                           // move A to C
  B_IN_A(n,m);                           // move B to A
  C_IN_B(m,n);                           // move C to B
  Matmult(n,m,n);                        // multiply A and B
  C_IN_A(n,n);                           // move C to A
  Matinv();                              // B=Inverse(A) 
  m = m4;                                // restore m
  B_IN_A(n,n);                           // move B to A
  for (i = 1; i < m+1; i++)
    for (j = 1; j < n+1; j++)
      B[j][i] = Z[i][j];
  Matmult(n,n,m);                        // multiply A and B
  C_IN_A(n,m);                           // move C to A
  for (i = 1; i < m+1; i++) B[i][1] = Y[i];
  Matmult(n,m,1);                        // multiply A and B

  // Product C is n by 1 - Regression coefficients are in C(i,1)
  // and put for convenience into vector D(i).
  for (i = 1; i < n+1; i++)  D[i] = C[i][1];

} // LstSqrN()

// Standard deviation calculation for a polynomial fit
void Standard_deviation() {
  int i, j; double b, yy;
  dd = 0;
  for (i = 1; i < m+1; i++) {
    yy = 0; b = 1;
    for (j = 1; j < n+2; j++) {
      yy = yy + D[j] * b;
      b = b * X[i];
    }
    dd += (yy - Y[i]) * (yy - Y[i]);
  }
  if (m - n - 1 > 0)
    dd /= m - n - 1;
  else
    dd = 0;
  dd = sqrt(dd);
}


void main() {
  for (i=0; i < NMAX+2; i++) D[i]=0;
  printf(" MULTI-NONLINEAR REGRESSION\n\n");
  printf(" Number of data points (mini 3, maxi 25)........: "); scanf("%d",&m);
  if (m<3) m=3; if (m>25) m=25;
  printf(" Degree of the polynomial to be fitted (maxi 10): "); scanf("%d",&n);
  if (n<1) n=1; if (n>10) n=10;
  if (m<n+2) m=n+2;
  printf("\n\n Input the data in (x,y) pairs as prompted\n");
  for (i = 1; i < m+1; i++) {
    printf(" %d  X  Y = ",i); scanf("%lf %lf",&X[i],&Y[i]);
  }
  printf("\n");
  // Call coefficients generation subroutine
  Gene_Coeffs();
  // Call regression subroutine
  n++;
  LstSqrN();
  printf(" The calculated coefficients are:\n");
  for (i = 1; i < n+1; i++) {
    printf(" %d  %f\n",i-1,D[i]);
  }
  n--;
  // Call standard deviation subroutine
  Standard_deviation();
  printf("\n\n Standard deviation: %f\n\n",dd);
  
} // end of main()

// End leastsqr.cpp