/******************************************************
*    Program to demonstrate the Parafit subroutine    *
* --------------------------------------------------- *
*  Reference: BASIC Scientific Subroutines, Vol. II   *
*  By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
*                                                     *
*               C++ Release 1.0 By J-P Moreau, Paris. *
*                        (www.jpmoreau.fr)            *
* --------------------------------------------------- *
* SAMPLE RUN:                                         *
*                                                     *
* Parametric least squares fit                        *
*                                                     *
* The input data are:                                 *
*                                                     *
* X( 1) =  1.000000   Y( 1) =  0.033702               *
* X( 2) =  2.000000   Y( 2) =  0.249029               *
* X( 3) =  3.000000   Y( 3) =  0.944733               *
* X( 4) =  4.000000   Y( 4) =  1.840089               *
* X( 5) =  5.000000   Y( 5) =  1.840089               *
* X( 6) =  6.000000   Y( 6) =  0.944733               *
* X( 7) =  7.000000   Y( 7) =  0.249029               *
* X( 8) =  8.000000   Y( 8) =  0.033702               *
* X( 9) =  9.000000   Y( 9) =  0.002342               *
* X(10) = 10.000000   Y(10) =  0.000084               *
*                                                     *
* The coefficients are:                               *
*  2.000000                                           *
*  4.500000                                           *
*  3.000000                                           *
*                                                     *
* The standard deviation of the fit is 0.000000       *
*                                                     *
* The number of iterations was 49                     *
******************************************************/
#include <stdio.h>
#include <math.h>

#define  SIZE  25

  double A[SIZE],E1[SIZE],X[SIZE],Y[SIZE];

  double a0,d,e,ee1,l1,l2,m0,m1,xx,yy;
  int    i,l,m,n;


  void S500()  {  //Function subroutine
    yy = A[1] * exp(-(xx - A[2]) * (xx - A[2]) / A[3]);
  }

  // Residual generation subroutine
  void S200()  {
    int j;
    l2 = 0;
    for (j = 1; j < n+1; j++) {
      xx = X[j];
      // Obtain function
      S500();
      l2 = l2 + (Y[j] - yy) * (Y[j] - yy);
    }
    d = sqrt(l2 / (n - l));
  }


/**************************************************************
*      Parametric least squares curve fit subroutine          *
* ----------------------------------------------------------- *
* This program least squares fits a function to a set of data *
* values by successively reducing the variance.  Convergence  *
* depends on the initial values and is not assured.           *
* n pairs of data values, X(i), Y(i), are given. There are l  *
* parameters, A(j), to be optimized across.                   *
* Required are initial values for the A(l) and e. Another     *
* important parameter which affects stability is e1, which is *
* initially converted to e1(l) for the first intervals.       *
* The parameters are multiplied by (1 - e1(i)) on each pass.  *
**************************************************************/
void Param_LS()  {
  // Labels: e50,e100, fin
  int i;
  for (i = 1; i < l+1; i++) E1[i] = ee1;
  m = 0;
  //Set up test residual
  l1 = 1e6;
  //Make sweep through all parameters
e50:
for (i = 1; i < l+1; i++)  {
    a0 = A[i];
    //Get value of residual
    A[i] = a0;
e100: S200();
    //Store result in m0
    m0 = l2;
    //Repeat for m1
    A[i] = a0 * (1.0 - E1[i]);
    S200();
    m1 = l2;
    //Change interval size if called for
    //If variance was increased, halve E1(i) 
    if (m1 > m0)  E1[i] = -E1[i] / 2.0;
    //If variance was reduced, increase step size by increasing E1(i)
    if (m1 < m0)  E1[i] = 1.2 * E1[i];
    //If variance was increased, try to reduce it
    if (m1 > m0)  A[i] = a0;
    if (m1 > m0)  goto e100;
  } // i loop
  //End of a complete pass
  //Test for convergence 
  m++;
  if (l2==0) goto fin;
  if (fabs((l1 - l2) / l2) < e)  goto fin;
  //If this point is reached, another pass is called for
  l1 = l2;
  goto e50;

fin: ;
} // Param_LS()


void main()  {
  n = 10; l = 3;
  printf(" Parametric least squares fit\n\n");
  
  printf(" The input data are:\n\n");
  
  for (i = 1; i < n+1; i++)  {
    X[i] = (double) i;
    Y[i] = 2.0 * exp(-(X[i] - 4.5) * (X[i] - 4.5) / 3.0);
    printf(" X(%2d) = %9.6f",i,X[i]); 
    printf("   Y(%2d) = %9.6f\n",i,Y[i]); 
  }
  printf("\n");
  e = 0.1; ee1 = 0.5; A[1] = 10.0; A[2] = 10.0; A[3] = 10.0;

  Param_LS();

  printf(" The coefficients are:\n");
  printf("  %f\n",A[1]); 
  printf("  %f\n",A[2]); 
  printf("  %f\n",A[3]); 
  
  printf("\n The standard deviation of the fit is %f\n\n",d);
  
  printf(" The number of iterations was %d\n\n",m);
}

// End parafit.cpp