Vahid Moosavi

04 April 2017

# Introduction to Representation Learning 1¶

### To be discussed¶

• Signal Processing is a main part of Data Driven Modeling cases
• Harmonic Analysis, Fourier Series and Idealized Filters
• Fourier transfrom in 1D and 2D
• Introduction to Representation Learning
In [1]:
import warnings
warnings.filterwarnings("ignore")
import datetime
import pandas as pd
# import pandas.io.data
import numpy as np
import matplotlib
from matplotlib import pyplot as plt
import sys
import sompylib.sompy as SOM# from pandas import Series, DataFrame

from ipywidgets import interact, HTML, FloatSlider

%matplotlib inline
# from bokeh.models import CustomJS, ColumnDataSource, Slider,TextInput
# from bokeh.models import TapTool, CustomJS, ColumnDataSource
# from bokeh.layouts import column
# from bokeh.layouts import gridplot
# from bokeh.models import ColumnDataSource
# from bokeh.layouts import gridplot
# from bokeh.plotting import figure
# from bokeh.io import output_notebook, show
# from bokeh.plotting import Figure, output_file, show
# from bokeh.layouts import widgetbox
# from bokeh.models.widgets import Slider
# from bokeh.io import push_notebook
# from bokeh.io import show
# from bokeh.models import (
#     ColumnDataSource,
#     HoverTool,
#     LinearColorMapper,
#     BasicTicker,
#     PrintfTickFormatter,
#     ColorBar,
# )
# output_notebook()
# from bokeh.layouts import column
# from bokeh.models import CustomJS, ColumnDataSource, Slider
# from bokeh.plotting import Figure, output_file, show


# Signal Processing is a main part of Data Driven Modeling cases¶

• Usually observations are homogeneous
• Pixels of image
• Sequence of values in sound signal or any other time series
• Usually, we are looking for something on top of the observations
• Message of the sentence
• A certain pattern in an image such as a face in the picture
• All the methods are developed in a way to capture "invariances" within categories of interest
• translation invariance
• rotation invariance
• scale invariance
• Deformation
• Hierarchical Representation and Compositionality
• Text
• Video
• A Building
• Cities

## An example of translation¶

### If we naively compare two signals elementwise we get nothing!¶

In [2]:
N = 256
t = np.arange(N)
fi = .1*t
r = 1
z = np.cos(.1*t)
z1 = np.cos(.1*(t-30))
# .0*np.random.rand(N)
plt.plot(z)
plt.plot(z1)

Out[2]:
[<matplotlib.lines.Line2D at 0x11265ed10>]

# A more interactive version¶

https://bl.ocks.org/jinroh/7524988

# Discrete Fourier Transform (DFT)¶

https://en.wikipedia.org/wiki/DFT_matrix

### Euler's formula¶

It manages to encapsulate both cosine and sine function in a single complex exponentials. https://en.wikipedia.org/wiki/Euler%27s_formula#Using_power_series

In [3]:
#Complex Exponentials
N = 256
t = np.arange(N)
fi = .1*t
r = 1
cos = np.cos(.1*t)
sin = np.sin(.1*t)
ex = cos + 1j*sin
plt.plot(cos,linewidth=3)
plt.plot(np.real(ex),'.r')
plt.plot(np.imag(ex),'.g')

Out[3]:
[<matplotlib.lines.Line2D at 0x117b7aad0>]

# DFT¶

## $$X= Wx\$$¶

In [4]:
#Base Filters 1D
def base_filter_1D(N):
dd = -2*np.pi/N
w = np.exp(dd*1j)
#     w = np.exp(-2*np.pi*1j/N)

W = np.ones((N,N),dtype=complex)

for i in range(N):
for j in range(N):
W[i,j] = np.power(w,i*j)

#     W = W/np.sqrt(N)
return W

N = 4
W = base_filter_1D(N)
print np.around(W,decimals=3)

[[ 1.+0.j  1.+0.j  1.+0.j  1.+0.j]
[ 1.+0.j  0.-1.j -1.-0.j -0.+1.j]
[ 1.+0.j -1.-0.j  1.+0.j -1.-0.j]
[ 1.+0.j -0.+1.j -1.-0.j  0.-1.j]]

In [5]:
def DFT_1D_basis_vis(i):
fig = plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
#Real Part of the vector
plt.plot(np.real(W)[i],'r')
#Imaginary Part of the vector
plt.title('real part')
plt.axis('off')
plt.subplot(1,2,2)
plt.title('imaginary part')
plt.plot(np.imag(W)[i],'b')
plt.axis('off')
N = 256
W = base_filter_1D(N)
interact(DFT_1D_basis_vis,i=(0,N-1,1));

In [6]:
fig = plt.figure(figsize=(10,10))
for i in range(16):
plt.subplot(4,4,i+1)
plt.plot(np.real(W)[i],'r')
plt.plot(np.imag(W)[i],'b')
plt.axis('off')

In [7]:
# They are orthogonal basis
#Each row of W is a basis

#Orthogonal Basis
a = np.imag(W)[2].dot(np.imag(W)[1])
print np.around(a)

#Sin and Cos
a = np.real(W)[2].dot(np.imag(W)[1])
print np.around(a)

-0.0
-0.0


## Now from Time domain to Frequency domain¶

In [8]:
N = 128
t = np.arange(N)
x = 1*np.sin(.05*t) + 1*np.cos(5*t+.1)
plt.plot(x)

Out[8]:
[<matplotlib.lines.Line2D at 0x1189a0fd0>]
In [9]:
#DFT 1D
# Now we can transform this original signals to a new vector space based on fourier base patterns
# How to transform the original signal (DATA) to a new basis system
W = base_filter_1D(N)
X1 = W.dot(x)
# plt.plot(X,'.r')
plt.subplot(3,1,1)
plt.plot(x,'.-');
plt.title('time domain');
plt.subplot(3,1,2)
plt.plot(np.real(X1),'.');
plt.title('frequencey domain real part')
plt.subplot(3,1,3)
plt.plot(np.imag(X1),'.');
plt.title('frequencey domain imaginary part')
plt.tight_layout()


## Fast Fourier Transform (FFT) in Numpy¶

In [10]:
plt.plot(np.abs(np.real(X1)),'or');

X = np.fft.fft(x)
plt.plot(np.abs(np.real(X)),'*b',markersize=3);


## How does it deal with translation?¶

### we have two similar signals with a different in time¶

In [11]:
# Two signals have the same frequencies but with a shift in time
N = 256
t = np.arange(N)
x1 = .4*np.sin(1*t+.1) + .6*np.cos(-15*t+.1) + .3*np.random.rand(N)
x2 = .4*np.sin(1*(t-2)+.1) + .6*np.cos(-15*(t-16)+.1) + .1*np.random.rand(N)
plt.plot(x1)
plt.plot(x2)

Out[11]:
[<matplotlib.lines.Line2D at 0x119273cd0>]
In [12]:
W = base_filter_1D(N)
X1 = W.dot(x1)
X2 = W.dot(x2)

In [13]:
plt.plot(np.abs(np.absolute(X1)),'b');
plt.plot(np.abs(np.absolute(X2)),'g');


## Reconstruction¶

### From freq to time¶

In [14]:
N = 128
t = np.arange(N)
x = 1*np.sin(5*t)

W = base_filter_1D(N)
X = W.dot(x)
W_inv = np.linalg.inv(W)

def reconstruct_signal(howmanybasis=1):
x_r = W_inv[:,:howmanybasis].dot(X[:howmanybasis])
plt.plot(x[:],'b');
plt.plot(x_r[:],'r--');
plt.ylim(-1,1)
interact(reconstruct_signal,howmanybasis=(0,N-1,1));


## Denoising¶

### Note that there of course many other methods for signal smoothing and denoising¶

In [15]:
N = 128
t = np.arange(N)
x = 1*np.sin(.05*t) + 1*np.cos(.15*t+.1)  + .83*np.random.rand(N)

plt.plot(x[:],'b');

In [16]:
X  = np.fft.fft(x)
# plt.hist(X,bins=100);
plt.plot(np.absolute(X),'.');

In [17]:
def denoise(howmany=1):
X  = np.fft.fft(x)
X[np.argsort(-1*np.absolute(X))[howmany:]]=0
plt.plot(np.fft.ifft(X),'r',linewidth=3)
plt.plot(x)
interact(denoise,howmany=(0,N,1));


# 2D Fourier Transform¶

http://www.robots.ox.ac.uk/~az/lectures/ia/lect2.pdf

In [18]:
N = 32
W1D = base_filter_1D(N)
def base_filter2d_vis(u=1,v=1):

r = W1D[u][np.newaxis,:].T
c = W1D[v][np.newaxis,:]
W2 = r.dot(c)
fig = plt.figure(figsize=(15,7))
plt.subplot(1,2,1)
plt.title('Real Part(CoSine Wave)')
plt.imshow(np.real(W2),cmap=plt.cm.gray)
plt.axis('off')
plt.subplot(1,2,2)
plt.title('Imaginary Part (Sine Wave)')
plt.axis('off')
plt.imshow(np.imag(W2),cmap=plt.cm.gray)

In [19]:
interact(base_filter2d_vis,u=(0,N-1,1),v=(0,N-1,1));

In [20]:
#Base Filters 2D
def base_filter_2D(N):
W1D = base_filter_1D(N)*np.sqrt(N)
W2D = np.ones((N,N,N,N),dtype=complex)

for u in range(0,N):
for v in range(0,N):
r = W1D[u][np.newaxis,:].T
c = W1D[v][np.newaxis,:]
W2D[u,v,:,:] = r.dot(c)
W2D = W2D/(np.sqrt(N*N))
return W2D


## Base Filters Dictionary¶

• Note that these base filters are defined independent of the data!
In [21]:
N = 8
W2D = base_filter_2D(N)
print W2D.shape
fig = plt.figure(figsize=(7,7))
k =1
for u in range(0,N):
for v in range(0,N):
W2 = W2D[u,v,:,:]
plt.subplot(N,N,k)
plt.imshow(np.real(W2),cmap=plt.cm.gray)
#         plt.imshow(np.imag(W2),cmap=plt.cm.gray)
k = k +1
plt.axis('off')

(8, 8, 8, 8)

In [22]:
#DFT  in 2D is much more time consuming.
# Each cell is calculate as the sum of the dot product of base filter corresponding of that cell and the whole image
## Very slow: Don't use it :)
def DFT2D(W2D,img):
N = img.shape[0]
M = img.shape[1]
img_fr = np.ones((N,N),dtype=complex)
for i in range(N):
for j in range(M):
img_fr[i,j]= W2D[i,j,:,:].dot(img).sum()
return img_fr

def InverseDFT2D(howmanybasis=1):
img_rec = img
for u in range(howmanybasis):
for v in range(howmanybasis):
W = W2D[u,v,:]
W_inv = np.linalg.inv(W)
for i in range(N):
for j in range(M):
img_rec[i,j]= W_inv.dot(IMG[:howmanybasis,:howmanybasis]).sum()


## Low Pass High Pass filters¶

In [35]:
from skimage import data
from skimage.color import rgb2gray

fig = plt.figure(figsize=(9,6))

img = data.astronaut()
img = rgb2gray(img)
img = img
# + np.random.rand(img.shape[0],img.shape[1])

N = img.shape[0]
img = np.zeros((N, N))

x = 60
y = 60

for i in range(20):
indx = np.random.randint(x,N-x)
indy = np.random.randint(y,N-y)
img[indx:indx+x,indy:indy+y] = 1

plt.subplot(2,3,1)

plt.imshow(img,cmap=plt.cm.gray);
plt.axis('off')

F_img = np.fft.fft2(img)
F_img = np.fft.fftshift(F_img)
plt.subplot(2,3,4);
plt.imshow(np.log(np.absolute(F_img)));
plt.axis('off');

Nx = img.shape[0]
Ny = img.shape[1]

x  =30
y = 30
High_pass_freq =  F_img.copy()
High_pass_freq[Nx/2-x:Nx/2+x,Ny/2-y:Ny/2+y] = 0

High_pass_img = np.fft.ifft2(High_pass_freq)
plt.subplot(2,3,2);
plt.title('High Pass Filter');
plt.imshow(np.absolute(High_pass_img),cmap=plt.cm.gray);
plt.axis('off');

plt.subplot(2,3,5);
plt.imshow(np.log(np.absolute(High_pass_freq)));
plt.axis('off');

Low_pass_freq =  F_img.copy()
Low_pass_freq[Nx/2+x:,:] = 0
Low_pass_freq[:,:Ny/2-y] = 0
Low_pass_freq[:Nx/2-x,:] = 0
Low_pass_freq[:,Ny/2+y:] = 0

Low_pass_img = np.fft.ifft2(Low_pass_freq)
plt.subplot(2,3,3);
plt.title('Low Pass Filter');
plt.imshow(np.absolute(Low_pass_img),cmap=plt.cm.gray);
plt.axis('off');

plt.subplot(2,3,6);

plt.imshow(np.log(np.absolute(Low_pass_freq)));
plt.axis('off');