#####    Some simple 2-D plots of geometric objects in R  #####

library(MASS)    # eqscplot(), lda(), etc.
library(mvtnorm) # rmvnorm(), dmvnorm(), etc.
oldpar<-par(no.readonly=T)

# First of all, we establish an empty plot that is big enough
# to hold all of the data.

# For this example, we pretend we have data in a matrix called X,
# so first we make a dummy matrix:

X <- matrix(c(0,5,0,10),ncol=2)

# The first question in plotting "geometric objects" is whether
# we want to try to preserve an isometric geometry, or whether
# we want a more efficient use of the display area.

# We have

plot(X,,type='n',
        xlim=c(min(X[,1]),max(X[,1])),
        ylim=c(min(X[,2]),max(X[,2])),
        xlab=expression(x[1]),ylab=expression(x[2]))

# or, to put the coordinate axes on the same scale,

plot(X,,type='n',
        xlim=c(min(X),max(X)),
        ylim=c(min(X),max(X)),
        xlab=expression(x[1]),ylab=expression(x[2]))
        
####### Lines

# To draw a line on the existing plot, we can use 
#   lines() for two points, p1 and p2, or  
#   abline() for a line given in intercept-slope form, (a,b).


# lines() draws a line through two points (or lines segments through points).
# abline() draws a line across the full plot surface.

# We can distinguish 4 tasks:
#   given two points,
#    1. draw a line to connect them: 
#         use lines() directly
#    2. draw a line across the full plot surface:
#         determine the slope and intercept and use abline();
#         if the slope is infinite, use abline(v= )
#   given slope and intercept,
#    3. draw a line across the full plot surface:
#         use abline() directly
#    4. draw a line segment (limited in some way):
#         determine relevant points and use lines()

# For the point setup, let p1=(x1,y1) and p2=(x2,y2).
# For task 1., 
#   the values are passed to lines() in the form (x1,x2), (y1,y2).
# For task 2.,
#   the slope is 
#     b = (y2-y1)/(x2-x1)
#   and the intercept is
#     a = y1 - b*x1

# Given the slope and intercept, a and b, and the limits 
# on x, x1<x2, the line segment goes from 
# p1=(x1,a+b*x1) to 
# p2=(x2,a+b*x2)

# Example: an equation in x, such as equation (63) in Duda et al.
# This equation is given in the form t(w)%*%(x-x0)=0.
# For x of order 2, in this notation, for 2-D plots, x2 takes the role
# of y above -- complicating the notation somewhat.  In any event, if
# w2 != 0, we have
# x2 = t(w)%*%x0/w2 - (w1/w2)*x1
# so 
#    a = t(w)%*%x0/w2
# and
#    b = w1/w2

# Examples
# A line segment:
lines(c(3,1),c(6,2))  # slope = 2

# An orthogonal line segment:
lines(c(3,7),c(6,4))  # slope = -1/2

# To try to make them look more orthogonal:

eqscplot(X,,type='n',
        xlim=c(min(X),max(X)),
        ylim=c(min(X),max(X)),
        xlab=expression(x[1]),ylab=expression(x[2]))
lines(c(3,1),c(6,2))
lines(c(3,7),c(6,4))

# To do the same thing with abline()
eqscplot(X,,type='n',
        xlim=c(min(X),max(X)),
        ylim=c(min(X),max(X)),
        xlab=expression(x[1]),ylab=expression(x[2]))
abline(6-2*3,2)
abline(6+.5*3,-.5)

# To do a line and a line segment
eqscplot(X,,type='n',
        xlim=c(min(X),max(X)),
        ylim=c(min(X),max(X)),
        xlab=expression(x[1]),ylab=expression(x[2]))
abline(6-2*3,2)
text(5.1,9.9,pos=4,'a:abline(6-2*3,2)',cex=.75)
lines(c(3,7),c(6,4))
text(6.1,4.9,pos=4,'b:lines(c(3,7),c(6,4))',cex=.75)


####### Vectors

# Vectors and matrices are our basic mathematical objects.

# As a graphical object in the ordinary coordinate system, a mathematical 
# vector can be thought of as a directed line segment from the origin to
# some specified point.  A 2-D object of this type complete with an 
# arrowhead can be plotted by arrows(). 
# The size of the arrowhead is usually too large by default.  
# It is controlled by length

# Example, for x = (6,2):
arrows(0,0,6,2,length=0.1)  
text(5.8,2.1,pos=4,'c:arrows(0,0,6,2)',cex=.75)

# Note the difference in the ordering of the arguments in
# arrows() and lines()

# arrows() is not really very useful for plotting vectors, except just
# for the arrowhead.
# arrows() is best thought of just as a function of convenience for
# that specific purpose. 

# Often we want to plot a line or line segment through a given point
# in the direction of a given vector.
# Let the point be p=(x1,y1), and the vector be w=(w1,w2).
# If w1!=0, we can use abline(y1-(w2/w1)*x1), w2/w1)

# Example, for p=(3,6) and w = (6,2):
abline(5,1/3)
text(8,7.2,pos=4,'d:abline(5,1/3)',cex=.75)

# Through the origin, of course, it's just abline(0,w[2]/w[1])

####### Angles

# The cosine of the angle between a given vector x1=(x11,x12) and the vector 
# x2=(x21,0), that is, between the given vector and the horizontal axis, 
# is   x11/sqrt(x11^2+x12^2).
# This is simple trignometry.

# Using the formulas for the sine and cosine of the sum of two angles, we find
# that the cosine of the angle between a vector x1=(x11,x12) and an arbitrary
# vector x2=(x21,x22) is
# (x11*x21 + x12*x22)/sqrt((x11^2+x12^2)*(x21^2+x22^2))

# From this, we can also get very easily the cosine of the angle between a
# line with slope b1 and a line with slope b2:
# (1 + b1*b2)/sqrt((1+b1^2)(1+b2^2))

####### Rotations

# A (column) vector x is rotated through a counterclockwise angle t 
# by multiplication by the matrix
# [ cos(t)  -sin(t) ]
# [ sin(t)   cos(t) ]

# If this rotation matrix is Q, then y=Qx is the rotated vector.

# The rotation of a vector is about the origin.

# A line segment from p1=(x1,y1) to p2=(x2,y2) can be rotated counterclockwise
# through the angle t by 
# Q*(p2-p1)+p1

# Example:
# Rotate a line segment through 60 degrees.
c <- 0.5
s <- sqrt(1-c^2)
Q <- matrix(c(c,s,-s,c),nrow=2)
p1 <- c(2,2)
p2 <- c(6,4)
lines(c(p1[1],p2[1]),c(p1[2],p2[2]))
text(6,4,pos=4,'e',cex=.75)
p2new <- Q%*%(p2-p1)+p1
lines(c(p1[1],p2new[1]),c(p1[2],p2new[2]))
text(2.2,6.3,pos=4,'f',cex=.75)

# The rotation of a line with slope b through a counterclockwise angle t 
# yields a line with slope
# (2*b - sqrt(4*b^2-4*(cos^2(t)-b^2*sin^2(t))*(b^2*cos^2(t)-sin^2(t))))/
#   (2*(cos^2(t)-b^2*sin^2(t)))

####### Rotation of a data matrix

# The observation vectors in a data matrix are the rows.
# Therefore, the rotation of the data matrix X by the rotation matrix Q is
# X%*%t(Q)

# Example:
# Rotate a data set through 60 degrees.
set.seed(3) ###################
n <- 100
d <- 2     ### this would work for any dimension; we just couldn't plot it easily
Sigma <- matrix(c(1,1,1,4),nrow=2)
X <- rmvnorm(n,rep(0,d),Sigma)
plot(X,,type='n',
        xlim=c(min(X),max(X)),
        ylim=c(min(X),max(X)),
        xlab=expression(x[1]),ylab=expression(x[2]))
points(X,col='blue',pch=16,cex=0.5)
Xrot <- X%*%t(Q)
points(Xrot,col='magenta',pch=16,cex=0.5)

# A rotation of a geometric object such as a line or a vector is equivalent
# to rotating the coordinate system in the opposite direction.

####### Projections

# A projection carries a geometric object onto a similar object in a 
# lower-dimensional manifold.

# In a vector space, by definition, a projection is a transformation by a
# symmetric idempotent matrix.  (We might exclude the identity transformation,
# so that the dimension of the projection is actually lower.)

# The rank of the projection matrix is the dimension of the space being projected onto.

# The simplest projection is of an arbitrary n-vector onto a given vector.
# An arbitrary n-vector has a dimension of n; a given vector has a dimension of 1.

# The picture that represents the projection is a right triangle in which the
# hypotenuse is the vector being projected.

# To project the vector y onto the vector x, the projection matrix is
#     x%*%t(x)/t(x)%*%x

# Example:
X <- matrix(c(0,5,0,10),ncol=2)
eqscplot(X,,type='n',
        xlim=c(min(X),max(X)),
        ylim=c(min(X),max(X)),
        xlab=expression(x[1]),ylab=expression(x[2]))

x<-c(7,4)
y<-c(2,7)
arrows(0,0,x[1],x[2],length=0.1)  
text(x[1]+.1,x[2]+.1,pos=4,'x',cex=.75)
arrows(0,0,y[1],y[2],length=0.1)  
text(y[1]+.1,y[2]+.1,pos=4,'y',cex=.75)

# R does not demote t(x)%*%x in the statement
#   yproj<-(x%*%t(x)/t(x)%*%x)%*%y
# so we will rewrite it:
yproj<-(x%*%t(x)/sum(x*x))%*%y
arrows(0,0,yproj[1],yproj[2],length=0.1)  
lines(c(y[1],yproj[1]),c(y[2],yproj[2]))  
text(yproj[1],yproj[2]-.2,pos=4,expression(tilde(y)),cex=.75)

# This projection is the regression of y on x (with no intercept).

# The mean of y is the projection of y onto the 1 vector.

# Question: why is the projection a rank 1 transformation, while the rotation is
# a rank n transformation?

# The projection of a vector onto a vector space defined by the columns of a
# full-column-rank matrix X is defined similarly by the projection matrix:
#     X%*%inv(t(X)%*%X)%*%t(X)
# What is this???

####### Projection of a data matrix onto a line through the origin.

# In either the rotation or the projection of a data matrix, we first recognize
# that the vectors that we rotate or project are the rows of the data matrix.
# All this means is that we form the transpose of the vector-matrix product.
# This of course means postmultiplication by the transpose.  Because a projection
# matrix, unlike a rotation matrix, is symmetric, it's just postmultiplication.

# If v is a vector in the direction of the line, the projection of the data
# matrix X onto the line is just
#         X%*%(v%*%t(v))/t(v)%*%v
        
# Example.  Project X onto its first principle component direction:
set.seed(3) ###################
n <- 100
d <- 2     ### this would work for any dimension; we just couldn't plot it easily
Sigma <- matrix(c(1,1,1,4),nrow=2)
X <- rmvnorm(n,rep(0,d),Sigma)
plot(X,,type='n',
        xlim=c(min(X),max(X)),
        ylim=c(min(X),max(X)),
        xlab=expression(x[1]),ylab=expression(x[2]))
points(X,col='blue',pch=16,cex=0.5)

v1<-princomp(X)$loadings[,1]
# if(v1[1]!=0)
# draw line through mean
b <- v1[2]/v1[1]
a <- mean(X[,2]) - b*mean(X[,1])
abline(a,b)

##### run

z1 <- X%*%(v1%*%t(v1))/sum(v1*v1)
hist(z1,xlab=expression(z[1]))

v2<-princomp(X)$loadings[,2]
# if(v1[1]!=0)
# draw line through mean
b <- v2[2]/v2[1]
a <- mean(X[,2]) - b*mean(X[,1])
abline(a,b)

##### run

z2 <- X%*%(v2%*%t(v2))/sum(v2*v2)
hist(z2,xlab=expression(z[2]))