/******************************************************
* Program to demonstrate multi-dimensional operation  *
* of the multi-nonlinear regression subroutine.       *
* --------------------------------------------------- *
*  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                       *
*                                                     *
******************************************************/
#include <stdio.h>
#include <math.h>

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


double D[NMAX+1];
double Y[MMAX+1];

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

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()


  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()


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 coefficients generation subroutine
  Gene_Coeffs_SD();

  // Call regression subroutine
  LstSqrN();

  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);

} // end main()

// End MLTNLREG.CPP