# hw13MattHohertz.txt

Help:=proc():
	PrintHelp([G]);
end:

G:=proc(f,var::list,pt::list,tol::numeric,MaxIters::posint)
	description "Runs G1, with sz = min(10, |Grad(f)| / 1000)\n\n",
	" >> INPUTS: A function <f> in variables <var>, an initial point <pt>, a tolerance <tol>, and a maximum number <MaxIters> of iterations.",
	" << OUTPUTS: Local minimum, or result after <MaxIters> iterations.";
	local i, cpt, g, avGrad, j, gg:
	
	if not (tol > 0) then
		printf("G: tol must be strictly positive");
		pskip();
		return(FAIL);
	fi:
	
	cpt:=pt:
	g:=Grad(f,var):
	
	for i from 1 to MaxIters do
		gg:=subs({seq(var[i]=cpt[i],i=1..nops(var))},g):
		avGrad:=evalf(sqrt(add(gg[j]^2,j=1..nops(gg)))):
		if (avGrad < tol) then
			printf("Approximate minimum attained (i = %d, tol = %f, avGrad = %f):", i, tol, avGrad);
			pskip();
			return(cpt);
		fi:
		
		cpt:=GD1(f,var,cpt,min(10,avGrad/100)):
		# print(cpt);
	od:
	
	printf("Max attempts exhausted; best guess is (tol = %f, avGrad = %f)", tol, avGrad);
	nl();
	cpt:
end:







































## BEGIN MRHSTDLIB (helper functions)
insertLines:=proc(n::integer) 
	description "Appends {n} lines.",
	"See also: nl, pskip":
	local i:

	for i from 1 to n do:
		printf("\n");
	od:
end:

nl:=proc(): insertLines(1): end:
pskip:=proc(): insertLines(2): end:

PrintHelp:=proc(LofPs::list(procedure)) 
	description "Prints descriptions for the procedures in {LofPs}."
	:
	local P, i, NProcs:

	NProcs:=nops(LofPs):

	for i from 1 to NProcs do
		P:=LofPs[i]:
		printf("### %s ###", P);
		nl();
		Describe(P);
		pskip();
	od:
end:
## END MRHSTDLIB

# hw13
#C13.txt: March 7, 2019 Dr. Z.'s Math640(RU)
#Implementing the gradient descent algorithm
#Help:=proc(): print(`OptF(f,var), Grad(f,var) , GD1(f,var,pt,sz), GD2(f,var,pt,sz) `): 
#print(`G1(f,var,pt,sz,MaxI,tol), G2(f,var,pt,sz,MaxI,tol) `):
#end:

#OptF(f,var): inputs a nice function f (expression)
# of a list of variables, and finds the local max and min locations
#and values. For example (x^2+y^2+1,[x,y]); should output {[[0,0],1]}.
OptF:=proc(f,var) local eq,i,sol,L:
eq:={seq(diff(f,var[i])=0,i=1..nops(var))}:
sol:=solve(eq,var):

if sol=NULL then
	RETURN(FAIL):
fi:

L:=[seq(subs(sol[i],var),i=1..nops(sol))]: 
L:
end:

Grad:=proc(f,var) local i: option remember:
 [seq(diff(f,var[i]),i=1..nops(var))]: end:

#GD1(f,var,pt,sz): inputs a function (expression) f in the list of variables var, a point pt
#a step-size sz, and a small parameter e indicating when to stop
#Performs one step, returns the next point NORMALIZED
GD1:=proc(f,var,pt,sz) local g,i,Ng:
if not ( type(var,list) and type(pt,list) and nops(var)=nops(pt) and type(sz,numeric) and sz>0) then
	print(`bad input`):
	RETURN(FAIL):
fi:
	
 g:=Grad(f,var):

 g:=subs({seq(var[i]=pt[i],i=1..nops(var))},g):
 Ng:=evalf(sqrt(add(g[i]^2,i=1..nops(g)))):
 g:=g/Ng:

 pt-sz*g:

 end:

#GD2(f,var,pt,sz): inputs a function (expression) f in the list of variables var, a point pt
#a step-size sz, and a small parameter e indicating when to stop
#Performs one step, returns the next point NOT Normalized
GD2:=proc(f,var,pt,sz) local g,i:
if not ( type(var,list) and type(pt,list) and nops(var)=nops(pt) and type(sz,numeric) and sz>0) then
	print(`bad input`):
	RETURN(FAIL):
fi:
	
 g:=Grad(f,var):

 g:=subs({seq(var[i]=pt[i],i=1..nops(var))},g):
 

 pt-sz*g:

 end:

		#G1(f,var,pt,sz,MaxI,tol): inputs a function f given as an expression, in a list of variables var
		#a starting point pt, MaxI: number of iterations, and tol a small parameter, we stop when the
		#norm of the gradient is less than tol, or hit MaxI, iterations. NORMALIZED
	 G1:=proc(f,var,pt,sz,MaxI,tol) local g, g1, pt1,i,N1:
	 pt1:=pt:
		 g:=Grad(f,var):
			 g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
			 g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
		 for i from 1 to MaxI while g1>tol do
			 pt1:=GD1(f,var,pt1,sz):
			 
			 g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
			 g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
	 od:
	 pt1:    
					end:
 #G2(f,var,pt,sz,MaxI,tol): inputs a function f given as an expression, in a list of variables var
		#a starting point pt, MaxI: number of iterations, and tol a small parameter, we stop when the
		#norm of the gradient is less than tol, or hit MaxI, iterations. NOT NORMALIZED
	 G2:=proc(f,var,pt,sz,MaxI,tol) local g, g1, pt1,i,N1:
	 pt1:=pt:
		 g:=Grad(f,var):
			 g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
			 g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
		 for i from 1 to MaxI while g1>tol do
			 pt1:=GD2(f,var,pt1,sz):
			 
			 g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
			 g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
	 od:
	 pt1:    
					end: