#OK to post homework
#TianHao Zhang,03/09/19

#################################Help######################################

Help := proc():print(`G1(f,var,pt,sz,MaxI,tol) G2(f,var,pt,sz,MaxI,tol)`):end:

#################################Part1#####################################
#First, in my opinion, I find GD1 works better when the curvature is large, as you can
#see, GD2 decreases more rapidly then GD1 which means when curvature is large, it may miss 
#the ture answer. Also means that GD1 can deal with the problem in a more precise way! 

#Second, the size of the step 'sz' can be changed during application. I take the value 
#norm(pt1(former one)-pt1(present one)) (in L2 space). Which is inspired by Newton method!

#So I design an function which have a cutoff, when the 'sz' is larger then the cut-off 
#I use GD2 and when the 'sz'  is less then this cutoffI use GD1. 




G1andG2:=proc(f,var,pt,sz,MaxI,cutoff) local g, g1, pt1,i,N1,flag,ptintial,sz1,it:
   pt1:=pt:
   sz1:=sz:
   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)))):
   it:=0:
   for i from 1 to MaxI while sz1 > 0.0000000000001 do
     
     ptintial:=pt1:     

     if sz1>cutoff then 
       pt1:=GD2(f,var,pt1,sz1):
       
       g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
       g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
       
     else
       pt1:=GD1(f,var,pt1,sz1):
       
       g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
       g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
     fi:
     
     sz1:=max(evalf(sqrt(add((pt1[i]-ptintial[i])^2,i=1..nops(pt1)))),sz):

     it:=it+1:
   od:
pt1:    

end:
   
       
#################################From C13##################################
#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 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,g1]:    
          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 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,g1]:    
          end: