# Initializing procedure for simplex method
proc simplex_init { } {
    global ndim p psub y iter size_simplex xopt num_dp opt_ite opt_loop
    set opt_loop 1
    set ndim $num_dp

    OpenOptLogWindow
    InitializeOptLog

    source opt_initial_values.tcl

    global opt_flag
    set opt_flag 2

    simplex_main
}


# Main routine for simplex method
proc simplex_main { } {
    global ndim p psub y iter size_simplex xopt num_dp opt_ite \
           opt_loop ilo

    set iter 0
    set opt_ite 1

    # first, compute initial "p"
    for {set n 1} {$n<=$ndim} {incr n} {
        if {$opt_loop==1} {
            # set initial vector for the first loop
            set p(1,$n) $xopt($n)
        } else {
            # On the nxet loop, the lowest (best) vector has been 
            # already set as p(1,$n), so do nothing here.
        }
    }

    # pn & qn
    set pn [ expr ( pow(double($ndim)+1.0, 0.5)-1.0+double($ndim) ) \
                / ( double($ndim)*pow(2.0,0.5) ) * $size_simplex ]
    set qn [ expr ( pow(double($ndim)+1.0, 0.5)-1.0 ) \
                / ( double($ndim)*pow(2.0,0.5) ) * $size_simplex ]

    for {set n 1} {$n<=$ndim} {incr n} {
        for {set m 1} {$m<=$ndim} {incr m} {
            set p([expr $m+1],$n) $p(1,$n)
            if {$m==$n} {
                set p([expr $m+1],$n) [expr $p(1,$n)+ $pn]
            } else {
                set p([expr $m+1],$n) [expr $p(1,$n)+ $qn]
            }
        }
    }

    # compute y for initial p
    for {set m 1} {$m<=[expr $ndim+1]} {incr m} {
        for {set n 1} {$n<=$ndim} {incr n} {
            #set ptry($n) $p($m,$n)
            #set psub($n) $ptry($n)
            set psub($n) $p($m,$n)
        }
        set y($m) [ func ]
    }

    # next, enter simplex method
    for {set loop 1} {$loop<=1000000} {incr loop} {
        simplex1
        simplex2
    }
    set opt_flag 0
}


proc simplex1 { } {
    global ndim p psum
    # compute psum
    for {set n 1} {$n<=$ndim} {incr n} {
        set sum 0.0
        for {set m 1} {$m<=[expr $ndim+1]} {incr m} {
            set sum [expr $sum + $p($m,$n)]
        }
        set psum($n) $sum
    }
}


proc simplex2 { } {
    global y ftol tiny p iter ndim ihi psub psum f ilo

    for {set loop 1} {$loop<=100000} {incr loop} {

        # Determine which point is the highest, next highest, and lowest
        set ilo 1
        if {$y(1) > $y(2)} {
            set ihi  1
            set inhi 2
        } else {
            set ihi  2
            set inhi 1
        }
        for {set i 1} {$i<=[expr $ndim+1]} {incr i} {
            if {$y($i) <= $y($ilo)} {set ilo $i}
            if {$y($i) >  $y($ihi)} {
                set inhi $ihi
                set ihi $i
            } elseif {$y($i) > $y($inhi) } {
                if {$i != $ihi} {set inhi $i}
            }
        }


        # Judge for termination
        set rtol [expr 2.0*abs($y($ihi)-$y($ilo)) \
                      / (abs($y($ihi))+abs($y($ilo))+$tiny)]
        if {$rtol < $ftol} {
            set swap $y(1)
            set y(1) $y($ilo)
            set y($ilo) $swap
            for {set n 1} {$n<=$ndim} {incr n} {
                set swap $p(1,$n)
                set p(1,$n) $p($ilo,$n)
                set p($ilo,$n) $swap
            }

            global opt_loop opt_lastloop
            if {$opt_loop==$opt_lastloop} {
                .opt.ft.t configure -state normal
                .opt.ft.t insert end \
                       "\nOptimization successfully completed!"
                .opt.ft.t configure -state disabled
                set f [open "optimization_log.dat" w]
                puts -nonewline $f [.opt.ft.t get 1.0 end]
                close $f
            }

            global opt_flag
            set opt_flag 0

            # enter next loop
            if {$opt_loop<$opt_lastloop} {
                incr opt_loop
                set opt_flag 1
                simplex_main
            }

            # end (wait forever)
            tkwait variable nonexistent_variable
        }

        set iter [expr $iter +2]

        # 1st try
        #puts $f "simptry -1.0"
        set ytry [ simptry -1.0 ]

        if {$ytry <= $y($ilo)} {
            # expansion
            #puts $f "simptry 2.0"
            set ytry [ simptry 2.0 ]
        } elseif {$ytry >= $y($inhi)} {
            # one-dimensional contraction
            set ysave $y($ihi)
            #puts $f "simptry 0.5"
            set ytry [ simptry 0.5 ]
            if {$ytry >= $ysave} {
                #puts $f "multi dim contraction"
                for {set i 1} {$i<=[expr $ndim+1]} {incr i} {
                    if {$i != $ilo} {
                        for {set j 1} {$j<=$ndim} {incr j} {
                            set psum($j) [expr 0.5*($p($i,$j)+$p($ilo,$j))]
                            set psub($j) $psum($j)
                            set p($i,$j) $psum($j)
                        }
                        set y($i) [ func ]
                    }
                }
                set iter [expr $iter + $ndim]
                return
            }
        } else {
            set iter [expr $iter -1]
        }

    }
}


proc simptry { fac } {
    global ndim psum ihi p y psub
    global f

    set fac1 [expr (1.0-$fac)/double($ndim)]
    set fac2 [expr $fac1 - $fac]
    for {set j 1} {$j<=$ndim} {incr j} {
        set ptry($j) [expr $psum($j)*$fac1-$p($ihi,$j)*$fac2]
        set psub($j) $ptry($j)
    }

    set ytry [ func ]

    if {$ytry < $y($ihi)} {
        #puts $f "try so-so succeed"
        set y($ihi) $ytry
        for {set j 1} {$j<=$ndim} {incr j} {
            set psum($j) [expr $psum($j)-$p($ihi,$j)+$ptry($j)]
            set p($ihi,$j) $ptry($j)
        }
    } else {
        #puts $f "try failed"
    }
    update
    return $ytry
}


proc func { } {
    global psub ndim xopt opt_ite opt_loop

    for {set j 1} {$j<=$ndim} {incr j} {
        set xopt($j) $psub($j)
    }

    set str0 [string range "\n$opt_loop               " 0 6]
    set str1  [string range "$opt_ite                " 0 9]
    set str "$str0$str1"
    for {set i 1} {$i<=$ndim} {incr i} {
        set str2 [string range "$xopt($i)                       " 0 16]
        #set str2 $xopt($i)
        set str "$str  $str2"
    }
    .opt.ft.t configure -state normal
    .opt.ft.t insert end $str
    .opt.ft.t configure -state disabled
    set f [open "optimization_log.dat" w]
    puts -nonewline $f [.opt.ft.t get 1.0 end]
    close $f

    update

    # load opt_design_parameters.tcl
    # (produce input data files from xopt)
    source opt_design_parameters.tcl

    global end_flag
    set end_flag 0
    if {[file exists "end_flag.dat"]=="1"} {
        file delete "end_flag.dat"
    }

    # execute analysis engine
    global EngineExeCommand EngineExeCommand2 silent_flag
    if {$silent_flag=="1"} {
        eval $EngineExeCommand2
    } else {
        eval $EngineExeCommand
    }

    global opt_flag
    set opt_flag 1

    # wait until calculation done
    WaitEngineEnd
    tkwait variable end_flag

    # calculate objective function from output values
    source opt_objective_function.tcl

    set str "  $yopt"
    .opt.ft.t configure -state normal
    .opt.ft.t insert end $str
    .opt.ft.t configure -state disabled
    set f [open "optimization_log.dat" w]
    puts -nonewline $f [.opt.ft.t get 1.0 end]
    close $f

    update

    incr opt_ite

    global itmax
    if {$opt_ite>$itmax} {
        .opt.ft.t configure -state normal
        .opt.ft.t insert end "\nIteration number exceeds limit!"
        .opt.ft.t configure -state disabled

        set f [open "optimization_log.dat" w]
        puts -nonewline $f [.opt.ft.t get 1.0 end]
        close $f

        global opt_flag
        set opt_flag 0
        tkwait variable nonexistent_variable
    }

    return [expr -$yopt]
}


proc WaitEngineEnd { } {
    if {[file exists end_flag.dat]=="1"} {
        global end_flag 1
        set end_flag 1
    } else {
        after 500 { WaitEngineEnd }
    }
}

# Open optimization log window
proc OpenOptLogWindow { } {
    global ndim xopt yopt
    if {[winfo exists .opt]==0} {
        toplevel .opt 
        wm title .opt "Swumsuit: Optimization Log"

        # log (text widget)
        frame .opt.ft
        text  .opt.ft.t -xscrollcommand ".opt.ft.h set" \
                        -yscrollcommand ".opt.ft.v set" -wrap none
        scrollbar .opt.ft.h -orient horizontal \
                            -command ".opt.ft.t xview" 
        scrollbar .opt.ft.v -orient vertical \
                            -command ".opt.ft.t yview"
        grid .opt.ft.t -row 0 -column 0 -sticky nsew
        grid .opt.ft.h -row 1 -column 0 -sticky ew
        grid .opt.ft.v -row 0 -column 1 -sticky ns
        grid rowconfigure    .opt.ft  0 -weight 100
        grid columnconfigure .opt.ft  0 -weight 100

        # panel
        frame .opt.p
        global Close Abort
        button .opt.p.a -text $Abort -command AbortOptimization
        button .opt.p.c -text $Close -command "wm withdraw .opt"
        pack .opt.p.a .opt.p.c -side left -padx 20 -pady 10 -expand 1

        pack .opt.ft -fill both -expand 1
        pack .opt.p

        ChangeIcon .opt
    } else {
        wm deiconify .opt
    }

    .opt.ft.t configure -state normal
    .opt.ft.t delete 1.0 end
    # if log file exists, load it
    if {[file exists "optimization_log.dat"]=="1"} {
        set f [open "optimization_log.dat" r]
        .opt.ft.t insert end [read $f]
        close $f
    }
    .opt.ft.t configure -state disabled
}

proc InitializeOptLog { } {
    .opt.ft.t configure -state normal

    # first, all clear
    .opt.ft.t delete 1.0 end

    # insert legends
    global num_dp
    set str "loop  iteration   "
    for {set i 1} {$i<=$num_dp} {incr i} {
        set str2 [string range "xopt($i)                        " 0 18]
        set str $str$str2
    }
        set str2 "yopt             "
        set str $str$str2

        .opt.ft.t insert end $str
        .opt.ft.t configure -state disabled
        set f [open "optimization_log.dat" w]
        puts -nonewline $f [.opt.ft.t get 1.0 end]
        close $f
}


proc AbortOptimization { } {
    global opt_flag
    if {$opt_flag=="0"} { return }
    .opt.ft.t configure -state normal
    .opt.ft.t insert end "\nOptimization Aborted!"
    .opt.ft.t configure -state disabled

    set f [open "optimization_log.dat" w]
    puts -nonewline $f [.opt.ft.t get 1.0 end]
    close $f

    set opt_flag 0

    # stop here (wait forever...)
    tkwait variable nonexistent_variable
}


# initial values for optimization
set tiny 1.0e-11
set itmax 10000
set ftol 1.0e-8
set opt_flag 0
set silent_flag 0
set opt_lastloop 1
