{******************************************************
* 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].                                 *
*                                                     *
*              Pascal 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    1.95e-14                                      *
*  2    1.000000                                      *
*  3    1.83e-14                                      *
*  4    1.000000                                      *
*  5    1.69e-14                                      *
*  6    -2.80e-15                                     *
*                                                     *
* Standard deviation:  6.49e-13                       *
*                                                     *
******************************************************}
USES WinCrtMy,Type_def,Graph_2D;

CONST   LMAX = 9; MMAX = 25; NMAX = 25;

VAR
        D, Y : RV;     {pointers to real vectors size 2048}
        i,j,l,m,n : INTEGER;
        bb,cc,dd : real_ar;

        X : Array[1..MMAX,1..LMAX] of real_ar;     {maxi l:=9  number of dimensions}
        Z : Array[1..MMAX,1..NMAX] of real_ar;     {maxi m:=25 number of data points}
        A : Array[1..MMAX,1..MMAX] of real_ar;     {maxi n:=25 order of regression}
        B : Array[1..MMAX,1..2*MMAX] of real_ar;
        C : Array[1..MMAX,1..MMAX] of real_ar;
        M1: Array[1..LMAX] of integer;

        s1,s2 : STRING[12];
        buf   : STRING;
        titre : Array[0..40] of char;


{Matrix transpose}
PROCEDURE Transpose;
Var i,j : integer;
Begin
  FOR i := 1 TO n DO
    FOR j := 1 TO m DO
      B[i, j] := A[j, i]
End;

{Matrix save (A in C) }
PROCEDURE A_IN_C(n1,n2:INTEGER);
Var i1,i2 : integer;
Begin
  IF n1 * n2 <> 0 THEN 
    FOR i1 := 1 TO n1 DO
      FOR i2 := 1 TO n2 DO
        C[i1, i2] := A[i1, i2]
End;

{Matrix save (B in A) }
PROCEDURE B_IN_A(n1,n2:INTEGER);
Var i1,i2 : integer;
Begin
  IF n1 * n2 <> 0 THEN 
    FOR i1 := 1 TO n1 DO
      FOR i2 := 1 TO n2 DO
        A[i1, i2] := B[i1, i2]
End;

{Matrix save (C in B) }
PROCEDURE C_IN_B(n1,n2:INTEGER);
Var i1,i2 : integer;
Begin
  IF n1 * n2 <> 0 THEN 
    FOR i1 := 1 TO n1 DO
      FOR i2 := 1 TO n2 DO
        B[i1, i2] := C[i1, i2]
End;

{Matrix save (C in A) }
PROCEDURE C_IN_A(n1,n2:INTEGER);
Var i1,i2 : integer;
Begin
  IF n1 * n2 <> 0 THEN 
    FOR i1 := 1 TO n1 DO
      FOR i2 := 1 TO n2 DO
        A[i1, i2] := C[i1, i2]
End;

{Matrix multiplication}
PROCEDURE Matmult(m1,n1,n2:INTEGER);
Var i,j,k : integer;
Begin
  FOR i := 1 TO m1 DO
    FOR j := 1 TO n2 DO
    begin
      C[i, j] := 0;
      FOR k := 1 TO n1 DO
        C[i, j] := C[i, j] + A[i, k] * B[k, j]
    end
End;

{Matrix inversion}
PROCEDURE Matinv;
Label 10,20,30;
Var i,j,k : integer; bb:real_ar;
Begin
  FOR i := 1 TO n DO
  begin
    FOR j := 1 TO n DO
    begin
      B[i, j + n] := 0;
      B[i, j] := a[i, j]
    end;
    B[i, i + n] := 1
  end;
  FOR k := 1 TO n DO
  begin
    IF k = n THEN GOTO 10;
    m := k;
    FOR i := k + 1 TO n DO
      IF ABS(B[i, k]) > ABS(B[m, k]) THEN m := i;
    IF m = k THEN GOTO 10;
    FOR j := k TO 2 * n DO
    begin
      bb := B[k, j];
      B[k, j] := B[m, j];
      B[m, j] := bb
    end;
10: FOR j := k + 1 TO 2 * n DO
      B[k, j] := B[k, j] / B[k, k];
    IF k = 1 THEN GOTO 20;
    FOR i := 1 TO k - 1 DO
      FOR j := k + 1 TO 2 * n DO
        B[i, j] := B[i, j] - B[i, k] * B[k, j];
    IF k = n THEN GOTO 30;
20: FOR i := k + 1 TO n DO
      FOR j := k + 1 TO 2 * n DO
        B[i, j] := B[i, j] - B[i, k] * B[k, j]
  end; {k loop}
30: FOR i := 1 TO n DO
      FOR j := 1 TO n DO
        B[i, j] := B[i, j + n]
End; {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.  *
*********************************************************}
PROCEDURE LstSqrN;
Var  i,j,m4,n4 : integer;

Begin
  m4 := m;
  n4 := n;
  FOR i := 1 TO m DO
    FOR j := 1 TO n DO
      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 TO m DO
    FOR j := 1 TO n DO
      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 TO m DO 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) }
  FOR i := 1 TO n DO D^[i] := C[i, 1]
End;

{*************************************************************
* 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 M(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.           *
*************************************************************}
PROCEDURE Gene_Coeffs_SD;
Label fin;
Var i,j,k : integer; yy : real_ar;

  Procedure S40;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[1] do
    begin
      j := j + 1; Z[i, j] := bb; bb := bb * X[i, 1]
    end;
    bb := cc
  end;

  Procedure S50;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[2] do
    begin
      S40; bb := bb * X[i, 2]
    end;
    bb := cc
  end;

  Procedure S60;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[3] do
    begin
      S50; bb := bb * X[i, 3]
    end;
    bb := cc
  end;

  Procedure S70;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[4] do
    begin
      S60; bb := bb * X[i, 4]
    end;
    bb := cc
  end;

  Procedure S80;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[5] do
    begin
      S70; bb := bb * X[i, 5]
    end;
    bb := cc
  end;

  Procedure S90;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[6] do
    begin
      S80; bb := bb * X[i, 6]
    end;
    bb := cc
  end;

  Procedure S100;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[7] do
    begin
      S90; bb := bb * X[i, 7]
    end;
    bb := cc
  end;

  Procedure S110;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[8] do
    begin
      S100; bb := bb * X[i, 8]
    end;
    bb := cc
  end;

  Procedure S120;
  var i1:integer;
  begin
    cc := bb;
    FOR i1 := 0 TO M1[9] do
    begin
      S110; bb := bb * X[i, 9]
    end;
    bb := cc
  end;


Begin
  {Calculate the total number of dimensions}
  n := 1;
  FOR i := 1 TO l DO
    n := n * (M1[i] + 1);
  dd := 0;
  FOR i := 1 TO m DO
  begin
    {Branch according to dimension l (return if l > 9) }
    IF l < 1 THEN
    begin
      l := 0; goto fin
    end;
    IF l > 9 THEN
    begin
      l := 0; goto fin
    end;
    j := 0;
    IF l = 1 THEN begin bb:=1; S40 end;
    IF l = 2 THEN begin bb:=1; S50 end;
    IF l = 3 THEN begin bb:=1; S60 end;
    IF l = 4 THEN begin bb:=1; S70 end;
    IF l = 5 THEN begin bb:=1; S80 end;
    IF l = 6 THEN begin bb:=1; S90 end;
    IF l = 7 THEN begin bb:=1; S100 end;
    IF l = 8 THEN begin bb:=1; S110 end;
    IF l = 9 THEN begin bb:=1; S120 end;

    yy := 0;
    FOR k := 1 TO n DO
      yy := yy + D^[k] * Z[i, k];
    dd := dd + (Y^[i] - yy) * (Y^[i] - yy)
  end; {i loop}
  {Calculate standard deviation (if m > n) }
  IF m - n < 1 THEN
    dd := 0
  else
  begin
    dd := dd / (m - n);
    dd := SQRT(dd)
  end;

fin: End; {Gene_Coeffs_SD}

{main program}
BEGIN
  WinCrtInit('MLTNLREG');
  New(D); New(Y);
  Writeln(' MULTI-DIMENSIONAL AND MULTI-NONLINEAR REGRESSION');
  Writeln;
  Write(' How many data points ? '); read(m);
  Writeln;
  Write(' How many dimensions ? '); read(l);
  Writeln;
  FOR i := 1 TO l DO
  begin
    Write(' What is the fit for dimension ',i,' ? '); read(M1[i])
  end;
  n := 1;
  FOR i := 1 TO l DO n := n * (M1[i] + 1);
  IF m < n + 2 THEN m := n + 2;

  Writeln;
  Writeln(' Input the data as prompted:');
  Writeln;
  FOR i := 1 TO m DO
  begin
    Write(' Y(',i,') = ? '); read(y^[i]);
    FOR j := 1 TO l DO
    begin
      Write(' X(',i,',',j,') =  ? '); read(X[i, j])
    end;
    Writeln
  end;

  {Call coefficients generation subroutine}
  Gene_Coeffs_SD;

  {Call regression subroutine}
  LstSqrN;

  Writeln(' The calculated coefficients are:');
  Writeln;
  FOR i := 1 TO n DO
    Writeln('  ',i,'  ',D^[i]:9:6); 
  {Get standard deviation}
  n := n - 1;
  Gene_Coeffs_SD;
  Writeln;
  Writeln(' Standard deviation: ',dd:9:6);
  Writeln;
  Dispose(D); Dispose(Y);
  Readkey; DoneWinCrt
END.

{End of file mltnlreg.pas}