Reinforcement Learning with Pygame Tutorial Part - 2

(Part 1: http://anixtech.blogspot.com/2016/12/reinforcement-learning-with-pygame-part.html)

I have created the second part of the tutorial as a video. It could be seen here:
https://www.youtube.com/watch?v=o5GiQkClbAY

The full code is also available in github:
https://github.com/hasanIqbalAnik/q-learning-python-example

Reinforcement Learning With Pygame Part - 1

Following is an example of a simple game which could be used to train agents. 

import pygame as pg
from pygame.locals import *
import sys
import random


def new_rect_after_action(newrct, act):
    if act == 'right':
        if newrct.right + newrct.width > windowWidth:
            return newrct
        else:
            return pg.Rect(newrct.left + newrct.width, newrct.top, newrct.width,
                           newrct.height)  # Rect(left, top, width, height)
    else:  # action is left
        if newrct.left - newrct.width < 0:
            return newrct
        else:
            return pg.Rect(newrct.left - newrct.width, newrct.top, newrct.width,
                           newrct.height)  # Rect(left, top, width, height)


def circle_falling(crclradius):
    newx = random.randint(0 + crclradius, 200 - crclradius)
    multiplier = random.randint(1, 4)
    newx *= multiplier
    return newx


windowWidth = 800
windowHeight = 400
FPS = 3  # frames per second setting
fpsClock = pg.time.Clock()

pg.init()

window = pg.display.set_mode((windowWidth, windowHeight))  # width, height
pg.display.set_caption('Catch the ball!')

# setup colors
RED = (255, 0, 0)
GREEN = (0, 255, 0)
WHITE = (255, 255, 255)
BLACK = (0, 0, 0)

# specify circle properties
crclCentreX = 50
crclCentreY = 50
crclRadius = 20

# specify rectangle properties
rctLeft = 400
rctTop = 350
rctWidth = 200
rctHeight = 50
rct = pg.Rect(rctLeft, rctTop, rctWidth, rctHeight)  # Rect(left, top, width, height)

action = 'left'
score = 0
font = pg.font.Font(None, 30)

while True:

    for event in pg.event.get():
        if event.type == QUIT:
            pg.quit()
            sys.exit()
        elif event.type == pg.KEYDOWN:
            if event.key == pg.K_LEFT:
                action = 'left'
                rct = new_rect_after_action(rct, action)
            elif event.key == pg.K_RIGHT:
                action = 'right'
                rct = new_rect_after_action(rct, action)

    window.fill(BLACK)
    # at this position, the rectangle should be here. else loses
    if crclCentreY >= windowHeight - rctHeight - crclRadius:
        if rct.left <= crclCentreX <= rct.right:
            score += 1
        else:
            score -= 1
        crclCentreX = circle_falling(crclRadius)
        crclCentreY = 50
    else:
        crclCentreY += 20

    pg.draw.circle(window, RED, (crclCentreX, crclCentreY),
                   crclRadius)  # circle(Surface, color, pos(x, y), radius, width=0)
    pg.draw.rect(window, GREEN, rct)  # rect(Surface, color, Rect, width=0)

    text = font.render('score: ' + str(score), True, (238, 58, 140))
    window.blit(text, (windowWidth - 100, 10))
    pg.display.update()
    fpsClock.tick(FPS)

This will create a game window that looks like the following screenshot. The green rectangle can be moved with left and right arrow of the keyboard. If the player can place it under the falling red circle, he gets a point, else loses one. That's all there is in this game. In the second part of this tutorial, we will add intelligence and train an agent to play this game.

Update: Second part is now available: http://anixtech.blogspot.com/2016/12/reinforcement-learning-with-pygame.html

Simplified Neural Network Backpropagation Example Python

This is an even simplified version of this fantastic tutorial: https://mattmazur.com/2015/03/17/a-step-by-step-backpropagation-example/
Let's consider a very simple neural network shown in the following image:

All we want is that given the input .1, the network would output the value .9.
We will be using $\frac{1}{1+e^{-x}}$ as our activation function. Initially both weight w1 and w2 are initialized to .5, though this can be any arbitrary value. However, different initialization shall converge to different final values.

First let's consider the net input in the hidden layer(neth). It's simply the input multiplied by w1. That is:(Here, $outh$ means the output of the hidden layer, $neto$ means the net input in the output layer, $outo$ is the output of the output layer. $E$ is the error. )

$neth = i1 * w1 = .1 * .5 = .05$
$outh = \frac{1}{1+e^{-neth}} = \frac{1}{1+e^{-.05}} = .51249$

$neto = outh * w2 = .51249 * .5 = .25624$
$outo = \frac{1}{1+e^{-neto}} = \frac{1}{1+e^{-.25624}} = .563711$

$E = \frac{1}{2}*(target - outo)^2 =  \frac{1}{2}*(.9 - .563711)^2 = .056545$

This concludes our forward pass of the network. We've selected that particular error function for our convenience, as it's convex and differentiable. Now, we can see that our target value was $.9$ but we've got $.5637$ . We are going to update our weights $w1, w2$ so that the network produces that desired output. How would we do that is the reason for using the famous Backpropagation algorithm. The basic idea is, we will calculate how much responsibility should each of the weights take for producing that error. Then the weight would be modified accordingly and the next iteration would begin.

First let's see what happens in case of $w2$. The partial differentiation expression $\frac{\partial E}{\partial w2}$ says the rate of change in the error with respect to $w2$. As it can't be calculated directly, we would use the chain rule to find it's constituting parts.

$\frac{\partial E}{\partial w2} = \frac{\partial E}{\partial outo} * \frac{\partial outo}{\partial neto} * \frac{\partial neto}{\partial w2}$


Let's find out each of those parts separately:

$\frac{\partial E}{\partial outo} = \frac{\partial}{\partial outo}(\frac{1}{2}*(target - outo)^2)  = \frac{1}{2} * 2 * ( target - outo) * (-1) $
$= (.9 - .563711) * (-1) = -.336289$

$\frac{\partial outo}{\partial neto} = \frac{\partial }{\partial neto} (\frac{1}{1+e^{-neto}}) = \frac{-1}{(1+e^{-neto})^2} * e^{-neto} * (-1) = \frac{e^{-neto}}{(1+e^{-neto})^2}$
$  = outo * (1-outo) = .563711 * ( 1- .563711) = .24594$

$\frac{\partial neto}{\partial w2} = \frac{\partial (w2*outh)}{\partial w2} = outh = .51249$

So finally we have found the responsibility of $w2$ for producing the error $E$:

$\frac{\partial E}{\partial w2} = -.33628 * .24594 * .51249 = -.042386$

Now let's turn our attention to $w1$, that is we want to calculate $\frac{\partial E}{\partial w1}$. Again, chain rule comes to our rescue.

$\frac{\partial E}{\partial w1} = \frac{\partial E}{\partial outo} * \frac{\partial outo}{\partial neto} *\frac{\partial neto}{\partial outh} * \frac{\partial outh}{\partial neth} * \frac{\partial neth}{\partial w1}$

We already have calculated $\frac{\partial E}{\partial outo}$ and $\frac{\partial outo}{\partial neto}$. Let's calculate the rest.

$\frac{\partial neto}{\partial outh} = \frac{\partial (w2*outh)}{\partial outh} = w2 = .5$
$\frac{\partial outh}{\partial neth} = outh * ( 1- outh) = .51249 * ( 1- .51249) = .24984$
$\frac{\partial neth}{\partial w1} = \frac{\partial (w1 * input)}{\partial w1} = input = .1 $

So,

$\frac{\partial E}{\partial w1} = -.336289 * .24594 * .5 * .24984 * .1 = -.001033$

Now it's time to update our initial guesses of $w1$ and $w2$. Let's say our learning rate $\eta = .9$
$w1 = w1 - (\eta * \frac{\partial E}{\partial w1}) =  .5 + .9 * 001033 = 0.5009298654$
Similarly,
$w2 = w2 - ( \eta * \frac{\partial E}{\partial w2}) = .5 + .9 * 042386 = 0.538148$

We can continue this process until our network's output is very close to our target. The following is a python implementation of that process:

from __future__ import division
import math

ctr = 0
eta = .9
w1 = .5
w2 = .5

i1 = .1
target = .9

while True:
    ctr += 1
    neth = i1 * w1
    outh = 1 / (1 + math.exp(-neth))

    neto = w2 * outh
    outo = 1 / (1 + math.exp(-neto))

    e = .5 * math.pow((target - outo), 2)
    print 'error now is: ', e
    if e < .00005:
        print '*' * 100
        print 'target is: ', target
        print 'output is now: ', outo
        print 'w1 and w2 is: ', w1, ',', w2
        print 'number of iterations: ', ctr
        break
    # update w2
    de_by_douto = (-1) * (target - outo)
    douto_by_dneto = outo * (1 - outo)
    dneto_by_dw2 = outh

    de_by_dw2 = de_by_douto * douto_by_dneto * dneto_by_dw2

    # update w1
    # de_by_dw1 = de_by_douto * douto_by_dneto * dneto_by_douth * douth_by_dneth * dneth_by_dw1
    dneto_by_douth = w2
    de_by_outh = de_by_douto * douto_by_dneto * dneto_by_douth
    douth_by_dneth = outh * (1 - outh)
    dneth_by_dw1 = i1

    de_by_dw1 = de_by_outh * douth_by_dneth * dneth_by_dw1

    # updated w1 and w2
    w1 = w1 - (eta * de_by_dw1)
    w2 = w2 - (eta * de_by_dw2)

Sample output:

Who Liked Your Last N Number of Facebook Posts Most

Previously it was possible to find out this information for any friend's profile too. But Facebook has recently revoked that access and apps now need to have this access directly from that friend. However, you can find out this info about yourself in the following way:

1. Go to this link: https://developers.facebook.com/tools/explorer/.
2. Click on 'Get Token' > 'Get User Access Token' > check user_posts > 'Get Access token'.
3. This would ask you for the permission to read your timeline. Click on 'Continue as your-name'.
4. You would get an access token like the following image. Copy the whole string.

5. Replace the 'your-token' part in the following code with that string.

from __future__ import division
import operator
import requests
from facepy import GraphAPI

large_dic = {}
token = 'your-token'
graph = GraphAPI(token)
n = 50  # number of posts in the history that we want to consider
ctr = 0
pcount = 0
res = graph.get('me/feed')

def update_progress(progress):
    print '\r{0}% [{1}]'.format(int(progress), '#'*int(progress/10)),

def merge_two_dicts(dict1, dict2):
    """ merges and sums up dict1 with the values from dict2. if key2 from dict2 is present in dict1, then dict1[key2] = dict1[key2] + dict2[key2]
    else it appends key2 with the value as it is
    example: if
    d1 = {'a': 1, 'b': 1, 'c': 1} and
    d2 = {'a': 4, 'd': 1}
    then it returns {'a': 5, 'c': 1, 'b': 1, 'd': 1}
    :param dict1, dict2
    :return merged dict1 and dict2
    """
    for key2 in dict2:
        if key2 in dict1:
            dict1[key2] = dict1[key2] + dict2[key2]
        else:
            dict1[key2] = dict2[key2]
    return dict1


def get_likes_list(post_id):
    """given a post id, this method queries the graphapi to get details of the post, gets the likes list
    and adds them into a dictionary by using the person's name as the key
    :param post_id
    :return dictionary containing likers
    """
    global graph
    global pcount
    d = {}  # to hold the liker's name with a count of 1

    likes_obj = graph.get(post_id + '/likes')
    while True:
        try:
            for liker in likes_obj['data']:
                d[liker['name']] = 1
            likes_obj = requests.get(likes_obj['paging']['next']).json()
        except KeyError:
            pcount += 1
            update_progress((pcount/n)*100)
            break
    return d


for item in res['data']:
    merge_two_dicts(large_dic, get_likes_list(item['id']))

while True:
    if ctr == (n - 25) / 25:  # as 25 is the default limit of /me/feed, to get the last n posts, set the value of ctr to (n-25)/25
        break
    try:
        res = requests.get(res['paging']['next']).json()
        for item in res['data']:
            merge_two_dicts(large_dic, get_likes_list(item['id']))
        ctr += 1
    except KeyError:
        print 'end of calculating posts'


sorted_x = sorted(large_dic.items(), key=operator.itemgetter(1))

print '\n'
print '*'*100
print 'total posts considered', pcount
print 'total likes generated', sum(large_dic.values())
print 'average likes per post', sum(large_dic.values()) / pcount
print '*'*100

for item in reversed(sorted_x):
    print item[0], 'liked your posts', item[1], 'times'

6. Change the value of n in the program to your preference, this is the number of posts that we want to consider while counting the likes.

7. Run that program, it might take a while depending on the response time of facebook server. The output would contain the list of likes in a sorted fashion. With the most liker at the top.

Sample output:

Comparing Radix, Merge, Selection and Bubble Sort

Implementation of Radix Sort:

bucket = {}

# creating bucket for each digits
def reset_bucket():
    global bucket
    for i in xrange(10): bucket[i] = []
# to handle variable length items`
def radix_group_numbers(numbers):
    joined = []
    max_length = max(len(str(x)) for x in numbers)
    bkt = {}
    for i in xrange(max_length + 1):
        bkt[i] = []
    for item in numbers:
        bkt[len(str(item))].append(item)
    for k, v in bkt.items():
        if v: joined.extend(radix_sort(v))
    return joined

def put_in_buckets(numbers, index):  # assuming that the numbers are of equal lenght
    global bucket
    try:
        for item in numbers:
            c = str(item)[index]
            if c is '0':
                bucket[0].append(item)
            elif c is '1':
                bucket[1].append(item)
            elif c is '2':
                bucket[2].append(item)
            elif c is '3':
                bucket[3].append(item)
            elif c is '4':
                bucket[4].append(item)
            elif c is '5':
                bucket[5].append(item)
            elif c is '6':
                bucket[6].append(item)
            elif c is '7':
                bucket[7].append(item)
            elif c is '8':
                bucket[8].append(item)
            elif c is '9':
                bucket[9].append(item)
    except IndexError:
        print 'error occured'
        print numbers, index

def radix_sort(lst):
    global bucket
    max_length = max(len(str(x)) for x in lst)
    for i in xrange(max_length):
        reset_bucket()
        put_in_buckets(lst, max_length - i - 1)
        lst = []
        for key, value in bucket.iteritems():
            if value is not None:
                for item in value:
                    lst.append(item)
        bucket.clear()
    return lst

Implementation of Merge Sort:

# merge two lists
class stack:
    def __init__(self, lst=None):
        if lst is None:
            self.items = []
        else:
            self.items = lst

    def pop(self):
        return self.items.pop()
    def peek(self):
        return self.items[len(self.items) - 1]
    def __str__(self):
        return str(self.items)
    def empty(self):
        return len(self.items) is 0
    def reverse(self):
        l = []
        while self.empty() is False:
            l.append(self.items.pop())
        return l

def merge(l1, l2):  # takes two sorted lists and merges them while keeping the order
    ll = []  # long list consists of l1 + l2
    l1.reverse()
    l2.reverse()
    st1 = stack(l1)
    st2 = stack(l2)
    while st1.empty() is False and st2.empty() is False:
        i1 = st1.peek()
        i2 = st2.peek()
        if i1 <= i2:
            ll.append(i1)
            l1.pop()
        else:
            ll.append(i2)
            l2.pop()
    if st1.empty() is True and st2.empty() is False:
        ll = ll + st2.reverse()
    elif st2.empty() is True and st1.empty() is False:
        ll = ll + st1.reverse()

    return ll


def merge_sort(lst):
    if len(lst) is 1:
        return lst
    else:
        mid = len(lst) // 2
        left = merge_sort(lst[:mid])
        right = merge_sort(lst[mid:])
        return merge(left, right)

Implementation of Selection Sort:

from utils import get_numbers_list

def selection_sort(unsorted_lst):
    for k in xrange(len(unsorted_lst) - 1):  # at max this much iteration required to fully sort

        current_minimum = unsorted_lst[k]  # first elem as the current minimum
        current_minimum_idx = k

        for i in xrange(k, len(unsorted_lst), 1):
            if current_minimum > unsorted_lst[i]:
                current_minimum = unsorted_lst[i]
                current_minimum_idx = i
        # swap the current minimum with the kth element
        tmp = unsorted_lst[k]
        unsorted_lst[k] = current_minimum
        unsorted_lst[current_minimum_idx] = tmp
    return unsorted_lst

Bubble Sort:

from utils import get_numbers_list

def bubble_sort(lst):
    max_iter = len(lst)
    counter = 0
    for m in xrange(max_iter):
        for i in xrange(max_iter - 1 - counter):  # -counter to execute faster
            if lst[i] > lst[i + 1]:  # should swap
                tmp = lst[i + 1]
                lst[i + 1] = lst[i]
                lst[i] = tmp
        counter += 1
    return lst

Now let's generate some sample data:

import sys
import random

# n is the number of samples generated
def generate_samples(n):
    lst = []
    while(n > 0):
        lst.append(random.randint(0, 1000000)) #random integer between 0 and a million
        n -= 1
    f = open("data.txt", "w+")
    for item in lst:
        f.write(str(item) + " ")
    f.close()
    print 'samples written at data.txt...'

def get_numbers_list(fname):
    num_arr = []
    f = open(fname, "r")
    for line in f:
        num_arr.append(list(int(x) for x in line.strip().split(" ")))
    return num_arr[0]

Then we measure the time required for their executions:

import time

from utils import generate_samples, get_numbers_list
from bubble_sort import bubble_sort
from merge_sort import merge_sort
from radix_sort import radix_sort, radix_group_numbers
from selection_sort import selection_sort

#first we need to generate a hugh amount of data

generate_samples(10000)

lst = get_numbers_list('data.txt')

start_time = time.time()
radix_group_numbers(lst)
end_time = time.time()

print 'radix sort took', end_time - start_time, 'seconds'

start_time = time.time()
merge_sort(lst)
end_time = time.time()

print 'merge sort took', end_time - start_time, 'seconds'

start_time = time.time()
bubble_sort(lst)
end_time = time.time()

print 'bubble sort took', end_time - start_time, 'seconds'

start_time = time.time()
selection_sort(lst)
end_time = time.time()

print 'selection sort took', end_time - start_time, 'seconds'

Sample outputs: for 1000 data points:

radix sort took 0.00797295570374 seconds
merge sort took 0.0257380008698 seconds
bubble sort took 0.0774970054626 seconds
selection sort took 0.0248939990997 seconds

for 10000 data points:

radix sort took 0.0430598258972 seconds
merge sort took 0.187104940414 seconds
bubble sort took 8.21462607384 seconds
selection sort took 2.59621596336 seconds

Fining all possible subsets of a given set

Three possible way to create all subsets of a given set is shown below:

1. Using Binary:

lst = ['a', 'b', 'c', 'd', 'e']
def interpret_combination(s):
 s = s[::-1] #reversing the string to match the index positions with the given list
 if not '1' in list(s): # all zeros
  return []
 else:
  indices_of_ones = [i for i, x in enumerate(list(s)) if x == '1'] # find the indices of 1s only
  l = list(s) #convert to list as you can't operate on strings
  for item in indices_of_ones:
   l[item] = lst[item] # replace the one's with the characters on the list
  l = [x for x in l if x != '0'] # replace the zeros
  return ''.join(map(str, l)) # join them to make a string again
all_subsets = []

for i in xrange(2 ** len(lst)): # until 2^(lenght of list)
 l1 = len(str(bin(i))[2:]) # binary representation without the preceding 0b
 l2 = len(lst)
 s = (l2 - l1)*'0' + str(bin(i))[2:] # prepending the zeros in front to make it more readable
 all_subsets.append(interpret_combination(s)) 

print all_subsets

2. Recursive:

def adder(elem, elem_set): # copy elem with each element of elem_set
 exp_lst = []
 for item in elem_set:
  exp_lst.append(elem + item)
 return exp_lst

def all_subsets(s): 
 if len(s) is 0:
  return ''
 elif len(s) is 1:
  return ['', s[0]]
 else:
  return all_subsets(s[:-1]) + adder(s[len(s)-1], all_subsets(s[:-1]))  

print all_subsets(['a','b', 'c', 'd'])  

3. Iterative:

lst = ['a', 'b', 'c', 'd', 'e']
sub_arr = []
sub_arr_tmp = []

for item in lst:
 for elem in sub_arr:
  sub_arr_tmp.append([item + elem])
 sub_arr.append(item)
 if len(sub_arr_tmp) is 0:
  sub_arr.append('')
 else:
   for e in sub_arr_tmp:
    if(isinstance(e, list)):
     sub_arr.append(e[0])
 del sub_arr_tmp[:]

print set(sub_arr)
print len(set(sub_arr))

Finding Gaussian Distribution Parameters Using Expectation Maximization Matlab

Expecation Maximization is a powerful approach to estimate parameters of statistical distributions. It consists of two steps, i) Expectation step, where the likelihood of the dataset is estimated using the current set of parameters. ii) Maximization step, where the parameters are recalculated based on the likelihood found in the expectation step.

Let's see an example: Suppose we have a dataset which has only one feature, namely, age. Suppose the data was taken from a country where most of the people are young and from another country where most of them are old. We are interested in finding out the mean and variance of these two distributions, assuming that they were normally distributed. 

Say our dataset looks like this: 

$[25, 27, 27, 28, 28, 30, 100, 105, 105, 107, 107, 108]$ . We shall implement the EM algorithm to find out the means and variances.

clc;clear all;
% sample dataset
d = [25 27 27 28 28 30 100 105 105 107 107 108];

%initial guesses of the two means and variances
mu1 = 10;
mu2 = 50;
sigma1 = 10;
sigma2 = 10;

% initial responsibilities
p1 = .5;
p2 = .5;

%until converges

for k = 1:100
    
% start expectation step    

v1 = p1 * normpdf(d, mu1, sqrt(sigma1)); 
v2 = p2 * normpdf(d, mu2, sqrt(sigma2));

r1 = v1 ./ (v1 + v2); %divide ith element of v1 by the sum of ith element of both v1 and v2
r2 = v2 ./ (v1 + v2) ;%divide ith element of v2 by the sum of ith element of both v1 and v2

% end expectation step    

% start maximization step

mu1 = sum(r1 .* d) / sum(r1); %recalculate mu1, weighted average of the datapoints
mu2 = sum(r2 .* d) / sum(r2); %recalculate mu2, weighted average of the datapoints

sigma1 = sum((d-mu1) .* (d-mu1) .* r1) / sum(r1); %recalculate weighted variance first distribution
sigma2 = sum((d-mu2) .* (d-mu2) .* r2) / sum(r2); %recalculate weighted variance second distribution

p1 = sum(r1) / length(d); % recalculate responsibilities, total weights normalized by number of data points
p2 = sum(r2) / length(d);

% end maximization step

end

mu1 % 27.5000, mean of the first distribution
mu2 % 105.3333, mean of the second distribution
sigma1 % 2.2500, variance of the first distribution
sigma2 % 6.8889, variance of the second distribution