# coding: utf-8

# # Orthogonal Projection

# In[1]:

import numpy as np
import numpy.linalg as la


# In[2]:


# for in-line plots
get_ipython().magic(u'matplotlib inline')

# for plots in a window
# %matplotlib qt

import matplotlib.pyplot as pt
from mpl_toolkits.mplot3d import Axes3D


# Make two random 3D vectors:

# In[3]:

np.random.seed(13)
x = np.random.randn(3)
y = np.random.randn(3)


# Make them orthonormal:

# In[4]:

y = y - y.dot(x)/x.dot(x)*x
x = x/la.norm(x)
y = y/la.norm(y)


# Check:

# In[5]:

print(y.dot(x))
print(la.norm(x))
print(la.norm(y))


# Plot the two vectors:

# In[6]:

fig = pt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.set_xlim3d([-3, 3])
ax.set_ylim3d([-3, 3])
ax.set_zlim3d([-3, 3])

xy = np.array([x, y]).T
ax.quiver(
   0, 0, 0,
   xy[0], xy[1], xy[2],)


# Make an array with the cornerpoints of a cube:

# In[7]:

points = np.array([
    [-1,-1,-1],
    [-1,-1,1],
    [-1,1,-1],
    [-1,1,1],

    [1,-1,-1],
    [1,-1,1],
    [1,1,-1],
    [1,1,1],
])


# Plot them:

# In[8]:

fig = pt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:, 1], points[:, 2])


# Construct the projection matrix:

# In[9]:

Q = np.array([
   x,y,np.zeros(3)
]).T
print(Q)


# In[10]:

P = Q.dot(Q.T)
print(P)


# Check that $P^2=P$:

# In[11]:

la.norm(P.dot(P)-P)


# Project the points, assign to `proj_points`:

# In[12]:

proj_points = np.einsum("ij,nj->ni", P, points)


# In[13]:

fig = pt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0], points[:, 1], points[:, 2])
ax.scatter(proj_points[:,0], proj_points[:, 1], proj_points[:, 2], color="red")

xy = np.array([x, y]).T
ax.quiver(
   0, 0, 0,
   xy[0], xy[1], xy[2],)