/********************************************************************/
/**                                                                **/
/**   mode.c                                                       **/
/**                                                                **/
/**   A universal generator for multivariate T-concave             **/
/**   distributions                                                **/
/**                                                                **/
/**   find mode of unimodal distribution                           **/
/**                                                                **/
/**   author: Josef Leydold                                        **/
/**   last modified: Mon Mar  9 13:38:00 MET 1998                  **/
/**                                                                **/
/********************************************************************/              

/*------------------------------------------------------------------*/

/* the header file */
#include "tdrmv.h"

/*------------------------------------------------------------------*/

#if MODE == 1

/*------------------------------------------------------------------*/

/** global variables                                               **/
/* pointer to parameters                                            */
static C_CONSTR *lhc;                      

/** functions                                                      **/
/* function to be minimized by Hooke's algorithm (see hooke.c)      */
double fkt_for_hooke(double *x, int n);

/* density function                                                 */
static double (*fdensity)();

/* write error message                                              */
extern void warning(char *msg, char *filename, int line);
/* fatal error, exit                                                */
extern void fatal(char *msg, char *filename, int line);


/*------------------------------------------------------------------*/

int get_mode(C_CONSTR *hc, double (*density)())
/*
  find mode of given density
*/
{
  int i;
  int nr_iter;            /* numbers of iterations */
  double startpt[N];      /* starting point for Hooke's algorithm */
  extern int hooke();     /* algorithm for finding minimum */
  
  /* set global pointer to density function and parameters */
  /* we need this information in fkt_for_hooke() */
  fdensity = density;
  lhc = hc;


#if RECTANGLE == 1
  /* first check bounds of rectangle */
  for( i=0; i<N; i++ )
    if( hc->rl[i] > hc->ru[i] )
      FATAL("rectangle: lower bound > upper bound");

  /* check if starting point is in domain (rectangle) */
  for( i=0; i<N; i++ )
    if( (hc->mode)[i] < (hc->rl)[i] || (hc->mode)[i] > (hc->ru)[i] ) {
      WARNING( "starting point for finding mode not in given domain (rectangle)" );
      (hc->mode)[i] = 0.5 * ((hc->rl)[i]+(hc->ru)[i]);
    }
#endif

  /* now check, if the starting point is not in the plain outside of the support of density */
  if( density(hc->mode) < TOLERANCE )
    WARNING( "density == 0 at starting point.\nProbably \"mode\" found in plain outside of support ==> algorithm fails!!" );


  /* find mode. first run */

  /* set starting point for Hooke's algorithm */
  /* (we use the default value for the mode */
  for( i=0; i<N; i++ )
    startpt[i] = hc->mode[i];
  /* run */
  nr_iter = hooke(fkt_for_hooke, N, startpt, hc->mode, hc->hooke_rho_begin, hc->hooke_epsilon, hc->hooke_imax);
  /*   printf("iter = %d\n",nr_iter); */

  /* find mode. second run */
  if( hc->hooke_rho_second > 0. ) {
    for( i=0; i<N; i++ )
      startpt[i] = hc->mode[i];
    nr_iter = hooke(fkt_for_hooke, N, startpt, hc->mode, hc->hooke_rho_second, hc->hooke_epsilon, hc->hooke_imax);
    /* printf("iter = %d\n",nr_iter); */
  }

#if RECTANGLE == 1
  /* check if mode is domain */
  for( i=0; i<N; i++ )
    if( hc->mode[i] < lhc->rl[i] || hc->mode[i] > lhc->ru[i] )
      FATAL( "mode not in given domain (rectangle)" );

  /* move mode to near boundary */
  for( i=0; i<N; i++ ) {
    if( (hc->mode[i] - hc->rl[i]) / (hc->ru[i] - hc->rl[i]) < hc->mode_to_boundary )
      hc->mode[i] = hc->rl[i];
    if( (lhc->ru[i] - hc->mode[i]) / (hc->ru[i] - hc->rl[i]) < hc->mode_to_boundary )
      hc->mode[i] = hc->ru[i];
  }
#endif

  return nr_iter;
} /* end of get_mode() */

/*-----------------------------------------------------------------*/

double fkt_for_hooke(double *x, int n)
/*
  function to be minimized by Hooke's algorithm
  = - density
*/
{

#if RECTANGLE == 1
  /* check if x is in the rectangle */
  {
    int i;
    for( i=0; i<N; i++ )
      if( x[i] < lhc->rl[i] || x[i] > lhc->ru[i] )
	return 0.;     /* not in rectangle --> f is set to 0 */
  }
#endif

  return( -fdensity(x) );

} /* end of fkt_for_hooke() */

/*-----------------------------------------------------------------*/
#endif
/*-----------------------------------------------------------------*/