Learning to speak to computers!

Good to be back! Forgive me of my sins for I was occupied with interests of my own. 🙂

So apparently my girlfriend is now studying Linguistics and programming can become handy for large dataset analysis and automating some stuffs, so she (who never programmed anything ever) asked me to teach “programming”. I thought for a while and thought it would be fun to teach someone programming from scratch, think about it, it is a “language” too, albeit a “written” one, which is quite similar to the set of instructions that we are accustomed to process.

Someone writes to you,

Find the vowels in the sentence:   “Programming is simple conversation with the machine!”

The word “Find” is a trigger your brain processes to “narrow” down your “attention” to your prior knowledge of vowels, namely “AEIOU”. You can do the same thing to a computer, using its language! Now as humans have developed different language, you can talk to computer using various “languages” which it can understand.

BTW computers mother toungue is just True(1) and Flase(0).

I have been a great fan of the language “Python”, because of its simplicity and user friendliness, so I will stick to that. The main aim of these set of blog tutorials will be 2 fold. Recently I have been interested in NLP, which is Natural Language Processing, and its applications, these blog posts will apart from doing the basics of programming, move on to some basic language processing tools, which will be my own study journal. I might think to introduce some maths later on in this course for understanding some statistical properties of language.

So are you ready to talk to computer and explore the fascinating world?

Lets begin with the analogies, so that it isnt too “unfamiliar”. Your friend who speaks only German says to you :

Schreib “back”

Now I cant understand what Schreib means! because I dont know the language! Similarly a computer would not know any language, but you can teach him by installing a “software” ! I will write about Installing python in your computer in the next post, but till then hang on. So the basic thing to know speak to computer is to know its language. Specifically forget the adjectives there are only Names(nouns) and Commands(verbs) :p

BTW did you know ?

Guido van Rossum, the creator of the Python language, named the language after the BBC show “Monty Python’s Flying Circus”. He doesn’t particularly like snakes that kill animals for food by winding their long bodies around them and crushing them.

Now it is important that you learn a language by practicing it! So it applies same here, the next set of tutorials will be on my view about how one should see programming, well since i have said that, I believe that we are fast evolving into a world where humans and technology will be inseparable, and programming will be an intrinsic part of everyone’s life, there will be pet cyborgs and humanoids to make your life easier, so you better learn how to control them.

Ok, enough digression, so lets get back to our second example, Schreib in german means write, to do the same thing to a computer you can just say

print “HelloWorld”

and as you might have imagined, it does exactly as it is told 🙂 We will get to these details later. To have a formal understanding on how programming works, you must be really clear on what you want to do, and how exactly you want to do it!

This is the first and last thing one must worry about. Like if I said, “you Good  are ” (well its a odd fact that you could still make out the meaning, but believe me “computers” take things quite “Literally”), this doesnt make sense because  what exactly I want to ask is “Are you Good?” So structure is important in program too. Again in the first example, how would you count all vowels in a sentence?

1. recall what vowels are! {A,E,I,O,U} , wow !

2. Read the sentence “Programming is simple conversation with the machine!” letter by letter and note down vowels from it!

3.Naievely, you can also take a copy and write 1+1+1….. as soon as you encounter a vowel in the sentence.

4. To be sane enough, start reading from left to right!

Simple isnt it? You see how I broke a “Task” to smaller steps of what needs to be done. This is an important concept.

Now first time my teacher would teach me how to find vowels like this, but the next time she says me to find vowels, its a “Task” we should already know. This is the concept of a function. It is always a good idea to learn to define new short “verbs” for actions which you do often, say you work in a brick Kline, and the owner asks you “to collect mud, shape it and then bake it”  , would it be easier for him to “invent” a single word  for the whole process say “comisab” or something! and everyone just understands what the owner is saying.

So I will end this topic by summarizing 2 things what i think is important for writing a program:

1. Manually first do what you want to instruct to your computer in pen and paper. Its like a translation, first your need to know your mother language!

Did you Know? Programme comes from word Pro meaning before, and graphe meaning write, which probably sums up,what you need to do before writing your first code!Think about it before writing!

2.Try to break the process to as simple as possible and group them as small “tasks”.

Next Tutorial:

Installing python in Linux machine/Windows. Editors

Fast Fourier Transform (FFT) examples

The Discrete Fourier Transform(DFT) is defined by X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i 2 \pi \frac{k}{N} n} latex and the inverse fourier transform is defined as x_n = \frac{1}{N} \sum_{k=0}^{N-1} X_k \cdot e^{+i 2 \pi \frac{k}{N} n} latex

In [2]:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import time

'''
	Define two functions to compute fourier transforms using two different methods
	1.	dft(x,inverse=False,timer=False): for INVERSE: inverse=true
						  for timed: timer=true
	2.	timedfft(x,inverse=False,timer):for INVERSE: inverse=true
						  for timed: timer=true
	3.	func(x): defines function at x
'''

def func(x):
    f=10.0 
    val=np.sin(2*np.pi*f*x)
    return val

'''
   	Computes the discrete fourier transform /inverse fourier transform
	using basic summation convention.
	OPTIONS:
		inverse: computes inverse fourier transform
		timer:	 prints time required for computation
'''
def dft(x, inverse = False, timer = False) :

    t = time.clock()							#get current time to calculate the computation time
    N = len(x)
    if inverse:
        inv=1
    else:
        inv=-1

    X =[0] * N								#make an array of size N

    for k in xrange(N) :						#compute the required sum
        for n in xrange(N) :
            X[k] += x[n] * np.e**(inv * 2j * np.pi * k * n / N)
        if inverse :
            X[k] /= N

    t = time.clock() - t						#get delta time
    if timer :
        print "Computed DFT of size",N,"in",t,"sec."
    return X

'''
	fourier transform using built in fast fourier method
'''

def timedfft(x,inverse=False,timer=False):

	t=time.clock()
	if not inverse:X=np.fft.fft(x)
	else:	X=np.fft.ifft(x)
	t=time.clock()-t
	if timer: print "Computed DFT using fft of size",len(x),"in",t,"sec."
	return X

'''
	Computing part: 1. Define a function over a interval
			2. Compute fourier transform using the two methods
			3. Compare time
'''

#fix sample size, this gives the the rate of sampling of the function
sample_size=2**8
L=5.0
x=np.linspace(0,L,sample_size)

#there are four plots 1.original signal 2.fourier transformed signal 3.fast fourier transformed signal 4.Reconstructed signal

plt.subplot(411)
plt.xlabel('t')
plt.ylabel('f($t$)')
plt.xlim([0,L/10])
signal=func(x)
plt.plot(x,signal,label='original signal')

'''
				In [20]: run fourier_basic.py
				Computed DFT using fft of size 256 in 0.0 sec.
				Computed DFT of size 256 in 0.48 sec.
				Computed DFT of size 256 in 0.16 sec.

				In [21]: run fourier_basic.py
				Computed DFT using fft of size 4096 in 0.0 sec.
				Computed DFT of size 4096 in 123.94 sec.
				Computed DFT of size 4096 in 43.18 sec.
'''
#compute timed fft
fourierfft=timedfft(signal,timer=True)

#caluclate discrete fourier transform using normal method
fourier =dft(signal,timer=True)

#get the frequencies
n=x.size
d=x[1]-x[0]
freqsfft=np.fft.fftfreq(n,d)

#get positive index
indices=np.where(freqsfft>=0)

#plot
plt.subplot(412)
plt.xlabel('freq')
plt.ylabel('f($\omega$)')
plt.plot(freqsfft,fourierfft,label='FFT')

plt.subplot(413)
plt.xlabel('freq')
plt.ylabel('f($\omega$)')
plt.plot(freqsfft,fourier,color='green',label='Without FFT')

#inverse fourier transform using dft
invfourier=(dft(fourier,inverse=True,timer=True))

plt.subplot(414)
plt.xlabel('t')
plt.ylabel('f($t$)')
plt.xlim([0,L/10])
plt.plot(x,invfourier,'+-',color='red')
display(gcf())
Computed DFT using fft of size 256 in 0.0 sec.
Computed DFT of size 256 in 0.48 sec.
Computed DFT of size 256 in 0.16 sec.

/numpy/core/numeric.py:235: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

N.B.: The fourier transform of the function is symmertric about 0.
Thus it suffices to display only the positive half,
as the negetive half gives us no more information

In [3]:
clf()
plt.subplot(111)
pi=np.pi
sampled=5*np.sin(2*pi*x)+np.sin(2*pi*10*x)
fourier=np.fft.fft(sampled)
#if you want to print only the positive frequencies
index=np.where(freqsfft>=0)
plt.plot(freqsfft[index],fourier[index])
display(gcf())

In [39]:
clf()
sampled=np.exp(-x**2)
fourier=np.fft.fft(sampled)
plt.plot(freqsfft,fourier)
display(gcf())

In [40]:
clf()
from scipy.special import jn
sampled=jn(2,x)
fourier=np.fft.fft(sampled)
plt.plot(freqsfft,fourier)
display(gcf())

2D Fast Poisson Solver

2D Fast Poisson Solver:

Introduction:

Here we will construct poisson solver for a 2D square boundary using a variation
fast fourier transform , dicrete sine transform.There are 4 types of discrete sine transform,you can read more about them in
wikipedia http://en.wikipedia.org/wiki/Discrete_sine_transform

We will solve the 2D electrostatic equation for square boundary,in

\epsilon_0=1 units

The poisson equation solver in 2D basically uses the matrix representation of the discrete poisson equation.
The poisson equation \nabla^2 u= \rho can be solved using fourier transform technique.The algorithm is basically
from
http://www.math.nus.edu.sg/~bao/teach/ma5233/lect10.ppt

This uses dst function(Type -I) from either fft.fftpack or a seperate
dst package I downloaded from internet called transforms

CODE:

In [2]:
from __future__ import division
import matplotlib.pyplot as plt
import numpy as np
from numpy.fft import fft,fft2,ifft2,ifft,irfft2,rfft2
import random as random
from mpl_toolkits.mplot3d import Axes3D

Define a function for discrete sin transform using fft algorithm,by symmetrizing the function x

In [3]:
def dst(x,axis=-1):
    """Discrete Sine Transform (DST-I)

    Implemented using 2(N+1)-point FFT
    xsym = r_[0,x,0,-x[::-1]]
    DST = (-imag(fft(xsym))/2)[1:(N+1)]

    adjusted to work over an arbitrary axis for entire n-dim array
    """
    n = len(x.shape)
    N = x.shape[axis]
    slices = [None]*3
    for k in range(3):
        slices[k] = []
        for j in range(n):
            slices[k].append(slice(None))
    newshape = list(x.shape)
    newshape[axis] = 2*(N+1)
    xsym = np.zeros(newshape,np.float)
    slices[0][axis] = slice(1,N+1)
    slices[1][axis] = slice(N+2,None)
    slices[2][axis] = slice(None,None,-1)
    for k in range(3):
        slices[k] = tuple(slices[k])
    xsym[slices[0]] = x
    xsym[slices[1]] = -x[slices[2]]
    DST = fft(xsym,axis=axis)
    #print xtilde
    return (-(DST.imag)/2)[slices[0]]

Define 2 dimensional DST, the idst is same as DST because of symmetry

In [4]:
def dst2(x,axes=(-1,-2)):
    return dst(dst(x,axis=axes[0]),axis=axes[1])

def idst2(x,axes=(-1,-2)):
    return dst(dst(x,axis=axes[0]),axis=axes[1])

Now define the poisson solver in 2D square boundary using the algorithm described above,this method can be
generalised to arbitary retangular boundary

In [5]:
def fft_poisson(f,h):

	m,n=f.shape
	f_bar=np.zeros([n,n])
	u_bar = f_bar			# make of the same shape
	u = u_bar

	f_bar=idst2(f)			# f_bar= fourier transform of f

	f_bar = f_bar * (2/n+1)**2  #Normalize
	#u_bar =np.zeros([n,n])
	pi=np.pi
	lam = np.arange(1,n+1)
	lam = -4/h**2 * (np.sin((lam*pi) / (2*(n + 1))))**2 #$compute $\lambda_x$
	#for rectangular domain add $lambda_y$
	for i in xrange(0,n):
		for j in xrange(0,n):
			u_bar[i,j] = (f_bar[i,j]) / (lam[i] + lam[j])

	u=dst2(u_bar)				#sine transform back
	u= u * (2/(n+1))**2			#normalize ,change for rectangular domain
	return u

Fix the initial parameters,domain of the equation etc.

In [6]:
#set bounds a,b,parameters
a = 0; b = 1;	
alpha=10				#alpha is grid points=2^alpha
n=2**alpha
L=b-a					#length of system

xe =np.linspace(a,b,n); 
ye =np.linspace(a,b,n);
x, y = np.meshgrid(xe,ye)
In [7]:
h=L/(n);			#size 
h2 = h**2;			#h squared
hx2=h2; 			#started with a cube,hx2=hy2
hy2=h2;
f=np.zeros([n,n]);

We put charges in the field.

Putting a random negetive charge in the field alog with a posituve charge

In [8]:
#initial conditions

#f[round(n/2):n,round(n/2):n]=1; #Initial condition
f[round((n+1)/2),round((n+1)/2-10)]=20
#f[round((n+1)/2),round((n+1)/2+10)]=-20
f[random.randint(0,n),random.randint(0,n)]=-10	#put a random charge
In [9]:
nx,ny=np.array(f.shape)-1		 #last index
U=np.zeros([n,n])

set boundary conditions

In [10]:
# BOUNDARY CONDITIONS
#set U(x,y)=g(x,y)
U[0,:]=0
U[nx,:]=0#5e-5
U[:,0]=0
U[:,ny]=0
In [11]:
##homogenize boundary condition

f[1,:]=f[1,:]+U[0,:]/hx2;
f[nx-1,:]=f[nx-1,:]+U[n-1,:]/hx2;

f[:,1]=f[:,1]+U[:,0]/hy2;
f[:,ny-1]=f[:,ny-1]+U[:,n-1]/hy2;

Solve the equation:

In [12]:
U=fft_poisson(f,h)
In [13]:
plt.figure()
#plt.imshow((U),cmap='hot')
plt.contour(U,50)
#plt.contour(f)

plt.show()
In [14]:
display(gcf())

In [15]:
imshow(U)
Out[15]:
<matplotlib.image.AxesImage at 0xb2b8a8c>
In [16]:
display(gcf())

 

In [ ]:

//