#include <newmatnl.h>

Inheritance diagram for FindMaximum2:

Public Member Functions
void	Fit (ColumnVector &, int)

virtual	~FindMaximum2 ()

Private Member Functions
virtual void	Value (const ColumnVector &, bool, Real &, bool &)=0

virtual bool	NextPoint (ColumnVector &, Real &)=0

virtual Real	LastDerivative (const ColumnVector &)=0

Detailed Description

Definition at line 171 of file newmatnl.h.

Constructor & Destructor Documentation

◆ ~FindMaximum2()

virtual FindMaximum2::~FindMaximum2 ( )

inlinevirtual

Definition at line 178 of file newmatnl.h.

178 {} // to keep gnu happy

Member Function Documentation

◆ Fit()

void FindMaximum2::Fit	(	ColumnVector &	Theta,
		int	n_it
	)

Definition at line 21 of file newmatnl.cpp.

References GeneralMatrix::Nrows(), Tracer::ReName(), and x.

Referenced by NonLinearLeastSquares::Fit(), and MLE_D_FI::Fit().

 {
    Tracer tr("FindMaximum2::Fit");
    enum State {Start, Restart, Continue, Interpolate, Extrapolate,
       Fail, Convergence};
    State TheState = Start;
    Real z,w,x,x2,g,l1,l2,l3,d1,d2=0,d3;
    ColumnVector Theta1, Theta2, Theta3;
    int np = Theta.Nrows();
    ColumnVector H1(np), H3, HP(np), K, K1(np);
    bool oorg, conv;
    int counter = 0;
    Theta1 = Theta; HP = 0.0; g = 0.0;
 
    // This is really a set of gotos and labels, but they do not work
    // correctly in AT&T C++ and Sun 4.01 C++.
 
    for(;;)
    {
       switch (TheState)
       {
       case Start:
      tr.ReName("FindMaximum2::Fit/Start");
      Value(Theta1, true, l1, oorg);
      if (oorg) Throw(ProgramException("invalid starting value\n"));
 
       case Restart:
      tr.ReName("FindMaximum2::Fit/ReStart");
      conv = NextPoint(H1, d1);
      if (conv) { TheState = Convergence; break; }
      if (counter++ > n_it) { TheState = Fail; break; }
 
      z = 1.0 / sqrt(d1);
      H3 = H1 * z; K = (H3 - HP) * g; HP = H3;
      g = 0.0;                     // de-activate to use curved projection
      if (g==0.0) K1 = 0.0; else K1 = K * 0.2 + K1 * 0.6;
      // (K - K1) * alpha + K1 * (1 - alpha)
      //     = K * alpha + K1 * (1 - 2 * alpha)
      K = K1 * d1; g = z;
 
       case Continue:
      tr.ReName("FindMaximum2::Fit/Continue");
      Theta2 = Theta1 + H1 + K;
      Value(Theta2, false, l2, oorg);
      if (counter++ > n_it) { TheState = Fail; break; }
      if (oorg)
      {
         H1 *= 0.5; K *= 0.25; d1 *= 0.5; g *= 2.0;
         TheState =  Continue; break;
      }
      d2 = LastDerivative(H1 + K * 2.0);
 
       case Interpolate:
      tr.ReName("FindMaximum2::Fit/Interpolate");
      z = d1 + d2 - 3.0 * (l2 - l1);
      w = z * z - d1 * d2;
      if (w < 0.0) { TheState = Extrapolate; break; }
      w = z + sqrt(w);
      if (1.5 * w + d1 < 0.0)
         { TheState = Extrapolate; break; }
      if (d2 > 0.0 && l2 > l1 && w > 0.0)
         { TheState = Extrapolate; break; }
      x = d1 / (w + d1); x2 = x * x; g /= x;
      Theta3 = Theta1 + H1 * x + K * x2;
      Value(Theta3, true, l3, oorg);
      if (counter++ > n_it) { TheState = Fail; break; }
      if (oorg)
      {
         if (x <= 1.0)
            { x *= 0.5; x2 = x*x; g *= 2.0; d1 *= x; H1 *= x; K *= x2; }
         else
         {
            x = 0.5 * (x-1.0); x2 = x*x; Theta1 = Theta2;
            H1 = (H1 + K * 2.0) * x;
            K *= x2; g = 0.0; d1 = x * d2; l1 = l2;
         }
         TheState = Continue; break;
      }
 
      if (l3 >= l1 && l3 >= l2)
         { Theta1 = Theta3; l1 = l3; TheState =  Restart; break; }
 
      d3 = LastDerivative(H1 + K * 2.0);
      if (l1 > l2)
         { H1 *= x; K *= x2; Theta2 = Theta3; d1 *= x; d2 = d3*x; }
      else
      {
         Theta1 = Theta2; Theta2 = Theta3;
         x -= 1.0; x2 = x*x; g = 0.0; H1 = (H1 + K * 2.0) * x;
         K *= x2; l1 = l2; l2 = l3; d1 = x*d2; d2 = x*d3;
         if (d1 <= 0.0) { TheState = Start; break; }
      }
      TheState =  Interpolate; break;
 
       case Extrapolate:
      tr.ReName("FindMaximum2::Fit/Extrapolate");
      Theta1 = Theta2; g = 0.0; K *= 4.0; H1 = (H1 * 2.0 + K);
      d1 = 2.0 * d2; l1 = l2;
      TheState = Continue; break;
 
       case Fail:
      Throw(ConvergenceException(Theta));
 
       case Convergence:
      Theta = Theta1; return;
       }
    }
 }