Robmin

            In the Pascal version the derivative of FA is discontinuous.  This means that for very small values of chi-square it is possible that no change is made.  This caused the code to go on and on.  The isame added here tells Robmin to quit if nothing is changing.

            uRobMin.pas

The isame is defined globally next to IAExt.

 

    if IAExt < 7 then begin

      if Chi <= Chl then begin

        if chi/chl -Fr >  0.95*(1-Fr) then begin  sets the flag if the minimum is less than 1/20th that requested.

          isame := isame + 1;

          if isame = 3 then Fr := 0;

          if isame = 5 then iend := 0;  Sets iend to end the code

          writeln(' isame = ',isame); 

          end

        else

          isame := 0;  Resets isame in the event that a new minimum is found.

        end;

        Fr := Fr*Fr*Fr;  In testing the cube worked much better than the square. ß

 

In the newest code this section is ..\RobminSA\uRobMin.pas

  if( (abs(Chi-Chl)/Chl)<1e-9 ) and ( abs(Chi-Chb)/Chb <1e-9 )or  ß This assumes that the new minimum will be larger than 0.999´(1-Fr)

    ( chi <= Chl ) and  ((Chi/Chl-Fr) > 0.999*(1-Fr)) then begin

    Ifin:=Ifin+1;

    alamda:=alamda;

    if BobRobMin then   writeln(' CHI=CHB=CHL IFIN ',Ifin,' Iend ',Iend);

    if Ifin > 5 then begin

      if BobRobMin then writeln(' Robmin finds no lower values ');

      Init:=-6;

      Iend:=0;

      exit;

      end

 

ROBMIN.FOR

        IF(IAEXT.LT.20)THEN

          IF((CHI/CHL-FR).GT.0.95D0*(1-FR).AND.FR.GT.0D0)THEN

            ISAME=ISAME+1

            IF(ISAME.EQ.3)FR=0

            IF(ISAME.EQ.5)IEND=0

            IF(IPRINT.GE.1)PRINT*,' ISAME = ',ISAME

            FR=DMIN1(1D0-1D-7,.1D0*(9*FR+1))

          ELSE

            ISAME=0

            FR=FR*FR

            IF(CHI.LE.CHB+0.5*(CHL-CHB)) FR=DMAX1(1.D-2,SQRT(FR*FR*FR)) ß also note line above this

            IF((CHI-CHB)/(CHL-CHI+1.E-6).GT..1)

     2        FR=DMIN1(1D0,.1*(9*FR+1))

          ENDIF   

        ENDIF 

 

The function robmfu has changed.  The old without the comments is

 

function RobMfu(

  var Chi: double;

  var Pc: MaxConArrE;

  var Ppcc: NDimAArr;

  var Sm: MaxConArr;

  var B: MaxConArrE;

  var NCon: integer;

  var Fr,ALamb: double): double;

 

var

  Ajp: double;

  itest: integer;

begin

  GSolve(Chi,Pc,Ppcc,Sm,Alamb,NCon,Ajp,B);  ß Ajp was zero for a non-inverted matrix

  Result := Ajp/Chi - Fr;

  end;

           

The new version is

 

function RobMfu(

  Chi: double;

  var Ajp:double;

  Pc: array of double;

  Ppcc: array of extended;

  Sm: array of double;

  var B: array of double;

  NCon: integer;

  Fr,ALamb: double): double;

var

  k,ki,kt,m,Iflt: integer;

begin

  Iflt:=0;

  GSolve(Chi,Pc,Ppcc,Sm,Alamb,NCon,Iflt,B);

  if Iflt < 0 then begin     ß This is the same flag as before for a non-inverted matrix.  The new code set Ajp/Chi to -1, rather than 0.

    Result:=-1-Fr;

    exit;

    end;

// the solution for Ajp has been moved to here to utilize the original Ppcc. ß

  ki:=1;

  Ajp:=Chi+B[1]*(Pc[1]+Ppcc[1]*B[1]/2);

  for k:=2 to NCon do begin

    ki:=ki+k;

    Ajp:=Ajp+B[k]*(Pc[k]+Ppcc[ki]*B[k]/2);

    kt:=ki-k;

    for m:=1 to k-1 do begin

      kt:=kt+1;

      Ajp:=Ajp+Ppcc[kt]*B[m]*B[k];

      end;

    end;

  if Ajp > Chi then Ajp := -Chi;  ß with better matrix inversion, I have learned that we may also predict a maximum details are in solving\BracketRM2.doc .htm

  Result := Ajp/Chi - Fr;

  end;