#Final version

#C24.txt: April 16, 2018, Math 640 (Spring 2018), Rutgers University, Dr. Z.
#Homemade linear algebra
#Rather than use Maple's linear algebra package, it is more instructive to write our own
#Today's virtual class will consist in understanding and testing the Maple code below
#we also use linalg[det] (the home-made version sometimes FAILs)
#At the top you ave Yukun Yao's homework hw22.txt, illustrating the amazing Iterated polygon phenomenon
#Try, e.g. IterPolygon(20,N), for N=10,20,30, ...,  several times (of course, the way it is programmed now
#the intial polygon changes every time, you may want to write a version that
#plots successive iterations, and even animate them).  Notice that sometimes it takes
#a fairly long time until it reaches the ellipse (i.e. oval) shape, and that slant is alway close to 45 degrees
#but sometimes pointing left and sometimes ponting right, and the "fatness" (eccentricity) changes
print(`Welcome to virtual class C24`):
print(`For a list of the new procedures, type: Help(); `):
print(`For a list of the procedures from hw22.txt (thanks to Yukun Yao)  type: HelpYY(); `):


Help:=proc(): print(` MtV(M,V),  LCdet(M), RandM(m,n,N), ChM(M,x), ChP(M,x), BigGuys(L), EV(M,eps)  , Evec(M,lambda,eps) , EVC(M) `): end:

HelpYY:=proc(): print(`IterPolygon(n,N)`): end:

#######From Yukun Yao's homework
#OK to post homework
#Yukun Yao, Apr 14, 2018, Assignment 22

##Problem 2## For the random polygon project, get a feel for it, by
#Writing a procedure
#RandomGoodPoly(n)
#that inputs a positive integer and outputs a random polygon with n vertices, whose centroid is the origin.
#Hint: generate two random vectors, x and y, and substract from them their averages, also call these new vectors x and y to make their sum to zero and let the polygon be [[x1,y1], ..., [xn,yn],[x1,y1]]

#Write a procedure
#AvePolygon(P)
#that inputs a polygon, (given by its list of vertices), and outputs the polygon whose vertices are the midpoints of the edges.

#Write a procedure
#IterPolygon(n,N)
#that inputs a RandomGoodPoly(n) and a large pos. integer N, and performs AvePolygon(P) N times, and also plots it (use the command plot(P)). Perform IterPolygon(20,200) many times and look at the picture.

#To generate a random polygon usually the x and y are rational. So by scaling let x,y be integers and give a range for them. 
RandomGoodPoly:=proc(n) local i,S,Ave,ra,L:
S:=[]:
ra:=rand(0..1000):
for i from 1 to n do
  S:=[op(S), [ra(),ra()]]:
od:
Ave:=convert(S,`+`)/n:
L:=[Ave$n]:
S:=S-L:
S:=[op(S), S[1]]:
end:


AvePolygon:=proc(P) local n,S,i:
n:=nops(P)-1:
S:=[]:
for i from 1 to n do
  S:=[op(S),(P[i]+P[i+1])/2]:
od:
S:=[op(S), S[1]]:
end:


IterPolygon:=proc(n,N) local P,i:
P:=RandomGoodPoly(n):
for i from 1 to N do 
  P:=AvePolygon(P):
od:
plot(P):
end:

#It seems that the picture will converge to an oval.
#######End From Yukun Yao's homework

#start today's class
#Our date structure for vectors will be lists;  For matrices, lists of lists, we will NOT use Maple's linear algebra
#packages EXCEPT the solve command

#MtV(M,V): inputs a matrix M and a column vector V (expressed as a list), outputs the vector MV
#For example, try
#MtV([[1,2],[3,4]],[5,6]); MtV([[a11,a12],[a21,a22]],[b1,b2]);
MtV:=proc(M,V) local i,j:

if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and type(V,list) and
         nops(M)>0 and nops(M[1])=nops(V)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

[seq(add(M[i][j]*V[j],j=1..nops(M[i])),i=1..nops(M))]:

end:

#LCdet(M): the determinant of a matrix using the Rev. Dodgson's method (as nicely proved in
# http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/twotime.html
LCdet:=proc(M) local M11,M12,M21,M22,Mid,i,n:

if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and nops(M[1])=nops(M)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

n:=nops(M):
if n=1 then
  RETURN(M[1][1]):
elif n=2 then
  RETURN(M[1][1]*M[2][2]-M[1][2]*M[2][1]):
else

 M11:=[seq([op(1..n-1,M[i])],i=1..n-1)]:
 M12:=[seq([op(2..n,M[i])],i=1..n-1)]:

 M21:=[seq([op(1..n-1,M[i])],i=2..n)]:
 M22:=[seq([op(2..n,M[i])],i=2..n)]:

 Mid:=[seq([op(2..n-1,M[i])],i=2..n-1)]:

 if LCdet(Mid)=0 then
  print(`Thou should not divide by zero`):
   RETURN(FAIL):
  else
 RETURN(normal((LCdet(M11)*LCdet(M22)-LCdet(M12)*LCdet(M21))/LCdet(Mid))):
 fi:
fi:
 
end:

#RandM(m,n,N): a random m by n matrix with entries from 1 to N
RandM:=proc(m,n,N) local ra,i,j:
ra:=rand(1..N):
[seq([seq(ra(),i=1..n)],j=1..m)]:
end:



#ChM(M,x): the characteristic matrix of the square matrix M, w.r.t the variable x, i.e. M-Ix
ChM:=proc(M,x) local n,i,j:
if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and nops(M[1])=nops(M)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

n:=nops(M):

[seq([seq(M[i][j],j=1..i-1),M[i][i]-x,seq(M[i][j],j=i+1..n)],i=1..n)]:
end:

#ChP(M,x): the characteristic polynomial of the square matrix M, w.r.t the variable x, i.e. M-Ix
ChP:=proc(M,x)
linalg[det](ChM(M,x)):
end:




#BigGuys(L,eps): inputs a set L and outputs the places of those that have largest absolute value, eps is very small
#Warning: very naive
#For example BigGuys([1+I,1,1,1+I], 10^(-10)); returns {1,4}
BigGuys:=proc(L,eps) local champs,rec,i:

if not (type(L,list) and nops(L)>0) then
  RETURN(FAIL):
fi:

champs:={1}:

rec:=evalf(  abs(L[1])  ) :

for i from 2 to nops(L) do
 if evalf(abs(L[i]))-rec>eps then
  champs:={i}:
  rec:=evalf(abs(L[i])):
 elif  evalf(abs(abs(L[i])-rec))<eps then
  champs:=champs union {i}:
 fi:
od:
champs:
end:

Digits:=20:
#EV(M,eps): the list of eigenvalues of the square matrix M, and decreasing order of their absolute values. Exact version
#Warning, sorting is very naive
EV:=proc(M,eps) local P,x,L,i,S,s,NewL, L1,L2:

if not (type(M,list) and {seq(type(M[i],list),i=1..nops(M))}={true} and nops({seq(nops(M[i]),i=1..nops(M))})=1 and nops(M[1])=nops(M)) then
   print(`Bad input`):
   RETURN(FAIL):
fi:
P:=ChP(M,x):

L:=[solve(P,x)]:

L1:=L:
NewL:=[]:

while nops(L1)>0 do
 S:=BigGuys(L1,eps):
 NewL:=[op(NewL),seq(L1[s], s in S)]:
 L2:=[]:
  for i from 1 to nops(L1) do
   if not member(i,S) then
       L2:=[op(L2),L1[i]]:
   fi:
  od:
 L1:=L2:
od:

NewL:

end:



#Evec(M,lambda,eps): the eigenspace of the matrix M corresponding to the eigenvalue lambda
Evec:=proc(M,lambda,eps) local eq,var,i,n,v,M1,var1,v1:

n:=nops(M):

if abs(ChP(M,lambda))>eps then
 RETURN(FAIL):
fi:

v:=[seq(v[i],i=1..nops(M))]:
var:={seq(v[i],i=1..n)}:

M1:=ChM(M,lambda):
eq:=convert(MtV(M1,v),set):
var:=solve(eq,var):

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

var1:={}:

for v1 in var do
 if op(1,v1)=op(2,v1) then
  var1:=var1 union {op(1,v1)}:
fi:
od:


if nops(var1)<>1 then
 print(`Eigenspace is more than 1 dimensions.`):
 RETURN(FAIL):
fi:


v:=subs(var,v):


v:=subs(var1[1]=1,v):

v:
end:


#EVC(M): inputs a square matrix M, and outputs the list of eigenvalues (in decreasing order of absolute value)  
#and the respective unit eigenvectors, all in floating point.
EVC:=proc(M) local eps,gu,v,nv,mu,i,i1:
 eps:=10^(-Digits+5):
 gu:=EV(M,eps):
 mu:=[]:

 for i from 1 to nops(gu) do
  v:=Evec(M,gu[i],eps):
   if v=FAIL then
     RETURN(FAIL):
   fi:
  v:=evalf(v):
  nv:=sqrt(add(v[i1]^2,i1=1..nops(v))):
  v:=v/nv:
  mu:=[op(mu),v]:
od:

evalf(gu,10),evalf(mu,10):
end:

#EVCm(M): Same as EVC(M) but using Maple's linalg package, except that the eignevalues are not sorted according to maximum value
EVCm:=proc(M) local n,lu,gu,mu,i:
n:=nops(M):
lu:=[linalg[eigenvectors](M)]:
if nops(lu)<>n then
 RETURN(FAIL):
fi:



if {seq(lu[i][2],i=1..n)}<>{1} then
  print(`Not all eigenspaces are 1-dim.`):
  RETURN(FAIL):
fi:

gu:=[seq(evalf(lu[i][1],10),i=1..n)]:
mu:= [seq(evalf(lu[i][3][1],10),i=1..nops(lu))]:

gu,mu:

end:


##