from __future__ import division
from matplotlib.patches import Circle
x1 = linspace(-1,1,10)
fig=figure()
ax=fig.add_subplot(111)
ax.plot(x1,(1-x1)/2)
ax.add_patch(Circle((0,0),1/sqrt(5),alpha=0.3))
ax.plot(1/5,2/5,'rs')
ax.axis('equal')
ax.set_xlabel('$x_1$',fontsize=24)
ax.set_ylabel('$x_2$',fontsize=24)
ax.grid()

from matplotlib.patches import Rectangle
import matplotlib.patches
import matplotlib.transforms

r=matplotlib.patches.RegularPolygon((0,0),4,1/2,pi/2,alpha=0.5)

fig=figure()
ax=fig.add_subplot(111)
ax.plot(x1,(1-x1)/2)
ax.plot(0,1/2,'rs')
ax.add_patch(r)
ax.grid()
ax.set_xlabel('$x_1$',fontsize=24)
ax.set_ylabel('$x_2$',fontsize=24)
ax.axis('equal')


from cvxopt import matrix as matrx # don't overrite numpy matrix class
from cvxopt import solvers

#t1,x1,t2,x2
c = matrx([1,0,1,0],(4,1),'d') 
G = matrx([   [-1,  -1, 0,  0],  #column-0
              [ 1,  -1, 0,  0],  #column-1
              [ 0,   0, -1,-1],  #column-2
              [ 0,   0,  1,-1],  #column-3
           ],(4,4),'d')

h = matrx([0,0,0,0],(4,1),'d') # (4,1) is 4-rows,1-column, 'd' is float type spec
A = matrx([0,1,0,2],(1,4),'d')
b = matrx([1],(1,1),'d')
sol = solvers.lp(c, G, h,A,b)
x1=sol['x'][1]
x2=sol['x'][3]
print 'x=%3.2f'% x1
print 'y=%3.2f'% x2

import scipy.linalg

def rearrange_G( x ): 
    'setup to put inequalities matrix with last 1/2 of elements as main variables'
    n=x.shape[0]
    return hstack([x[:,arange(0,n,2)+1], x[:,arange(0,n,2)]])

K=2 # components
Nf=128 # number of samples
M = 12 # > K log2(Nf/K); num of measurements
s=zeros((Nf,1)) # sparse vector we want to find
s[0]=1 # set the K nonzero entries
s[1]=0.5
np.random.seed(5489) # set random seed for reproducibility
Phi = matrix(randn(M,Nf)*sqrt(1/Nf)) # random Gaussian matrix
y=Phi*s # measurements

#-- setup L1 minimization problem -- 

# inequalities matrix with 
G=matrx(rearrange_G(scipy.linalg.block_diag(*[matrix([[-1,-1],[1,-1.0]]),]*Nf) ))
# objective function row-matrix
c=matrx(hstack([ones(Nf),zeros(Nf)]))
# RHS for inequalities
h = matrx([0.0,]*(Nf*2),(Nf*2,1),'d') 
# equality constraint matrix
A = matrx(hstack([Phi*0,Phi]))
# RHS for equality constraints 
b=matrx(y)

sol = solvers.lp(c, G, h,A,b)

#nonzero entries
nze= array(sol['x']).flatten()[:Nf].round(2).nonzero()
print array(sol['x'])[nze]

from cStringIO import StringIO
import sys

def L1_min(Phi,y,K):
    # inequalities matrix with 
    M,Nf = Phi.shape
    G=matrx(rearrange_G(scipy.linalg.block_diag(*[matrix([[-1,-1],[1,-1.0]]),]*Nf) ))
    # objective function row-matrix
    c=matrx(hstack([ones(Nf),zeros(Nf)]))
    # RHS for inequalities
    h = matrx([0.0,]*(Nf*2),(Nf*2,1),'d') 
    # equality constraint matrix
    A = matrx(hstack([Phi*0,Phi]))
    # RHS for equality constraints 
    b=matrx(y)
    # suppress standard output
    old_stdout = sys.stdout
    sys.stdout = mystdout = StringIO()
    sol = solvers.lp(c, G, h,A,b)
    # restore standard output
    sys.stdout = old_stdout
    sln = array(sol['x']).flatten()[:Nf].round(4)
    return sln



def dftmatrix(N=8): 
    'compute inverse DFT matrices'
    n = arange(N)
    U=matrix( exp(1j*2*pi/N*n*n[:,None] ))/sqrt(N)
    return matrix(U)

Nf=128
K=3 # components
M = 8 # > K log2(Nf/K); num of measurements
s=zeros((Nf,1)) # sparse vector we want to find
s[0]=1 # set the K nonzero entries
s[1]=0.5
s[-1] = 0.5 # symmetric to keep inverse Fourier transform real
Phi = dftmatrix(Nf)[:M,:] # take M-rows
y=Phi*s # measurements
# have to assert the type here on my hardware

sol=L1_min(Phi.real,y.real.astype(np.float64),K)

print np.allclose(s.flatten(),sol)

plot(sol)
plot(y.real)

def dftmatrix(N=8): 
    'compute inverse DFT matrices'
    n = arange(N)
    U=matrix( exp(1j*2*pi/N*n*n[:,None] ))/sqrt(N)
    return matrix(U)

def Q_rmatrix(Nf=8):
    'implements the reordering, adding, and stacking of the matrices above'
    Q_r=matrix(hstack([eye(Nf/2),eye(Nf/2)*0])
               +hstack([zeros((Nf/2,1)),fliplr(eye(Nf/2)),zeros((Nf/2,Nf/2-1))]))
    return Q_r

Nf=8
F=zeros((Nf,1)) # 8-point DFT
F[0]= 1 # DC-term, constant signal
n = arange(Nf/2)

ft = dftmatrix(Nf).H*F # this gives the constant signal

Q_r=Q_rmatrix(Nf)
U=dftmatrix(Nf/2) #half inverse DFT matrix
feven= U.real*Q_r*F # half the size
print np.allclose(feven,ft[::2]) # retrieved even-numbered samples



# let's try this with another sparse frequency-domain signal
F=zeros((Nf,1))  
F[1]=1
F[Nf-1]=1 # symmetric part
ft = dftmatrix(Nf).H*F # this gives the constant signal
feven= U.real*Q_r*F  # half the size
print np.allclose(feven,ft[::2]) # retrieved even-numbered samples

plot(arange(Nf),ft.real,arange(Nf)[::2],feven,'o')
xlabel('$t$',fontsize=22)
ylabel('$f(t)$',fontsize=22)
title('even-numbered samples')

Nf=32 # must be even
F=zeros((Nf,1))  

# set values and corresponding symmetry conditions

F[7]=1 
F[12]=0.5
F[9]=-0.25
F[Nf-9]=-0.25
F[Nf-12] = 0.5
F[Nf-7]=1 # symmetric part

Q_r=Q_rmatrix(Nf)
U=dftmatrix(Nf/2) #half inverse DFT matrix
ft = dftmatrix(Nf).H*F # this gives the constant signal
feven= U.real*Q_r*F  # half the size
print np.allclose(feven,ft[::2]) # retrieved even-numbered samples

plot(arange(Nf),ft.real,arange(Nf)[::2],feven,'o')
xlabel('$t$',fontsize=22)
ylabel('$f(t)$',fontsize=22)
title('even-numbered samples')

def rearrange_G( x ): 
    'setup to put inequalities matrix with first 1/2 of elements as main variables'
    n=x.shape[0]
    return hstack([x[:,arange(0,n,2)+1], x[:,arange(0,n,2)]])

K=2 # components
Nf=128 # number of samples
M = 18 # > K log(N); num of measurements

# setup signal DFT as F
F=zeros((Nf,1))  
F[1]=1
F[2]=0.5
F[Nf-1]=1 # symmetric parts
F[Nf-2]=0.5
ftime = dftmatrix(Nf).H*F # this gives the time-domain signal
ftime = ftime.real # it's real anyway
time_samples=[0, 2, 4, 12, 14, 16, 18, 24, 34, 36, 38, 40, 44, 46, 52, 56, 54,62]
half_indexed_time_samples = (array(time_samples)/2).astype(int)
Phi = dftmatrix(Nf/2).real*Q_rmatrix(Nf)
Phi_i = Phi[half_indexed_time_samples,:]

# inequalities matrix with 
G=matrx(rearrange_G(scipy.linalg.block_diag(*[matrix([[-1,-1],[1,-1.0]]),]*Nf) ))
# objective function row-matrix
c=matrx(hstack([zeros(Nf),ones(Nf)]))
# RHS for inequalities
h = matrx([0.0,]*(Nf*2),(Nf*2,1),'d') 
# equality constraint matrix
A = matrx(hstack([Phi_i,Phi_i*0]))
# RHS for equality constraints 
b=matrx(ftime[time_samples])

sol = solvers.lp(c, G, h,A,b)


import itertools as it

def dftmatrix(N=8): 
    'compute inverse DFT matrices'
    n = arange(N)
    U=matrix( exp(1j*2*pi/N*n*n[:,None] ))/sqrt(N)
    return matrix(U)

M = 3

np.random.seed(5489) # set random seed for reproducibility
Psi= dftmatrix(128)
Phi= randn(M,128)
s=zeros((128,1))
s[0]=1
s[10]=1

Theta = Phi*Psi
y = Theta*s

for i in it.combinations(range(128),2):
   sstar=zeros((128,1))
   sstar[array(i)]=1
   if np.allclose(Theta*sstar,y):
      break
else:
   print 'no solution'


%qtconsole