# Polygon Project Methods and Details

# For all parts of this help, a Polygon on n vertices is represented by a list of n+1 points [x,y] with the first and last entries the same.

read `PolygonSupport.txt`:
with(plots):

read `PolygonHelp.txt`:


### Polygon Generation

# RandomGoodPoly(n): Generates a random polygon on n vertices with centroid (0,0).
# Input: n, a number of vertices
# Output: P, a list of n+1 pairs of numbers [x,y], where the last entry is a repeat of the first to aid in plotting the polygon.

RandomGoodPoly:=proc(n) local x, y, xbar, ybar, poly, i, roll:
roll:=rand(-20..20):
x:=[]:
y:=[]:

#randomly generates x and y values
for i from 1 to n do
x:=[op(x), roll()]:
y:=[op(y), roll()]:
od:

#calculates the average of the x and y values
xbar:=1/n*add(x[i], i=1..n):
ybar:=1/n*add(y[i], i=1..n):

#creates the polygon
poly:=[]:
for i from 1 to n do
poly:=[op(poly), [x[i]-xbar, y[i]-ybar]]:
od:

poly:=[op(poly), poly[1]]:
end:

### Iteration and Rescaling Methods

# WeightedMatrix(n, alpha): Generates the iterating matrix (implemented as a list of lists) for the process defined by P_new(i) = alpha*P(i) + (1-alpha)*P(i+1)
# Input: n, a natural number, and alpha, a real number between 0 and 1.
# Output: An n by n matrix A corresponding to the iteration process.

WeightedMatrix:=proc(n, alpha) local i:
[seq([0$i, alpha, 1-alpha, 0$(n-2-i)], i=0..n-2), [1-alpha, 0$(n-2), alpha]]:
end:

# Vector Version of the Weighted Matrix Stuff. See section 4.1 in the paper.
# Input: n, a natural number, and eta, a list of at most n positive real numbers whose sum is at most 1.
# Output: An n by n matrix A corresponding to the iteration process with vector eta, where a new element is added (if needed) to make the total sum equal to 1.
WeightedMatrixVec := proc(n,eta) local i,j, eta_F:
if nops(eta) > n or evalf(add(eta[i], i=1..nops(eta))) > 1 then
	return FAIL:
fi:
if nops(eta) = n and evalf(add(eta[i], i=1..nops(eta))) < 1 then
	return FAIL:
fi:

eta_F := [op(eta), 1 - evalf(add(eta[i], i=1..nops(eta))), 0$(n - nops(eta)-1)]:
[seq([seq(eta_F[((j-i) mod n) + 1], j=1..n)], i=1..n)]:
end:

# WeightedAvePolygon(P, alpha): Computes the polygon generated by averaging the vertices of P according to the formula P_new(i) = alpha*P(i) + (1-alpha)*P(i+1)
# Input: P, a polygon on n vertices, and alpha, which can be either a number between 0 and 1 or eta, the weighting vector.
# Output: A new polygon P2 that is the average of these points according to the formula.

WeightedAvePolygon:=proc(P,alpha) local matrixProd, newPoly,i,j,n, x, y,M:
n := nops(P) - 1:
x := [seq(P[i][1], i=1..n)]:
y := [seq(P[i][2], i=1..n)]:
newPoly := []:
if not type(alpha,list) then
M := WeightedMatrix(nops(P)-1, alpha):
else
M := WeightedMatrixVec(n,alpha):
fi:
for i from 1 to n do
newPoly := [op(newPoly), [add(M[i][j]*x[j], j=1..nops(x)), add(M[i][j]*y[j], j=1..nops(y))]]:
od:
newPoly:=[op(newPoly), newPoly[1]]:
end:

# NormRescale(P): Takes a polygon P and rescales it so that the x and y vectors have 2-norm equal to 1.
# Input: P, a polygon on n vertices.
# Output: A new polygon P rescaled so that each of the x and y vectors have norm 1.

NormRescale:=proc(P) local xnorm, ynorm, i:
xnorm:=(add((P[i][1])^2, i=1..nops(P)-1))^(1/2):
ynorm:=(add((P[i][2])^2, i=1..nops(P)-1))^(1/2):

[seq([simplify(P[i][1]/xnorm), simplify(P[i][2]/ynorm)], i=1..nops(P))]:
end:

#### The next set of methods will initially only work for the strict averaging problem. If we can prove that they will still be ellipses, then they will all work.

# EValRescale(P,alpha): Takes a polygon P and rescales it *once* by the magnitude of the second largest eigenvalue of the iterating matrix A(n). 
# Input: P, a polygon on n vertices, and alpha, either a number between 0 and 1, or a vector corresponding to weights of the averaging procedure.
# Output: A new polygon P where each entry is scaled up by 1/|lambda_1|.

EValRescale := proc(P,alpha) local lambda, n,i, P_new:
n := nops(P)-1:
if not type(a,list) then
lambda := abs(alpha + (1-alpha)*exp(2*Pi*I/n)):
else
lambda := abs(EVC(WeightedMatrixVec(n, alpha))[2][1]):
fi:
P_new := []:
for i from 1 to nops(P) do
P_new := [op(P_new), P[i]*~(1/lambda)]:
od:
evalf(Recenter(P_new)):
end:

Recenter := proc(P) local xbar, ybar, n, i, poly:
#calculates the average of the x and y values
n := nops(P) - 1:
xbar:=1/n*add(P[i][1], i=1..n):
ybar:=1/n*add(P[i][2], i=1..n):

poly := []:
for i from 1 to n do
poly :=[op(poly), [P[i][1]-xbar, P[i][2]-ybar]]:
od:

poly:=[op(poly), poly[1]]:
end:

# IterPolygon(P, N, {alpha:=1/2}, {output:="plot"}, {scaling:="eval"})
# Performs N averaging iterations on the polygon P.
# INPUT:
# P: a polygon
# N: number of iterations
# alpha: optional keyword argument for the averaging parameter, 1/2 by default
# output: optional keyword argument, "plot" (default) or "list"
# scaling: optional keyword argument, "eval" (default), "none", or "norm"
# OUTPUT: 
# a plot of the final iteration or list of all iterations.
# For example, try this:
# IterPolygon(RandomGoodPoly(10),20);

IterPolygon:=proc(P, N, {alpha:=1/2}, {output:="plot"}, {scaling:="eval"}) local L,P1,i:
if not member(output,{"plot", "list"}) then 
	print(`ERROR in IterPolygon: output must be "plot" or "list"`):
	return FAIL:
fi:

if not member(scaling,{"eval", "none", "norm"}) then 
	print(`ERROR in IterPolygon: scaling must be "eval", "none", or "norm"`):
	return FAIL:
fi:

L:=[P]:
for i from 1 to N do
	P1:= WeightedAvePolygon(L[-1],alpha):
	if scaling="norm" then
		P1:=NormRescale(P1):
	elif scaling="eval" then
		P1:=EValRescale(P1, alpha):
	fi:
	L:=[op(L), P1]:
od:

if output="plot" then
	plot(L[-1], ':-scaling'=constrained):
else
	L:
fi:
end:

#### Ellipse Analysis Methods

# ParEllipse(P,t): Takes a polygon P and finds the parametric equation for the ellipse that it should converge to in the limit. 
# Input: P, a polygon on n vertices, and a variable t
# Output: [x(t), y(t)], the parametric equation for the ellipse given by the coefficient of the eigenvectors of the second largest absolute value for the iterating matrix 

ParEllipse:=proc(P,t, {alpha:=1/2}, {scaling:="eval"}) local n, omega, Px, Py, alph, beta, x, y,M,i:
n:=nops(P)-1:
omega := exp(2*Pi*I/n):
Px := [seq(P[i][1],i=1..n)]:
Py := [seq(P[i][2],i=1..n)]:

if not type(alpha,list) then
M := WeightedMatrix(n, alpha):
else
M := WeightedMatrixVec(n,alpha):
fi:

# alph:=CB(M,Px)[2]:
# beta:=CB(M,Py)[2]:

alph := add(Px[i]*omega^(i-1)/sqrt(n), i=1..n):
beta := add(Py[i]*omega^(i-1)/sqrt(n), i=1..n):

if scaling = "norm" then
	alph := alph/(sqrt(2)*abs(alph)):
	beta := beta/(sqrt(2)*abs(beta)):
fi:

x := 1/sqrt(n)*(2*Re(alph)*cos(t)-2*Im(alph)*sin(t)):
y := 1/sqrt(n)*(2*Re(beta)*cos(t) - 2*Im(beta)*sin(t)):
[x,y]:
end:

# NormalFormEllipse(P): Takes a polygon P and finds the normal form equation Ax^2 + Bxy + Cy^2 = 1 that the limiting ellipse will satisfy.
# Input: P, a polygon on n vertices.
# Output: [A,B,C] the coefficients in the Ax^2 + Bxy + Cy^2 = 1 normal form that the ellipse satisfies

NormalFormEllipse:=proc(P, {alpha:= 1/2}, {scaling:="eval"}) local n,omega, M, Px, Py, alph, beta, a,b,c,d,G,H,K:
n:=nops(P)-1:
Px := [seq(P[i][1],i=1..n)]:
Py := [seq(P[i][2],i=1..n)]:
omega := exp(2*Pi*I/n):

if not type(alpha,list) then
M := WeightedMatrix(n, alpha):
else
M := WeightedMatrixVec(n,alpha):
fi:

# alph:=CB(M,Px)[2]:
# beta:=CB(M,Py)[2]:

alph := add(Px[i]*omega^(i-1)/sqrt(n), i=1..n):
beta := add(Py[i]*omega^(i-1)/sqrt(n), i=1..n):

if scaling = "norm" then
	alph := alph/(sqrt(2)*abs(alph)):
	beta := beta/(sqrt(2)*abs(beta)):
fi:

a := 1/sqrt(n)*2*Re(alph):
b := -1/sqrt(n)*2*Im(alph):
c := 1/sqrt(n)*2*Re(beta):
d := -1/sqrt(n)*2*Im(beta):
G := (d^2+c^2)/(a*d-b*c)^2:
H := (2*b*d+2*a*c)/(a*d-b*c)^2:
K := (a^2+b^2)/(a*d-b*c)^2:
evalf([G,H,K]):
end:

# EllipseData([A,B,C]): Determines the orientation and eccentricity of the ellipse Ax^2 + Bxy + Cy^2 = 1
# Input: [A,B,C], the three coefficients of the equation of the ellipse
# Output: [ang, ecc, sMAx, smAx], where sMAx is the length of the semi-Major axis, smAx is the semi-minor axis, ecc is the eccentricity of the ellips, and ang is the angle of rotation assuming that the semi-major axis lies along the x-axis before the rotation

EllipseData:=proc(L) local a,b,c,e,orien,A1,B1,C1,minor,maj,truth,X,cand,can:
A1:=L[1]:
B1:=L[2]:
C1:=L[3]:

e:=evalf(sqrt(2*sqrt((A1-C1)^2+B1^2)/(A1+C1+sqrt((A1-C1)^2+B1^2)))):
if e >= 1 then
 print(`That's no ellipse!`):
 return(FAIL):
fi:

if A1 = C1 then
	orien := evalf(Pi/4*sign(B1)):
else
	orien := evalf(1/2*arctan(B1/(A1-C1))):
fi:

if orien < 0 then
	orien := orien + Pi/2:
fi:

if B1 > 0 then
	orien := orien + Pi/2:
fi:

orien := Pi - orien:

if orien = 0 or orien = Pi/2 then
	a := sqrt(1/A1):
	b := sqrt(1/B1):
else
	a := sqrt(abs(2/(A1+C1+B1/sin(2*orien)))):
	b := sqrt(abs(2/(A1+C1-B1/sin(2*orien)))):
fi:

maj := max(a,b):
minor := min(a,b):

e := sqrt(1 - minor^2/maj^2):

[orien,e,maj,minor]:
end:


EllipseProperties := proc(P, {alpha:=1/2}, {scaling:="eval"}) local yMaj, ang, ecc, sMAx, smAx, A, B, C:
A,B,C := op(NormalFormEllipse(P, ':-alpha' = alpha,':-scaling'=scaling)):
ang, ecc, sMAx, smAx := op(EllipseData([A,B,C])):

if scaling = "eval" then
	printf(`Under eigenvalue rescaling, `);
elif scaling = "norm" then
	printf(`Under norm rescaling, `);
else
	print("No rescaling given!");
	return(FAIL):
fi:

printf(`the limiting ellipse will have normal form equation:\n`);
printf(" %.2f x^2 + %.2f xy + %.2f y^2 = 1\n", A, B, C); 
printf(`which has:\n`); 
printf("semi-major axis length = %.4f\n", sMAx); 
printf("semi-minor axis length = %.4f\n", smAx); 
printf("eccentricity = %.5f\n", ecc); 
printf("with the semi-major axis tilted at an angle of %.2f degrees ", ang*180/Pi); 
printf(`counter-clockwise from the x-axis.\n`)

end:


# AnimateConvergence(P,N, {alpha:=1/2}, {scaling:="eval"})
# Displays N iterations of averaging on the polygon P.
# INPUT:
# P: a polygon
# N: number of iterations
# alpha: optional keyword argument for the averaging matrix, 1/2 by default
# scaling: optional keyword argument, "eval" (default), "none", or "norm"
# OUTPUT: 
# an animated plot showing the iterations, and, if scaling="eval", the limiting ellipse.
# For example, try this:
# AnimateConvergence(RandomGoodPoly(10),20);

AnimateConvergence:= proc(P,N, {alpha:=1/2}, {scaling:="eval"}) local t,q, plot2, anim, X, Y, L, LL:
L:=IterPolygon(P, N, output="list", ':-alpha'=alpha, ':-scaling'=scaling):
L:=L[2..nops(L)]: # throw away the first (not rescaled) polygon
LL:=ListTools[Flatten](L,1):
X:={seq(p[1], p in LL)}:
Y:={seq(p[2], p in LL)}:

anim := display([seq(plot(pl, view=[min(X)..max(X), min(Y)..max(Y)]), pl in L)], insequence = true);

display([anim, plot([op(ParEllipse(P, t, ':-alpha'=alpha,':-scaling'=scaling )), t = -2*Pi .. 2*Pi], color=blue, view=[min(X)..max(X), min(Y)..max(Y)])]):

end:

Compare:= proc(P,N,{alpha:=1/2}, {scaling:="eval"}) local EParPlot, PFinPlot:
PFinPlot := IterPolygon(P,N, ':-alpha'=alpha, ':-scaling'=scaling):
EParPlot := plot([op(ParEllipse(P, t, ':-alpha'=alpha,':-scaling'=scaling )), t = -2*Pi .. 2*Pi], color=blue):
display(PFinPlot,EParPlot):
end:
####