Quick and Dirty histogram in awk

awk '{histogram[int(($2*100)/100)]++;}END{for (i in histogram) { print i, histogram[i] }}' somerandomfilewithatagandcountdataincol2 | sort -n -k 1 > hist

Print every other line of a file

cat blah | awk 'BEGIN{n=0}{if(n%2==0) print $0;n++;}'

Getting started with Cython on Centos 5.4 64Bit

Centos is a nightmare. OK so first up I wanted to do with with Python 2.6 for which there is no official Centos package. So get the geekymedia repos which has a version of Python 2.6 which will install along side the existing 2.4 (and can be access as python26) here: here.

Then download the Cython tarball. Before you build Cython make sure you setup compile flags for 32bitness and make sure python 2.6 gets used i.e.

export CFLAGS="-m32 -march=i386"
export PYTHON="python26"

then run python26 setup.py install in the Cython directory. I had to do some other hacking of -m32 and -march=i386 on to things, but I think doing it this way should work. That should build and install Cython.

Put the following example in a file called demo.pyx:

cdef extern from "math.h":
    double sin(double)

cdef double f(double x):
    return sin(x)

print f(2.0)

Then make the following Makefile. Make sure you have tabs!

demo: demo.pyx
        cython demo.pyx --embed
        gcc -g -pthread -m32 -march=i386 demo.c -L/usr/lib -lm -lpython2.6 -o ./demo -I/usr/include/python2.6

Then just type make and it should build a binary called demo.

Digital Archaeology Tags

It's fairly clear that in the future archaeology will move in to the digital domain with researchers digging through hard discs, USB keys, floppy discs, tapes, EEPROMs etc. I've often though it would be nice to drop them a message every so often, or just leave a note apologizing for my appalling spelling and grammar. So I suggest the following tags, to be encoded in 8bit ASCII when possible:

DIGIARCHNOTESTART

DIGIARCHNOTEEND

Dear future archaeologists, I shall use this in my documents from now on when I want to give you some background, point you in the direction of more information, or just say Hi, so grep away!

Grab all hrefs from a html page which have text containing View associated with them

Quick and dirty perl program to grab all links from a webpage which have have anchors with the text "View" in them:

#!/usr/bin/perl

use constant false => 0;
use constant true  => 1;

use HTML::TreeBuilder;
use HTML::FormatText;

$html = HTML::TreeBuilder->new();
$html->parse_file($ARGV[0]);

my @stuff = $html->look_down( '_tag' , 'a' );

my $seqtag = "";
my $use_next = false;

for my $i (@stuff) {
  my @thing = $i->content();
  my $target = $i->attr('href');
  my $str = $thing[0][0];

  # string contains View
  if($str =~ m/View/) {
    print $target . "\n";
  }
}

samtools, bam cigar to string cigar

This code fragment converts a bam style binary cigar string in to a normal cigar string.

    CigarLen = b->core.n_cigar;  // Cigar length
    uint32_t *cigar    = bam1_cigar(b);

    Cigar.clear();
    for(size_t n=0;n<cigar_len;n++) {
      uint32_t opcode = cigar[n] & 0x0000000F;
      uint32_t len    = cigar[n] & 0xFFFFFFF0;

      len = len >> 4;

      if(opcode == 0)  { Cigar += stringify(len); Cigar += "M"; }
      if(opcode == 1)  { Cigar += stringify(len); Cigar += "I"; }
      if(opcode == 2)  { Cigar += stringify(len); Cigar += "D"; }
      if(opcode == 3)  { Cigar += stringify(len); Cigar += "N"; }
      if(opcode == 4)  { Cigar += stringify(len); Cigar += "S"; }
      if(opcode == 5)  { Cigar += stringify(len); Cigar += "H"; }
      if(opcode == 6)  { Cigar += stringify(len); Cigar += "P"; }
      if(opcode == 7)  { Cigar += stringify(len); Cigar += "a"; }
      if(opcode == 8)  { Cigar += stringify(len); Cigar += "b"; }
      if(opcode == 9)  { Cigar += stringify(len); Cigar += "c"; }
      if(opcode == 10) { Cigar += stringify(len); Cigar += "d"; }
      if(opcode == 11) { Cigar += stringify(len); Cigar += "e"; }
      if(opcode == 12) { Cigar += stringify(len); Cigar += "f"; }
      if(opcode == 13) { Cigar += stringify(len); Cigar += "g"; }
      if(opcode == 14) { Cigar += stringify(len); Cigar += "h"; }
    }
  }

string to int conversion

I was recently asked to write a string to int conversion function (without using library functions). I initially came up with a solution using the pow function (which is quite expensive). I had a think about it and found there were a surprising number of solutions. Briefly I came up with the following methods:

Method	Summary
Pow	my initial pow based solution (after converting a position to an in calculating 10^val)
Mul	Rather than using pow generating powers of 10 in the loop (multiplier = multipler*10)
Table	Just use a lookup table for the multipliers
Case	More of less the same as table, but encode the table in a switch statement (very ugly!)

Results are probably compiler/CPU/platform dependent. But on my Atom Z530 (1.6GHz) based netbook using GCC 4.3.3 I obtained the following results when performing the conversion 10 million times:

Method	User time
Pow	43.67s
Mul	28.21s
Table	28.22s
Case	29.13s

There's a big difference between the pow method and the others, but I was reasonably surprised that multiplier and table based methods performed similarly. It would be interesting to look at the assembler generated for these.

For reference, source code follows (note I sum and output the converted values to prevent the call to string_to_int from being optimised away). I was slightly concerned that something funky /might/ be going on in string::size() however benchmarked with this in and outside the loop and didn't observe any difference. Note: Following programs don't process signs, but in terms of benchmarking I don't believe this should be relevant.

Pow:

#include <string>
#include <iostream>
#include <math.h>
#include <stdlib.h>

using namespace std;

int string_to_int(string s) {

  int output=0;

  for(int n=0;n<s.size();n++) {
    int cval = s[n]-'0';
    output += cval*pow(10,s.size()-n-1);
  }

  return output;
}

int main() {

  // Simple tests
  cout << "1 is: " << string_to_int("1") << endl;
  cout << "10 is: " << string_to_int("10") << endl;
  cout << "14532 is: " << string_to_int("14532") << endl;

  int rsum=0;
  for(int i=0;i<10000000;i++) {
    string s;
    int numlen = rand()%11;
    for(int n=0;n<numlen;n++) {
      int rval;
      if((numlen==10) && (n==0)) { rval = rand()%2;  }
                           else  { rval = rand()%10; }
      s.push_back('0'+rval);
    }
    int v = string_to_int(s);
    rsum += v;
  }
  cout << rsum << endl;
}

Mul:

#include <string>
#include <iostream>
#include <math.h>
#include <stdlib.h>

using namespace std;

int string_to_int(string s) {

  int output=0;

  int mul=1;
  for(int n=s.size()-1;n>=0;n--) {
    int cval = s[n]-'0';
    output += cval*mul;
    mul = mul * 10;
  }

  return output;
}

int main() {

  // Simple tests
  cout << "1 is: " << string_to_int("1") << endl;
  cout << "10 is: " << string_to_int("10") << endl;
  cout << "14532 is: " << string_to_int("14532") << endl;

  int rsum=0;
  for(int i=0;i<10000000;i++) {
    string s;
    int numlen = rand()%11;
    for(int n=0;n<numlen;n++) {
      int rval;
      if((numlen==10) && (n==0)) { rval = rand()%2;  }
                           else  { rval = rand()%10; }
      s.push_back('0'+rval);
    }
    int v = string_to_int(s);
    rsum += v;
  }
  cout << rsum << endl;
}

Table:

#include <string>
#include <iostream>
#include <math.h>
#include <stdlib.h>

using namespace std;

const int powtable [] = { 1,
                          10,
                          100,
                          1000,
                          10000,
                          100000,
                          1000000,
                          10000000,
                          100000000,
                          1000000000
                        };

int string_to_int(string s) {

  int output=0;

  for(int n=s.size()-1;n>=0;n--) {
    int cval = s[n]-'0';
    output += cval*(powtable[s.size()-n-1]);
  }

  return output;
}

int main() {

  // Simple tests
  cout << "1 is: " << string_to_int("1") << endl;
  cout << "10 is: " << string_to_int("10") << endl;
  cout << "14532 is: " << string_to_int("14532") << endl;

  int rsum=0;
  for(int i=0;i<10000000;i++) {
    string s;
    int numlen = rand()%11;
    for(int n=0;n<numlen;n++) {
      int rval;
      if((numlen==10) && (n==0)) { rval = rand()%2;  }
                           else  { rval = rand()%10; }
      s.push_back('0'+rval);
    }
    int v = string_to_int(s);
    rsum += v;
  }
  cout << rsum << endl;
}

Case:

#include <string>
#include <iostream>
#include <math.h>
#include <stdlib.h>

using namespace std;

int string_to_int(string s) {

  int output=0;

  int pos=0;
  for(int n=s.size()-1;n>=0;n--) {
    int cval = s[n]-'0';

    switch(pos) {
      case 0:
        output += cval*1;
        break;
      case 1:
        output += cval*10;
        break;
      case 2:
        output += cval*100;
        break;
      case 3:
        output += cval*1000;
        break;
      case 4:
        output += cval*10000;
        break;
      case 5:
        output += cval*100000;
        break;
      case 6:
        output += cval*1000000;
        break;
      case 7:
        output += cval*10000000;
        break;
      case 9:
        output += cval*100000000;
        break;
      case 10:
        output += cval*1000000000;
        break;
    }

    pos++;
  }

  return output;
}

int main() {

  // Simple tests
  cout << "1 is: " << string_to_int("1") << endl;
  cout << "10 is: " << string_to_int("10") << endl;
  cout << "14532 is: " << string_to_int("14532") << endl;

  int rsum=0;
  for(int i=0;i<10000000;i++) {
    string s;
    int numlen = rand()%11;
    for(int n=0;n<numlen;n++) {
      int rval;
      if((numlen==10) && (n==0)) { rval = rand()%2;  }
                           else  { rval = rand()%10; }
      s.push_back('0'+rval);
    }
    int v = string_to_int(s);
    rsum += v;
  }
  cout << rsum << endl;
}