#!/bin/sh
# If running in a UNIX shell, restart wish on the next line \
exec wish "$0" ${1+"$@"}


####  generate.tcl  ####

# Copyright (c) 2003, 2009 Oscar Garcia
# MIT OSI License.  See License.txt file.


package require Tk  ;#  for tclkit

set progID {Generate TR subroutines  2.0}

source getValue.tcl

# GUI: Get a derivative of Expression wrt Variable
proc getDer {Expression Variable} {
    set result [getValue .der Derivative \
    "Derivative of $Expression with respect to $Variable :" 30 {} quit 0 0]
    if {$result == {quit}} exit
    return $result
}


#  Generate derivative lists
#  Input: model parameters list
#  Output: pair with first and second derivatives lists
#    Items in firsts:  (basic_param variable derivative)
#    Items in seconds:  (basic_param variable1 variable2 derivative),
#    with variable1 <= variable2
#  Uses getDer
proc ders {Modlist} {
	set Firsts [list]
	set Seconds [list]
	set Basic 0
	# first derivatives
	foreach Expression $Modlist {
		incr Basic
		if {[regexp {^(G|L)[1-8]$} $Expression]} {
			lappend Firsts "$Basic $Expression 1"
		} elseif {![string is double $Expression]} {
			set Substring $Expression
			while {[regexp {(G|L)[1-8]} $Substring Var]} {
				lappend Firsts "$Basic $Var [getDer $Expression $Var]"
				regsub -all $Var $Substring {} Substring
			}
		}
	}
	# second derivatives
	foreach Der $Firsts {
		set Expression [lindex $Der 2]
		if {[regexp {^(G|L)[1-8]$} $Expression]} {
			if {$Expression >= [lindex $Der 1]} {
				lappend Seconds "[lindex $Der 0] [lindex $Der 1] $Expression 1"
			}
		} elseif {![string is double $Expression]} {
			set Substring $Expression
			while {[regexp {(G|L)[1-8]} $Substring Var]} {
				if {$Var >= [lindex $Der 1]} {
					lappend Seconds "[lindex $Der 0] [lindex $Der 1]  $Var [getDer $Expression $Var]"
				}
				regsub -all $Var $Substring {} Substring
			}
		}
	}
	return [list $Firsts $Seconds]
}

#  Take list of model parameter clasification, return TRPARS code
proc trpars {Modlist} {
	set Code {}
	set I 0
	foreach X $Modlist {
		incr I
		regsub -all {G([1-8])} $X {T0(\1)} X
		regsub -all {L([1-8])} $X {TK(\1)} X
		append Code T($I)= $X \n
	}
	return $Code
}


#  Generate TRDERS code
proc trders {Modlist} {
	# Get number of globals and locals, check
	set X $Modlist
	while {[regexp {G([1-8])} $X Var N]} {
		set Glob($N) 1
		regsub $Var $X {} X
	}
	while {[regexp {L([1-8])} $X Var N]} {
		set Loc($N) 1
		regsub $Var $X {} X
	}
	set Nglobals [array size Glob]
	for {set I 1} {$I <= $Nglobals} {incr I} {
		if {[catch {set x $Glob($I)}]} {
			error "ERROR: There is no G$I"
		}
	}
	set Nlocals [array size Loc]
	for {set I 1} {$I <= $Nlocals} {incr I} {
		if {[catch {set x $Loc($I)}]} {
			error "ERROR: There is no L$I"
		}
	}
	# OK, continue
	set X [ders $Modlist]
	set Firsts [lindex $X 0]
	set Seconds [lindex $X 1]
	set Code {}
	# Intermediate variables
	set X [list]
	foreach First $Firsts {
		set Der [lindex  $First 2]
		if {[string length $Der] > 2} {
			set Name [string tolower "d[lindex $First 0][lindex $First 1]"]
			append Code $Name = $Der \n
			set First [lreplace $First 2 2 $Name]
		} elseif {[regexp {\-.} $Der]} {
		    set First [lreplace $First 2 2 ($Der)]
		}
		lappend X $First
	}
	set Firsts $X
	set X [list]
	foreach Second $Seconds {
		set Der [lindex  $Second 3]
		if {[string length $Der] > 2} {
			set Name [string tolower "d[lindex $Second 0][lindex $Second 1][lindex $Second 2]"]
			append Code $Name = $Der \n
			set Second [lreplace $Second 3 3 $Name]
		} elseif {[regexp {\-.} $Der]} {
		    set First [lreplace $Second 3 3 ($Der)]
		}
		lappend X $Second
	}
	set Seconds $X
	# Terms with first derivatives
	foreach First $Firsts {
		set Basic [lindex $First 0]
		set Var [lindex $First 1]
		set Der [lindex $First 2]
		if {$Der != {1}} {
			set Der $Der*
		} else {
			set Der {}
		}
		append G($Var) + $Der G($Basic)
		foreach First2 $Firsts {
			set Var2 [lindex $First2 1]
			if {$Var <= $Var2} {
				set Basic2 [lindex $First2 0]
				set Der2 [lindex $First2 2]
				if {$Der2 != {1}} {
					set Der2 $Der2*
				} else {
					set Der2 {}
				}
				if {$Basic <= $Basic2} {
					append H($Var$Var2) + $Der $Der2 H($Basic,$Basic2)
				} else {
					append H($Var$Var2) + $Der $Der2 H($Basic2,$Basic)
				}
			}
		}
	}
	# Terms with second derivatives
	foreach Second $Seconds {
	    set Der [lindex $Second 3]
		if {$Der != {1}} {
			set Der $Der*
		} else {
			set Der {}
		}
		append H([lindex $Second 1][lindex $Second 2]) + $Der G([lindex $Second 0])
	}
	# Collect gradient terms
	foreach I [array names G] {
		if {[regexp {^G([1-8])$} $I X N]} {
			append Code G0($N) = $G($I) \n
		} elseif {[regexp {^L([1-8])$} $I X N]} {
			append Code GK($N) = $G($I) \n
		} else {
			error "OOPS! Wrong gradient array index: $I"
		}
	}
	# Collect Hessian terms
	foreach I [array names H] {
		append Code [eval hlhs [split $I {}] $Nlocals] $H($I) \n
	}
	# Substitute global and local names
	regsub -all {G([1-8])} $Code {T0(\1)} Code
	regsub -all {L([1-8])} $Code {TK(\1)} Code
	regsub -all {=\+} $Code {=} Code
	return $Code
}


#  Work out LHS for Hessian
proc hlhs {GL1 N1 GL2 N2 Nlocals} {
	if {$GL1 == $GL2} {
		if {$GL1 == {G}} {
			set Var H00(
		} else {
			set Var HKK(
		}	     
		return "$Var[expr {$N1+($N2*($N2-1))/2}])="
	} else {
		if {$GL1 == {L}} {
			set X $N1
			set N1 $N2
			set N2 $X
		}
		return "H0K([expr {($N1-1)*$Nlocals+$N2}])="
	}
}


#  Sort lines of code, indent, make continuations
proc sortLines {Input} {
		set Input [lsort -dictionary [split $Input "\n"]]
		set Output {}
		foreach Line $Input {
		    set Line "      $Line"
		    set Length [string length $Line]
		    while {$Length} {
		        if {$Length <= 72} {
		            append Output $Line \n
		            set Length 0
		        } else {
		            append Output [string range $Line 0 71] \n
		            set Line "     +  [string range $Line 72 end]"
		            set Length [string length $Line]
		        }
		    }
		}
		return $Output
}

#  Main function: Generate TRPARS and TRDERS
proc generateTR {Model} {
	set i {      }
	set  Code "C  ---  Model:  $Model  ---\n\n"
	append Code $i {SUBROUTINE TRPARS(T0,TK,T)} \n
	append Code $i {DOUBLE PRECISION T0(*),TK(*),T(8)} \n
	append Code $i {double precision tScale, hScale} \n
    append Code $i {integer model} \n
    append Code $i {integer*4 plotID} \n
    append Code $i {common /ESDE/ tScale, hScale, plotID, model}
	append Code [sortLines [trpars $Model]]
	append Code $i RETURN \n $i END \n \n
	append Code $i {SUBROUTINE TRDERS(T,T0,TK,F,G,H,G0,GK,H00,H0K,HKK)} \n
	set Ders [trders $Model]
	if {[regexp d $Ders]} {
    	append Code $i {IMPLICIT DOUBLE PRECISION (d)} \n
	}
	append Code $i {DOUBLE PRECISION T(8),T0(*),TK(*),F,G(8),H(8,8),G0(*),GK(*),} \n
	append Code {     +  H00(*),H0K(*),HKK(*)} \n
	append Code $i {double precision tScale, hScale} \n
    append Code $i {integer model} \n
    append Code $i {integer*4 plotID} \n
    append Code $i {common /ESDE/ tScale, hScale, plotID, model}
	append Code [sortLines $Ders]
	return [append Code $i RETURN \n $i END \n]
}


########  GUI  #########

# GUI: Results log
proc resLog X {
    toplevel .subs
	wm title .subs "Subroutines"
    source resLog.ui.tcl
    resLog_ui .subs $X
}


# GUI: Main
wm withdraw .
set Model [getValue .model $progID Model: 66 {} quit 0 0]
if {$Model == {quit}} exit
resLog [generateTR $Model]

