Comparison of iterative and recursive way of Fibonacci Sequence calculation

The Fibonacci numbers are the sequence of numbers

$\{F_n\}_{n=1}^\infty$

defined by the linear recurrence equation

$$F_n = F_{n-1} + F_{n-2}$$

with $F_1 = F_2 = 1$, and conventionally defining $F_0 = 0$.

We can implement Fibonacci numbers by iteratively or by using recursion. These two functions are implemented in Matlab to compare the running time of the algorithms.

function [result] = fibonacciRecursive(n)
%calculates the fibonacci output in recursive way
    if n==0 || n==1
        result=1;
    else
        result=fibonacciRecursive(n-1)+ fibonacciRecursive(n-2);
    end
end

function [result] = fibonacciIterative( n )
%calculates the fibonacci output in iterative way
    y1=1;
    y2=2;
    for i=2:n
        result=y1+y2;
        y2=y1;
        y1=result;
    end
end

For testing:

N=10;

fibonacciIter=zeros(N,1);
fibonacciRecur=zeros(N,1);

for i=1:N
     tic
        for j=1:25
            fibonacciIterative(i);
        end
     
      fibonacciIter(i)=toc;
     
      tic
        for j=1:25
            fibonacciRecursive(i);
        end
      fibonacciRecur(i)=toc;
end
X=1:N;
plot(X,fibonacciIter,X,fibonacciRecur);
legend('Iterative', 'Recursive')
xlabel('input size');
ylabel('time(seconds)');

The plot is:

As shown in plot recursive implementation is increasing exponentially in time with the given input size

Implementing Principle Component Analysis (PCA) in Python

i take a look at PCA (principle component analysis). i'm not sure this is implemented somewhere else but a quick review of my collage notes (reference needed) lead me the code below, and  data is (reference needed):

x y
2.5 2.4
0.5 0.7
2.2 2.9
1.9 2.2
3.1 3.0
2.3 2.7
2 1.6
1 1.1
1.5 1.6
1.1 0.9

'''
 *@author beck 
 *@date Sep 14, 2012
 *PCA with Python 
 *bekoc.blogspot.com 
'''
import numpy as np
import matplotlib.pyplot as plt
import pylab


xs= np.loadtxt("pcaData",delimiter=" ", skiprows=1, usecols=(0,1)) # numpy array - similar to C array notation.
#get mean
meanx=np.average(xs[:,0])
meany=np.average(xs[:,1])

correctedX=[value-meanx for value in (xs[:,0])] #X data with the means subtracted
correctedY=[value-meany for value in (xs[:,1])] #Y data with the means subtracted
data= np.array([correctedX,correctedY])
print data.shape
covData=np.cov(data)#calculate covariance matrix

eigenvalues, eigenvectors = np.linalg.eig(covData)

print eigenvectors
print eigenvectors[0][0] #eigenvectors are both unit eigenvectors
print eigenvectors[1][0]
x= [n for n in range (-2,3)]
y=  [eigenvectors[1][0]*i/eigenvectors[0][0] for i in x ] 
y1=  [eigenvectors[1][1]*i/eigenvectors[0][1] for i in x ] 

print x
print y 
plt.plot(x, y,linestyle='--', label='eigenvector1')
plt.plot(x, y1, linestyle='--', label='eigenvector2')
plt.plot(data[0,:],data[1,:], marker='+', linestyle=' ',  label= "Normalized data" )

#plt.plot(xs[:,0],xs[:,1],marker='+',linestyle=' ')
pylab.ylim([-2,2])
pylab.xlim([-2,2])
plt.title('PCA example')
plt.legend()
plt.show()

The code includes step 1 to 5

 PCA summary :
1- Given a dataset calculate normalized data (mean substructed data), let's say n dimension (feature) data
2-calculate covariance matrix of normalized data
3-calculate eigenvalues and eigenvectors of the covariance matrix
4-eigenvector with the largest eigenvalue is the principal component
5-choose p eigenvectors and multiply with your data
6-now your data is p dimension.

The green dotted plot of the eigenvector shows the most significant relation between dimensions

Please refer to simple and consise tutorial at georgemdallas blog 

how to implement two stacks in one array

My friend came up with a new question  that he is asked in his interview,

Explain how to implement two stacks in one array A[1..n] in such a way that neither stack overflows unless the total number of elements in both stacks together is n. The PUSH and POP operations should run in O(1) time.

The code is implemented in C++.

/* 
 * File:   main.cpp
 * Author: bekoC
 */

#include <cstdlib>
#include<stdio.h>
#include <stdlib.h>
#include<iostream>

#define N 10 

struct Stacks {
    char charArray[N];
    int top1, top2;
    void olustur();
    bool Insert(int stackNumber, char veri);
    char deleteFromStack(int stackNumber);
    bool isEmpty(int stackNumber);
    void printStack(int stackNumber);
};

bool Stacks::Insert(int stackNumber, char newChar) {
    if (top1 == top2)
    {
        return false;
    }
    if (stackNumber == 1)
        charArray[top1++] = newChar;
    else
        charArray[top2--] = newChar;
    return true;
}

char Stacks::deleteFromStack(int stackNumber) {
    if (stackNumber == 1)
        return charArray[--top1];
    else
        return charArray[++top2];
}

bool Stacks::isEmpty(int stackNumber) {
    if (stackNumber == 1)
        return (top1 == 0);
    else
        return (top2 == N - 1);
}

void Stacks::create() {
    top1 = 0; 
    top2 = N - 1; 
}

void Stacks::printStack(int stackNumber) {

    if (stackNumber == 1 && !isEmpty(stackNumber)) {
        {
            for (int i = 0; i <= top1; ++i) {
                printf("%c", charArray[i]);

            }
        }
    } else if (isEmpty(stackNumber)) {
        printf("Stack 1 is empty");
    }
    if (stackNumber == 2 && !isEmpty(stackNumber)) {
        for (int i = N - 1; i >= top2; i--) {
            printf("%c", charArray[i]);
        }
    } else if (isEmpty(stackNumber)) {
        printf("Stack 2 is empty");
    }
}

int main() {
    Stacks myStack;
    myStack.create();
    myStack.Insert(1, 'A');
    myStack.Insert(1, 'B');
    myStack.Insert(1, 'C');
    myStack.Insert(1, 'D');
    myStack.Insert(1, 'E');
    myStack.Insert(1, 'F');

    myStack.Insert(2, 'W');
    myStack.Insert(2, 'X');
    myStack.Insert(2, 'Y');
    myStack.Insert(2, 'Z');
    myStack.printStack(1);
    myStack.printStack(2);
    if (!myStack.isEmpty(1))
        myStack.deleteFromStack(1);
    if (!myStack.isEmpty(2))
        myStack.deleteFromStack(2);
    return 0;
}

Linked List Interview Question

Recently, my friend asked me to write a basic program in C to add a indicator (e.g. any number or char) between the two nodes of the given linked list. He is asked this question in one of the job interviews.
The output of the program with the indicator of  99 is shown as below.

The code

/* 
 * File:   main.cpp
 * Author: bekoC
 */

//Some notes about my compilation errors:
//There are many ways to get a segfault, 
//at least in the lower-level languages such as C(++). 
//A common way to get a segfault is to dereference a null pointer:
//int *p = NULL;
//*p = 1;


#include <cstdlib>
#include <stdlib.h>
#include <stdio.h>
using namespace std;

struct Node {
    int data;
    struct Node *next;
};

struct Node* insert(struct Node *head, int data) {
    struct Node *temp = (struct Node*) malloc(sizeof (struct Node)); // for the new node address is taken by malloc in heap
    temp->data = data;
    temp->next = NULL;
    if (head == NULL) {
        head = temp; // first node give the head temp's address
    } else {
        Node *temp1 = head;
        while (temp1->next != NULL) {
            temp1 = temp1->next;
        }
        temp1->next = temp;
    }
    return head;
}

struct Node* addIndicator(struct Node *head, int indicator) {
    struct Node *newNode, *scan;
    struct Node *temp1 = head;
    while (temp1 != NULL && temp1->next != NULL) {
        newNode = (struct Node*) malloc(sizeof (struct Node));
        newNode->data = indicator;
        newNode->next = temp1->next;
        temp1->next = newNode;
        temp1 = newNode->next;
    }
    return head;
}

void printList(struct Node *head) {
    //using recursion
    if (head == NULL) {
        return;
    }
    printf("%d  ", head->data);
    printList(head->next);
}

int main(int argc, char** argv) {
    struct Node *head = NULL; // head point to no address
    head = insert(head, 2);
    head = insert(head, 4);
    head = insert(head, 6);
    head = insert(head, 8);
    printList(head);
    printf("\n");
    //Adding a indicator which is 99
    int indicator = 99;
    head = addIndicator(head, indicator);
    printList(head);
    return 0;
}

Yet another another language for scientific computing

i have written a post about which programming  language to use for the implementation of algorithms in machine learning, yet i have came across another scientific computing language while i was searching about recurrence. It's called Julia:

Julia is a high-level, high-performance dynamic programming language for technical computing, with syntax that is familiar to users of other technical computing environments. It provides a sophisticated compiler, distributed parallel execution, numerical accuracy, and an extensive mathematical function library. The library, largely written in Julia itself, also integrates mature, best-of-breed C and Fortran libraries for linear algebra, random number generation, signal processing, and string processing. In addition, the Julia developer community is contributing a number of external packages through Julia’s built-in package manager at a rapid pace. IJulia, a collaboration between the IPython and Julia communities, provides a powerful browser-based graphical notebook interface to Julia.

i'm not sure that i need to use it but it desires to look at it because it is stated in official release that Julia beats all other high-level systems (i.e. everything besides C and Fortran) on all micro-benchmarks.

you can read  more detailed post "MATLAB, R, and Julia: Languages for data analysis" if you'd like to give a decision of using it or not. In addition, i'd like to recommend reading the post "A Matlab Programmer's Take On Julia",which  criticizes Matlab. Finally, take a look at some benchmark results of fibonacci and matrix multiplication computations implemented with Julia.

Understanding thread basics with Python

Recently, I needed to use the threads in order to increase the CPU efficiency in terms of idle time in Python. I came a cross great examples from agiliq's blog. I have reimplemented his code and add some minor comments.
You can refer the original post with more details.

Examining thread order:

'''
Created on Sep 29, 2013
@author: Bekoc::algorithms
'''

from threading import Thread
import time
import urllib2

class GetUrlThread(Thread):
    def __init__(self, url):
        self.url = url 
        super(GetUrlThread, self).__init__()

    def run(self):
        resp = urllib2.urlopen(self.url)
        print self.url, resp.getcode()

def get_responses():
    urls = ['http://www.google.com', 'http://www.amazon.com', 'http://www.ebay.com', 'http://www.alibaba.com', 'http://www.reddit.com']
    start = time.time()
    threads = []
    for url in urls:
        t = GetUrlThread(url)
        threads.append(t)
        print ('Thread %s is calling %s' %(t.getName(), url))
        t.start()
    for t in threads:
        t.join()
    print "Elapsed time: %s" % (time.time()-start)

get_responses()

Race Condition example without use of lock.acquire and lock.release:

'''
Created on Sep 30, 2013
@author: Bekoc::algorithms
'''
from threading import Thread
import time
#define a global variable
some_var = 0

class IncrementThreadRaceCondition(Thread):
    def run(self):
        #we want to read a global variable
        #and then increment it
        global some_var
        read_value = some_var
        time.sleep(.001)
        print "some_var in %s is %d" % (self.name, read_value)
        some_var = read_value + 1 
        print "some_var in %s after increment is %d" % (self.name, some_var)

def use_increment_thread():
    threads2 = []
    for i in range(50):
        t = IncrementThreadRaceCondition()
        threads2.append(t)
        t.start()
       
    for t in threads2:
        t.join()
    print "After 50 modifications, some_var should have become 50"
    print "After 50 modifications, some_var is %d" % (some_var,)

use_increment_thread()

in order to prevent the race condition change the run function  and import Lock:

from threading import Lock

    def run(self):
        #we want to read a global variable
        #and then increment it
        global some_var
        lock.acquire()
        read_value = some_var
        time.sleep(.001)
        print "some_var in %s is %d" % (self.name, read_value)
        some_var = read_value + 1 
        print "some_var in %s after increment is %d" % (self.name, some_var)
        lock.release()

Khan Academy offers Python Programming course

There are lots of Python online courses going on around, but take a look at khan academy playlist. It includes the basic programming concepts with Python.

1. Introduction to Programs Data Types and Variables
2. Binary Numbers
3. Python Lists
4. For Loops in Python
5. While Loops in Python
6. Fun with Strings
7. Writing a Simple Factorial Program. (Python 2)
8. Stepping Through the Factorial Program
9. Flowchart for the Factorial Program
10. Python 3 Not Backwards Compatible with Python 2
11. Defining a Factorial Function
12. Diagramming What Happens with a Function Call

1. Recursive Factorial Function
2. Comparing Iterative and Recursive Factorial Functions
3. Exercise - Write a Fibonacci Function
4. Iterative Fibonacci Function Example
5. Stepping Through Iterative Fibonacci Function
6. Recursive Fibonacci Example
7. Stepping Through Recursive Fibonacci Function
8. Exercise - Write a Sorting Function
9. Insertion Sort Algorithm
10. Insertion Sort in Python
11. Stepping Through Insertion Sort Function
12. Simpler Insertion Sort Function

Linear Regression with Matlab Using Batch Gradient Descent Algorithm

i will implement linear regression which can be adapted classification easily,  i use Matlab by following the  Dr. Andrew Ng's class. You can watch the classes online from here.

While implementing i also came across a very nice blog post, actually only dataset differs, in this case i use the original dataset given by the Dr. Ng, more details of the code below can be reached from here (DSPlog )

download data

 the linear regression model is 

 and the batch gradient descent update rule is

theta_vec = [0 0]';
alpha = 0.007;
err = [0 0]';
for kk = 1:1500
 h_theta = (x*theta_vec);
 h_theta_v = h_theta*ones(1,n);
 y_v = y*ones(1,n);
 theta_vec = theta_vec - alpha*1/m*sum((h_theta_v - y_v).*x).';
 err(:,kk) = 1/m*sum((h_theta_v - y_v).*x).';
end

Cost Function:

For different values of theta, in this case theta0 and theta1, we can plot the cost function J(theta) in 3d space or as a contour.

j_theta = zeros(250, 250);   % initialize j_theta
theta0_vals = linspace(-5000, 5000, 250);
theta1_vals = linspace(-200, 200, 250);
for i = 1:length(theta0_vals)
   for j = 1:length(theta1_vals)
  theta_val_vec = [theta0_vals(i) theta1_vals(j)]';
  h_theta = (x*theta_val_vec);
  j_theta(i,j) = 1/(2*m)*sum((h_theta - y).^2);
    end
end
figure;
surf(theta0_vals, theta1_vals,10*log10(j_theta.'));
xlabel('theta_0'); ylabel('theta_1');zlabel('10*log10(Jtheta)');
title('Cost function J(theta)');
figure;
contour(theta0_vals,theta1_vals,10*log10(j_theta.'))
xlabel('theta_0'); ylabel('theta_1')
title('Cost function J(theta)');

%linear regression using gradient descent algorithm to find the coefficients, 
%there are many ways to find the coefficients but now we apply gradient descent
%one dimensional data, with an output 
x = load('ex2x.dat');
y = load('ex2y.dat');

figure % open a new figure window
plot(x, y, 'o');
ylabel('Height in meters')
xlabel('Age in years')


m = length(y); % store the number of training examples
x = [ones(m, 1), x]; % Add a column of ones to x
n = size(x,2);
%this part is for minimizing the Theta_vec: coefficients.
theta_vec = [0 0]';
alpha = 0.007;
err = [0 0]';
for kk = 1:10000
 h_theta = (x*theta_vec);
 h_theta_v = h_theta*ones(1,n);
 y_v = y*ones(1,n);
 theta_vec = theta_vec - alpha*1/m*sum((h_theta_v - y_v).*x).';
 err(:,kk) = 1/m*sum((h_theta_v - y_v).*x)';
end

figure;
plot(x(:,2),y,'bs-');
hold on
plot(x(:,2),x*theta_vec,'rp-');
legend('measured', 'predicted');
grid on;
xlabel('x');
ylabel('y');
title('Measured and predicted ');



j_theta = zeros(250, 250);   % initialize j_theta
theta0_vals = linspace(-5000, 5000, 250);
theta1_vals = linspace(-200, 200, 250);
for i = 1:length(theta0_vals)
   for j = 1:length(theta1_vals)
  theta_val_vec = [theta0_vals(i) theta1_vals(j)]';
  h_theta = (x*theta_val_vec);
  j_theta(i,j) = 1/(2*m)*sum((h_theta - y).^2);
    end
end
figure;
surf(theta0_vals, theta1_vals,10*log10(j_theta.'));
xlabel('theta_0'); ylabel('theta_1');zlabel('10*log10(Jtheta)');
title('Cost function J(theta)');
figure;
contour(theta0_vals,theta1_vals,10*log10(j_theta.'))
xlabel('theta_0'); ylabel('theta_1')
title('Cost function J(theta)');