optimize.tcl

Go to the documentation of this file.

/* ----------------------------------------------------------------------*/
/* */
/*  math/optimize.tcl --*/
/* */
/*  This file contains functions for optimization of a function*/
/*  or expression.*/
/* */
/*  Copyright (c) 2004, by Arjen Markus.*/
/*  Copyright (c) 2004, 2005 by Kevin B. Kenny.  All rights reserved.*/
/* */
/*  See the file "license.terms" for information on usage and redistribution*/
/*  of this file, and for a DISCLAIMER OF ALL WARRANTIES.*/
/* */
/*  RCS: @(#) $Id: optimize.tcl,v 1.11 2005/09/28 04:51:22 andreas_kupries Exp $*/
/* */
/* ----------------------------------------------------------------------*/

package require Tcl 8.4

/*  math::optimize --*/
/*     Namespace for the commands*/
/* */
namespace ::math::optimize {
   namespace export minimum  maximum solveLinearProgram linearProgramMaximum
   namespace export min_bound_1d min_unbound_1d

   /*  Possible extension: minimumExpr, maximumExpr*/
}

/*  minimum --*/
/*     Minimize a given function over a given interval*/
/* */
/*  Arguments:*/
/*     begin       Start of the interval*/
/*     end         End of the interval*/
/*     func        Name of the function to be minimized (takes one*/
/*                 argument)*/
/*     maxerr      Maximum relative error (defaults to 1.0e-4)*/
/*  Return value:*/
/*     Computed value for which the function is minimal*/
/*  Notes:*/
/*     The function needs not to be differentiable, but it is supposed*/
/*     to be continuous. There is no provision for sub-intervals where*/
/*     the function is constant (this might happen when the maximum*/
/*     error is very small, < 1.0e-15)*/
/* */
/*  Warning:*/
/*     This procedure is deprecated - use min_bound_1d instead*/
/* */
ret  ::math::optimize::minimum ( type begin , type end , type func , optional maxerr =1.0e-4 ) {

   set nosteps  [expr {3+int(-log($maxerr)/log(2.0))}]
   set delta    [expr {0.5*($end-$begin)*$maxerr}]

   for { set step 0 } { $step < $nosteps } { incr step } {
      set x1 [expr {($end+$begin)/2.0}]
      set x2 [expr {$x1+$delta}]

      set fx1 [uplevel 1 $func $x1]
      set fx2 [uplevel 1 $func $x2]

      if {$fx1 < $fx2} {
         set end   $x1
      } else {
         set begin $x1
      }
   }
   return $x1
}

/*  maximum --*/
/*     Maximize a given function over a given interval*/
/* */
/*  Arguments:*/
/*     begin       Start of the interval*/
/*     end         End of the interval*/
/*     func        Name of the function to be maximized (takes one*/
/*                 argument)*/
/*     maxerr      Maximum relative error (defaults to 1.0e-4)*/
/*  Return value:*/
/*     Computed value for which the function is maximal*/
/*  Notes:*/
/*     The function needs not to be differentiable, but it is supposed*/
/*     to be continuous. There is no provision for sub-intervals where*/
/*     the function is constant (this might happen when the maximum*/
/*     error is very small, < 1.0e-15)*/
/* */
/*  Warning:*/
/*     This procedure is deprecated - use max_bound_1d instead*/
/* */
ret  ::math::optimize::maximum ( type begin , type end , type func , optional maxerr =1.0e-4 ) {

   set nosteps  [expr {3+int(-log($maxerr)/log(2.0))}]
   set delta    [expr {0.5*($end-$begin)*$maxerr}]

   for { set step 0 } { $step < $nosteps } { incr step } {
      set x1 [expr {($end+$begin)/2.0}]
      set x2 [expr {$x1+$delta}]

      set fx1 [uplevel 1 $func $x1]
      set fx2 [uplevel 1 $func $x2]

      if {$fx1 > $fx2} {
         set end   $x1
      } else {
         set begin $x1
      }
   }
   return $x1
}

/* ----------------------------------------------------------------------*/
/* */
/*  min_bound_1d --*/
/* */
/*        Find a local minimum of a function between two given*/
/*        abscissae. Derivative of f is not required.*/
/* */
/*  Usage:*/
/*        min_bound_1d f x1 x2 ?-option value?,,,*/
/* */
/*  Parameters:*/
/*        f - Function to minimize.  Must be expressed as a Tcl*/
/*            command, to which will be appended the value at which*/
/*            to evaluate the function.*/
/*        x1 - Lower bound of the interval in which to search for a*/
/*             minimum*/
/*        x2 - Upper bound of the interval in which to search for a minimum*/
/* */
/*  Options:*/
/*        -relerror value*/
/*                Gives the tolerance desired for the returned*/
/*                abscissa.  Default is 1.0e-7.  Should never be less*/
/*                than the square root of the machine precision.*/
/*        -maxiter n*/
/*                Constrains minimize_bound_1d to evaluate the function*/
/*                no more than n times.  Default is 100.  If convergence*/
/*                is not achieved after the specified number of iterations,*/
/*                an error is thrown.*/
/*        -guess value*/
/*                Gives a point between x1 and x2 that is an initial guess*/
/*                for the minimum.  f(guess) must be at most f(x1) or*/
/*                f(x2).*/
/*         -fguess value*/
/*                 Gives the value of the ordinate at the value of '-guess'*/
/*                 if known.  Default is to evaluate the function*/
/*        -abserror value*/
/*                Gives the desired absolute error for the returned*/
/*                abscissa.  Default is 1.0e-10.*/
/*        -trace boolean*/
/*                A true value causes a trace to the standard output*/
/*                of the function evaluations. Default is 0.*/
/* */
/*  Results:*/
/*        Returns a two-element list comprising the abscissa at which*/
/*        the function reaches a local minimum within the interval,*/
/*        and the value of the function at that point.*/
/* */
/*  Side effects:*/
/*        Whatever side effects arise from evaluating the given function.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::optimize::min_bound_1d ( type f , type x1 , type x2 , type args ) {

    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]

    set phim1 0.6180339887498949
    set twomphi 0.3819660112501051

    array set params {
        -relerror 1.0e-7
        -abserror 1.0e-10
        -maxiter 100
        -trace 0
        -fguess {}
    }
    set params(-guess) [expr { $phim1 * $x1 + $twomphi * $x2 }]

    if { ( [llength $args] % 2 ) != 0 } {
        return -code error -errorcode [list min_bound_1d wrongNumArgs] \
            "wrong \# args, should be\
                 \"[lreplace [info level 0] 1 end f x1 x2 ?-option value?...]\""
    }
    foreach { key value } $args {
        if { ![info exists params($key)] } {
            return -code error -errorcode [list min_bound_1d badoption $key] \
                "unknown option \"$key\",\
                     should be -abserror,\
                     -fguess, -guess, -initial, -maxiter, -relerror,\
                     or -trace"
        }
    set params($key) $value
    }

    # a and b presumably bracket the minimum of the function.  Make sure
    # they're in ascending order.

    if { $x1 < $x2 } {
        set a $x1; set b $x2
    } else {
        set b $x1; set a $x2
    }

    set x $params(-guess);              # Best abscissa found so far
    set w $x;                           # Second best abscissa found so far
    set v $x;                           # Most recent earlier value of w

    set e 0.0;                          # Distance moved on the step before
                    # last.

    # Evaluate the function at the initial guess

    if { $params(-fguess) ne {} } {
        set fx $params(-fguess)
    } else {
        set s $f; lappend s $x; set fx [eval $s]
        if { $params(-trace) } {
            puts stdout "f($x) = $fx (initialisation)"
        }
    }
    set fw $fx
    set fv $fx

    for { set iter 0 } { $iter < $params(-maxiter) } { incr iter } {

        # Find the midpoint of the current interval

        set xm [expr { 0.5 * ( $a + $b ) }]

        # Compute the current tolerance for x, and twice its value

        set tol [expr { $params(-relerror) * abs($x) + $params(-abserror) }]
        set tol2 [expr { $tol + $tol }]
        if { abs( $x - $xm ) <= $tol2 - 0.5 * ($b - $a) } {
            return [list $x $fx]
        }
        set golden 1
        if { abs($e) > $tol } {

            # Use parabolic interpolation to find a minimum determined
            # by the evaluations at x, v, and w.  The size of the step
            # to take will be $p/$q.

            set r [expr { ( $x - $w ) * ( $fx - $fv ) }]
            set q [expr { ( $x - $v ) * ( $fx - $fw ) }]
            set p [expr { ( $x - $v ) * $q - ( $x - $w ) * $r }]
            set q [expr { 2. * ( $q - $r ) }]
            if { $q > 0 } {
                set p [expr { - $p }]
            } else {
                set q [expr { - $q }]
            }
            set olde $e
            set e $d

            # Test if parabolic interpolation results in less than half
            # the movement of the step two steps ago.

            if { abs($p) < abs( .5 * $q * $olde )
                 && $p > $q * ( $a - $x )
                 && $p < $q * ( $b - $x ) } {

                set d [expr { $p / $q }]
                set u [expr { $x + $d }]
                if { ( $u - $a ) < $tol2 || ( $b - $u ) < $tol2 } {
                    if { $xm-$x < 0 } {
                        set d [expr { - $tol }]
                    } else {
                        set d $tol
                    }
                }
                set golden 0
            }
        }

        # If parabolic interpolation didn't come up with an acceptable
        # result, use Golden Section instead.

        if { $golden } {
            if { $x >= $xm } {
                set e [expr { $a - $x }]
            } else {
                set e [expr { $b - $x }]
            }
            set d [expr { $twomphi * $e }]
        }

        # At this point, d is the size of the step to take.  Make sure
        # that it's at least $tol.

        if { abs($d) >= $tol } {
            set u [expr { $x + $d }]
        } elseif { $d < 0 } {
            set u [expr { $x - $tol }]
        } else {
            set u [expr { $x + $tol }]
        }

        # Evaluate the function

        set s $f; lappend s $u; set fu [eval $s]
        if { $params(-trace) } {
            if { $golden } {
                puts stdout "f($u)=$fu (golden section)"
            } else {
                puts stdout "f($u)=$fu (parabolic interpolation)"
            }
        }

        if { $fu <= $fx } {
            # We've the best abscissa so far.

            if { $u >= $x } {
                set a $x
            } else {
                set b $x
            }
            set v $w
            set fv $fw
            set w $x
            set fw $fx
            set x $u
            set fx $fu
        } else {

            if { $u < $x } {
                set a $u
            } else {
                set b $u
            }
            if { $fu <= $fw || $w == $x } {
                # We've the second-best abscissa so far
                set v $w
                set fv $fw
                set w $u
                set fw $fu
            } elseif { $fu <= $fv || $v == $x || $v == $w } {
                # We've the third-best so far
                set v $u
                set fv $fu
            }
        }
    }

    return -code error -errorcode [list min_bound_1d noconverge $iter] \
        "[lindex [info level 0] 0] failed to converge after $iter steps."

}

/* ----------------------------------------------------------------------*/
/* */
/*  brackmin --*/
/* */
/*        Find a place along the number line where a given function has*/
/*        a local minimum.*/
/* */
/*  Usage:*/
/*        brackmin f x1 x2 ?trace?*/
/* */
/*  Parameters:*/
/*        f - Function to minimize*/
/*        x1 - Abscissa thought to be near the minimum*/
/*        x2 - Additional abscissa thought to be near the minimum*/
/*  trace - Boolean variable that, if true,*/
/*                causes 'brackmin' to print a trace of its function*/
/*                evaluations to the standard output.  Default is 0.*/
/* */
/*  Results:*/
/*        Returns a three element list {x1 y1 x2 y2 x3 y3} where*/
/*        y1=f(x1), y2=f(x2), y3=f(x3).  x2 lies between x1 and x3, and*/
/*        y1>y2, y3>y2, proving that there is a local minimum somewhere*/
/*        in the interval (x1,x3).*/
/* */
/*  Side effects:*/
/*        Whatever effects the evaluation of f has.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::optimize::brackmin ( type f , type x1 , type x2 , optional trace =0 ) {

    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]

    set phi 1.6180339887498949
    set epsilon 1.0e-20
    set limit 50.

    # Choose a and b so that f(a) < f(b)

    set cmd $f; lappend cmd $x1; set fx1 [eval $cmd]
    if { $trace } {
        puts "f($x1) = $fx1 (initialisation)"
    }
    set cmd $f; lappend cmd $x2; set fx2 [eval $cmd]
    if { $trace } {
        puts "f($x2) = $fx2 (initialisation)"
    }
    if { $fx1 > $fx2 } {
        set a $x1; set fa $fx1
        set b $x2; set fb $fx2
    } else {
        set a $x2; set fa $fx2
        set b $x1; set fb $fx1
    }

    # Choose a c in the downhill direction

    set c [expr { $b + $phi * ($b - $a) }]
    set cmd $f; lappend cmd $c; set fc [eval $cmd]
    if { $trace } {
        puts "f($c) = $fc (initial dilatation by phi)"
    }

    while { $fb >= $fc } {

        # Try to do parabolic extrapolation to the minimum

        set r [expr { ($b - $a) * ($fb - $fc) }]
        set q [expr { ($b - $c) * ($fb - $fa) }]
        if { abs( $q - $r ) > $epsilon } {
            set denom [expr { $q - $r }]
        } elseif { $q > $r } {
            set denom $epsilon
        } else {
            set denom -$epsilon
        }
        set u [expr { $b - ( (($b - $c) * $q - ($b - $a) * $r)
                             / (2. * $denom) ) }]
        set ulimit [expr { $b + $limit * ( $c - $b ) }]

        # Test the extrapolated abscissa

        if { ($b - $u) * ($u - $c) > 0 } {

            # u lies between b and c.  Try to interpolate

            set cmd $f; lappend cmd $u; set fu [eval $cmd]
            if { $trace } {
                puts "f($u) = $fu (parabolic interpolation)"
            }

            if { $fu < $fc } {

                # fb > fu and fc > fu, so there is a minimum between b and c
                # with u as a starting guess.

                return [list $b $fb $u $fu $c $fc]

            }

            if { $fu > $fb } {

                # fb < fu, fb < fa, and u cannot lie between a and b
                # (because it lies between a and c).  There is a minimum
                # somewhere between a and u, with b a starting guess.

                return [list $a $fa $b $fb $u $fu]

            }

            # Parabolic interpolation was useless. Expand the
            # distance by a factor of phi and try again.

            set u [expr { $c + $phi * ($c - $b) }]
            set cmd $f; lappend cmd $u; set fu [eval $cmd]
            if { $trace } {
                puts "f($u) = $fu (parabolic interpolation failed)"
            }


        } elseif { ( $c - $u ) * ( $u - $ulimit ) > 0 } {

            # u lies between $c and $ulimit.

            set cmd $f; lappend cmd $u; set fu [eval $cmd]
            if { $trace } {
                puts "f($u) = $fu (parabolic extrapolation)"
            }

            if { $fu > $fc } {

                # minimum lies between b and u, with c an initial guess.

                return [list $b $fb $c $fc $u $fu]

            }

            # function is still decreasing fa > fb > fc > fu. Take
            # another factor-of-phi step.

            set b $c; set fb $fc
            set c $u; set fc $fu
            set u [expr { $c + $phi * ( $c - $b ) }]
            set cmd $f; lappend cmd $u; set fu [eval $cmd]
            if { $trace } {
                puts "f($u) = $fu (parabolic extrapolation ok)"
            }

        } elseif { ($u - $ulimit) * ( $ulimit - $c ) >= 0 } {

            # u went past ulimit.  Pull in to ulimit and evaluate there.

            set u $ulimit
            set cmd $f; lappend cmd $u; set fu [eval $cmd]
            if { $trace } {
                puts "f($u) = $fu (limited step)"
            }

        } else {

            # parabolic extrapolation gave a useless value.

            set u [expr { $c + $phi * ( $c - $b ) }]
            set cmd $f; lappend cmd $u; set fu [eval $cmd]
            if { $trace } {
                puts "f($u) = $fu (parabolic extrapolation failed)"
            }

        }

        set a $b; set fa $fb
        set b $c; set fb $fc
        set c $u; set fc $fu
    }

    return [list $a $fa $b $fb $c $fc]
}

/* ----------------------------------------------------------------------*/
/* */
/*  min_unbound_1d --*/
/* */
/*  Minimize a function of one variable, unconstrained, derivatives*/
/*  not required.*/
/* */
/*  Usage:*/
/*        min_bound_1d f x1 x2 ?-option value?,,,*/
/* */
/*  Parameters:*/
/*        f - Function to minimize.  Must be expressed as a Tcl*/
/*            command, to which will be appended the value at which*/
/*            to evaluate the function.*/
/*        x1 - Initial guess at the minimum*/
/*        x2 - Second initial guess at the minimum, used to set the*/
/*       initial length scale for the search.*/
/* */
/*  Options:*/
/*        -relerror value*/
/*                Gives the tolerance desired for the returned*/
/*                abscissa.  Default is 1.0e-7.  Should never be less*/
/*                than the square root of the machine precision.*/
/*        -maxiter n*/
/*                Constrains min_bound_1d to evaluate the function*/
/*                no more than n times.  Default is 100.  If convergence*/
/*                is not achieved after the specified number of iterations,*/
/*                an error is thrown.*/
/*        -abserror value*/
/*                Gives the desired absolute error for the returned*/
/*                abscissa.  Default is 1.0e-10.*/
/*        -trace boolean*/
/*                A true value causes a trace to the standard output*/
/*                of the function evaluations. Default is 0.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::optimize::min_unbound_1d ( type f , type x1 , type x2 , type args ) {

    set f [lreplace $f 0 0 [uplevel 1 [list namespace which [lindex $f 0]]]]

    array set params {
    -relerror 1.0e-7
    -abserror 1.0e-10
    -maxiter 100
        -trace 0
    }
    if { ( [llength $args] % 2 ) != 0 } {
        return -code error -errorcode [list min_unbound_1d wrongNumArgs] \
            "wrong \# args, should be\
                 \"[lreplace [info level 0] 1 end \
                         f x1 x2 ?-option value?...]\""
    }
    foreach { key value } $args {
        if { ![info exists params($key)] } {
            return -code error -errorcode [list min_unbound_1d badoption $key] \
                "unknown option \"$key\",\
                     should be -trace"
        }
        set params($key) $value
    }
    foreach { a fa b fb c fc } [brackmin $f $x1 $x2 $params(-trace)] {
    break
    }
    return [eval [linsert [array get params] 0 \
              min_bound_1d $f $a $c -guess $b -fguess $fb]]
}

/* ----------------------------------------------------------------------*/
/* */
/*  nelderMead --*/
/* */
/*  Attempt to minimize/maximize a function using the downhill*/
/*  simplex method of Nelder and Mead.*/
/* */
/*  Usage:*/
/*  nelderMead f x ?-keyword value?*/
/* */
/*  Parameters:*/
/*  f - The function to minimize.  The function must be an incomplete*/
/*      Tcl command, to which will be appended N parameters.*/
/*  x - The starting guess for the minimum; a vector of N parameters*/
/*      to be passed to the function f.*/
/* */
/*  Options:*/
/*  -scale xscale*/
/*      Initial guess as to the problem scale.  If '-scale' is*/
/*      supplied, then the parameters will be varied by the*/
/*          specified amounts.  The '-scale' parameter must of the*/
/*      same dimension as the 'x' vector, and all elements must*/
/*      be nonzero.  Default is 0.0001 times the 'x' vector,*/
/*      or 0.0001 for zero elements in the 'x' vector.*/
/* */
/*  -ftol epsilon*/
/*      Requested tolerance in the function value; nelderMead*/
/*      returns if N+1 consecutive iterates all differ by less*/
/*      than the -ftol value.  Default is 1.0e-7*/
/* */
/*  -maxiter N*/
/*      Maximum number of iterations to attempt.  Default is*/
/*      500.*/
/* */
/*  -trace flag*/
/*      If '-trace 1' is supplied, nelderMead writes a record*/
/*      of function evaluations to the standard output as it*/
/*      goes.  Default is 0.*/
/* */
/* ----------------------------------------------------------------------*/

ret  ::math::optimize::nelderMead ( type f , type startx , type args ) {
    array set params {
    -ftol 1.e-7
    -maxiter 500
    -scale {}
    -trace 0
    }

    # Check arguments

    if { ( [llength $args] % 2 ) != 0 } {
        return -code error -errorcode [list nelderMead wrongNumArgs] \
            "wrong \# args, should be\
                 \"[lreplace [info level 0] 1 end \
                         f x1 x2 ?-option value?...]\""
    }
    foreach { key value } $args {
        if { ![info exists params($key)] } {
            return -code error -errorcode [list nelderMead badoption $key] \
                "unknown option \"$key\",\
                     should be -ftol, -maxiter, -scale or -trace"
        }
        set params($key) $value
    }
    
    # Construct the initial simplex

    set vertices [list $startx]
    if { [llength $params(-scale)] == 0 } {
    set i 0
    foreach x0 $startx {
        if { $x0 == 0 } {
        set x1 0.0001
        } else {
        set x1 [expr {1.0001 * $x0}]
        }
        lappend vertices [lreplace $startx $i $i $x1]
        incr i
    }
    } elseif { [llength $params(-scale)] != [llength $startx] } {
    return -code error -errorcode [list nelderMead badOption -scale] \
        "-scale vector must be of same size as starting x vector"
    } else {
    set i 0
    foreach x0 $startx s $params(-scale) {
        lappend vertices [lreplace $startx $i $i [expr { $x0 + $s }]]
        incr i
    }
    }

    # Evaluate at the initial points

    set n [llength $startx]
    foreach x $vertices {
    set cmd $f
    foreach xx $x {
        lappend cmd $xx
    }
    set y [uplevel 1 $cmd]
    if {$params(-trace)} {
        puts "nelderMead: evaluating initial point: x=[list $x] y=$y"
    }
    lappend yvec $y
    }


    # Loop adjusting the simplex in the 'vertices' array.

    set nIter 0
    while { 1 } {

    # Find the highest, next highest, and lowest value in y,
    # and save the indices.

    set iBot 0
    set yBot [lindex $yvec 0]
    set iTop -1
    set yTop [lindex $yvec 0]
    set iNext -1
    set i 0
    foreach y $yvec {
        if { $y <= $yBot } {
        set yBot $y
        set iBot $i
        }
        if { $iTop < 0 || $y >= $yTop } {
        set iNext $iTop
        set yNext $yTop
        set iTop $i
        set yTop $y
        } elseif { $iNext < 0 || $y >= $yNext } {
        set iNext $i
        set yNext $y
        }
        incr i
    }

    # Return if the relative error is within an acceptable range

    set rerror [expr { 2. * abs( $yTop - $yBot )
               / ( abs( $yTop ) + abs( $yBot ) ) }]
    if { $rerror < $params(-ftol) } {
        set status ok
        break
    }

    # Count iterations

    if { [incr nIter] > $params(-maxiter) } {
        set status too-many-iterations
        break
    }
    incr nIter

    # Find the centroid of the face opposite the vertex that
    # maximizes the function value.

    set centroid {}
    for { set i 0 } { $i < $n } { incr i } {
        lappend centroid 0.0
    }
    set i 0
    foreach v $vertices {
        if { $i != $iTop } {
        set newCentroid {}
        foreach x0 $centroid x1 $v {
            lappend newCentroid [expr { $x0 + $x1 }]
        }
        set centroid $newCentroid
        }
        incr i
    }
    set newCentroid {}
    foreach x $centroid {
        lappend newCentroid [expr { $x / $n }]
    }
    set centroid $newCentroid

    # The first trial point is a reflection of the high point
    # around the centroid

    set trial {}
    foreach x0 [lindex $vertices $iTop] x1 $centroid {
        lappend trial [expr {$x1 + ($x1 - $x0)}]
    }
    set cmd $f
    foreach xx $trial {
        lappend cmd $xx
    } 
    set yTrial [uplevel 1 $cmd]
    if { $params(-trace) } {
        puts "nelderMead: trying reflection: x=[list $trial] y=$yTrial"
    }

    # If that reflection yields a new minimum, replace the high point,
    # and additionally try dilating in the same direction.

    if { $yTrial < $yBot } {
        set trial2 {}
        foreach x0 $centroid x1 $trial {
        lappend trial2 [expr { $x1 + ($x1 - $x0) }]
        }
        set cmd $f
        foreach xx $trial2 {
        lappend cmd $xx
        }
        set yTrial2 [uplevel 1 $cmd]
        if { $params(-trace) } {
        puts "nelderMead: trying dilated reflection:\
                      x=[list $trial2] y=$y"
        }
        if { $yTrial2 < $yBot } {

        # Additional dilation yields a new minimum

        lset vertices $iTop $trial2
        lset yvec $iTop $yTrial2
        } else {

        # Additional dilation failed, but we can still use
        # the first trial point.

        lset vertices $iTop $trial
        lset yvec $iTop $yTrial

        }
    } elseif { $yTrial < $yNext } {

        # The reflected point isn't a new minimum, but it's
        # better than the second-highest.  Replace the old high
        # point and try again.

        lset vertices $iTop $trial
        lset yvec $iTop $yTrial

    } else {

        # The reflected point is worse than the second-highest point.
        # If it's better than the highest, keep it... but in any case,
        # we want to try contracting the simplex, because a further
        # reflection will simply bring us back to the starting point.

        if { $yTrial < $yTop } {
        lset vertices $iTop $trial
        lset yvec $iTop $yTrial
        set yTop $yTrial
        }
        set trial {}
        foreach x0 [lindex $vertices $iTop] x1 $centroid {
        lappend trial [expr { ( $x0 + $x1 ) / 2. }]
        }
        set cmd $f
        foreach xx $trial {
        lappend cmd $xx
        }
        set yTrial [uplevel 1 $cmd]
        if { $params(-trace) } {
        puts "nelderMead: contracting from high point:\
                      x=[list $trial] y=$y"
        }
        if { $yTrial < $yTop } {

        # Contraction gave an improvement, so continue with
        # the smaller simplex

        lset vertices $iTop $trial
        lset yvec $iTop $yTrial

        } else {

        # Contraction gave no improvement either; we seem to
        # be in a valley of peculiar topology.  Contract the
        # simplex about the low point and try again.

        set newVertices {}
        set newYvec {}
        set i 0
        foreach v $vertices y $yvec {
            if { $i == $iBot } {
            lappend newVertices $v
            lappend newYvec $y
            } else {
            set newv {}
            foreach x0 $v x1 [lindex $vertices $iBot] {
                lappend newv [expr { ($x0 + $x1) / 2. }]
            }
            lappend newVertices $newv
            set cmd $f
            foreach xx $newv {
                lappend cmd $xx
            }
            lappend newYvec [uplevel 1 $cmd]
            if { $params(-trace) } {
                puts "nelderMead: contracting about low point:\
                                  x=[list $newv] y=$y"
            }
            }
            incr i
        }
        set vertices $newVertices
        set yvec $newYvec
        }

    }

    }
    return [list y $yBot x [lindex $vertices $iBot] vertices $vertices yvec $yvec nIter $nIter status $status]

}

/*  solveLinearProgram*/
/*     Solve a linear program in standard form*/
/* */
/*  Arguments:*/
/*     objective     Vector defining the objective function*/
/*     constraints   Matrix of constraints (as a list of lists)*/
/* */
/*  Return value:*/
/*     Computed values for the coordinates or "unbounded" or "infeasible"*/
/* */
ret  ::math::optimize::solveLinearProgram ( type objective , type constraints ) {
    #
    # Check the arguments first and then put them in a more convenient
    # form
    #

    foreach {nconst nvars matrix} \
        [SimplexPrepareMatrix $objective $constraints] {break}

    set solution [SimplexSolve $nconst nvars $matrix]

    if { [llength $solution] > 1 } {
        return [lrange $solution 0 [expr {$nvars-1}]]
    } else {
        return $solution
    }
}

/*  linearProgramMaximum --*/
/*     Compute the value attained at the optimum*/
/* */
/*  Arguments:*/
/*     objective     The coefficients of the objective function*/
/*     result        The coordinate values as obtained by solving the program*/
/* */
/*  Return value:*/
/*     Value at the maximum point*/
/* */
ret  ::math::optimize::linearProgramMaximum (type objective , type result) {

    set value    0.0

    foreach coeff $objective coord $result {
        set value [expr {$value+$coeff*$coord}]
    }

    return $value
}

/*  SimplexPrintMatrix*/
/*     Debugging routine: print the matrix in easy to read form*/
/* */
/*  Arguments:*/
/*     matrix        Matrix to be printed*/
/* */
/*  Return value:*/
/*     None*/
/* */
/*  Note:*/
/*     The tableau should be transposed ...*/
/* */
ret  ::math::optimize::SimplexPrintMatrix (type matrix) {
    puts "\nBasis:\t[join [lindex $matrix 0] \t]"
    foreach col [lrange $matrix 1 end] {
        puts "      \t[join $col \t]"
    }
}

/*  SimplexPrepareMatrix*/
/*     Prepare the standard tableau from all program data*/
/* */
/*  Arguments:*/
/*     objective     Vector defining the objective function*/
/*     constraints   Matrix of constraints (as a list of lists)*/
/* */
/*  Return value:*/
/*     List of values as a standard tableau and two values*/
/*     for the sizes*/
/* */
ret  ::math::optimize::SimplexPrepareMatrix (type objective , type constraints) {

    #
    # Check the arguments first
    #
    set nconst [llength $constraints]
    set ncols {}
    foreach row $constraints {
        if { $ncols == {} } {
            set ncols [llength $row]
        } else {
            if { $ncols != [llength $row] } {
                return -code error -errorcode ARGS "Incorrectly formed constraints matrix"
            }
        }
    }

    set nvars [expr {$ncols-1}]

    if { [llength $objective] != $nvars } {
        return -code error -errorcode ARGS "Incorrect length for objective vector"
    }

    #
    # Set up the tableau:
    # Easiest manipulations if we store the columns first
    # So:
    # - First column is the list of variable indices in the basis
    # - Second column is the list of maximum values
    # - "nvars" columns that follow: the coefficients for the actual
    #   variables
    # - last "nconst" columns: the slack variables
    #
    set matrix   [list]
    set lastrow  [concat $objective [list 0.0]]

    set newcol   [list]
    for {set idx 0} {$idx < $nconst} {incr idx} {
        lappend newcol [expr {$nvars+$idx}]
    }
    lappend newcol "?"
    lappend matrix $newcol

    set zvector [list]
    foreach row $constraints {
        lappend zvector [lindex $row end]
    }
    lappend zvector 0.0
    lappend matrix $zvector

    for {set idx 0} {$idx < $nvars} {incr idx} {
        set newcol [list]
        foreach row $constraints {
            lappend newcol [expr {double([lindex $row $idx])}]
        }
        lappend newcol [expr {-double([lindex $lastrow $idx])}]
         lappend matrix $newcol
    }

    #
    # Add the columns for the slack variables
    #
    set zeros {}
    for {set idx 0} {$idx <= $nconst} {incr idx} {
        lappend zeros 0.0
    }
    for {set idx 0} {$idx < $nconst} {incr idx} {
        lappend matrix [lreplace $zeros $idx $idx 1.0]
    }

    return [list $nconst $nvars $matrix]
}

/*  SimplexSolve --*/
/*     Solve the given linear program using the simplex method*/
/* */
/*  Arguments:*/
/*     nconst        Number of constraints*/
/*     nvars         Number of actual variables*/
/*     tableau       Standard tableau (as a list of columns)*/
/* */
/*  Return value:*/
/*     List of values for the actual variables*/
/* */
ret  ::math::optimize::SimplexSolve (type nconst , type nvars , type tableau) {
    set end 0
    while { !$end } {

        #
        # Find the new variable to put in the basis
        #
        set nextcol [SimplexFindNextColumn $tableau]
        if { $nextcol == -1 } {
            set end 1
            continue
        }

        #
        # Now determine which one should leave
        # TODO: is a lack of a proper row indeed an
        #       indication of the infeasibility?
        #
        set nextrow [SimplexFindNextRow $tableau $nextcol]
        if { $nextrow == -1 } {
            return "infeasible"
        }

        #
        # Make the vector for sweeping through the tableau
        #
        set vector [SimplexMakeVector $tableau $nextcol $nextrow]

        #
        # Sweep through the tableau
        #
        set tableau [SimplexNewTableau $tableau $nextcol $nextrow $vector]
    }

    #
    # Now we can return the result
    #
    SimplexResult $tableau
}

/*  SimplexResult --*/
/*     Reconstruct the result vector*/
/* */
/*  Arguments:*/
/*     tableau       Standard tableau (as a list of columns)*/
/* */
/*  Return value:*/
/*     Vector of values representing the maximum point*/
/* */
ret  ::math::optimize::SimplexResult (type tableau) {
    set result {}

    set firstcol  [lindex $tableau 0]
    set secondcol [lindex $tableau 1]
    set result    {}

    set nvars     [expr {[llength $tableau]-2}]
    for {set i 0} {$i < $nvars } { incr i } {
        lappend result 0.0
    }

    set idx 0
    foreach col [lrange $firstcol 0 end-1] {
        set result [lreplace $result $col $col [lindex $secondcol $idx]]
        incr idx
    }

    return $result
}

/*  SimplexFindNextColumn --*/
/*     Find the next column - the one with the largest negative*/
/*     coefficient*/
/* */
/*  Arguments:*/
/*     tableau       Standard tableau (as a list of columns)*/
/* */
/*  Return value:*/
/*     Index of the column*/
/* */
ret  ::math::optimize::SimplexFindNextColumn (type tableau) {
    set idx        0
    set minidx    -1
    set mincoeff   0.0

    foreach col [lrange $tableau 2 end] {
        set coeff [lindex $col end]
        if { $coeff < 0.0 } {
            if { $coeff < $mincoeff } {
                set minidx $idx
               set mincoeff $coeff
            }
        }
        incr idx
    }

    return $minidx
}

/*  SimplexFindNextRow --*/
/*     Find the next row - the one with the largest negative*/
/*     coefficient*/
/* */
/*  Arguments:*/
/*     tableau       Standard tableau (as a list of columns)*/
/*     nextcol       Index of the variable that will replace this one*/
/* */
/*  Return value:*/
/*     Index of the row*/
/* */
ret  ::math::optimize::SimplexFindNextRow (type tableau , type nextcol) {
    set idx        0
    set minidx    -1
    set mincoeff   {}

    set bvalues [lrange [lindex $tableau 1] 0 end-1]
    set yvalues [lrange [lindex $tableau [expr {2+$nextcol}]] 0 end-1]

    foreach rowcoeff $bvalues divcoeff $yvalues {
        if { $divcoeff > 0.0 } {
            set coeff [expr {$rowcoeff/$divcoeff}]

            if { $mincoeff == {} || $coeff < $mincoeff } {
                set minidx $idx
                set mincoeff $coeff
            }
        }
        incr idx
    }

    return $minidx
}

/*  SimplexMakeVector --*/
/*     Make the "sweep" vector*/
/* */
/*  Arguments:*/
/*     tableau       Standard tableau (as a list of columns)*/
/*     nextcol       Index of the variable that will replace this one*/
/*     nextrow       Index of the variable in the base that will be replaced*/
/* */
/*  Return value:*/
/*     Vector to be used to update the coefficients of the tableau*/
/* */
ret  ::math::optimize::SimplexMakeVector (type tableau , type nextcol , type nextrow) {

    set idx      0
    set vector   {}
    set column   [lindex $tableau [expr {2+$nextcol}]]
    set divcoeff [lindex $column $nextrow]

    foreach colcoeff $column {
        if { $idx != $nextrow } {
            set coeff [expr {-$colcoeff/$divcoeff}]
        } else {
            set coeff [expr {1.0/$divcoeff-1.0}]
        }
        lappend vector $coeff
        incr idx
    }

    return $vector
}

/*  SimplexNewTableau --*/
/*     Sweep through the tableau and create the new one*/
/* */
/*  Arguments:*/
/*     tableau       Standard tableau (as a list of columns)*/
/*     nextcol       Index of the variable that will replace this one*/
/*     nextrow       Index of the variable in the base that will be replaced*/
/*     vector        Vector to sweep with*/
/* */
/*  Return value:*/
/*     New tableau*/
/* */
ret  ::math::optimize::SimplexNewTableau (type tableau , type nextcol , type nextrow , type vector) {

    #
    # The first column: replace the nextrow-th element
    # The second column: replace the value at the nextrow-th element
    # For all the others: the same receipe
    #
    set firstcol   [lreplace [lindex $tableau 0] $nextrow $nextrow $nextcol]
    set newtableau [list $firstcol]

    #
    # The rest of the matrix
    #
    foreach column [lrange $tableau 1 end] {
        set yval   [lindex $column $nextrow]
        set newcol {}
        foreach c $column vcoeff $vector {
            set newval [expr {$c+$yval*$vcoeff}]
            lappend newcol $newval
        }
        lappend newtableau $newcol
    }

    return $newtableau
}

/*  Now we can announce our presence*/
package provide math::optimize 1.0

if { ![info exists ::argv0] || [string compare $::argv0 [info script]] } {
    return
}

namespace import math::optimize::min_bound_1d
namespace import math::optimize::maximum
namespace import math::optimize::nelderMead

ret  f (type x , type y) {
    set xx [expr { $x - 3.1415926535897932 / 2. }]
    set v1 [expr { 0.3 * exp( -$xx*$xx / 2. ) }]
    set d [expr { 10. * $y - sin(9. * $x) }]
    set v2 [expr { exp(-10.*$d*$d)}]
    set rv [expr { -$v1 - $v2 }]
    return $rv
}

ret  g (type a , type b) {
    set x1 [expr {0.1 - $a + $b}]
    set x2 [expr {$a + $b - 1.}]
    set x3 [expr {3.-8.*$a+8.*$a*$a-8.*$b+8.*$b*$b}]
    set x4 [expr {$a/10. + $b/10. + $x1*$x1/3. + $x2*$x2 - $x2 * exp(1-$x3*$x3)}]
    return $x4
}

if { ![package vsatisfies [package provide Tcl] 8.5] } {
     tcl = _precision 17
}
puts "f"
puts [math::optimize::nelderMead f {1. 0.} -scale {0.1 0.01} -trace 1]
puts "g"
puts [math::optimize::nelderMead g {0. 0.} -scale {1. 1.} -trace 1]