
x = np.mat(np.arange(-1.,1.,0.01)).T
N = len(x)

degree = 10
#A = np.hstack((np.power(x,0), np.power(x,1), np.power(x,2)))
B = np.hstack((np.power(x,i) for i in range(degree+1)))
B = B[:,np.random.permutation(degree+1)]
#B = np.hstack((np.power(x,i) for i in range(degree+1)))
#B = np.random.randn(N,degree+1)
A, R, = np.linalg.qr(B)
# Append an extra basis for outliers

K = A.shape[1]
plt.figure(figsize=(5,20))
#plt.show()
idx = range(0,N,20)
idx.append(N-1)
for i in range(K):
    plt.subplot(K,1,i+1)
    plt.stem(x[idx],A[idx,i])
    plt.plot(x,A[:,i])
    plt.gcf().gca().set_xlim([-1.2, 1.2])
    plt.gcf().gca().set_ylim([-0.3, 0.3])
    plt.gcf().gca().axis('off')

plt.show()

"""
Explanation: Given a $p \geq 1$, a $p$-norm ball is the set of points 
$$
B_p = { x : \| x \|_p \leq 1 }
$$
Design an algorithm to sample uniform points in a given norm-ball $B_p$, where $x \in \mathbb{R}^m$.
Let $C = AB$. Hence, $c_{i,j} = a_i^\top b_j$ where $a_i^\top$ is the $i$'th row of $A$ and $b_j$ is the $j$'th column of $B$
\begin{eqnarray}
\| A B \|F^2 & = & \sum{i} \sum_j |c_{i,j} |^2 \
& = & \sum_{i} \sum_j |a^\top_{i} b_j |^2 \
& \leq & \sum_{i} \sum_j \|a_i\|^2 \|b_j\|^2  \
& = & \left(\sum_{i} \|a_i\|^2 \right) \left(\sum_j  \|b_j\|^2 \right)\
& = & \| A \|_F^2 \| B \|_F^2 \
\end{eqnarray}
Hence 
$$
\| A B \|_F \leq  \| A \|_F \| B \|_F
$$
QR for polynomials
End of explanation
"""


import numpy as np
import numpy.linalg as la

np.set_printoptions(precision=3, suppress=True)

def house(x):
    '''Householder Reflection'''
    m = len(x)
    e = np.mat(np.zeros((m,1)))
    e[0] = 1
    v = la.norm(x)*e - x
    z = v/la.norm(v)
    H = np.eye(m) - 2*z*z.T
    return z,H

x = np.mat('[1;2;3;4]')
z,H = house(x)

print(x)

print(H*x)

print(x - 2*z*(z.T*x))


np.set_printoptions(precision=3, suppress=True)

m = 5
n = 7
A = np.mat(np.random.randn(m, n))


R = A.copy()
Q = np.eye(m)
print('A')
print(A)
for i in range(np.min((n,m))-1):
    z,H = house(R[i:,i])
    R[i:,i:] = H*R[i:,i:]
    Q[:,i:] = Q[:,i:]*H.T
    #print('q =')
    #print(q)
    #print('Q =')
    #print(Q)
    print('R =')
    print(R)

print('--'*4)
print(Q*R)
print('--'*4)
print(A)

"""
Explanation: Householder Reflection
Linearly map a vector $x$ onto the coordinate axis without changing its length.
\begin{eqnarray}
H \left(
\begin{array}{c}
x_0 \ x_1 \ \vdots \ x_{N-1} 
\end{array}
\right)
=
\left(
\begin{array}{c}
\pm \|x\|  \ 0  \ \vdots \ 0 
\end{array}
\right) 
\end{eqnarray}
A reflector is like a projection onto some hyperplane but we go twice further. 
$$
H = I - 2P
$$
The projection is onto the hyperplane orthogonal to $v$
$$
P = I - \frac{v}{\|v\|} \frac{v^\top}{\|v\|} = I - \frac{v v^\top}{v^\top v}
$$
where, $v$ is the difference vector 
$
v = \|x\|e_1 - x
$. 
$$
H = I - 2I + 2qq' = -I + 2qq' = -(I - 2qq')
$$
If $P$ is an orthogonal projector, the reflector $I-2P$ is an orthonormal transformation as we have $P^\top P = P^2 = P$ for an orthogonal projector.
$(I-2P)^\top (I-2P) = I - 2P^\top -2P + 4P^\top P = I$
End of explanation
"""


def random_house(m, k=None):
    '''Generate a random Householder Reflection'''
    if k is None:
        k = m
        
    e = np.mat(np.zeros((k,1)))
    e[0] = 1
    x = np.mat(np.random.randn(k,1))
    v = la.norm(x)*e - x
    q = v/la.norm(v)
    q = np.vstack((np.zeros((m-k,1)), q))
    H = np.eye(m) - 2*q*q.T
    return q,H

q, H = random_house(5)

print(q)
np.linalg.matrix_rank(H)

N = 5
A = [random_house(N,k)[1] for k in range(N,1,-1)]
qq = np.hstack([random_house(N,k)[0] for k in range(N,1,-1)])
qq

A[:,0].T*A[:,1]

"""
Explanation: We have explicitely constructed $H$ (or $P$) but this is not needed. It is sufficient to store the Householder vectors $v$ to construct/implement transformations by $H$.
End of explanation
"""


