# -*- tcl -*-
# Tests for 1-d optimisation functions in math library  -*- tcl -*-
#
# This file contains a collection of tests for one or more of the Tcllib
# procedures.  Sourcing this file into Tcl runs the tests and
# generates output for errors.  No output means no errors were found.
#
# $Id: optimize.test,v 1.16 2006/10/09 21:41:41 andreas_kupries Exp $
#
# Copyright (c) 2004 by Arjen Markus
# Copyright (c) 2004, 2005 by Kevin B. Kenny
# All rights reserved.
#
# Note:
#    By evaluating the tests in a different namespace than global,
#    we assure that the namespace issue (Bug #...) is checked.
#

# -------------------------------------------------------------------------

source [file join \
	[file dirname [file dirname [file join [pwd] [info script]]]] \
	devtools testutilities.tcl]

testsNeedTcl     8.4
testsNeedTcltest 2.1

support {
    useLocal math.tcl math
}
testing {
    useLocal optimize.tcl math::optimize
}

# -------------------------------------------------------------------------

namespace eval optimizetest {

namespace import ::math::optimize::*

if {![package vsatisfies [package present Tcl] 8.5]} {
    set old_precision $::tcl_precision
    set ::tcl_precision 17
}

#
# Simple test functions
#
proc const_func { x } {
   return 1.0
}
proc ffunc { x } {
   expr {$x*(1.0-$x*$x)}
}
proc minfunc { x } {
   expr {-$x*(1.0-$x*$x)}
}
proc absfunc { x } {
   expr {abs($x*(1.0-$x*$x))}
}

proc within_range { result min max } {
   #puts "Within range? $result $min $max"
   #puts "[expr {2.0*abs($result-$min)/abs($max+$min)}]"
   if { $result >= $min && $result <= $max } {
      set ok 1
   } else {
      set ok 0
   }
   return $ok
}

#
# Test the minimum procedure
#
# Note about the uneven and even functions:
# the initial interval is chosen symmetrical, so that the
# three function values are equal.
#
test optimize-1.1 "Minimum of constant function" {
   set result [minimum -1.0 1.0 ::optimizetest::const_func]
   within_range $result -1.0 1.0
} 1

test optimize-1.2 "Minimum of odd function, case 1" {
   set result [minimum -1.0 1.0 ::optimizetest::ffunc]
   set xmin   [expr {-sqrt(1.0/3.0)-0.0001}]
   set xmax   [expr {-sqrt(1.0/3.0)+0.0001}]
   within_range $result $xmin $xmax
} 1

test optimize-1.3 "Minimum of odd function, asymmetric interval" {
   set result [minimum -0.8 1.2 ::optimizetest::ffunc]
   set xmin   [expr {-sqrt(1.0/3.0)-0.0001}]
   set xmax   [expr {-sqrt(1.0/3.0)+0.0001}]
   within_range $result $xmin $xmax
} 1

test optimize-1.4 "Minimum of odd function, case 2" {
   set result [minimum -1.0 1.0 ::optimizetest::minfunc]
   set xmin   [expr {sqrt(1.0/3.0)-0.0001}]
   set xmax   [expr {sqrt(1.0/3.0)+0.0001}]
   within_range $result $xmin $xmax
} 1

test optimize-1.5 "Minimum of even function" {
   set result [minimum -1.0 1.0 ::optimizetest::absfunc]
   set xmin   -0.0001
   set xmax    0.0001
   within_range $result $xmin $xmax
} 1

#
# Test the maximum procedure
#
# Note about the uneven and even functions:
# the initial interval is chosen symmetrical, so that the
# three function values are equal.
#
test optimize-2.1 "Maximum of constant function" {
   set result [maximum -1.0 1.0 ::optimizetest::const_func]
   within_range $result -1.0 1.0
} 1

test optimize-2.2 "Maximum of odd function, case 1" {
   set result [maximum -1.0 1.0 ::optimizetest::ffunc]
   set xmin   [expr {sqrt(1.0/3.0)-0.0001}]
   set xmax   [expr {sqrt(1.0/3.0)+0.0001}]
   within_range $result $xmin $xmax
} 1

test optimize-2.3 "Maximum of odd function, case 2" {
   set result [maximum -1.0 1.0 ::optimizetest::minfunc]
   set xmin   [expr {-sqrt(1.0/3.0)-0.0001}]
   set xmax   [expr {-sqrt(1.0/3.0)+0.0001}]
   within_range $result $xmin $xmax
} 1

#
# Either of the two maxima will do
#
test optimize-2.4 "Maximum of even function" {
   set result [maximum -1.0 1.0 ::optimizetest::absfunc]
   set xmin   [expr {-sqrt(1.0/3.0)-0.0001}]
   set xmax   [expr {-sqrt(1.0/3.0)+0.0001}]
   set ok     [within_range $result $xmin $xmax]
   set xmin   [expr {sqrt(1.0/3.0)-0.0001}]
   set xmax   [expr {sqrt(1.0/3.0)+0.0001}]
   incr ok    [within_range $result $xmin $xmax]
} 1


# Custom match procedure for approximate results

proc withinEpsilon { shouldBe is } {
    expr { [string is double $is]
	   && abs( $is - $shouldBe ) < 1.e-07 * abs($shouldBe) }
}

::tcltest::customMatch withinEpsilon [namespace code withinEpsilon]

test linmin-1.1 {find minimum of a parabola - constrained} \
    -setup {
	proc f x { expr { ($x + 3.) * ($x - 1.) } }
    } \
    -body {
	foreach {x y} [min_bound_1d f 10. -10.] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result -1. \
    -match withinEpsilon

test linmin-1.2 {find minimum of cosine} \
    -setup {
	proc f x { expr { cos($x) } }
    } \
    -body {
	foreach { x y } [min_bound_1d f 0. 6.28318] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 3.1415926535897932 \
    -match withinEpsilon

test linmin-1.3 {find minimum of a bell-shaped function} \
    -setup {
	proc f x {
	    set t [expr { $x - 3. }]
	    return [expr { -exp ( -$t * $t / 2 ) }]
	}
    } \
    -body {
	foreach { x y } [min_bound_1d f 0 30.] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 3. \
    -match withinEpsilon

test linmin-1.4 {function where parabolic extrapolation never works} \
    -setup {
	proc f x { expr { -1. / ( 0.01 + abs( $x - 5.) ) } }
    } \
    -body {
	foreach {x y} [min_bound_1d f 0 20.] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 5. \
    -match withinEpsilon

test linmin-2.1 {wrong \# args} \
    -body {
	min_bound_1d f
    } \
    -returnCodes 1 \
    -result [tcltest::wrongNumArgs min_bound_1d {f x1 x2 args} 1]

test linmin-2.2 {wrong \# args} \
    -body {
	min_bound_1d f 0 1 -bad
    } \
    -returnCodes 1 \
    -result "wrong # args, should be \"min_bound_1d f x1 x2 ?-option value?...\""

test linmin-2.3 {bad arg} \
    -body {
	min_bound_1d f 0 1 -bad option
    } \
    -returnCodes 1 \
    -result "unknown option \"-bad\", should be -abserror,\
             -fguess, -guess, -initial,\
             -maxiter, -relerror, or -trace"

test linmin-2.4 {iteration limit} \
    -setup {
	proc f x { expr { -1. / ( 0.01 + abs( $x - 5.) ) } }
    } \
    -body {
	min_bound_1d f 20. 0 -maxiter 10
    } \
    -cleanup {
	rename f {}
    } \
    -returnCodes 1 \
    -result "min_bound_1d failed to converge after \\d* steps" \
    -match regexp

test linmin-3.1 {minimise cos(x), unbounded} \
    -setup {
	proc f x { expr { cos($x) } }
    } -body {
	foreach { x y } [min_unbound_1d f 3. 3.01] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 3.1415926535897932 \
    -match withinEpsilon

test linmin-3.2 {minimise cos(x), unbounded, too eager} \
    -setup {
	proc f x { expr { cos($x) } }
    } -body {
	foreach { x y } [min_unbound_1d f 0.1 0.15] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result [expr { 3. * 3.1415926535897932 }] \
    -match withinEpsilon

test linmin-3.3 {near underflow in parabolic extrapolation} \
    -setup {
	proc f x {
	    expr { ( 1.12712e-22 * $x * $x * $x - 1e-15 ) * $x + 1e-15 }
	}
    } \
    -body {
	foreach { x y } [min_unbound_1d f 1. 0.] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 130.41372 \
    -match withinEpsilon

test linmin-3.4 {near underflow in parabolic extrapolation} \
    -setup {
	proc f x {
	    expr { ( ( 1e-30 * $x * $x - 1.12712e-22 )
		     * $x * $x * $x - 1e-15 )
		   * $x + 1e-15 }
	}
    } \
    -body {
	foreach { x y } [min_unbound_1d f 1. 0. -relerror 1e-08] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 8668.4248 \
    -match withinEpsilon

test linmin-3.5 {parabolic interpolation finds a minimum - case 1} \
    -setup {
	proc f x {
	    expr { ( ( ( 1e-5 * $x - 2.69672 )
		       * $x + 10.0902 )
		     * $x - 8.39345 )
		   * $x + 1. }
	}
    } \
    -body {
	foreach { x y } [min_unbound_1d f 1. 0. -relerror 1e-08] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 0.527450252 \
    -match withinEpsilon

test linmin-3.6 {parabolic interpolation finds a minimum - case 2} \
    -setup {
	proc f x {
	    expr { ( ( 0.125669 * $x * $x - 0.982687 )
		     * $x - 0.142982 )
		   * $x + 1 }
	}
    } \
    -body {
	foreach { x y } [min_unbound_1d f 1. 0. -relerror 1e-08] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 2.0127451 \
    -match withinEpsilon

test linmin-3.7 {parabolic interpolation is useless} \
    -setup {
	proc f x {
	    expr { ( ( ( 1e-5 * $x - 6.79171 )
		       * $x + 24.8107 )
		     * $x - 19.019 )
		   * $x + 1. }
	}
    } \
    -body {
	foreach { x y } [min_unbound_1d f 1 0 -relerror 1e-8] break
	set x
    } \
    -cleanup {
	rename f {}
    } \
    -result 509375.81 \
    -match withinEpsilon

test linmin-4.1 {wrong \# args} \
    -body {
	min_unbound_1d f
    } \
    -returnCodes 1 \
    -result [tcltest::wrongNumArgs min_unbound_1d {f x1 x2 args} 1]

test linmin-4.2 {wrong \# args} \
    -body {
	min_unbound_1d f 0 1 -bad
    } \
    -returnCodes 1 \
    -result "wrong # args, should be \"min_unbound_1d f x1 x2 ?-option value?...\""

test linmin-4.3 {bad arg} \
    -body {
	min_unbound_1d f 0 1 -bad option
    } \
    -returnCodes 1 \
    -result "unknown option \"-bad\", should be -trace"

#
# Test the solveLinearProgram procedure
#

set ::symm_constraints {
       { 1.0   2.0  1.0 }
       { 2.0   1.0  1.0 } }

test linprog-1.0 "Symmetric constraints, case 1" \
     -body {
	set result [solveLinearProgram {1.0 1.0} $::symm_constraints]
	set ok 1
	if { ! [within_range [lindex $result 0]  0.333300 0.333360] ||
	     ! [within_range [lindex $result 1]  0.333300 0.333360] } {
	   set ok 0
	}
	set ok
    } \
    -result 1

test linprog-1.1 "Symmetric constraints, case 2" \
     -body {
	set result [solveLinearProgram {1.0 0.0} $::symm_constraints]
	set ok 1
	if { ! [within_range [lindex $result 0]  0.49900 0.50100] ||
	     ! [within_range [lindex $result 1] -0.00100 0.00100] } {
	    set ok 0
	}
	set ok
    } \
    -result 1

test linprog-1.2 "Symmetric constraints, case 3" \
    -body {
	set result [solveLinearProgram {0.0 1.0} $::symm_constraints]
	set ok 1
	if { ! [within_range [lindex $result 1]  0.499900 0.500100] ||
	     ! [within_range [lindex $result 0] -0.000100 0.000100] } {
	      set ok 0
	}
	set ok
    } \
    -result 1

test linprog-1.3 "Symmetric constraints, case 4" \
    -body {
	set result [solveLinearProgram {3.0 4.0} $::symm_constraints]
	set ok 1
	if { ! [within_range [lindex $result 0]  0.333300 0.333360] ||
	     ! [within_range [lindex $result 1]  0.333300 0.333360] } {
	      set ok 0
	}
	set ok
    } \
    -result 1

#
# TODO: Current algorithm makes no difference between infeasible
#       and unbounded
#
test linprog-2.1 "Unbounded program" \
    -body {
	set result [solveLinearProgram {3.0 4.0} {{1.0 -2.0 1.0} {-2.0 1.0 1.0}} ]
    } \
    -result "infeasible"

test linprog-3.1 "Simple 3D program" \
   -body {
	set result [solveLinearProgram \
	   {1.0 1.0 1.0} \
	   {{1.0  1.0  2.0  1.0}
	    {1.0  2.0  1.0  1.0}
	    {2.0  1.0  1.0  1.0}}]
	set ok 1
	if { ! [within_range [lindex $result 0]  0.249900 0.250100] ||
	     ! [within_range [lindex $result 1]  0.249900 0.250100] ||
	     ! [within_range [lindex $result 2]  0.249900 0.250100] } {
	      set ok 0
	}
	set ok
   } \
   -result 1

test nelderMead-1.1 "Nelder-Mead - wrong \# args" \
    -body {
	::math::optimize::nelderMead f {0.0 0.0} -bogus
    } \
    -returnCodes error \
    -match glob \
    -result "wrong \# args*"
test nelderMead-1.2 "Nelder-Mead - bad param" \
    -body {
	::math::optimize::nelderMead f {0.0 0.0} -bogus 1
    } \
    -returnCodes error \
    -match glob \
    -result {unknown option "-bogus"*}
test nelderMead-1.3 "Nelder-Mead - bad size of scale" \
    -body {
	::math::optimize::nelderMead f {0.0 0.0} -scale {0 0 0}
    } \
    -returnCodes error \
    -result {-scale vector must be of same size as starting x vector}

# Easy case - minimize in a paraboloid

test nelderMead-2.1 "Nelder-Mead - easy" \
    -setup {
	proc f {x y} {
	    expr {($x-3.)*($x-3.) + ($y-2.)*($y-2.) + 1.}
	}
    } \
    -body {
	array set dd [::math::optimize::nelderMead f {1. 1.}]
	foreach {x y} $dd(x) break
	expr { abs($x-3.) < 0.001 && abs($y-2.) < 0.001 }
    } \
    -cleanup {
	rename f {}; unset dd
    } \
    -result 1

test nelderMead-2.2 "Nelder-Mead - easy" \
    -setup {
	proc f {x y} {
	    expr {($x-3.)*($x-3.) + ($y-2.)*($y-2.) + 1.}
	}
    } \
    -body {
	array set dd [::math::optimize::nelderMead f {0. 0.}]
	foreach {x y} $dd(x) break
	expr { abs($x-3.) < 0.001 && abs($y-2.) < 0.001 }
    } \
    -cleanup {
	rename f {}; unset dd
    } \
    -result 1

# Slalom down a sinuous valley - exercises most of the code

test nelderMead-2.3 "Nelder-Mead - sinuous valley" \
    -setup {
	set pi 3.1415926535897932
	proc f {x 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
	}
    } \
    -body {
	array set dd [::math::optimize::nelderMead f {1. 0.} -scale {0.1 0.01}]
	foreach {x y} $dd(x) break
	expr { abs($x-$pi/2) < 0.001 && abs($y-0.1) < 0.001 }
    } \
    -cleanup {rename f {}; unset dd} \
    -result 1

# Exercise the difficult case where the simplex has to contract about the
# low point because all else has failed.

test nelderMead-2.4 "Nelder-Mead - simplex contracts about the minimum" \
    -setup {
	proc g {a 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
	}
    } \
    -body {
	array set dd [::math::optimize::nelderMead g {0. 0.} \
			 -scale {1. 1.} -ftol 1e-10]
	foreach {x y} $dd(x) break
	expr { abs($x-0.774561) < 0.00005 && abs($y-0.755644) < 0.00005 }
    } \
    -cleanup {
	rename g {}; unset dd
    } \
    -result 1

testsuiteCleanup

#  Restore precision
if {![package vsatisfies [package present Tcl] 8.5]} {
    set ::tcl_precision $old_precision
}

# Local Variables:
# mode: tcl
# End:

} ;# End of optimizetest namespace