/*----------------------------------------------------------------------------*
TEST FOR CONVERGENCE

This routine returns a measure of convergence in parameters being fit. It
carries a running average of each parameter and its rate of change, calculating
a discrete representation of 'w=(1/x)dx/dt' for each parameter. When that
measure sufficiently approaches zero for each parameter, the system has
converged.

For a single number to serve as the measure of convergence, this routine
computes and returns the euclidean distance from the origin of all 'n'
individual 'w's, scaled to the length of the diagonal of the unit hypercube for
that number of parameters. In other words, it returns the square root of the sum
of the squares of the individual 'w's, divided by the square root of 'n'.

There is a potential problem if any of the parameters actually approach zero,
since convergence will then be tested by the ratio of two random numbers
approaching zero. The distribution for that ratio can easily be without mean or
variance (e.g., the Cauchy results if the two numbers are normally distributed).
In such cases, the parameter should be biased before being passed to this
routine, such as by adding one to it.

Note: The routine drops the oldest set of parameters by subtracting it from the
totals, rather than re-summing all the values. This could introduce some
floating-point inaccuracies if it were done hundred of thousands of times, but
in convergence testing it will happen at most thousands of times, so
inaccuracies will not become an issue.

ENTRY: 't' contains the step number, which must be 0 on the very first call
         and advance by one on each call thereafter.
       'n' contains the number of parameters, which does not exceed 'NPAR'.
       'cpar' contains the current set of parameters.

EXIT:  'converge' contains the absolute-value rate of change of parameters,
         averaged over the last 'NSTEP' steps, or over how many steps are
         available. The closer this is to zero, the closer the fitting system is
         to equilibrium. At the very first call, when no differences can be
         calculated, 'converge' is negative.

WORK:  'vpar', 'dpar', 'vtot', 'dtot', and 'nstep' are updated on each call.
*/

#include "common.h"

#define NPAR  10               //Maximum number of parameters.
#define NSTEP 20               //Maximum number of steps in the running average.

static dec vpar[NSTEP][NPAR];  //Most recent parameters for running average.
static dec dpar[NSTEP][NPAR];  //Most recent differences for running average.
static dec vtot[NPAR];         //Most recent totals.
static dec dtot[NPAR];         //Most recent total absolute differences.
static dec nstep;              //Number of steps in the array.

dec converge(int t, int n, dec cpar[])
{ int i,j, h,k; dec d, e;

  if(t==0)                                   //If this is the first call, clear
  { for(i=0; i<NSTEP; i++)                   //the arrays to make the routine
    for(j=0; j<NPAR;  j++)                   //serially reuseable.
      vpar[i][j] = dpar[i][j] = 0;
    for(j=0; j<NPAR;  j++)
      dtot[j] = vtot[j] = 0;
    nstep = 0; }

  if(n>NPAR) { Error(477.0); n = NPAR; }     //Check the number of parameters.
  h = (t-1) % NSTEP; k = t % NSTEP;          //Locate the current and previous.

  for(j=0; j<n; j++)                         //Advance to the next parameter.
  {
    if(t>NSTEP)                              //If a new derivative is available,
    { dtot[j] -= abs(dpar[k][j]);            //subtract the old derivative from
      vtot[j] -= abs(vpar[k][j]); }          //its running total, also the value.

    vpar[k][j] = cpar[j];                    //Store the new value.

    if(t>0)                                  //If a new derivative is available,
    { dpar[k][j] = vpar[h][j]-vpar[k][j];    //add to its running total, and
      dtot[j] += abs(dpar[k][j]);            //also the value.
      vtot[j] += abs(vpar[k][j]); }

    if(nstep<NSTEP) nstep += 1; }            //Advance the number of items.

  if(t==0) return -1;                        //On the second and subsequent
  for(e=0, j=0; j<n; j++)                    //calls, compute the average
  { if(vtot[j]==0) { e = -1; break; }        //absolute derivative relative to
    d = dtot[j]/vtot[j];                     //the average parameter value, and
    e += d*d; }                              //accumulate the sum of squares.

  if(e<0) return e;                          //Return the euclidean distance to
  return sqrt(e/n);                          //the origin, relative to the
}                                            //diagonal.

// CLARENCE LEHMAN, JUNE 2011