/******************************************************
*  Program to demonstrate multidimensional operation  *
*  of the multi-nonlinear regression subroutine with  *
*  iterative error reduction                          *
* --------------------------------------------------- *
*  Ref.:  BASIC Scientific Subroutines Vol. II, By    *
*         F.R. Ruckdeschel, Byte/McGRAW-HILL, 1981    *
*         [BIBLI 01].                                 *
*                                                     *
*                C++ version by J-P Moreau, Paris     *
*                        (www.jpmoreau.fr)            *
* --------------------------------------------------- *
* SAMPLE RUN:                                         *
*                                                     *
* MULTI-DIMENSIONAL AND MULTI-NONLINEAR REGRESSION    *
*                                                     *
* How many data points ? 10                           *
*                                                     *
* How many dimensions ? 2                             *
*                                                     *
* What is the fit for dimension 1 ? 2                 *
* What is the fit for dimension 2 ? 1                 *
*                                                     *
* Input the data as prompted:                         *
*                                                     *
* Y( 1) = ? 7                                         *
* X( 1, 1) = ? 1                                      *
* X( 1, 2) = ? 6                                      *
*                                                     *
* Y( 2) = ? 7                                         *
* X( 2, 1) = ? 6                                      *
* X( 2, 2) = ? 1                                      *
*                                                     *
* Y( 3) = ? 6                                         *
* X( 3, 1) = ? 3                                      *
* X( 3, 2) = ? 3                                      *
*                                                     *
* Y( 4) = ? 8                                         *
* X( 4, 1) = ? 2                                      *
* X( 4, 2) = ? 6                                      *
*                                                     *
* Y( 5) = ? 9                                         *
* X( 5, 1) = ? 1                                      *
* X( 5, 2) = ? 8                                      *
*                                                     *
* Y( 6) = ? 9                                         *
* X( 6, 1) = ? 7                                      *
* X( 6, 2) = ? 2                                      *
*                                                     *
* Y( 7) = ? 6                                         *
* X( 7, 1) = ? 3                                      *
* X( 7, 2) = ? 3                                      *
*                                                     *
* Y( 8) = ? 7                                         *
* X( 8, 1) = ? 3                                      *
* X( 8, 2) = ? 4                                      *
*                                                     *
* Y( 9) = ? 7                                         *
* X( 9, 1) = ? 4                                      *
* X( 9, 2) = ? 3                                      *
*                                                     *
* Y( 10) = ? 2                                        *
* X( 10, 1) = ? 0                                     *
* X( 10, 2) = ? 2                                     *
*                                                     *
*                                                     *
* The calculated coefficients are:                    *
*                                                     *
*  1    -0.000000                                     *
*  2    1.000000                                      *
*  3    -0.000000                                     *
*  4    1.000000                                      *
*  5    -0.000000                                     *
*  6    0.000000                                      *
*                                                     *
* Standard deviation:  0.000000                       *
*                                                     *
* Number of iterations: 2                             *
*                                                     *
******************************************************/
#include <stdio.h>
#include <math.h>

#define  LMAX   9
#define  MMAX  25
#define  NMAX  25


double D[NMAX+1], D1[NMAX+1];
double Y[MMAX+1], Y1[MMAX+1];

int    i,j,l,l1,m,n;
double bb,cc,dd,dd1;

double X[MMAX+1][LMAX+1];              // maxi l:=9  number of dimensions
double Z[MMAX+1][NMAX+1];              // maxi m:=25 number of data points
double A[MMAX+1][MMAX+1];              // maxi n:=25 order of regression
double B[MMAX+1][2*MMAX+1];
double C[MMAX+1][MMAX+1];
int    M1[LMAX+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()


/*********************************************************
*   Least squares fitting subroutine, general purpose    *
* subroutine for multidimensional, nonlinear regression. *
* ------------------------------------------------------ *
* The equation fitted has the form:                      *
*     Y = D(1)X1 + D(2)X2 + ... + D(n)Xn                 *
* The coefficients are returned by the program in D(i).  *  
* The X[i) can be simple powers of x, or functions.      *      
* Note that the X[i) are assumed to be independent.      *   
* The measured responses are Y[i), there are m of them.  * 
* Y is a m row column vector, Z(i,j) is a m by n matrix. *
* m must be >= n+2. The subroutine inputs are m, n, Y[i) *
* and Z(i,j) previously calculated. The subroutine calls *
* several other matrix routines during the calculation.  *
*********************************************************/
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()

  // Functions used by Gene_Coeffs_SD()
  void S40(int i, int *j) {
    int i1,j1; 
    cc = bb; j1 = *j;
    for (i1 = 0; i1 < M1[1]+1; i1++) {
      j1++;  
      Z[i][j1] = bb; bb *= X[i][1];
    }
    bb = cc; *j = j1; 
  }

  void S50(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[2]+1; i1++) {
      S40(i,j); bb *= X[i][2];
    }
    bb = cc;
  }

  void S60(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[3]+1; i1++) {
      S50(i,j); bb *= X[i][3];
    }
    bb = cc;
  }

  void S70(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[4]+1; i1++) {
      S60(i,j); bb *= X[i][4];
    }
    bb = cc;
  }

  void S80(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[5]+1; i1++) {
      S70(i,j); bb *= X[i][5];
    }
    bb = cc;
  }

  void S90(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[6]+1; i1++) {
      S80(i,j); bb *= X[i][6];
    }
    bb = cc;
  }

  void S100(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[7]+1; i1++) {
      S90(i,j); bb *= X[i][7];
    }
    bb = cc;
  }

  void S110(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[8]+1; i1++) {
      S100(i,j); bb *= X[i][8];
    }
    bb = cc;
  }

  void S120(int i, int *j) {
    int i1;
    cc = bb;
    for (i1 = 0; i1 < M1[9]+1; i1++) {
      S110(i,j); bb *= X[i][9];
    }
    bb = cc;
  }

/***********************************************************
*       Coefficient matrix generation subroutine for       *
*              multiple non-linear regression.             *
* -------------------------------------------------------- *
* Also calculates the standard deviation d, even though    *
* there is some redundant computing.                       *
* The maximum number of dimensions is 9.                   *
* The input data set consists of m data sets of the form:  *
*   Y(i),X(i,1),X(i,2) ... X(i,l)                          *
* The number of dimensions is l.                           *
* The order of the fit to each dimension is M1(j).         *
* The result is an (m1+1)(m2+1)...(ml+1)+1 column by m row *
* matrix, Z. This matrix is arranged as follows            *
* (Ex.: l=2,M(1)=2,M(2)=2):                                *
* 1 X1 X1*X1 X2 X2*X1 X2*X1*X1 X2*X2 X2*X2*X1 X2*X2*X1*X1  *
* This matrix should be dimensioned in the calling program *
* as should also the X(i,j) matrix of data values.         *
***********************************************************/
void  Gene_Coeffs_SD()  {
  // Label: fin
  int i,j,k; double yy;
  // Calculate the total number of dimensions
  n = 1;
  for (i=1; i<l+1; i++) 
    n = n * (M1[i] + 1);
  dd = 0;
  for (i=1; i<m+1; i++) {
    //Branch according to dimension l (return if l < 1 or l > 9) 
    if (l < 1) {
      l = 0; goto fin;
    }
    if (l > 9) {
      l = 0; goto fin;
    }
    j = 0;
    if (l == 1) {bb=1; S40(i,&j);}
    if (l == 2) {bb=1; S50(i,&j);}
    if (l == 3) {bb=1; S60(i,&j);}
    if (l == 4) {bb=1; S70(i,&j);}
    if (l == 5) {bb=1; S80(i,&j);}
    if (l == 6) {bb=1; S90(i,&j);}
    if (l == 7) {bb=1; S100(i,&j);}
    if (l == 8) {bb=1; S110(i,&j);}
    if (l == 9) {bb=1; S120(i,&j);}

    yy = 0;
    for (k=1; k<n+1; k++)
      yy = yy + D[k] * Z[i][k];
    dd = dd + (Y[i] - yy) * (Y[i] - yy);
  } //i loop
  
  //Calculate standard deviation (if m > n) 
  if (m - n < 1) 
    dd = 0;
  else {
    dd = dd / (m - n);
    dd = sqrt(dd);
  }

fin:; } // Gene_Coeffs_SD()

  // begin internal procedures of Iter_Supervisor()
  void S160(int i, int *j) {
  int i1, j1;
    cc = bb; j1 = *j;
    for (i1 = 0; i1 < M1[1]+1; i1++) {
      j1++;
      Z[i][j1] = bb; bb = bb * X[i][1];
    }
    bb = cc; *j = j1;
  }

  void S170(int i, int *j) {
    int i2;
    cc = bb;
    for (i2 = 0; i2 < M1[2]+1; i2++) {
      S160(i,j); bb = bb * X[i][2];
    }
    bb = cc;
  }

  void S180(int i, int *j) {
    int i3;
    cc = bb;
    for (i3 = 0; i3 < M1[3]+1; i3++) {
      S170(i,j); bb = bb * X[i][3];
    }
    bb = cc;
  }

  void S190(int i, int *j) {
    int i4;
    cc = bb;
    for (i4 = 0; i4 < M1[4]+1; i4++) {
      S180(i,j); bb = bb * X[i][4];
    }
    bb = cc;
  }
 
  void S200(int i, int *j) {
    int i5;
    cc = bb;
    for (i5 = 0; i5 < M1[5]+1; i5++) {
      S190(i,j); bb = bb * X[i][5];
    }
    bb = cc;
  }
 
  void S210(int i, int *j) {
    int i6;
    cc = bb;
    for (i6 = 0; i6 < M1[6]+1; i6++) {
      S200(i,j); bb = bb * X[i][6];
    }
    bb = cc;
  }

  void S220(int i, int *j) {
    int i7;
    cc = bb;
    for (i7 = 0; i7 < M1[7]+1; i7++) {
      S210(i,j); bb = bb * X[i][7];
    }
    bb = cc;
  }

  void S230(int i, int *j) {
    int i8;
    cc = bb;
    for (i8 = 0; i8 < M1[8]+1; i8++) {
      S220(i,j); bb = bb * X[i][8];
    }
    bb = cc;
  }

  void S240(int i, int *j) {
    int i9;
    for (i9 = 0; i9 < M1[9]+1; i9++) {
      S230(i,j); bb = bb * X[i][9];
    }
    bb = cc;
  }

  void S250() {
    int i,j,k; double yy;
    for (i = 1; i < m+1; i++) {
      j = 0;
      if (l == 1) {bb=1; S160(i,&j);}
      if (l == 2) {bb=1; S170(i,&j);}
      if (l == 3) {bb=1; S180(i,&j);}
      if (l == 4) {bb=1; S190(i,&j);}
      if (l == 5) {bb=1; S200(i,&j);}
      if (l == 6) {bb=1; S210(i,&j);}
      if (l == 7) {bb=1; S220(i,&j);}
      if (l == 8) {bb=1; S230(i,&j);}
      if (l == 9) {bb=1; S240(i,&j);}
      //Array generated for row i
      yy = 0;
      for (k = 1; k < n+1; k++)  yy = yy + D[k] * Z[i][k];
      Y[i] = Y[i] - yy;
    }
  } //S250()


/*******************************************************************
*   Multi-dimensional polynomial regression iteration subroutine   *
* ---------------------------------------------------------------- *
* This routine supervises the calling of several other subroutines *
* in order to iteratively fit least squares polynomials in more    * 
* than one dimension.                                              *
* The routine repeatedly calculates improved coefficients until    *
* the standard deviation is no longer reduced. The inputs to the   *
* subroutine are the number of dimensions l, the degree of fit for *
* each dimension m(i), and the input data, X[i) and Y[i).          *
* The coefficients are returned in d(i), with the standard devia-  *
* tion in d. Also returned is the number of iterations tried, l1.  *
* y1(i), d1(i) and d1 are used respectively to save the original   *
* values of Y[i) and the current values of d(i) and d.             *
*******************************************************************/
void Iter_Supervisor()  {
  //Labels: e50, e100, fin
  l1 = 0;
  //Save the Y(i)
  for (i = 1; i < m+1; i++)  Y1[i] = Y[i];
  //Zero d1(i) 
  for (i = 1; i < n+1; i++)  D1[i] = 0;
  //Set the initial standard deviation high
  dd1 = 1e7;
  //Call coefficients subroutine
e50: Gene_Coeffs_SD();

  //Call regression subroutine
  LstSqrN();
  
  //Get standard deviation
  Gene_Coeffs_SD();

  //if standard deviation is decreasing, continue
  if (dd1 > dd)  goto e100;
  //Terminate iteration
  for (i = 1; i < n+1; i++)  D[i] = D1[i];
  //Restore the Y[i)
  for (i = 1; i < m+1; i++)  Y[i] = Y1[i];
  //Get the final standard deviation
  Gene_Coeffs_SD();
  goto fin;
  //Save the standard deviation
e100: dd1 = dd; l1++;
  //Augment coefficient matrix
  for (i = 1; i < n+1; i++) {
    D[i]  = D1[i] + D[i];
    D1[i] = D[i];
  }
  //Restore Y[i)
  for (i = 1; i < m+1; i++) Y[i] = Y1[i];
  //Reduce Y[i) according to the d(i) 
  S250();
  //We now have a set of error values
  goto e50;

fin:; } //Iter_Supervisor()


void main()  {

  printf(" MULTI-DIMENSIONAL AND MULTI-NONLINEAR REGRESSION\n\n");
  
  printf(" How many data points ? "); scanf("%d",&m);
  
  printf("\n How many dimensions ? "); scanf("%d",&l);

  printf("\n");

  for (i=1; i<l+1; i++) {
    printf(" What is the fit for dimension %d ? ",i); scanf("%d",&M1[i]);
  } 

  n = 1;
  for (i=1; i<l+1; i++)  n *= M1[i] + 1;
  if (m < n+2) m = n+2;  // m must be >= n+2 for LstSqrN()

  printf("\n Input the data as prompted:\n\n");
  
  for (i=1; i<m+1; i++) {
    printf(" Y[%d) = ? ",i); scanf("%lf",&Y[i]);
    for (j=1; j<l+1; j++) {
      printf(" X[%d,%d) =  ? ",i,j); scanf("%lf",&X[i][j]);
    }
    printf("\n");
  } 

  //Call iterations supervisor subroutine
  Iter_Supervisor();

  printf("\n The calculated coefficients are:\n\n");
  
  for (i=1; i<n+1; i++) {
    printf("  %d   %12.6f\n",i, D[i]); 
  }
  // Get standard deviation
  n--;
  Gene_Coeffs_SD();
  
  printf("\n Standard deviation: %12.6f\n\n",dd);

  printf(" Number of iterations: %d\n\n",l1);

}

// End REGITER.cpp