## Raspberry Pi & Arduino: a laser pointer communication and a LDR voltage sigmoid

So I finally got some more time to play with my Raspberry Pi GPIOs and Arduino, this post will explain how to use a LDR (Photoresistor, Light Dependent Resistor) on the Raspberry Pi to detect a laser light emitted by an Arduino.

# A small intro on the rationale

So I’m working on a Symbolic Regression Machine written in C/C++ called Shine, which is intended to be a JIT for Genetic Programming libraries (like Pyevolve for instance). The main rationale behind Shine is that we have today a lot of research on speeding Genetic Programming using GPUs (the GPU fever !) or any other special hardware, etc, however we don’t have many papers talking about optimizing GP using the state of art compilers optimizations like we have on clang, gcc, etc.

The “hot spot” or the component that consumes a lot of CPU resources today on Genetic Programming is the evaluation of each individual in order to calculate the fitness of the program tree. This evaluation is often executed on each set of parameters of the “training” set. Suppose you want to make a symbolic regression of a single expression like the Pythagoras Theorem and you have a linear space of parameters from 1.0 to 1000.0 with a step of 0.1 you have 10.000 evaluations for each individual (program tree) of your population !

What Shine does is described on the image below:

It takes the individual of the Genetic Programming engine and then converts it to LLVM Intermediate Representation (LLVM assembly language), after that it runs the transformation passes of the LLVM (here is where the true power of modern compilers enter on the GP context) and then the LLVM JIT converts the optimized LLVM IR into native code for the specified target (X86, PowerPC, etc).

You can see below the Shine architecture:

This architecture brings a lot of flexibility for Genetic Programming, you can for instance write functions that could be used later on your individuals on any language supported by the LLVM, what matters to Shine is the LLVM IR, you can use any language that LLVM supports and then use the IR generated by LLVM, you can mix code from C, C++, Ada, Fortran, D, etc and use your functions as non-terminal nodes of your Genetic Programming trees.

Shine is still on its earlier development, it looks a simple idea but I still have a lot of problems to solve, things like how to JIT the evaluation process itself instead of doing calls from Python using ctypes bindings of the JITed trees.

# Doing Genetic Programming on the Python AST itself

During the development of Shine, an idea happened to me, that I could use a restricted Python Abstract Syntax Tree (AST) as the representation of individuals on a Genetic Programming engine, the main advantage of this is the flexibility and the possibility to reuse a lot of things. Of course that a shared library written in C/C++ would be useful for a lot of Genetic Programming engines that doesn’t uses Python, but since my spare time to work on this is becoming more and more rare I started to rethink the approach and use Python and the LLVM bindings for LLVM (LLVMPY) and I just discovered that is pretty easy to JIT a restricted set of the Python AST to native code using LLVM, and this is what this post is going to show.

# JIT’ing a restricted Python AST

The most amazing part of LLVM is obviously the amount of transformation passes, the JIT and of course the ability to use the entire framework through a simple API (ok, not so simple sometimes). To simplify this example, I’m going to use an arbitrary restricted AST set of the Python AST that supports only subtraction (-), addition (+), multiplication (*) and division (/).

To understand the Python AST, you can use the Python parser that converts source into AST:

 12345678910 >>> import ast >>> astp = ast.parse("2*7") >>> ast.dump(astp) 'Module(     body=[Expr(         value=BinOp(             left=Num(n=2), op=Mult(), right=Num(n=7)         )     )] )'

What the parse created was an Abstract Syntax Tree containing the BinOp (Binary Operation) with the left operator as the number 2, the right operator as the number 7 and the operation itself as Multiplication(Mult), very easy to understand. What we are going to do to create the LLVM IR is to create a visitor that is going to visit each node of the tree. To do that, we can subclass the Python NodeVisitor class from the ast module. What the NodeVisitor does is to visit each node of the tree and then call the method ‘visit_OPERATOR’ if it exists, when the NodeVisitor is going to visit the node for the BinOp for example, it will call the method ‘visit_BinOp’ passing as parameter the BinOp node itself.

The structure of the class for for the JIT visitor will look like the code below:

 12345678910 # Import the ast and the llvm Python bindings import ast from llvm import * from llvm.core import * from llvm.ee import * import llvm.passes as lp class AstJit(ast.NodeVisitor):     def __init__(self):         pass

What we need to do now is to create an initialization method to keep the last state of the JIT visitor, this is needed because we are going to JIT the content of the Python AST into a function and the last instruction of the function needs to return what was the result of the last instruction visited by the JIT. We also need to receive a LLVM Module object in which our function will be created as well the closure type, for the sake of simplicity I’m not type any object, I’m just assuming that all numbers from the expression are integers, so the closure type will be the LLVM integer type.

 1234567891011121314151617181920212223242526272829 def __init__(self, module, parameters):         self.last_state = None         self.module = module         # Parameters that will be created on the IR function         self.parameters = parameters         self.closure_type = Type.int()         # An attribute to hold a link to the created function         # so we can use it to JIT later         self.func_obj = None         self._create_builder()     def _create_builder(self):         # How many parameters of integer type         params = [self.closure_type] * len(self.parameters)         # The prototype of the function, returning a integer         # and receiving the integer parameters         ty_func = Type.function(self.closure_type, params)         # Add the function to the module with the name 'func_ast_jit'         self.func_obj = self.module.add_function(ty_func, 'func_ast_jit')         # Create an argument in the function for each parameter specified         for index, pname in enumerate(self.parameters):             self.func_obj.args[index].name = pname         # Create a basic block and the builder         bb = self.func_obj.append_basic_block("entry")         self.builder = Builder.new(bb)

Now what we need to implement on our visitor is the ‘visit_OPERATOR’ methods for the BinOp and for the Numand Name operators. We will also implement the method to create the return instruction that will return the last state.

 12345678910111213141516171819202122232425262728293031323334353637383940 # A 'Name' is a node produced in the AST when you     # access a variable, like '2+x+y', 'x' and 'y' are     # the two names created on the AST for the expression.     def visit_Name(self, node):         # This variable is what function argument ?         index = self.parameters.index(node.id)         self.last_state = self.func_obj.args[index]         return self.last_state     # Here we create a LLVM IR integer constant using the     # Num node, on the expression '2+3' you'll have two     # Num nodes, the Num(n=2) and the Num(n=3).     def visit_Num(self, node):         self.last_state = Constant.int(self.closure_type, node.n)         return self.last_state       # The visitor for the binary operation     def visit_BinOp(self, node):         # Get the operation, left and right arguments         lhs = self.visit(node.left)         rhs = self.visit(node.right)         op = node.op         # Convert each operation (Sub, Add, Mult, Div) to their         # LLVM IR integer instruction equivalent         if isinstance(op, ast.Sub):             op = self.builder.sub(lhs, rhs, 'sub_t')         elif isinstance(op, ast.Add):             op = self.builder.add(lhs, rhs, 'add_t')         elif isinstance(op, ast.Mult):             op = self.builder.mul(lhs, rhs, 'mul_t')         elif isinstance(op, ast.Div):             op = self.builder.sdiv(lhs, rhs, 'sdiv_t')                 self.last_state = op         return self.last_state     # Build the return (ret) statement with the last state     def build_return(self):         self.builder.ret(self.last_state)

And that is it, our visitor is ready to convert a Python AST to a LLVM IR assembly language, to run it we’ll first create a LLVM module and an expression:

 12345 module = Module.new('ast_jit_module') # Note that I'm using two variables 'a' and 'b' expr = "(2+3*b+33*(10/2)+1+3/3+a)/2" node = ast.parse(expr) print ast.dump(node)

Will output:

Module(body=[Expr(value=BinOp(left=BinOp(left=BinOp(left=BinOp(


Now we can finally run our visitor on that generated AST the check the LLVM IR output:

 1234 visitor = AstJit(module, ['a', 'b']) visitor.visit(node) visitor.build_return() print module

Will output the LLVM IR:

; ModuleID = 'ast_jit_module'

define i32 @func_ast_jit(i32 %a, i32 %b) {
entry:
%mul_t = mul i32 3, %b
%sdiv_t = sdiv i32 %add_t4, 2
ret i32 %sdiv_t
}


Now is when the real fun begins, we want to run LLVM optimization passes to optimize our code with an equivalent GCC -O2 optimization level, to do that we create a PassManagerBuilder and a PassManager, the PassManagerBuilder is the component that adds the passes to the PassManager, you can also manually add arbitrary transformations like dead code elimination, function inlining, etc:

 12345678910 pmb = lp.PassManagerBuilder.new() # Optimization level pmb.opt_level = 2 pm = lp.PassManager.new() pmb.populate(pm) # Run the passes into the module pm.run(module) print module

Will output:

; ModuleID = 'ast_jit_module'

define i32 @func_ast_jit(i32 %a, i32 %b) nounwind readnone {
entry:
%mul_t = mul i32 %b, 3
%sdiv_t = sdiv i32 %add_t4, 2
ret i32 %sdiv_t
}


And here we have the optimized LLVM IR of the Python AST expression. The next step is to JIT that IR into native code and then execute it with some parameters:

 123456 ee = ExecutionEngine.new(module)     arg_a = GenericValue.int(Type.int(), 100)     arg_b = GenericValue.int(Type.int(), 42)         retval = ee.run_function(visitor.func_obj, [arg_a, arg_b])     print "Return: %d" % retval.as_int()

Will output:

Return: 197


And that’s it, you have created a AST->LLVM IR converter, optimized the LLVM IR with the transformation passes and then converted it to native code using the LLVM execution engine. I hope you liked =)

## Review Board – Review Requests plotting using Python

Review Board is one of these projects that Python community is always proud of, I really think that it became the de facto standard for code reviews nowadays. Review Board comes with an interesting an very useful REST API which you can use to retrieve information about the code reviews, comments, diffs and every single information stored on its database. So, I created a small Python script called rbstats that retrieves information about the Review Requests done by the users and then plots a heat map using matplotlib. I’ll show an example of the use on the Review Board system of the Apache foundation.

To use the tool, just call it with the API URL of the Review Boars system, i.e.:

python rb_stats.py
--max-results 80 https://reviews.apache.org/api

an then you’ll get a graphical plot like this (click to enlarge):

Where the “hottest” points are weighted according to the number of the code reviews that the user have created to the other axis user. You can also plot the statistics by user, for instance of the user benjaminhindman using the autumn colormap and with 400 max results:

python rb_stats.py
--max-results 400
--from-user benjaminhindman
--colormap autumn https://reviews.apache.org/api

Click to enlarge the image:

## Accessing HP Cloud OpenStack Nova using Python and Requests

So, my request to enter on the free and private beta season of the new HP Cloud Services was gently accepted by the HP Cloud team, and today I finally got some time to play with the OpenStack API at HP Cloud. I’ll start with the first impressions I had with the service:

The user interface of the management is very user-friendly, the design is much like of the Twitter Bootstrap, see the screenshot below of the “Compute” page from the “Manage” section:

As you can see, they have a set of 4 Ubuntu images and a CentOS, I think that since they are still in the beta period, soon we’ll have more default images to use.

Here is a screenshot of the instance size set:

Since they are using OpenStack, I really think that they should have imported the vocabulary of the OpenStack into the user interface, and instead of calling it “Size”, it would be more sensible to use “Flavour“.

The user interface still doesn’t have many features, something that I would really like to have is a “Stop” or something like that for the instances, only the “Terminate” function is present on the Manage interface, but those are details that they should be still working on since they’re only in beta.

Another important info to cite is that the access to the instances are done through SSH using a generated RSA key that they provide to you.

Let’s dig into the OpenStack API now.

## OpenStack API

To access the OpenStack API you’ll need the credentials for the authentication, HP Cloud services provide these keys on the Manage interface for each zone/service you have, see the screenshot below (with keys anonymized of course):

Now, OpenStack authentication could be done in different schemes, the scheme that I know that HP supports is the token authentication. I know that there is a lot of clients already supporting the OpenStack API (some have no documentation, some have weird API design, etc.), but the aim of this post is to show how easy would be to create a simple interface to access the OpenStack API using Python and Requests (HTTP for Humans !).

Let’s start defining our authentication scheme by sub-classing Requests AuthBase:

 123456789 class OpenStackAuth(AuthBase):     def __init__(self, auth_user, auth_key):         self.auth_key = auth_key         self.auth_user = auth_user def __call__(self, r):     r.headers['X-Auth-User'] = self.auth_user     r.headers['X-Auth-Key'] = self.auth_key     return r

As you can see, we’re defining the X-Auth-User and the X-Auth-Key in the header of the request with the parameters. These parameters are respectively your Account ID and  Access Key we cited earlier. Now, all you have to do is to make the request itself using the authentication scheme, which is pretty easy using Requests:

 1234 ENDPOINT_URL = 'https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/' ACCESS_KEY = 'Your Access Key' ACCOUNT_ID = 'Your Account ID' response = requests.get(ENDPOINT_URL, auth=OpenStackAuth(ACCOUNT_ID, ACCESS_KEY))

And that is it, you’re done with the authentication mechanism using just a few lines of code, and this is how the request is going to be sent to the HP Cloud service server:

This request is sent to the HP Cloud Endpoint URL (https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/). Let’s see now how the server answered this authentication request:

You can show this authentication response using Requests by printing the header attribute of the request Response object. You can see that the server answered our request with two important header items: X-Server-Management-URL and the X-Auth-Token. The management URL is now our new endpoint, is the URL we should use to do further requests to the HP Cloud services and the X-Auth-Token is the authentication Token that the server generated based on our credentials, these tokens are usually valid for 24 hours, although I haven’t tested it.

What we need to do now is to sub-class the Requests AuthBase class again but this time defining only the authentication token that we need to use on each new request we’re going to make to the management URL:

 1234567 class OpenStackAuthToken(AuthBase):     def __init__(self, request):         self.auth_token = request.headers['x-auth-token'] def __call__(self, r):     r.headers['X-Auth-Token'] = self.auth_token     return r

Note that the OpenStackAuthToken is receiving now a response request as parameter, copying the X-Auth-Token and setting it on the request.

Let’s consume a service from the OpenStack API v.1.1, I’m going to call the List Servers API function, parse the results using JSON and then show the results on the screen:

 12345678910 # Get the management URL from the response header mgmt_url = response.headers['x-server-management-url'] # Create a new request to the management URL using the /servers path # and the OpenStackAuthToken scheme we created r_server = requests.get(mgmt_url + '/servers', auth=OpenStackAuthToken(response)) # Parse the response and show it to the screen json_parse = json.loads(r_server.text) print json.dumps(json_parse, indent=4)

And this is what we get in response to this request:

 12345678910111213141516171819202122232425262728293031323334 {     "servers": [         {             "id": 22378,             "uuid": "e2964d51-fe98-48f3-9428-f3083aa0318e",             "links": [                 {                     "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/20817201684751/servers/22378",                     "rel": "self"                 },                 {                     "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/20817201684751/servers/22378",                     "rel": "bookmark"                 }             ],             "name": "Server 22378"         },         {             "id": 11921,             "uuid": "312ff473-3d5d-433e-b7ee-e46e4efa0e5e",             "links": [                 {                     "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/20817201684751/servers/11921",                     "rel": "self"                 },                 {                     "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/20817201684751/servers/11921",                     "rel": "bookmark"                 }             ],             "name": "Server 11921"         }     ] }

And that is it, now you know how to use Requests and Python to consume OpenStack API. If you wish to read more information about the API and how does it works, you can read the documentation here.

- Christian S. Perone

## Announce: Stallion v0.2 released !

I just tagged and released the v0.2 version of the Stallion. In the change log (Github project page), you can see that a lot of bugs were fixed and some new features were introduced in this release. I added compatibility with almost all Python 2.x versions, PyPy 1.7+ (and probably older versions too), I also fixed the compatibility with the Internet Explorer browser, now you should be able to use Stallion with Chrome, Firefox and IE.

The most important feature introduced is the global checking for updates (a lot of people requested it):

The new checking is under the menu “PyPI Repository”. Another new feature is the refactoring on the visual appearance of the package classifiers:

Some small visual enhancements were also introduced, like the little gray marker next to the selected package:

I hope you liked, I’m looking forward to implement more features as soon as possible, but a new version shouldn’t be released until next year.

Visit the project page at Github to get instructions on how to update or install Stallion.

- Christian S. Perone

## Announce: ‘Stallion’ – Python Package Manager

I’m happy to announce the first release v.0.1 of the Stallion project. Stallion is a visual Python package manager compatible with Python 2.6 and 2.7 (I still haven’t tested it with Python 2.5).

The motivation behind Stallion is to provide an user friendly visualization with some management features (most of them are still under development) for Python packages installed on your local Python distribution. Stallion is intended to be used specially for Python newcomers.

The project is currently hosted at Github, so feel free to fork, contribute, make suggestion, report bugs, etc.

### Installation

All you need to do to install Stallion is to use your favorite Python distribution system, examples:

 123 user@machine:~/$pip install stallion or user@machine:~/$ easy_install stallion

After installing Stallion, you need to start the local server by using:

 1 user@machine:~/\$ python -m stallion.main

And if it’s all ok, Stallion will start the server on localhost only at the port 5000, so all you need to do now is to browse into the URL http://localhost:5000

### See some screenshots (click to enlarge)

Click on the screenshots below to enlarge.

Home

Installed package information

PyPI version mismatch diagnosis

## Hacking into Python objects internals

You know, Python represents every object using the low-level C API PyObject (or PyVarObject for variable-size objects) structure, so, concretely, you can cast any Python object pointer to this type; this inheritance is built by hand, every new object must have a leading macro called PyObject_HEAD which defines the PyObject header for the object. The PyObject structure is declared in Include/object.h as:

 123 typedef struct _object {     PyObject_HEAD } PyObject;

and the PyObject_HEAD macro is defined as:

 1234 #define PyObject_HEAD                   \     _PyObject_HEAD_EXTRA                \     Py_ssize_t ob_refcnt;               \     struct _typeobject *ob_type;

… with two fields (forget the _PyObject_HEAD_EXTRA, it’s only used for a tracing debug feature) called ob_refcnt and ob_type, representing the reference counting for the object and the type of the object. I know you can use sys.getrefcount to get the reference counting of an object, but hacking the object memory using ctypes is by far more powerful, since you can get the contents of any field of the object (in cases where you don’t have a native API for that), I’ll show more examples later, but lets focus on the reference counting field of the object.

### Getting the reference count (ob_refcnt)

So, in Python, we have the built-in function id(), this function returns the identity of the object, but, looking at its definition on CPython implementation, you’ll notice that id() returns the memory address of the object, see the source in Python/bltinmodule.c:

 12345 static PyObject * builtin_id(PyObject *self, PyObject *v) {     return PyLong_FromVoidPtr(v); }

… the function PyLong_FromVoidPtr returns a Python long object from a void pointer. So, in CPython, this value is the address of the object in the memory as shown below:

 123 >>> value = 666 >>> hex(id(value)) '0x8998e50' # memory address of the 'value' object

Now that we have the memory address of the object, we can use the Python ctypes module to get the reference counting by accessing the attribute ob_refcnt, here is the code needed to do that:

 123456 >>> value = 666 >>> value_address = id(value) >>> >>> ob_refcnt = ctypes.c_long.from_address(value_address) >>> ob_refcnt c_long(1)

What I’m doing here is getting the integer value from the ob_refcnt attribute of the PyObject in memory.  Let’s add a new reference for the object ‘value’ we created, and then check the reference count again:

 12345 >>> value_ref = value >>> id(value_ref) == id(value) True >>> ob_refcnt c_long(2)

Note that the reference counting was increased by 1 due to the new reference variable called ‘value_ref’.

### Interned strings state (ob_sstate)

Now, getting the reference count wasn’t even funny, we already had the sys.getrefcount API for that, but what about the interned state of the strings ? In order to avoid the creation of different allocations for the same string (and to speed comparisons), Python uses a dictionary that works like a “cache” for strings, this dictionary is defined in Objects/stringobject.c:

 123456789 /* This dictionary holds all interned strings.  Note that references to strings in this dictionary are *not* counted in the string's ob_refcnt. When the interned string reaches a refcnt of 0 the string deallocation function will delete the reference from this dictionary. Another way to look at this is that to say that the actual reference count of a string is:  s->ob_refcnt + (s->ob_sstate?2:0) */ static PyObject *interned;

I also copied here the comment about the dictionary, because is interesting to note that the strings in the dictionary aren’t counted in the string’s ob_refcnt.

So, the interned state of a string object is hold in the attribute ob_sstate of the string object, let’s see the definition of the Python string object:

 123456789101112131415 typedef struct {     PyObject_VAR_HEAD     long ob_shash;     int ob_sstate;     char ob_sval[1];     /* Invariants:     *     ob_sval contains space for 'ob_size+1' elements.     *     ob_sval[ob_size] == 0.     *     ob_shash is the hash of the string or -1 if not computed yet.     *     ob_sstate != 0 iff the string object is in stringobject.c's     *       'interned' dictionary; in this case the two references     *       from 'interned' to this object are *not counted* in ob_refcnt.     */ } PyStringObject;

As you can note, strings objects inherit from the PyObject_VAR_HEAD macro, which defines another header attribute, let’s see the definition to get the complete idea of the structure:

 123 #define PyObject_VAR_HEAD               \     PyObject_HEAD                       \     Py_ssize_t ob_size; /* Number of items in variable part */

The PyObject_VAR_HEAD macro adds another field called ob_size, which is the number of items on the variable part of the Python object (i.e. the number of items on a list object). So, before getting to the ob_sstate field, we need to shift our offset to skip the fields ob_refcnt (long), ob_type (void*) (from PyObject_HEAD), the field ob_size (long) (from PyObject_VAR_HEAD) and the field ob_shash (long) from the PyStringObject. Concretely, we need to skip this offset (3 fields with size long and one field with size void*) of bytes:

 123 >>> ob_sstate_offset = ctypes.sizeof(ctypes.c_long)*3 + ctypes.sizeof(ctypes.c_voidp) >>> ob_sstate_offset 16

Now, let’s prepare two cases, one that we know that isn’t interned and another that is surely interned, then we’ll force the interned state of the other non-interned string to check the result of the ob_sstate attribute:

 12345678 >>> a = "lero" >>> b = "".join(["l", "e", "r", "o"]) >>> ctypes.c_long.from_address(id(a) + ob_sstate_offset) c_long(1) >>> ctypes.c_long.from_address(id(b) + ob_sstate_offset) c_long(0) >>> ctypes.c_long.from_address(id(intern(b)) + ob_sstate_offset) c_long(1)

Note that the interned state for the object “a” is 1 and for the object “b” is 0. After forcing the interned state of the variable “b”, we can see that the field ob_sstate has changed to 1.

### Changing internal states (evil mode)

Now, let’s suppose we want to change some internal state of a Python object through the interpreter. Let’s try to change the value of an int object. Int objects are defined in Include/intobject.h:

 1234 typedef struct {     PyObject_HEAD     long ob_ival; } PyIntObject;

As you can see, the internal value of an int is stored in the field ob_ival, to change it, we just need to skip the ob_refcnt (long) and the ob_type (void*) from the PyObject_HEAD:

 12345678 >>> value = 666 >>> ob_ival_offset = ctypes.sizeof(ctypes.c_long) + ctypes.sizeof(ctypes.c_voidp) >>> ob_ival = ctypes.c_int.from_address(id(value)+ob_ival_offset) >>> ob_ival c_long(666) >>> ob_ival.value = 8 >>> value 8

And that is it, we have changed the value of the int value directly in the memory.

I hope you liked it, you can play with lots of other Python objects like lists and dicts, note that this method is just intended to show how the Python objects are structured in the memory and how you can change them using the native API, but obviously, you’re not supposed to use this to change the value of ints lol.

Update 11/29/11: you’re not supposed to do such things on your production code or something like that, in this post I’m doing lazy assumptions about arch details like sizes of primitives, etc. Be warned.

## Machine Learning :: Text feature extraction (tf-idf) – Part II

Read the first part of this tutorial: Text feature extraction (tf-idf) – Part I.

This post is a continuation of the first part where we started to learn the theory and practice about text feature extraction and vector space model representation. I really recommend you to read the first part of the post series in order to follow this second post.

Since a lot of people liked the first part of this tutorial, this second part is a little longer than the first.

### Introduction

In the first post, we learned how to use the term-frequency to represent textual information in the vector space. However, the main problem with the term-frequency approach is that it scales up frequent terms and scales down rare terms which are empirically more informative than the high frequency terms. The basic intuition is that a term that occurs frequently in many documents is not a good discriminator, and really makes sense (at least in many experimental tests); the important question here is: why would you, in a classification problem for instance, emphasize a term which is almost present in the entire corpus of your documents ?

The tf-idf weight comes to solve this problem. What tf-idf gives is how important is a word to a document in a collection, and that’s why tf-idf incorporates local and global parameters, because it takes in consideration not only the isolated term but also the term within the document collection. What tf-idf then does to solve that problem, is to scale down the frequent terms while scaling up the rare terms; a term that occurs 10 times more than another isn’t 10 times more important than it, that’s why tf-idf uses the logarithmic scale to do that.

But let’s go back to our definition of the $\mathrm{tf}(t,d)$ which is actually the term count of the term $t$ in the document $d$. The use of this simple term frequency could lead us to problems like keyword spamming, which is when we have a repeated term in a document with the purpose of improving its ranking on an IR (Information Retrieval) system or even create a bias towards long documents, making them look more important than they are just because of the high frequency of the term in the document.

To overcome this problem, the term frequency $\mathrm{tf}(t,d)$ of a document on a vector space is usually also normalized. Let’s see how we normalize this vector.

### Vector normalization

Suppose we are going to normalize the term-frequency vector $\vec{v_{d_4}}$ that we have calculated in the first part of this tutorial. The document $d4$ from the first part of this tutorial had this textual representation:

 1 d4: We can see the shining sun, the bright sun.

And the vector space representation using the non-normalized term-frequency of that document was:

$\vec{v_{d_4}} = (0,2,1,0)$

To normalize the vector, is the same as calculating the Unit Vector of the vector, and they are denoted using the “hat” notation: $\hat{v}$. The definition of the unit vector $\hat{v}$ of a vector $\vec{v}$ is:

$\displaystyle \hat{v} = \frac{\vec{v}}{\|\vec{v}\|_p}$

Where the $\hat{v}$ is the unit vector, or the normalized vector, the $\vec{v}$ is the vector going to be normalized and the $\|\vec{v}\|_p$ is the norm (magnitude, length) of the vector $\vec{v}$ in the $L^p$ space (don’t worry, I’m going to explain it all).

The unit vector is actually nothing more than a normalized version of the vector, is a vector which the length is 1.

Source: http://processing.org/learning/pvector/

But the important question here is how the length of the vector is calculated and to understand this, you must understand the motivation of the $L^p$ spaces, also called Lebesgue spaces.

### Lebesgue spaces

Source: http://processing.org/learning/pvector/

Usually, the length of a vector $\vec{u} = (u_1, u_2, u_3, \ldots, u_n)$ is calculated using the Euclidean norma norm is a function that assigns a strictly positive length or size to all vectors in a vector space -, which is defined by:

Source: http://processing.org/learning/pvector/

$\|\vec{u}\| = \sqrt{u^2_1 + u^2_2 + u^2_3 + \ldots + u^2_n}$

But this isn’t the only way to define length, and that’s why you see (sometimes) a number $p$ together with the norm notation, like in $\|\vec{u}\|_p$. That’s because it could be generalized as:

$\displaystyle \|\vec{u}\|_p = ( \left|u_1\right|^p + \left|u_2\right|^p + \left|u_3\right|^p + \ldots + \left|u_n\right|^p )^\frac{1}{p}$

and simplified as:

$\displaystyle \|\vec{u}\|_p = (\sum\limits_{i=1}^{n}\left|\vec{u}_i\right|^p)^\frac{1}{p}$

So when you read about a L2-norm, you’re reading about the Euclidean norm, a norm with $p=2$, the most common norm used to measure the length of a vector, typically called “magnitude”; actually, when you have an unqualified length measure (without the $p$ number), you have the L2-norm (Euclidean norm).

When you read about a L1-norm, you’re reading about the norm with $p=1$, defined as:

$\displaystyle \|\vec{u}\|_1 = ( \left|u_1\right| + \left|u_2\right| + \left|u_3\right| + \ldots + \left|u_n\right|)$

Which is nothing more than a simple sum of the components of the vector, also known as Taxicab distance, also called Manhattan distance.

Taxicab geometry versus Euclidean distance: In taxicab geometry all three pictured lines have the same length (12) for the same route. In Euclidean geometry, the green line has length $6 \times \sqrt{2} \approx 8.48$, and is the unique shortest path.
Source: Wikipedia :: Taxicab Geometry

Note that you can also use any norm to normalize the vector, but we’re going to use the most common norm, the L2-Norm, which is also the default in the 0.9 release of the scikits.learn. You can also find papers comparing the performance of the two approaches among other methods to normalize the document vector, actually you can use any other method, but you have to be concise, once you’ve used a norm, you have to use it for the whole process directly involving the norm (a unit vector that used a L1-norm isn’t going to have the length 1 if you’re going to take its L2-norm later).

### Back to vector normalization

Now that you know what the vector normalization process is, we can try a concrete example, the process of using the L2-norm (we’ll use the right terms now) to normalize our vector $\vec{v_{d_4}} = (0,2,1,0)$ in order to get its unit vector $\hat{v_{d_4}}$. To do that, we’ll simple plug it into the definition of the unit vector to evaluate it:

$\hat{v} = \frac{\vec{v}}{\|\vec{v}\|_p} \\ \\ \hat{v_{d_4}} = \frac{\vec{v_{d_4}}}{||\vec{v_{d_4}}||_2} \\ \\ \\ \hat{v_{d_4}} = \frac{(0,2,1,0)}{\sqrt{0^2 + 2^2 + 1^2 + 0^2}} \\ \\ \hat{v_{d_4}} = \frac{(0,2,1,0)}{\sqrt{5}} \\ \\ \small \hat{v_{d_4}} = (0.0, 0.89442719, 0.4472136, 0.0)$

And that is it ! Our normalized vector $\hat{v_{d_4}}$ has now a L2-norm $\|\hat{v_{d_4}}\|_2 = 1.0$.

Note that here we have normalized our term frequency document vector, but later we’re going to do that after the calculation of the tf-idf.

### The term frequency – inverse document frequency (tf-idf) weight

Now you have understood how the vector normalization works in theory and practice, let’s continue our tutorial. Suppose you have the following documents in your collection (taken from the first part of tutorial):

 123456789 Train Document Set: d1: The sky is blue. d2: The sun is bright. Test Document Set: d3: The sun in the sky is bright. d4: We can see the shining sun, the bright sun.

Your document space can be defined then as $D = \{ d_1, d_2, \ldots, d_n \}$ where $n$ is the number of documents in your corpus, and in our case as $D_{train} = \{d_1, d_2\}$ and $D_{test} = \{d_3, d_4\}$. The cardinality of our document space is defined by $\left|{D_{train}}\right| = 2$ and $\left|{D_{test}}\right| = 2$, since we have only 2 two documents for training and testing, but they obviously don’t need to have the same cardinality.

Let’s see now, how idf (inverse document frequency) is then defined:

$\displaystyle \mathrm{idf}(t) = \log{\frac{\left|D\right|}{1+\left|\{d : t \in d\}\right|}}$

where $\left|\{d : t \in d\}\right|$ is the number of documents where the term $t$ appears, when the term-frequency function satisfies $\mathrm{tf}(t,d) \neq 0$, we’re only adding 1 into the formula to avoid zero-division.

The formula for the tf-idf is then:

$\mathrm{tf\mbox{-}idf}(t) = \mathrm{tf}(t, d) \times \mathrm{idf}(t)$

and this formula has an important consequence: a high weight of the tf-idf calculation is reached when you have a high term frequency (tf) in the given document (local parameter) and a low document frequency of the term in the whole collection (global parameter).

Now let’s calculate the idf for each feature present in the feature matrix with the term frequency we have calculated in the first tutorial:

$M_{train} = \begin{bmatrix} 0 & 1 & 1 & 1\\ 0 & 2 & 1 & 0 \end{bmatrix}$

Since we have 4 features, we have to calculate $\mathrm{idf}(t_1)$, $\mathrm{idf}(t_2)$, $\mathrm{idf}(t_3)$, $\mathrm{idf}(t_4)$:

$\mathrm{idf}(t_1) = \log{\frac{\left|D\right|}{1+\left|\{d : t_1 \in d\}\right|}} = \log{\frac{2}{1}} = 0.69314718$

$\mathrm{idf}(t_2) = \log{\frac{\left|D\right|}{1+\left|\{d : t_2 \in d\}\right|}} = \log{\frac{2}{3}} = -0.40546511$

$\mathrm{idf}(t_3) = \log{\frac{\left|D\right|}{1+\left|\{d : t_3 \in d\}\right|}} = \log{\frac{2}{3}} = -0.40546511$

$\mathrm{idf}(t_4) = \log{\frac{\left|D\right|}{1+\left|\{d : t_4 \in d\}\right|}} = \log{\frac{2}{2}} = 0.0$

These idf weights can be represented by a vector as:

$\vec{idf_{train}} = (0.69314718, -0.40546511, -0.40546511, 0.0)$

Now that we have our matrix with the term frequency ($M_{train}$) and the vector representing the idf for each feature of our matrix ($\vec{idf_{train}}$), we can calculate our tf-idf weights. What we have to do is a simple multiplication of each column of the matrix $M_{train}$ with the respective $\vec{idf_{train}}$ vector dimension. To do that, we can create a square diagonal matrix called $M_{idf}$ with both the vertical and horizontal dimensions equal to the vector $\vec{idf_{train}}$ dimension:

$M_{idf} = \begin{bmatrix} 0.69314718 & 0 & 0 & 0\\ 0 & -0.40546511 & 0 & 0\\ 0 & 0 & -0.40546511 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}$

and then multiply it to the term frequency matrix, so the final result can be defined then as:

$M_{tf\mbox{-}idf} = M_{train} \times M_{idf}$

Please note that the matrix multiplication isn’t commutative, the result of $A \times B$ will be different than the result of the $B \times A$, and this is why the $M_{idf}$ is on the right side of the multiplication, to accomplish the desired effect of multiplying each idf value to its corresponding feature:

$\begin{bmatrix} \mathrm{tf}(t_1, d_1) & \mathrm{tf}(t_2, d_1) & \mathrm{tf}(t_3, d_1) & \mathrm{tf}(t_4, d_1)\\ \mathrm{tf}(t_1, d_2) & \mathrm{tf}(t_2, d_2) & \mathrm{tf}(t_3, d_2) & \mathrm{tf}(t_4, d_2) \end{bmatrix} \times \begin{bmatrix} \mathrm{idf}(t_1) & 0 & 0 & 0\\ 0 & \mathrm{idf}(t_2) & 0 & 0\\ 0 & 0 & \mathrm{idf}(t_3) & 0\\ 0 & 0 & 0 & \mathrm{idf}(t_4) \end{bmatrix} \\ = \begin{bmatrix} \mathrm{tf}(t_1, d_1) \times \mathrm{idf}(t_1) & \mathrm{tf}(t_2, d_1) \times \mathrm{idf}(t_2) & \mathrm{tf}(t_3, d_1) \times \mathrm{idf}(t_3) & \mathrm{tf}(t_4, d_1) \times \mathrm{idf}(t_4)\\ \mathrm{tf}(t_1, d_2) \times \mathrm{idf}(t_1) & \mathrm{tf}(t_2, d_2) \times \mathrm{idf}(t_2) & \mathrm{tf}(t_3, d_2) \times \mathrm{idf}(t_3) & \mathrm{tf}(t_4, d_2) \times \mathrm{idf}(t_4) \end{bmatrix}$

Let’s see now a concrete example of this multiplication:

$M_{tf\mbox{-}idf} = M_{train} \times M_{idf} = \\ \begin{bmatrix} 0 & 1 & 1 & 1\\ 0 & 2 & 1 & 0 \end{bmatrix} \times \begin{bmatrix} 0.69314718 & 0 & 0 & 0\\ 0 & -0.40546511 & 0 & 0\\ 0 & 0 & -0.40546511 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix} \\ = \begin{bmatrix} 0 & -0.40546511 & -0.40546511 & 0\\ 0 & -0.81093022 & -0.40546511 & 0 \end{bmatrix}$

And finally, we can apply our L2 normalization process to the $M_{tf\mbox{-}idf}$ matrix. Please note that this normalization is “row-wise” because we’re going to handle each row of the matrix as a separated vector to be normalized, and not the matrix as a whole:

$M_{tf\mbox{-}idf} = \frac{M_{tf\mbox{-}idf}}{\|M_{tf\mbox{-}idf}\|_2}$ $= \begin{bmatrix} 0 & -0.70710678 & -0.70710678 & 0\\ 0 & -0.89442719 & -0.4472136 & 0 \end{bmatrix}$

And that is our pretty normalized tf-idf weight of our testing document set, which is actually a collection of unit vectors. If you take the L2-norm of each row of the matrix, you’ll see that they all have a L2-norm of 1.

### Python practice

Environment Used: Python v.2.7.2, Numpy 1.6.1, Scipy v.0.9.0, Sklearn (Scikits.learn) v.0.9.

Now the section you were waiting for ! In this section I’ll use Python to show each step of the tf-idf calculation using the Scikit.learn feature extraction module.

The first step is to create our training and testing document set and computing the term frequency matrix:

 1234567891011121314151617 from sklearn.feature_extraction.text import CountVectorizer train_set = ("The sky is blue.", "The sun is bright.") test_set = ("The sun in the sky is bright.", "We can see the shining sun, the bright sun.") count_vectorizer = CountVectorizer() count_vectorizer.fit_transform(train_set) print "Vocabulary:", count_vectorizer.vocabulary # Vocabulary: {'blue': 0, 'sun': 1, 'bright': 2, 'sky': 3} freq_term_matrix = count_vectorizer.transform(test_set) print freq_term_matrix.todense() #[[0 1 1 1] #[0 2 1 0]]

Now that we have the frequency term matrix (called freq_term_matrix), we can instantiate the TfidfTransformer, which is going to be responsible to calculate the tf-idf weights for our term frequency matrix:

 12345678 from sklearn.feature_extraction.text import TfidfTransformer tfidf = TfidfTransformer(norm="l2") tfidf.fit(freq_term_matrix) print "IDF:", tfidf.idf_ # IDF: [ 0.69314718 -0.40546511 -0.40546511  0.        ]

Note that I’ve specified the norm as L2, this is optional (actually the default is L2-norm), but I’ve added the parameter to make it explicit to you that it it’s going to use the L2-norm. Also note that you can see the calculated idf weight by accessing the internal attribute called idf_. Now that fit() method has calculated the idf for the matrix, let’s transform the freq_term_matrix to the tf-idf weight matrix:

 12345 tf_idf_matrix = tfidf.transform(freq_term_matrix) print tf_idf_matrix.todense() # [[ 0.         -0.70710678 -0.70710678  0.        ] # [ 0.         -0.89442719 -0.4472136   0.        ]]

And that is it, the tf_idf_matrix is actually our previous $M_{tf\mbox{-}idf}$ matrix. You can accomplish the same effect by using the Vectorizer class of the Scikit.learn which is a vectorizer that automatically combines the CountVectorizer and the TfidfTransformer to you. See this example to know how to use it for the text classification process.

I really hope you liked the post, I tried to make it simple as possible even for people without the required mathematical background of linear algebra, etc. In the next Machine Learning post I’m expecting to show how you can use the tf-idf to calculate the cosine similarity.

If you liked it, feel free to comment and make suggestions, corrections, etc.

### References

Understanding Inverse Document Frequency: on theoretical arguments for IDF

Wikipedia :: tf-idf

The classic Vector Space Model

Sklearn text feature extraction code

03 Oct 2011 - Added the info about the environment used for Python examples

## Machine Learning :: Text feature extraction (tf-idf) – Part I

### Short introduction to Vector Space Model (VSM)

In information retrieval or text mining, the term frequency – inverse document frequency also called tf-idf, is a well know method to evaluate how important is a word in a document. tf-idf are also a very interesting way to convert the textual representation of information into a Vector Space Model (VSM), or into sparse features, we’ll discuss more about it later, but first, let’s try to understand what is tf-idf and the VSM.

VSM has a very confusing past, see for example the paper The most influential paper Gerard Salton Never Wrote that explains the history behind the ghost cited paper which in fact never existed; in sum, VSM is an algebraic model representing textual information as a vector, the components of this vector could represent the importance of a term (tf–idf) or even the absence or presence (Bag of Words) of it in a document; it is important to note that the classical VSM proposed by Salton incorporates local and global parameters/information (in a sense that it uses both the isolated term being analyzed as well the entire collection of documents). VSM, interpreted in a lato sensu, is a space where text is represented as a vector of numbers instead of its original string textual representation; the VSM represents the features extracted from the document.

Let’s try to mathematically define the VSM and tf-idf together with concrete examples, for the concrete examples I’ll be using Python (as well the amazing scikits.learn Python module).

### Going to the vector space

The first step in modeling the document into a vector space is to create a dictionary of terms present in documents. To do that, you can simple select all terms from the document and convert it to a dimension in the vector space, but we know that there are some kind of words (stop words) that are present in almost all documents, and what we’re doing is extracting important features from documents, features do identify them among other similar documents, so using terms like “the, is, at, on”, etc.. isn’t going to help us, so in the information extraction, we’ll just ignore them.

Let’s take the documents below to define our (stupid) document space:

 123456789 Train Document Set: d1: The sky is blue. d2: The sun is bright. Test Document Set: d3: The sun in the sky is bright. d4: We can see the shining sun, the bright sun.

Now, what we have to do is to create a index vocabulary (dictionary) of the words of the train document set, using the documents $d1$ and $d2$ from the document set, we’ll have the following index vocabulary denoted as $\mathrm{E}(t)$ where the $t$ is the term:

$\mathrm{E}(t) = \begin{cases} 1, & \mbox{if } t\mbox{ is blue''} \\ 2, & \mbox{if } t\mbox{ is sun''} \\ 3, & \mbox{if } t\mbox{ is bright''} \\ 4, & \mbox{if } t\mbox{ is sky''} \\ \end{cases}$

Note that the terms like “is” and “the” were ignored as cited before. Now that we have an index vocabulary, we can convert the test document set into a vector space where each term of the vector is indexed as our index vocabulary, so the first term of the vector represents the “blue” term of our vocabulary, the second represents “sun” and so on. Now, we’re going to use the term-frequency to represent each term in our vector space; the term-frequency is nothing more than a measure of how many times the terms present in our vocabulary $\mathrm{E}(t)$ are present in the documents $d3$ or $d4$, we define the term-frequency as a couting function:

$\mathrm{tf}(t,d) = \sum\limits_{x\in d} \mathrm{fr}(x, t)$

where the $\mathrm{fr}(x, t)$ is a simple function defined as:

$\mathrm{fr}(x,t) = \begin{cases} 1, & \mbox{if } x = t \\ 0, & \mbox{otherwise} \\ \end{cases}$

So, what the $tf(t,d)$ returns is how many times is the term $t$ is present in the document $d$. An example of this, could be $tf(sun'', d4) = 2$ since we have only two occurrences of the term “sun” in the document $d4$. Now you understood how the term-frequency works, we can go on into the creation of the document vector, which is represented by:

$\displaystyle \vec{v_{d_n}} =(\mathrm{tf}(t_1,d_n), \mathrm{tf}(t_2,d_n), \mathrm{tf}(t_3,d_n), \ldots, \mathrm{tf}(t_n,d_n))$

Each dimension of the document vector is represented by the term of the vocabulary, for example, the $\mathrm{tf}(t_1,d_2)$ represents the frequency-term of the term 1 or $t_1$ (which is our “blue” term of the vocabulary) in the document $d_2$.

Let’s now show a concrete example of how the documents $d_3$ and $d_4$ are represented as vectors:

$\vec{v_{d_3}} = (\mathrm{tf}(t_1,d_3), \mathrm{tf}(t_2,d_3), \mathrm{tf}(t_3,d_3), \ldots, \mathrm{tf}(t_n,d_3)) \\ \vec{v_{d_4}} = (\mathrm{tf}(t_1,d_4), \mathrm{tf}(t_2,d_4), \mathrm{tf}(t_3,d_4), \ldots, \mathrm{tf}(t_n,d_4))$

which evaluates to:

$\vec{v_{d_3}} = (0, 1, 1, 1) \\ \vec{v_{d_4}} = (0, 2, 1, 0)$

As you can see, since the documents $d_3$ and $d_4$ are:

 12 d3: The sun in the sky is bright. d4: We can see the shining sun, the bright sun.

The resulting vector $\vec{v_{d_3}}$ shows that we have, in order, 0 occurrences of the term “blue”, 1 occurrence of the term “sun”, and so on. In the $\vec{v_{d_3}}$, we have 0 occurences of the term “blue”, 2 occurrences of the term “sun”, etc.

But wait, since we have a collection of documents, now represented by vectors, we can represent them as a matrix with $|D| \times F$ shape, where $|D|$ is the cardinality of the document space, or how many documents we have and the $F$ is the number of features, in our case represented by the vocabulary size. An example of the matrix representation of the vectors described above is:

$M_{|D| \times F} = \begin{bmatrix} 0 & 1 & 1 & 1\\ 0 & 2 & 1 & 0 \end{bmatrix}$

As you may have noted, these matrices representing the term frequencies tend to be very sparse (with majority of terms zeroed), and that’s why you’ll see a common representation of these matrix as sparse matrices.

### Python practice

Environment Used: Python v.2.7.2, Numpy 1.6.1, Scipy v.0.9.0, Sklearn (Scikits.learn) v.0.9.

Since we know the  theory behind the term frequency and the vector space conversion, let’s show how easy is to do that using the amazing scikit.learn Python module.

Scikit.learn comes with lots of examples as well real-life interesting datasets you can use and also some helper functions to download 18k newsgroups posts for instance.

Since we already defined our small train/test dataset before, let’s use them to define the dataset in a way that scikit.learn can use:

 123 train_set = ("The sky is blue.", "The sun is bright.") test_set = ("The sun in the sky is bright.", "We can see the shining sun, the bright sun.")

In scikit.learn, what we have presented as the term-frequency, is called CountVectorizer, so we need to import it and create a news instance:

 12 from sklearn.feature_extraction.text import CountVectorizer vectorizer = CountVectorizer()

The CountVectorizer already uses as default “analyzer” called WordNGramAnalyzer, which is responsible to convert the text to lowercase, accents removal, token extraction, filter stop words, etc… you can see more information by printing the class information:

 123 print vectorizer CountVectorizer(analyzer__min_n=1, analyzer__stop_words=set(['all', 'six', 'less', 'being', 'indeed', 'over', 'move', 'anyway', 'four', 'not', 'own', 'through', 'yourselves', (...)

Let’s create now the vocabulary index:

 123 vectorizer.fit_transform(train_set) print vectorizer.vocabulary {'blue': 0, 'sun': 1, 'bright': 2, 'sky': 3}

See that the vocabulary created is the same as $E(t)$ (except because it is zero-indexed).

Let’s use the same vectorizer now to create the sparse matrix of our test_set documents:

 123456789 smatrix = vectorizer.transform(test_set) print smatrix (0, 1)        1 (0, 2)        1 (0, 3)        1 (1, 1)        2 (1, 2)        1

Note that the sparse matrix created called smatrix is a Scipy sparse matrix with elements stored in a Coordinate format. But you can convert it into a dense format:

 1234 smatrix.todense() matrix([[0, 1, 1, 1], ........[0, 2, 1, 0]], dtype=int64)

Note that the sparse matrix created is the same matrix $M_{|D| \times F}$ we cited earlier in this post, which represents the two document vectors $\vec{v_{d_3}}$ and $\vec{v_{d_4}}$.

We’ll see in the next post how we define the idf (inverse document frequency) instead of the simple term-frequency, as well how logarithmic scale is used to adjust the measurement of term frequencies according to its importance, and how we can use it to classify documents using some of the well-know machine learning approaches.

I hope you liked this post, and if you really liked, leave a comment so I’ll able to know if there are enough people interested in these series of posts in Machine Learning topics.

As promised, here is the second part of this tutorial series.

### References

The classic Vector Space Model

The most influential paper Gerard Salton never wrote

Wikipedia: tf-idf

Wikipedia: Vector space model

Scikits.learn Examples