unit ParmSrch;

interface

uses  SysUtils, LBRMTypes, MyUtils;

const NumParms = 9;

type
   TParmSrch = class
   private
      StartParm     : integer;          // index number (1..NumParms) -- param to change first
      CurrParm      : integer;          // index number (1..NumParms) -- param currently changing
      PInfo         : TParameterInfo;
      ChangedParm   : array[1..NumParms] of integer;
      OkToReverse   : boolean;      // OK to reverse direction from up to down?
      VaryDirection : integer;
      ParamStartVal : single;
      LastParamChanged : integer;

      function ParamsAreFullyOptimized : boolean;

   public
      constructor Create;
      destructor  Destroy; override;

      procedure SetStartParm(SP : integer);
      procedure SetCurrentParm(CP : integer);
      procedure SetParmNow(Parms : TParmArray);
      procedure SetParmMin(Parms : TParmArray);
      procedure SetParmMax(Parms : TParmArray);
      procedure SetParmOpt(Parms : TParmArray);
      procedure SetParmDigitNow(Digits : TDigitArray);
      procedure SetParmDigitMax(Digits : TDigitArray);
      procedure SetAllParmInfo(PI : TParameterInfo);

      procedure GetStartParm(var SP : integer);
      procedure GetCurrentParm(var CP : integer);
      procedure GetParmNow(var Parms : TParmArray);
      procedure GetParmMin(var Parms : TParmArray);
      procedure GetParmMax(var Parms : TParmArray);
      procedure GetParmOpt(var Parms : TParmArray);
      procedure GetParmDigitNow(var Digits : TDigitArray);
      procedure GetParmDigitMax(var Digits : TDigitArray);
      procedure GetAllParmInfo(var PI : TParameterInfo);

      procedure ChangeDirection;
      procedure GoToNextParameter(var Done : boolean);

      //
      // overloaded procedure
      //
      procedure GetNextParmSet(RMSE_Better  : boolean;
                               var SearchIsDone : boolean); overload;
      procedure GetNextParmSet(var SearchIsDone : boolean); overload;

   end;


function  RoundToDigit(X : single; D : integer) : single;
procedure BumpIt(var X : single; D, Dir : integer);


implementation

{$RANGECHECKS OFF}

//----------------------------------
constructor TParmSrch.Create;
var
   I : integer;
begin
   StartParm := 0;
   CurrParm := 0;
   for I:=1 to NumParms do begin
      PInfo.ParmNow[I] := -99.99;
      PInfo.ParmMin[I] := -99.99;
      PInfo.ParmMax[I] := -99.99;
      PInfo.ParmOpt[I] := -99.99;
      PInfo.DigitNow[I]  := -9;
      PInfo.DigitMax[I]  := -9;
      ChangedParm[I] := 1;     // default is to try to change it again
   end;

   // We don't use the boundary condition info for the parameter search
   // so we don't need to initialize it -- just ignore it.

   OkToReverse      := TRUE;
   VaryDirection    := 1;  // up, to start
   LastParamChanged := 1;
   ParamStartVal    := -1;
end;

//----------------------------------
destructor TParmSrch.Destroy;
begin
   // anything to do here?

   //
   inherited Destroy;
end;


//---------------------------------------------------------------------------
//  Parameter search path is always from parameter 1 to NumParms and for each
//  parameter up then down.  (Of course, if we start with parameter 5 it
//  is 5,6,7,8,...NumParms,1,2,....)
//
procedure TParmSrch.GetNextParmSet(RMSE_Better  : boolean;
                               var SearchIsDone : boolean);
var
   CP, I : integer;
begin
   //
   // First, verify that the parameter info is valid.
   //
   SearchIsDone := TRUE;
   if (StartParm < 1) then exit;
   if (CurrParm < 1)  then exit;
   for I:=1 to NumParms do begin
      if (abs(PInfo.ParmNow[I] + 99.99) < 0.1) then exit;
      if (abs(PInfo.ParmMin[I] + 99.99) < 0.1) then exit;
      if (abs(PInfo.ParmMax[I] + 99.99) < 0.1) then exit;
      if (abs(PInfo.ParmOpt[I] + 99.99) < 0.1) then exit;
      if (PInfo.DigitNow[I] = -9) then exit;
   end;

   //
   // Set up for this iteration
   //
   CP := CurrParm;

   //
   // If the last set of parameters yielded an improved RMSE, then we want
   // to try to just keep on going as we were.  We will leave VaryDirection
   // the way it was.  The only thing we will do is set OkToReverse=FALSE.
   // This is because we've now established that we get better RMSE by going
   // this direction with this parameter and we'll not be going back the other
   // way, at least in this iteration.  Next time around (after optimizing the
   // rest of the parameters) we may want to go down, but that's another story.
   //
   // Note that if we're currently heading DOWN, this is redundant.  We never
   // switch from DOWN to UP on the same parameter anyway.  It doesn't hurt
   // anything, though, and reduces code complexity to just always do it when
   // the RMSE is better.
   //
   if (RMSE_Better) then begin
      OkToReverse := FALSE;
      ChangedParm[CP] := 1;    // it has been changed
   end;

   //
   // If the last set of parameters did NOT yield a better RMSE, then we
   // restore the last value and change either the direction or the parameter.
   // Note the assumption that the last parameter set was the optimal set
   // found so far.  This should always be the case, because we should never
   // go "2 away" from the optimum.
   //
   if (not RMSE_Better) then begin
      for I:=1 to NumParms do PInfo.ParmNow[I] := PInfo.ParmOpt[I];
      ChangeDirection;
   end;

   // If VaryDirection=0, we need to move to the next parameter and set
   // VaryDirection=1.
   //
   if (VaryDirection = 0) then begin
      GoToNextParameter(SearchIsDone);
      if (SearchIsDone) then exit;
   end;

   //
   // Get the new candidate parameter set.  We will use the
   // "Don't Care If It's Better" version of this procedure
   // to calculate it.  We've already set everything that we needed to.
   //
   GetNextParmSet(SearchIsDone);
end;


//---------------------------------------------------------------------------
//  Parameter search path is always from parameter 1 to NumParms and for each
//  parameter up then down.  (Of course, if we start with parameter 5 it
//  is 5,6,7,8,...NumParms,1,2,....)
//
//  This version of the procedure doesn't care if the RMSE was better or not.
//  It just goes to the next valid parameter set.
//
procedure TParmSrch.GetNextParmSet(var SearchIsDone : boolean);
var
   CP, D, I : integer;
   P, Pmax, Pmin : single;
   ValidParms : boolean;
begin
   //
   // First, verify that the parameter info is valid.
   //
   SearchIsDone := TRUE;
   if (StartParm < 1) then exit;
   if (CurrParm < 1)  then exit;
   for I:=1 to NumParms do begin
      if (abs(PInfo.ParmNow[I] + 99.99) < 0.1) then exit;
      if (abs(PInfo.ParmMin[I] + 99.99) < 0.1) then exit;
      if (abs(PInfo.ParmMax[I] + 99.99) < 0.1) then exit;
      if (abs(PInfo.ParmOpt[I] + 99.99) < 0.1) then exit;
      if (PInfo.DigitNow[I] = -9) then exit;
   end;

   //
   //  Note that if the Digits value for this parameter is 0, then we
   //  are not allowed to change it.
   //
   CP := CurrParm;
   D  := PInfo.DigitNow[CP];
   P  := RoundToDigit(PInfo.ParmNow[CP], D);
   if (D <> 0) then BumpIt(P, D, VaryDirection);
   PInfo.ParmNow[CP] := P;

   //
   // Now verify that this is a valid parameter set.
   // Set ValidParms to FALSE in order to force at least one iteration
   // through the checking process.
   //
   ValidParms := FALSE;
   while (not ValidParms) do begin
      ValidParms := TRUE;    // now assume it is a good one until proven otherwise.
      CP   := CurrParm;
      D    := PInfo.DigitNow[CP];
      P    := RoundToDigit(PInfo.ParmNow[CP], D);
      Pmin := PInfo.ParmMin[CP];
      Pmax := PInfo.ParmMax[CP];

      //
      //  If the DigitToChange was 0 then we didn't actually do ANYTHING, thus
      //  this is NOT a valid "NextParameterSet" (in general, it's the same as
      //  the optimal set calculated 2 iterations ago).  We need to go to the
      //  next parameter, so set VaryDirection to 0.
      //
      if (D = 0) then begin
         ValidParms := FALSE;
         VaryDirection := 0;
      end;

      //
      //  If the new value of the parameter is outside the max/min established,
      //  then it won't work either.  Handle this just as if the RMSE had been
      //  worse.
      //
      if ((P > Pmax) OR (P < Pmin)) then begin
         ValidParms := FALSE;
         // restore that best parameter set from last time
         for I:=1 to NumParms do PInfo.ParmNow[I] := PInfo.ParmOpt[I];
         ChangeDirection;
      end;

      if (not ValidParms) then begin
        //
        // Since VaryDirection=0, we need to move to the next parameter and
        // set VaryDirection=1.
        //
        if (VaryDirection = 0) then begin
           GoToNextParameter(SearchIsDone);
           if (SearchIsDone) then exit;
        end;

        CP   := CurrParm;
        D    := PInfo.DigitNow[CP];
        P    := RoundToDigit(PInfo.ParmNow[CP], D);
        BumpIt(P, D, VaryDirection);
        PInfo.ParmNow[CP] := P;
      end;

   end;

   //
   // Now we have a valid parameter set that needs to be run through the LBRM.
   // SearchIsDone = FALSE  because we are not done searching the parameters.
   // We will set it to TRUE when all of the parameters have been fully optimized.
   //
   SearchIsDone := FALSE;
end;


//---------------------------------------------------------------------
// Change the search direction.  In general we follow the Up-Down-Next
// pattern.
//
// Once we've gone UP once with better results we will never go DOWN
// (for this parameter, this iteration).  When we reach the max allowed
// value or get a worse RMSE we will go to the next parameter.
//
//
// There are 3 cases here.
//
//  1)  VaryDirection = +1 (up) and this was the first attempt to go
//      up with this parameter (this time).  In this case the variable
//      OkToReverse will be TRUE.  We can simply go back to the original
//      value then go down 1 (VaryDirection = -1).
//
//  2)  VaryDirection = +1 (up) and we had already established an "up" trend
//      (previous "ups" had improved the RMSE).  In this case the variable
//      OkToReverse will be FALSE.  We don't want to start going down, because
//      we already know that the "1 down" value of this parameter is the one
//      to keep.  We need to go back to it, but then jump to the next parameter.
//      So we set VaryDirection = 0.
//
//  3)  VaryDirection = -1 (down).  This one is simple -- go back up 1 then
//      start with the next parameter.  So we set VaryDirection = 0.
//
//  To simplify the code:
//     Cases 2 and 3 both have us setting VaryDirection to 0, so just test for
//     case 1.  If it is not case 1 then set VaryDirection to 0.
//
procedure TParmSrch.ChangeDirection;
var
   I, CP : integer;
begin
   CP := CurrParm;
   I := 0;    // for cases 2 & 3
   if ((VaryDirection = 1) AND (OkToReverse)) then begin  // if case 1
      I := -1;
      ChangedParm[CP] := -1;
   end;
   VaryDirection := I;
end;

//--------------------------------------------------------------------------
//  This routine is where we will move to the next parameter.
//  Before doing so, though, we will check to see if we have completely
//  optimized the parameters.  If so then Done will be set to TRUE.
//
procedure TParmSrch.GoToNextParameter(var Done : boolean);
var
   CP, D : integer;
   P     : single;
begin
   CP := CurrParm;
   D  := PInfo.DigitNow[CP];
   P := RoundToDigit(PInfo.ParmNow[CP], D);
   ChangedParm[CP] := -1;                       // assume it was not changed
   if (not REqual(P,ParamStartVal)) then begin  // it actually was changed
      ChangedParm[CP] := 1;
      LastParamChanged := CP;
   end;

   Done := FALSE;
   if (ParamsAreFullyOptimized) then begin
      Done := TRUE;
      exit;
   end;

   inc(CurrParm);
   if (CurrParm > NumParms) then CurrParm := 1;
   OkToReverse := TRUE;
   VaryDirection := 1;
   CP := CurrParm;
   D  := PInfo.DigitNow[CP];
   ParamStartVal := RoundToDigit(PInfo.ParmNow[CP], D);
end;


//--------------------------------------------------------------------------
function TParmSrch.ParamsAreFullyOptimized : boolean;
var
   I, D1, D2, MinD   : integer;
   Pmax, Pmin        : single;
   Optimum, MoreToDo : boolean;
begin
   if (CurrParm <> LastParamChanged) then begin
      Result := FALSE;
      exit;
   end;

   //
   // If the Digits setting is 0 or the max/min values are the same, then it
   // is impossible to change that parameter.  Force all such parameters to
   // indicate "not changed", since we have set them to "changed" by default
   // when we created the object or when we increment the digits.
   //
   for I:=1 to NumParms do begin
      D1   := PInfo.DigitNow[I];
      Pmin := RoundToDigit(PInfo.ParmMin[I], D1);
      Pmax := RoundToDigit(PInfo.ParmMax[I], D1);
      if (D1=0) then ChangedParm[I] := -1;
      if (Pmax = Pmin) then ChangedParm[I] := -1;
   end;


   //
   // Has the parameter set been completely optimized for the current
   // settings of the DigitToChange array?  This really means just checking to
   // see if all of the parameters have been "touched" without any being changed.
   //
   Optimum := TRUE;
   for I:=1 to NumParms do begin
      if (ChangedParm[I] = 1) then Optimum := FALSE;
   end;

   //
   // If the param set IS fully optimized for the current DigitToChange, verify
   // that the DigitToChange is set to all DigitMax (as far as we go) or 0 (user
   // doesn't want it changed).  If it is NOT, then increment as necessary then
   // return to the calling procedure (Result = FALSE).
   //
   if (Optimum) then begin
      MoreToDo := FALSE;
      MinD := 999999;
      for I:=1 to NumParms do begin
         D1 := PInfo.DigitNow[I];
         D2 := PInfo.DigitMax[I];
         if ((D1<>0) AND (D1<D2)) then begin
            MoreToDo := TRUE;
            if (D1 < MinD) then MinD := D1;
         end;
      end;
      if (MoreToDo) then begin
         Optimum := FALSE;      // more processing needs to be done
         for I:=1 to NumParms do begin
            D1 := PInfo.DigitNow[I];
            D2 := PInfo.DigitMax[I];
            if ((D1<>0) AND (D1<D2) AND (D1=MinD)) then inc(PInfo.DigitNow[I]);
            ChangedParm[I] := 1;   // try to change each one again
         end;
      end;
   end;

   Result := Optimum;
end;


//-------------------------------------
procedure TParmSrch.SetStartParm(SP : integer);
begin
   StartParm := SP;
end;

//-------------------------------------
procedure TParmSrch.SetCurrentParm(CP : integer);
begin
   CurrParm := CP;
end;

//-------------------------------------
procedure TParmSrch.SetParmNow(Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do PInfo.ParmNow[I] := Parms[I];
end;

//-------------------------------------
procedure TParmSrch.SetParmMin(Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do PInfo.ParmMin[I] := Parms[I];
end;

//-------------------------------------
procedure TParmSrch.SetParmMax(Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do PInfo.ParmMax[I] := Parms[I];
end;

//-------------------------------------
procedure TParmSrch.SetParmOpt(Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do PInfo.ParmOpt[I] := Parms[I];
end;

//-------------------------------------
procedure TParmSrch.SetParmDigitNow(Digits : TDigitArray);
var
   I : integer;
begin
   for I:=1 to NumParms do PInfo.DigitNow[I] := Digits[I];
end;

//-------------------------------------
procedure TParmSrch.SetParmDigitMax(Digits : TDigitArray);
var
   I : integer;
begin
   for I:=1 to NumParms do PInfo.DigitMax[I] := Digits[I];
end;

//-------------------------------------
procedure TParmSrch.SetAllParmInfo(PI : TParameterInfo);
var
   I : integer;
begin
   PInfo.IDString := PI.IDString;                     // actually unused
   PInfo.Area     := PI.Area;                         // actually unused
   PInfo.SDate    := PI.SDate;                        // actually unused
   PInfo.CDate    := PI.CDate;                        // actually unused
   PInfo.EDate    := PI.EDate;                        // actually unused
   PInfo.KCons    := PI.KCons;                        // actually unused
   PInfo.B1       := PI.B1;                           // actually unused
   PInfo.B2       := PI.B2;                           // actually unused
   PInfo.USZC     := PI.USZC;                         // actually unused
   PInfo.CurrentParm := PI.CurrentParm;               // actually unused
   for I:=1 to NumParms do begin
      PInfo.ParmNow[I]  := PI.ParmNow[I];
      PInfo.ParmMin[I]  := PI.ParmMin[I];
      PInfo.ParmMax[I]  := PI.ParmMax[I];
      PInfo.ParmOpt[I]  := PI.ParmOpt[I];
      PInfo.DigitNow[I] := PI.DigitNow[I];
      PInfo.DigitMax[I] := PI.DigitMax[I];
   end;
   for I:=1 to 12 do begin
      PInfo.Insolation[I] := PI.Insolation[I];
   end;
   for I:=1 to 50 do begin
      PInfo.Comments[I] := PI.Comments[I];
   end;
end;

//-------------------------------------
procedure TParmSrch.GetStartParm(var SP : integer);
begin
   SP := StartParm;
end;

//-------------------------------------
procedure TParmSrch.GetCurrentParm(var CP : integer);
begin
   CP := CurrParm;
end;

//-------------------------------------
procedure TParmSrch.GetParmNow(var Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do Parms[I] := PInfo.ParmNow[I];
end;

//-------------------------------------
procedure TParmSrch.GetParmMin(var Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do Parms[I] := PInfo.ParmMin[I];
end;

//-------------------------------------
procedure TParmSrch.GetParmMax(var Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do Parms[I] := PInfo.ParmMax[I];
end;

//-------------------------------------
procedure TParmSrch.GetParmOpt(var Parms : TParmArray);
var
   I : integer;
begin
   for I:=1 to NumParms do Parms[I] := PInfo.ParmOpt[I];
end;

//-------------------------------------
procedure TParmSrch.GetParmDigitNow(var Digits : TDigitArray);
var
   I : integer;
begin
   for I:=1 to NumParms do Digits[I] := PInfo.DigitNow[I];
end;

//-------------------------------------
procedure TParmSrch.GetParmDigitMax(var Digits : TDigitArray);
var
   I : integer;
begin
   for I:=1 to NumParms do Digits[I] := PInfo.DigitMax[I];
end;

//-------------------------------------
procedure TParmSrch.GetAllParmInfo(var PI : TParameterInfo);
var
   I : integer;
begin
   PI.IDString := PInfo.IDString;
   PI.Area     := PInfo.Area;
   PI.SDate    := PInfo.SDate;
   PI.EDate    := PInfo.EDate;
   PI.CDate    := PInfo.CDate;
   PI.KCons    := PInfo.KCons;
   PI.B1       := PInfo.B1;
   PI.B2       := PInfo.B2;
   PI.USZC     := PInfo.USZC;
   PI.CurrentParm := PInfo.CurrentParm;
   for I:=1 to NumParms do begin
      PI.ParmNow[I]  := PInfo.ParmNow[I];
      PI.ParmMin[I]  := PInfo.ParmMin[I];
      PI.ParmMax[I]  := PInfo.ParmMax[I];
      PI.ParmOpt[I]  := PInfo.ParmOpt[I];
      PI.DigitNow[I] := PInfo.DigitNow[I];
      PI.DigitMax[I] := PInfo.DigitMax[I];
   end;
   for I:=1 to 12 do begin
      PI.Insolation[I] := PInfo.Insolation[I];
   end;
   for I:=1 to 50 do begin
      PI.Comments[I] := PInfo.Comments[I];
   end;
end;


//--------------------------------------------------------------------------
function RoundToDigit(X : single; D : integer) : single;
var
   S1, S2 : string;
begin
   if ((D<1) OR (D>4)) then begin
      Result := X;
      exit;
   end;
   S1 := FormatFloat(' 0.000E+00;0.000E+00', X);    // string of form '-1.234E+12'

   // brute force technique .... ugly but effective
   if (D = 1) then S2 := Copy(S1, 1, 3) + '000' + Copy(S1, 7, 4);
   if (D = 2) then S2 := Copy(S1, 1, 4) + '00'  + Copy(S1, 7, 4);
   if (D = 3) then S2 := Copy(S1, 1, 5) + '0'   + Copy(S1, 7, 4);
   if (D = 4) then S2 := Copy(S1, 1, 6) +         Copy(S1, 7, 4);

   Result := StrToFloat(S2);
end;

//--------------------------------------------------------------------------
//
//  D   = digit to change (1..4)
//  Dir = direction to change (+1 or -1)
//
procedure BumpIt(var X : single; D, Dir : integer);
var
   Y, Xt  : single;
   S1, S2 : string;
begin
   if ((D<1) OR (D>4)) then exit;

   Xt := X;                                          // save for possible later use
   S1 := FormatFloat(' 0.000E+00;-0.000E+00', X);    // string of form '-1.234E+12'

   // brute force technique .... ugly but effective
   if (D = 1) then S2 := '1.000' + Copy(S1, 7, 4);
   if (D = 2) then S2 := '0.100' + Copy(S1, 7, 4);
   if (D = 3) then S2 := '0.010' + Copy(S1, 7, 4);
   if (D = 4) then S2 := '0.001' + Copy(S1, 7, 4);

   Y := StrToFloat(S2) * Dir;
   X := X + Y;

   //
   // If we just bumped down X to zero, then that's not what we wanted to do.
   // We need to adjust the exponent so that we can change the mantissa to 9.
   //
   if (abs(X) < 1.0e-20) then begin
      if (D = 1) then S2 := '.10000' + Copy(S1, 7, 4);
      if (D = 2) then S2 := '0.0100' + Copy(S1, 7, 4);
      if (D = 3) then S2 := '0.0010' + Copy(S1, 7, 4);
      if (D = 4) then S2 := '0.0001' + Copy(S1, 7, 4);
      Y := StrToFloat(S2) * Dir;
      X := Xt + Y;
   end;
end;


end.
