from numpy import*

n = 100
s = zeros((2,n))

u = array([1.0,0.01])# seed for uniformly distributed pseudorandom generation
for i in range (0,n):
    u = ((8121*u + 28411)%134456)      
    s[:,i] = u/134456    

# Mix the data with a square matrix A, and plot
A = array([[1,2],[2,-1]])
x = dot(A,s)

from pylab import*

figure(1)
title('All three plots')

subplot(311)
scatter(s[0,:],s[1,:], c ='b')
title('source (s) points')
subplot(312)
scatter(x[0,:],x[1,:], c ='r')
title('mixed (x) points')

# find the averages of the x[0] and x[1] values
avg = sum(x, axis = 1)/n
for i in range(0,n):
    # subtract average from actual values to give zero mean
    x[:,i] -= avg

from numpy import* 
from numpy import linalg as LA

# perform PCA to diagonalize
Cx = cov(x, rowvar = 1)
eig, phi = LA.eig(Cx)
print '\neigenvalues of the covariance matrix of x:\n' + str(eig)

# Scale the eigenvectors so that the transformed data has unit variance
phi[:,0] /= eig[0]**(1.0/len(eig))
phi[:,1] /= eig[1]**(1.0/len(eig))

# M contains the eigenvectors of Cx (corresponding to nonzero eigenvalues) as columns
M = transpose(phi[:,fabs(eig)>1e-8])
# Transfrom the data points x
xt = dot(M,x)
Cxt = cov(xt)
# xt should have a diagonal covariance matrix, with no zero eigenvalues
print '\nCovariance matrix of transformed x:\n' + str(Cxt)

# plot transformed data:
subplot(313)
scatter(x[0,:],x[1,:], c ='r')
title('trasformed (zero mean and unit var) x')
show()

# Now, find independent components using log cosh contrast function

#STILL TO DO....