-- optm2.adb   find minimum of an arbitrary function
--             no assumption that the function is continuous
--             this is reduced version z:=f(x,y)
--             given xmin, xmax, ymin, ymax, dxmin, dymin
--             and of course  f(x,y)
with TEXT_IO ; use TEXT_IO ;
with GENERIC_DIGITS_TYPES ;

procedure OPTM2 is
  package FLT_IO is new FLOAT_IO(LONG_FLOAT) ;
  use FLT_IO ;
  package INT_IO is new INTEGER_IO(INTEGER) ;
  use INT_IO ;

  type SIZE_100 is range 0 .. 100 ;

  package DIGITS_100 is new GENERIC_DIGITS_TYPES(SIZE_100, LONG_FLOAT) ;
  use DIGITS_100 ;

  function F(X, Y: D_NUMBER) return D_NUMBER is
    D10 : D_NUMBER := INTEGER_TO_D_NUMBER(10);
    D50 : D_NUMBER := INTEGER_TO_D_NUMBER(50);
    D60 : D_NUMBER := INTEGER_TO_D_NUMBER(60);
    D70 : D_NUMBER := INTEGER_TO_D_NUMBER(70);
    D80 : D_NUMBER := INTEGER_TO_D_NUMBER(80);
    Dfourth : D_NUMBER := INTEGER_TO_D_NUMBER(1) / INTEGER_TO_D_NUMBER(4);
  begin
    return Exp(Sin(D50*X))+ Sin(D60*Exp(Y)) +Sin(D70*Sin(X))+
           Sin(Sin(D80*Y))-Sin(D10*(X+Y)) +Dfourth*(X*X+Y*Y);
  end F;

  function STRING_TO_D_NUMBER(S : string) return D_NUMBER is
    result : D_NUMBER;
    last : integer;
  begin
    GET(S, result, last);
    return result;
  end STRING_TO_D_NUMBER;


  xmax  : D_NUMBER := STRING_TO_D_NUMBER(
          "-0.244030796943751719036133055699193898194358369381263E-1");
  xmin  : D_NUMBER := STRING_TO_D_NUMBER(
          "-0.244030796943751719036133083733373202681760784075469E-1");
  ymax  : D_NUMBER := STRING_TO_D_NUMBER(
          "0.210612427155355770591591103322758299678106470385244E0");
  ymin  : D_NUMBER := STRING_TO_D_NUMBER(
          "0.210612427155355770591591094208881035390578010456610E0");
  dxmin : D_NUMBER := REAL_TO_D_NUMBER( 1.0E-80);
  dymin : D_NUMBER := REAL_TO_D_NUMBER( 1.0E-80);

  D2 : D_NUMBER := INTEGER_TO_D_NUMBER(2);
  dx : D_NUMBER ;
  dy : D_NUMBER ;
  x, y, zxy : D_NUMBER ;
  xm, zxmy, zxmym : D_NUMBER ;
  xp, zxpy, zxmyp : D_NUMBER ;
  ym, zxym, zxpym : D_NUMBER ;
  yp, zxyp, zxpyp : D_NUMBER ;
  hxm, zhxmym, zhxmy, zhxmyp : D_NUMBER ;
  hxp, zhxpym, zhxpy, zhxpyp : D_NUMBER ;
  hym, zxmhym, zxhym, zxphym : D_NUMBER ;
  hyp, zxmhyp, zxhyp, zxphyp : D_NUMBER ;
  xbest, ybest, zbest : D_NUMBER ;
  xhbox, yhbox : D_NUMBER ;
  improve : integer :=1;
  iterations : integer := 400;

begin

  PUT("optm2  xmin="); PUT(xmin); NEW_LINE;
  PUT("       xmax="); PUT(xmax); NEW_LINE;
  PUT("       ymin="); PUT(ymin); NEW_LINE;
  PUT("       ymax="); PUT(ymax); NEW_LINE;
  PUT("      dxmin="); PUT(dxmin); NEW_LINE;
  PUT("      dymin="); PUT(dymin); NEW_LINE;
  PUT(" iterations="); PUT(iterations); NEW_LINE;
  PUT("     digits="); PUT("SIZE_100"); NEW_LINE;

  dx:=dxmin;
  dy:=dymin;
  x:=(xmin+xmax)/D2;
  y:=(ymin+ymax)/D2;

  -- initial consistency test
  if Xmax <= xmin then PUT("xmax<=xmin error"); return; end if;
  if Ymax <= ymin then PUT("ymax<=ymin error"); return; end if;
  zxy:=f(x,y);
  PUT("f(x,y)    z="); PUT(zxy); NEW_LINE;
  PUT("          x="); PUT(x); NEW_LINE;
  PUT("          y="); PUT(y); NEW_LINE;
  zbest:=zxy;
  xbest:=x;
  ybest:=y;

  for I in 0..iterations loop -- new zxy,x,y  new xmin,ymin  from below

    if improve=0 then PUT("no improvement"); exit; end if;
    improve:=0;
    NEW_LINE;

    xhbox:=(xmax-xmin)/D2;
    if xhbox<dxmin then  PUT("solved dxmin"); exit; end if;
    yhbox:=(ymax-ymin)/D2;
    if yhbox<dymin then  PUT("solved dymin"); exit; end if;
    xmin:=x-xhbox;
    xmax:=x+xhbox;
    ymin:=y-yhbox;
    ymax:=y+yhbox;

    xm:=xmin;
    zxmy:=f(xm,y);
    PUT("f(xm,y)   z="); PUT(zxmy); NEW_LINE;
    PUT("         xm="); PUT(xm); NEW_LINE;
    PUT("          y="); PUT(y); NEW_LINE;
    if zxmy < zbest then
      zbest:=zxmy;
      xbest:=xm;
      ybest:=y;
      improve:=1;
      PUT("bad xmin");
      goto continue;
    end if;

    xp:=xmax;
    zxpy:=f(xp,y);
    PUT("f(xp,y)   z="); PUT(zxpy); NEW_LINE;
    PUT("         xp="); PUT(xp); NEW_LINE;
    PUT("          y="); PUT(y); NEW_LINE;
    if zxpy < zbest then
      zbest:=zxpy;
      xbest:=xp;
      ybest:=y;
      improve:=1;
      PUT("bad xmax");
      goto continue;
    end if;

    ym:=ymin;
    zxym:=f(x,ym);
    PUT("f(x,ym)   z="); PUT(zxym); NEW_LINE;
    PUT("          x="); PUT(x); NEW_LINE;
    PUT("         ym="); PUT(ym); NEW_LINE;
    if zxym < zbest then
      zbest:=zxym;
      xbest:=x;
      ybest:=ym;
      improve:=1;
      PUT("bad ymin");
      goto continue;
    end if;

    yp:=ymax;
    zxyp:=f(x,yp);
    PUT("f(x,yp)   z="); PUT(zxyp); NEW_LINE;
    PUT("          x="); PUT(x); NEW_LINE;
    PUT("         yp="); PUT(yp); NEW_LINE;
    if zxyp < zbest then
      zbest:=zxyp;
      xbest:=x;
      ybest:=yp;
      improve:=1;
      PUT("bad ymax");
      goto continue;
    end if;

    zxmym:=f(xm,ym);
    PUT("f(xm,ym)  z="); PUT(zxmym); NEW_LINE;
    PUT("         xm="); PUT(xm); NEW_LINE;
    PUT("         ym="); PUT(ym); NEW_LINE;
    if zxmym < zbest then
      zbest:=zxmym;
      xbest:=xm;
      ybest:=ym;
      improve:=1;
      PUT("bad xmin ymin");
      goto continue;
    end if;

    zxmyp:=f(xm,yp);
    PUT("f(xm,yp)  z="); PUT(zxmyp); NEW_LINE;
    PUT("         xm="); PUT(xm); NEW_LINE;
    PUT("         yp="); PUT(yp); NEW_LINE;
    if zxmyp < zbest then
      zbest:=zxmyp;
      xbest:=xm;
      ybest:=yp;
      improve:=1;
      PUT("bad xmin ymax");
      goto continue;
    end if;

    zxpym:=f(xp,ym);
    PUT("f(xp,ym)  z="); PUT(zxpym); NEW_LINE;
    PUT("         xp="); PUT(xp); NEW_LINE;
    PUT("         ym="); PUT(ym); NEW_LINE;
    if zxpym < zbest then
      zbest:=zxpym;
      xbest:=xp;
      ybest:=ym;
      improve:=1;
      PUT("bad xmax ymin");
      goto continue;
    end if;

    zxpyp:=f(xp,yp);
    PUT("f(xp,yp)  z="); PUT(zxpyp); NEW_LINE;
    PUT("         xp="); PUT(xp); NEW_LINE;
    PUT("         yp="); PUT(yp); NEW_LINE;
    if zxpyp < zbest then
      zbest:=zxpyp;
      xbest:=xp;
      ybest:=yp;
      improve:=1;
      PUT("bad xmax ymax");
      goto continue;
    end if;

    -- contract xm
    hxm:=(xmin+x)/D2;
    zhxmym:=f(hxm,ym);
    PUT("f(hxm,ym) z="); PUT(zhxmym); NEW_LINE;
    PUT("        hxm="); PUT(hxm); NEW_LINE;
    PUT("         ym="); PUT(ym); NEW_LINE;
    if zhxmym < zbest then
      zbest:=zhxmym;
      xbest:=hxm;
      ybest:=ym;
      improve:=1;
    end if;
    zhxmy:=f(hxm,y);
    PUT("f(hxm,y)  z="); PUT(zhxmy); NEW_LINE;
    PUT("        hxm="); PUT(hxm); NEW_LINE;
    PUT("          y="); PUT(y); NEW_LINE;
    if zhxmy < zbest then
      zbest:=zhxmy;
      xbest:=hxm;
      ybest:=y;
      improve:=1;
    end if;
    zhxmyp:=f(hxm,yp);
    PUT("f(hxm,yp) z="); PUT(zhxmyp); NEW_LINE;
    PUT("        hxm="); PUT(hxm); NEW_LINE;
    PUT("         yp="); PUT(yp); NEW_LINE;
    if zhxmyp < zbest then
      zbest:=zhxmyp;
      xbest:=hxm;
      ybest:=yp;
      improve:=1;
    end if;
    if zhxmym<zxmym and then zhxmy<zxmy and then zhxmyp<zxmyp then
        xmin:=hxm;
        xm:=hxm;
        zxmym:=zhxmym;
        zxmy:=zhxmy;
        zxmyp:=zhxmyp;
        PUT("to xm  xmin="); PUT(xmin); NEW_LINE;
        PUT("best      z="); PUT(zbest); NEW_LINE;
        PUT("          x="); PUT(xbest); NEW_LINE;
        PUT("          y="); PUT(ybest); NEW_LINE;
        improve:=1;
    end if;

    -- contract to xp
    hxp:=(xp+x)/D2;
    zhxpym:=f(hxp,ym);
    PUT("f(hxp,ym) z="); PUT(zhxpym); NEW_LINE;
    PUT("        hxp="); PUT(hxp); NEW_LINE;
    PUT("         ym="); PUT(ym); NEW_LINE;
    if zhxpym < zbest then
      zbest:=zhxpym;
      xbest:=hxp;
      ybest:=ym;
      improve:=1;
    end if;
    zhxpy:=f(hxp,y);
    PUT("f(hxp,y)  z="); PUT(zhxpy); NEW_LINE;
    PUT("        hxp="); PUT(hxp); NEW_LINE;
    PUT("          y="); PUT(y); NEW_LINE;
    if zhxpy < zbest then
      zbest:=zhxpy;
      xbest:=hxp;
      ybest:=y;
      improve:=1;
    end if;
    zhxpyp:=f(hxp,yp);
    PUT("f(hxp,yp) z="); PUT(zhxpyp); NEW_LINE;
    PUT("        hxp="); PUT(hxp); NEW_LINE;
    PUT("         yp="); PUT(yp); NEW_LINE;
    if zhxpyp < zbest then
      zbest:=zhxpyp;
      xbest:=hxp;
      ybest:=yp;
      improve:=1;
    end if;
    if zhxpym<zxpym and then zhxpy<zxpy and then zhxpyp<zxpyp then
        xmax:=hxp;
        xp:=hxp;
        zxpym:=zhxpym;
        zxpy:=zhxpy;
        zxpyp:=zhxpyp;
        PUT("to xp  xmax="); PUT(xmax); NEW_LINE;
        PUT("best      z="); PUT(zbest); NEW_LINE;
        PUT("          x="); PUT(xbest); NEW_LINE;
        PUT("          y="); PUT(ybest); NEW_LINE;
        improve:=1;
    end if;

    -- contract ym
    hym:=(ymin+y)/D2;
    zxmhym:=f(xm,hym);
    PUT("f(xm,hym) z="); PUT(zxmhym); NEW_LINE;
    PUT("         xm="); PUT(xm); NEW_LINE;
    PUT("        hym="); PUT(hxm); NEW_LINE;
    if zxmhym < zbest then
      zbest:=zxmhym;
      xbest:=xm;
      ybest:=hym;
      improve:=1;
    end if;
    zxhym:=f(x,hym);
    PUT("f(x,hym)  z="); PUT(zxhym); NEW_LINE;
    PUT("          x="); PUT(x); NEW_LINE;
    PUT("        hym="); PUT(hxm); NEW_LINE;
    if zxhym < zbest then
      zbest:=zxhym;
      xbest:=x;
      ybest:=hym;
      improve:=1;
    end if;
    zxphym:=f(xp,hym);
    PUT("f(xp,hym) z="); PUT(zxphym); NEW_LINE;
    PUT("         xp="); PUT(xp); NEW_LINE;
    PUT("        hym="); PUT(hxm); NEW_LINE;
    if zxphym < zbest then
      zbest:=zxphym;
      xbest:=xp;
      ybest:=hym;
      improve:=1;
    end if;
    if zxmhym<zxmym and then zxhym<zxym and then zxphym<zxpym then
        ymin:=hym;
        ym:=hym;
        zxmym:=zxmhym;
        zxym:=zxhym;
        zxpym:=zxphym;
        PUT("to ym  ymin="); PUT(ymin); NEW_LINE;
        PUT("best      z="); PUT(zbest); NEW_LINE;
        PUT("          x="); PUT(xbest); NEW_LINE;
        PUT("          y="); PUT(ybest); NEW_LINE;
        improve:=1;
    end if;

    -- contract to yp
    hyp:=(yp+y)/D2;
    zxmhyp:=f(xm,hyp);
    PUT("f(xm,hyp) z="); PUT(zxmhyp); NEW_LINE;
    PUT("         xm="); PUT(xm); NEW_LINE;
    PUT("        hyp="); PUT(hxp); NEW_LINE;
    if zxmhyp < zbest then
      zbest:=zxmhyp;
      xbest:=xm;
      ybest:=hyp;
      improve:=1;
    end if;
    zxhyp:=f(x,hyp);
    PUT("f(x,hyp)  z="); PUT(zxhyp); NEW_LINE;
    PUT("          x="); PUT(x); NEW_LINE;
    PUT("        hyp="); PUT(hxp); NEW_LINE;
    if zxhyp < zbest then
      zbest:=zxhyp;
      xbest:=x;
      ybest:=hyp;
      improve:=1;
    end if;
    zxphyp:=f(xp,hyp);
    PUT("f(xp,hyp) z="); PUT(zxphyp); NEW_LINE;
    PUT("         xp="); PUT(xp); NEW_LINE;
    PUT("        hyp="); PUT(hxp); NEW_LINE;
    if zxphyp < zbest then
      zbest:=zxphyp;
      xbest:=xp;
      ybest:=hyp;
      improve:=1;
    end if;
    if zxmhyp<zxmyp and then zxhyp<zxpy and then zxphyp<Zxpyp then
        ymax:=hyp;
        yp:=hyp;
        zxmyp:=zxmhyp;
        zxyp:=zxhyp;
        zxpyp:=zxphyp;
        PUT("to yp  ymax="); PUT(ymax); NEW_LINE;
        PUT("best      z="); PUT(zbest); NEW_LINE;
        PUT("          x="); PUT(xbest); NEW_LINE;
        PUT("          y="); PUT(ybest); NEW_LINE;
        improve:=1;
    end if;

    <<continue>>
    zxy:=zbest; -- re-center
    x:=xbest;
    y:=ybest;
  end loop; -- on i
  NEW_LINE;
  PUT("best      z="); PUT(zbest); NEW_LINE;
  PUT("          x="); PUT(xbest); NEW_LINE;
  PUT("          y="); PUT(ybest); NEW_LINE;
  PUT("       xmax="); PUT(xmax); NEW_LINE;
  PUT("       xmin="); PUT(xmin); NEW_LINE;
  PUT("       ymax="); PUT(ymax); NEW_LINE;
  PUT("       ymin="); PUT(ymin); NEW_LINE;
  PUT_LINE("end of optm2");
end OPTM2 ;