new; 
@output file = optdemo.out reset; format /rd 7,4; output on; @
"Demonstration of Gauss Optimization Library";"";

library optmum;             @   Load the Optimization Library   @
c0 = 1;                     @   Initial values                  @
a = {0.5,2};                @   Fixed parameters                @

{c, val, gradi, retcode} = optmum(&tstr,c0);
val = -val; gradi = -gradi;         @   Note change in sign     @
gosub prit; "";

@   c = argmin (-argmax)    @
@   val = minimum (-max)    @
@   gradi = gradient at c   @
@   retcode = tells if c satisfies Gauss's criteria (1 = OK)    @
@   tstr = name of function to minimize:  "&" denotes functions @
@   c0 = initial values     @

@   Now optimize over constrained set: cmin < c < cmax.         @

cmin = 1e-3; cmax = 10; c0 = 0;
{c, val, gradi, retcode} = optmum(&tstr2,c0);

@   Find transformed parameter and associated gradient.         @

c0 = trans(c0,cmin,cmax); c = trans(c,cmin,cmax);
val = -val;
gradi = -gradp(&tstr,c);            @   Note change in sign     @
hss = -hessp(&tstr,c);              @   2nd Order Condition!    @
gosub prit; "Hessian =       ";; hss;

output off; end;

PRIT:
    format /rd 2,1;
    "Function = ";; a[2,1];;"c^0.5 - ";; a[1,1];;"c.";
    format /rd 7,4;
    "Initial Value = ";; c0;
    "Argmax =        ";; c;
    "Maximum =       ";; val;
    "Gradient =      ";; gradi;
return;


@   Optmum requires a procedure with one argument--a vector of      @
@   parameters to be minimized over--and one return--a value        @

proc(1) = tstr(c);          @   c is an argument, not the same as   @
                            @   c in the main program.              @
    local val;              @   val defined locally, not the same   @
                            @   as val in main program.             @

    val = a[2,1]*(c^.5) - a[1,1]*c;

@   Note that a is defined in the main program and not explicitly   @
@   passed in.  This is VERY bad style, but necessary here.  Always @
@   identify variables that are handled like this.                  @

    val = -val;             @   Optmum minimizes.  We want to max.  @
retp(val); endp;

proc(1) = tstr2(c);
    local val, c2;
    c2 = trans(c,cmin,cmax);            @   cmin < c2 < cmax        @
                                        @   cmin, cmax are globals  @
    val = a[2,1]*(c2^.5) - a[1,1]*c2;   @   a is a global           @
    val = -val;
retp(val); endp;

proc(1) = trans(c,cmin,cmax);
    local c2;
    c2 = cmin + (cmax-cmin)*exp(c)/(1+exp(c));
retp(c2); endp;