Ed's Big Plans

Computing for Science and Awesome

Archive for the ‘linkedin’ tag

SOM Indexing Logic

without comments

Update: Added some code too.

Brief: A classmate and I started talking about how we implemented the Kohonen Self Organizing Maps (SOM)s. I used C and indexed first the rows and the columns of the SOM before the index corresponding to the weight vectors (same as the index for the input vectors); he used C++ and indexed the weight vectors first before the columns and the rows.

Either way, we could use a three-deep array like this (switching the indexers as appropriate) …

const double low = 0.0; // minimum allowed random value to initialize weights
const double high = 1.0; // maximum allowed random value to initialize weights
const int nrows = 4; // number of rows in the map
const int ncolumns = 4; // number of columns in the map
const int ninputs = 3; // number of elements in an input vector, each weight vector

double*** weight; // the weight array

weight = calloc(nrows, sizeof(double**));
for(int r = 0; r < nrows; r ++) {
    weight[r] = calloc(ncolumns, sizeof(double*));
    for(int c = 0; c < ncolumns; c ++) {
        weight[r] = calloc(ninputs, sizeof(double));
        for(int i = 0; i < ninputs; i ++) {
            weight[r][i] =
            (((double)random() / (double)INT_MAX) * (high - low)) + low;

In the below diagram, the left side is a schematic of his approach and the right side is a schematic of my approach.

SOM Indexing

Figure above: SOM Indexing — Left (his): SOM indexed as input, row, column; Right (mine): SOM indexed as row, column, input.

Both schematics in the above diagram have four rows and four columns in the map where each weight (and input) vector has three elements.

I think my logic is better since we’ll often be using some distance function to evaluate how similar a weight vector is to a given input — to me, it’s natural to thus index these at the inner most nesting while looping over the rows and columns of the map.

The opposing approach was apparently used because my classmate had previously developed a matrix manipulation library. I’m actually kind of curious to take a look at it later.

Eddie Ma

March 31st, 2011 at 3:04 pm

Posted in Machine Learning

Tagged with , , ,

SOM in Regression Data (animation)

with 2 comments

The below animation is what happens when you train a Kohonen Self-Organizing Map (SOM) using data meant for regression instead of clustering or classification. A total of 417300 training cycles are performed — a single training cycle is the presentation of one exemplar to the SOM. The SOM has dimensions 65×65. Exemplars are chosen at random and 242 snapshots are taken in total (once every 1738 training cycles).

(animation in full post here …)

Eddie Ma

March 17th, 2011 at 7:52 pm

Python’s List Multiply — Use List Comprehension Instead

without comments

Brief: I ran into a snag parsing a CSV today — I have lines upon lines of 27 integers each. Each line represents a cube of values — each cube has dimensions 3×3×3. Here’s an example line and the corresponding cube it represents.

A line …

0, 0, 0, 0, 0, 0, 0, 50, 0, 2, 22, 0, 0, 4, 0, 5, 0, 17, 0, 26, 24, 0, 0, 0, 0, 0, 0,

The corresponding cube …

0 0 0
0 0 0
0 50 0
2 22 0
0 4 0
5 0 17
0 26 24
0 0 0
0 0 0

… Where the three major squares represents three levels of depth; each major square has three rows and three columns — the values in each depth are organized row-major.

So let’s say the lines in a file are read from stdin thanks to the magic of posix pipes.

We might be able to use a construct like this to store all of our lines in the variable cases as below …

import sys
cases = []
for line in sys.stdin:
    entry = [int(i) for i in line[:-1].split(",") if i != '']
    current_cases = [[[[0] * 3] * 3] * 3] # Doesn't work ...
    for ii, i in enumerate(entry):
        current_cases[ii/9%3][ii/3%3][ii%3] = i
        # ii is the index (of enumeration), i is the value (from the 'entry')

Let’s break apart the line that assigns entry first — it’s equal to the below …

entry = line[:-1]
# get rid of the trailing newline character

entry = entry.split(",")
# break apart the line by comma character

entry = [int(i) for i in entry if i != '']
# convert each element into an integer
# except for the trailing empty string

Now let’s break apart the line that assigns current_cases — here’s the logic behind getting each element of the 27-element cube …

[ i x x i x x i x x i x x i x x i x x i x x i x x i x x ] - column index
[ i i i x x x x x x i i i x x x x x x i i i x x x x x x ] - row index
[ i i i i i i i i i x x x x x x x x x x x x x x x x x x ] - depth index

The columns are indexed by [index (mod 3)] — every third item falls in the same column. The rows are indexed by [index /3 (mod 3)] — every run of three items out of nine are part of the same row. The levels of depth are indexed by [index /9 (mod 3)] — every run of nine items out of the 27 elements is part of the same depth.

So here’s the line that doesn’t work

current_cases = [[[[0] * 3] * 3] * 3]

In current_cases as defined above, there are only three integers — not 27 that are allocated in memory. The above saves cubes where the depth index and the row index don’t matter, only the inner-most column index has any meaning. The strange thing is, this wasn’t immediately intuitive to me — I expected the list multiply operation to create nested lists — instead, each of the two enclosing levels of lists just creates additional references to the same list.

This is the line we must use instead to make the nested lists correctly save the cubes …

current_cases = [[[0 for i in xrange(3)] for i in xrange(3)] for i in xrange(3)]
# or ...
current_cases = [[[0] * 3] for i in xrange(3)] for i in xrange(3)]
# or ...
current_cases = [[[0, 0, 0] for i in xrange(3)] for i in xrange(3)]

Here, the comprehension iteratively creates the correct 27 memory locations needed for each cube. Each nested comprehension is responsible for creating a unique list at the correct depth.

Putting it together, the correct listing for this task is …

import sys
cases = []
for line in sys.stdin:
entry = [int(i) for i in line[:-1].split(",") if i != '']
    current_cases = [[[0 for i in xrange(3)] for i in xrange(3)] for i in xrange(3)]
    for ii, i in enumerate(entry):
        current_cases[ii/9%3][ii/3%3][ii%3] = i

I’ve written this solution down this time because I seem to rediscover it every time I encounter it. Hopefully, this will save you some work too 😀

Eddie Ma

March 5th, 2011 at 11:15 pm

C & Bioinformatics: ASCII Nucleotide Comparison Circuit

with 4 comments

Here’s a function I developed for Andre about a week ago. The C function takes two arguments. Both arguments are C characters. The first argument corresponds to a degenerate nucleotide as in the below table. The second argument corresponds to a non-degenerate nucleotide {‘A’, ‘C’, ‘G’, ‘T’} or any nucleotide ‘N’. The function returns zero if the logical intersection between the two arguments is zero. The function returns a one-hot binary encoding for the logical intersection if it exists so that {‘A’ = 1, ‘C’ = 2, ‘G’ = 4, ‘T’ = 8} and {‘N’ = 15}. All of this is done treating the lower 5-bits of the ASCII encoding for each character as wires of a circuit.

Character ASCII
(Low 5-Bits)
Represents One Hot Equals
A 0001
B 00010 CGT 1110
C 00011 C 0010
D 00100 AGT 1101
G 00111 G 0100
H 01000 ACT 1011
01011 GT 1100
M 01101 AC 0011
N 01110 ACGT 1111
R 10010 AG 0101 5
S 10011 GC 0110
T 10100 T
V 10110 ACG
W 10111 AT
Y 11001 CT

The premise is that removing all of the logical branching and using only binary operators would make things a bit faster — I’m actually not sure if the following solution is faster because there are twelve variables local to the function scope — we can be assured that at least half of these variables will be stored outside of cache and will have to live in memory. We’d get a moderate speed boost if at all.

        Bitwise comparison circuit that treats nucleotide and degenerate
        nucleotide ascii characters as mathematical sets.
        The operation performed is set i intersect j.
    arg char i:
        Primer -- accepts ascii |ABCDG HKMNR STVWY| = 15.
    arg char j:
        Sequence -- accepts ascii |ACGTN| = 5.
    return char (k = 0):
        false -- i intersect j is empty.
    return char (k > 0):
        1 -- the intersection is 'A'
        2 -- the intersection is 'C'
        4 -- the intersection is 'G'
        8 -- the intersection is 'T'
       15 -- the intersection is 'N'
    return char (undefined value):
        ? -- if any other characters are placed in i or j.
char f(char i, char j) {

    // break apart Primer into bits ...
    char p = (i >> 4) &1;
    char q = (i >> 3) &1;
    char r = (i >> 2) &1;
    char s = (i >> 1) &1;
    char t =  i       &1;

    // break apart Sequence into bits ...
    char a = (j >> 4) &1;
    char b = (j >> 3) &1;
    char c = (j >> 2) &1;
    char d = (j >> 1) &1;
    char e =  j       &1;


    ( // == PRIMER CIRCUIT ==
        ( // -- A --
            ((p|q|r|s  )^1) &         t |
            ((p|q|  s|t)^1) &     r     |
            ((p|  r|s|t)^1) &   q       |
            ((p|    s  )^1) &   q&r&  t |
            ((p|      t)^1) &   q&r&s   |
            ((  q|    t)^1) & p&s       |
            ((  q      )^1) & p&r&s
        ( // -- C --
            ((p|q|r    )^1) &       s   |
            ((p|  r|s|t)^1) &   q       |
            ((p|    s  )^1) &   q&r&  t |
            ((p|      t)^1) &   q&r&s   |
            ((  q|r    )^1) &       s&t |
            ((  q|    t)^1) & p  &r&s   |
            ((    r|s  )^1) & p&q&    t
        ) << 1
        ( // -- G --
            ((  q|r|  t)^1) &       s   |
            ((p|q|  s|t)^1) &     r     |
            ((p|q      )^1) &     r&s&t |
            ((p|  r    )^1) &   q&  s&t |
            ((p|      t)^1) &   q&r&s   |
            ((  q|r    )^1) & p&    s   |
            ((  q|    t)^1) & p&    s
        ) << 2
        ( // -- T --
            ((p|q|r|  t)^1) &       s   |
            ((  q|  s|t)^1) &     r     |
            ((p|  r|s|t)^1) &   q       |
            ((p|  r    )^1) &   q&  s&t |
            ((p|      t)^1) &   q&r&s   |
            ((  q      )^1) & p&  r&s&t |
            ((    r|s  )^1) & p&q&    t
        ) << 3
    ( // == SEQUENCE CIRCUIT ==
        ( // -- A --
            ((a|b|c|d  )^1) &         e |
            ((a|      e)^1) &   b&c&d
        ( // -- C --
            ((a|b|c    )^1) &       d&e |
            ((a|      e)^1) &   b&c&d
        ) << 1
        ( // -- G --
            ((a|b      )^1) &     c&d&e |
            ((a|      e)^1) &   b&c&d
        ) << 2
        ( // -- T --
            ((a|      e)^1) &   b&c&d   |
            ((  b|  d|e)^1) & a&  c
        ) << 3

Andre’s eventual solution was to use a look-up table which very likely proves faster in practice. At the very least, this was a nice refresher and practical example for circuit logic using four sets of minterms (one for each one-hot output wire).

Should you need this logic to build a fast physical circuit or have some magical architecture with a dozen registers (accessible to the compiler), be my guest 😀

Eddie Ma

February 27th, 2011 at 11:49 pm

The Right Tool (Why I chose Java to do RSA)

without comments

Brief: I learned something valuable last week when working on this RSA encryption/decryption assignment for my Computer Security class. It’s important to be versatile when doing computer science — we must ensure we always use the most efficient tool. If we aren’t versatile, we risk taking tremendous amounts of time trying to reimplement something that already exists elsewhere.

So, what tool is the right tool to quickly throw together RSA encryption?

It turns out that Java does an excellent job. Its BigInteger class has all the ingredients you’ll ever need.

// This function generates a new probable prime ...
BigInteger p = BigInteger.probablePrime(bits, random);

// This function performs modulus inverse ...
BigInteger d = e.modInverse(phi_n);

// These functions can be used to check your work ...
BigInteger one = d.multiply(e).mod(phi_n);
BigInteger one = p.gcd(q);

// This function is at the heart of RSA ...
BigInteger C = M.modPow(d, n);

Before I looked at the Java documentation, I had plans to do this with Python and some of my classmates had plans with MATLAB. It’s not that these are inherently bad technologies — they’re just not the right tool.

As much as I have gripes with Java, I’ve got to say that this made me and a lot of my class very happy 😀

Eddie Ma

February 24th, 2011 at 6:26 pm

Posted in Pure Programming

Tagged with , ,

Move a Subdirectory to a Subdomain (Apache mod_alias)

without comments

In this post, I’m going to explain how to move resources located at a subdirectory of your domain to its own subdomain. The reason why I came across this was because I decided to move my MediaWiki installation to better separate the URLs corresponding to my blog and my notebook — this is paired with better separation of the actual files and directories corresponding to my blog and notebook as they exist on my server. I’ll assume that you’re running Apache; I’ll further assume that you’ve already gotten your A-RECORD set up with your domain registrar or name server provider (DNS) — if you don’t know what that means, ask me.

I’ll use my MediaWiki move as my running example for this post.

Here are some properties we want a redirect to satisfy — I’ll lead you through to making a redirect such that each of these conditions are met.

  1. Moving the subdirectory to a subdomain corresponds to moving files on the server filesystem.
    • (URL) http://eddiema.ca/wiki = ~/Sites/eddiema/notebook (server filesystem)
    • (URL) http://wiki.eddiema.ca = ~/Sites/notebook (server filesystem)
  2. Requests to the subdomain must work.
    • http://wiki.eddiema.ca
  3. Requests to the subdirectory must be forwarded to the subdomain.
    • http://eddiema.ca/wiki →
    • http://wiki.eddiema.ca
  4. Requests to subdirectories of the former subdirectory must be forwarded to subdirectories of the subdomain.
    • http://eddiema.ca/wiki/index.php/Special:AllPages
    • http://wiki.eddiema.ca/index.php/Special:AllPages
  5. URLs with GET requests to the subdirectory must be forwarded to the subdomain without losing the GET requests.
    • http://eddiema.ca/wiki/?search=Notes+CIS+6050
    • http://wiki.eddiema.ca/?search=Notes+CIS+6050

Let’s meet each of our above properties.

Changing the URL corresponds to a physical move (Property 1)

Move your subdirectory on the server to its target location on the server.

For my example, I moved my notebook out of my main site directory.

mv ~/Sites/eddiema/notebook ~/Sites/notebook

(If you’re moving a MediaWiki installation like me, there’s one more thing you have to do — this is explained at the end.)

Requests to the subdomain must work (Property 2)

Update your Apache httpd.conf configuration file with a new VirtualHost entry. This new entry stipulates a new document root for your target subdomain.

#<VirtualHost *:80>
#    ServerName {full URL to your new subdomain}
#    DocumentRoot {full path to your new document root}

# Example ...
<VirtualHost *:80>
    ServerName wiki.eddiema.ca
    DocumentRoot /Users/eddiema/Sites/notebook

Requests to the subdirectory must be forwarded to the subdomain (Properties 3~5)

There are two ways to accomplish this task — you can either modify Apache’s httpd.conf file, or you can create a new .htaccess file. I chose the latter because I wanted to keep all of the changes related to this move localized.

To accomplish this task, you must perform two steps.

First, you must recreate the subdirectory on your server filesystem that you have just moved elsewhere. We do this so that a request for that directory causes Apache to look there for a .htaccess file in the correct place. This .htaccess file will tell Apache to redirect to your new subdomain instead.

For my example, I recreated the subdirectory with…

mkdir ~/Sites/eddiema/notebook

Second, Create a new plaintext .htaccess file in the recreated subdirectory with a RedirectMatch rule.

#RedirectMatch 301 ^/{subdirectory}/(.*)$ http://{full subdomain URL}/$1

# Example ...
RedirectMatch 301 ^/wiki/(.*)$ http://wiki.eddiema.ca/$1

RedirectMatch is part of mod_alias. If you’re doing something more sophisticated than just forwarding URLs, then you will need to use mod_redirect which has the ability to manipulate query strings.

The number 301 tells browsers and search engines that the URL pointing at your subdirectory is no longer valid and has permanently moved to the one pointing at your new subdomain. The next part is a PERL style regular expression. Let’s break it apart.

  • ^ and $ match the start of, and end of the string respectively.
  • .* means we want to match zero or more characters in a string.
  • (.*) means we want to save the matched string as an autovariable $1.

The last string is appended with $1 to push the tail of any URL of your subdirectory onto the tail of your subdomain.

A final note for moving a MediaWiki installation

In order to stop MediaWiki from appending the subdirectory “wiki” to the tail of your new subdomain (“http://wiki.{your domain}.*/wiki“), you can change the variable $wgScriptPath to the empty string in LocalSettings.php found in the root directory of your MediaWiki installation.

#$wgScriptPath       = "/wiki";
$wgScriptPath       = "";

Note that MediaWiki doesn’t seem to like this idea — but they don’t really say explain why.

After this step, you can use Short URLs (same as above link) to get rid of index.php appearing in the URL, but I won’t go into that today.

That’s everything 😀

Eddie Ma

February 23rd, 2011 at 12:36 pm