/******************************************************
* Program to demonstrate the use of multi-dimensional *
*     Steepest Descent Optimization subroutine        *
* --------------------------------------------------- *
*  Reference: BASIC Scientific Subroutines, Vol. II   *
*  By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].  *
*                                                     *
*                 C++ version by J-P Moreau, Paris.   *
*                        (www.jpmoreau.fr)            *
* --------------------------------------------------- *
* Example:   Find a local maximum of                  *
*            F(x,y,z) = sin(x)+2*cos(y)-sin(z)        *
*                                                     *
* SAMPLE RUN:                                         *
*                                                     *
* How many dimensions: 3                              *
*                                                     *
* Convergence criterion: .000000001                   *
*                                                     *
* Maximum number of iterations: 50                    *
*                                                     *
* Starting constant: 1                                *
*                                                     *
* Input the starting point:                           *
*                                                     *
*     X[1] = 1                                        *
*     X[2] = 1                                        *
*     X[3] = 1                                        *
*                                                     *
* The results are:                                    *
*                                                     *
*     X[1] = 1.5707963                                *
*     X[2] = -0.0000435                               *
*     X[3] = -1.5707963                               *
*                                                     *
* Local maximum = 4.0000000                           *
*                                                     *
* The number of steps was: 30                         *
*                                                     *
******************************************************/
#include <stdio.h>
#include <math.h>

#define  MACHEPS  1e-15

double  D[4], Y[4];
double  X[11],X1[11];
double  e,xk;
int i,l,m,n;


/********************************************************
  Function subroutine with derivatives [L=3]           */   
double Eval() {
  D[1]=cos(X[1]); D[2]=-2.0*sin(X[2]); D[3]=-cos(X[3]);
  return sin(X[1])+2.0*cos(X[2])-sin(X[3]);
}
/*******************************************************/

  // Function called by Steepds()
  void Utilit() {
    int i; double dd;
    // Find the magnitude of the gradient
    dd=0.0;
    for (i=1; i<l+1; i++) dd += D[i]*D[i];
    dd=sqrt(dd);
    // Update the X[i] 
    for (i=1; i<l+1; i++) {
      // Save old values
      X1[i]=X[i];
      X[i] += xk*D[i]/dd;
    }
    Y[3]=Eval();
  }

/**************************************************
*    Steepest descent optimization subroutine     *
* ----------------------------------------------- *
* This routine finds the local maximum or minimum *
* of an L-dimensional function using the method   *
* of steepest decent, or the gradient.            *
* The function and the L derivatives must be      *
* available in subroutine Eval.                   *
* ----------------------------------------------- *
* INPUTS:                                         *
*   l - The dimension of function to study        *
*   e - The convergence criteria                  *
*   m - The maximum number of iterations          *
*   xk - A starting constant                      *
*   X[i] - Initial values of variables            *
* OUTPUTS:                                        *
*   X[i] - The locally optimum set                *
*   Eval - The value of local maximum             *
*   n - The number of iterations performed,       *
**************************************************/
  void Steepds()  {
  // Labels: e50,e51,e100,e200
  int i;
  n=0;
  // Start initial probe
  for (i=1; i<l+1; i++) {
    // Obtain yy and D[i]
    Y[i]=Eval();
    // Update X[i]
    Utilit();
  }
  // We now have a history to base the subsequent search on
  // Accelerate search if approach is monotonic 
e50: if (fabs(Y[2]-Y[1])<MACHEPS) goto e51;
  if ((Y[3]-Y[2])/(Y[2]-Y[1])>0.0) xk=xk*1.2;
  // Decelerate if heading the wrong way
e51: if (Y[3]<Y[2]) xk=xk/2.0;
  // Update the Y[i] if value has decreased
  if (Y[3]>Y[2]) goto e100;
  // Restore the X[i]
  for (i=1; i<l+1; i++) {
    X[i]=X1[i];
  }
  goto e200;
e100: Y[1]=Y[2]; Y[2]=Y[3];
  // Obtain new values
e200: Y[3]=Eval();
  // Update X[i]
  Utilit();
  // Check for maximum iterations and convergence
  n++;
  if (n>=m) return;
  if (fabs(Y[3]-Y[2])<e) return;
  // Try another iteration
  goto e50;
} // Steepds()


void main() {
  printf("\n How many dimensions: "); scanf("%d",&l);
  printf("\n Convergence criterion: "); scanf("%lf",&e);
  printf("\n Maximum number of iterations: "); scanf("%d",&m);
  printf("\n Starting constant: "); scanf("%lf",&xk);
  printf("\n Input the starting point:\n\n");
  for (i=1; i<l+1; i++) {
    printf("     X(%d) = ",i); scanf("%lf",&X[i]);
  }

  Steepds();   // Call steepest descent optimization subroutine

  printf("\n\n The results are:\n\n");
  for (i=1; i<l+1; i++) {
    printf("     X(%d) = %1.7f\n",i,X[i]);   
  }
  printf("\n Local maximum found = %10.7f\n",Eval());
  printf("\n The number of iterations was %d\n\n",n);
}

// End of file Steepds.cpp