Sum weights

I have a file with a bunch of sequences and some weights at the top of the file:

>WEIGHTS 0.926434 1.000000 1.000000 0.926434 1.000000 0.892712 1.000000 1.000000 1.000000 1.000000 1.000000 0.892712 
>CRTC_EUGGR__Q9ZNY3 Calreticulin precursor.
XRKELWXXXXXXXXXXXXXXXXXXXXXXXXXTRWTHSTXXSDYXKFKLTSGKFYGDKAKDAGIQTSQDAKFYAISSPIASXXSXEXXXLVLQFSVKHXXXXXXGXGXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXEPRCEXDTLSHTYXAXXXXDXXXEVLVDQVKKESGTLEEDWEILKPKTIPDPEDKKPADWVDEPDMVDPEDKKPEDWDKEPAQIPDPDATQPDDWDEEEDGKWEAPMISNPKYKGEWKAKKIPNPAYKGVWKPRDIPNPEYEADDKVXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXFYDQTNGATKDAEKKAFDSAEADKRKKEEDERKKQEEEEKKTAEEDEXXXDEXXXEDDKKDEL
>HSP47_RAT__P29457 47 kDa heat shock protein precursor (Collagen-binding protein 1) (GP46).
XRSLXXXXXXXXXXXXXXXXXXXXEAAAPGTAEKLSSKATTLAEXSXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXSQAKAVLSAEKLRDEEVHTGLGELVRSLSNXTARNVTWKLGSXXXXXXXXSFADDFVRSSKQHYNCEHSKINFRDKRSALQSINEWASQTTDGKLPEVTKDVERTDXXLLXXAMXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXYXYXXXXXXXXQXVEMXXXXXXXXXXXXXXXXXXXXXRLEKXXTKXQLKTWMGKMQKKAXXISLPXGVVXVTHDLQKXXAGLGLTEAIXKNKADLSXXSGXXXXXXXXXXXXXXXEWDTEGNPFDQDIYGRXXXRSXXXXXXXXXXXXXXXXXXXXXXXXIGRLXXXXGDKMRDEL
>ENPL_PIG__Q29092 Endoplasmin precursor (94 kDa glucose-regulated protein) (GRP94) (GP96 homolog) (98 kDa protein kinase) (PPK 98) (ppk98).
XRAXXXXXXXXXXXXXXXXXXXXEVDVDGTVEEDLGKSREGSRTDDEVVQREEEAIQLDGLNASQIRELREKSEKXAFXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXELTVKXKCDKEKNLLHVTDTGVGMTREELVKNLGTIAKSGTSEFLNKMAEAQEDGQSTSELIGXXXXXXXXXXXXXXXXXVTXXHNNDTQHIWESDSNEFSVIADPRGNTLGRGTTITLXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXSXKTETVEXPMXXXXAAKXEKEESDDEAAVXXXXXEKXPXTXKVEKTVWDWELMNDIKPIWQRPSKEVEDDEYKAFYKSFSXXXXXPMAYIHFTXXXXXXXXXILXXXXXXXXXLFDEYXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXLNVSREXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXGVIXDHXXXXRLAKLLRFQSSHHPSDITSLDQYVERMKEKQDKIYFMAGSSRKEAESSPFVEXXXXXXXXXXXXXXXXXXXXXQALPXXXXKRFQNVAKEGVKFDESEKSKENREAVEKEFEPLLNWMKDKALXDKIEKAVVSQRXXEXXXXLVASQYGWXGNXERIMKAQAYQTGKDISTNYYASQKKTFEINPRHPLIRDMLRRVKEDEDDKXXSDLXXXXXXXXXXXXXXLLPDTKAYXXRIERMLRLSLNIDPDAKVEXXPXXXPXXTTEDTTEDTEQDDDEEMDAGADEXXQXTSETSTAEKDEL
>CRTC_HUMAN__P27797 Calreticulin precursor (CRP55) (Calregulin) (HACBP) (ERp60) (grp60).
XLXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXTSXXIESXXXSDFXXFVLSSGKFYGDEEKDKGLQTSQDARFYALSASFEXXSXXXXXLVVXFXXKHXXXXXXGGGYVKLFPNSLDQTDMHGDSEYNIMFGPDIXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXTYEVKIDNSQVESGSLEDDWDFLPPKKIKDPDASKPEDWDERAKIDDPTDSKPEDWDKPEHIPDPDAKKPEDWDEEMDGEWEPPVIQNPEYKGEWKPRQIDNPDYKGTWIHPEIDNPEYSPDPSIYXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXAYAEEFGNETWGVTKAAEKQMKDKQDEEQRLKEEEEDKKRKEXXXAXDKEDXEXKXEDXXDXXDKXXDXXEDVPGQAKDEL

I want to sum the weights which is fine for this example with like 12 sequences. However, some of the files have a couple of hundred entries. Step in Bash:

 head -n 1 filename.txt | awk '{for (i=1; i<=NF; i++) s=s+$i}; END{print s}'

Giving me the sum of the weights! Woohoo.

find a random webserver using libcurl

This post has moved here

Metasploit simple example (3.5.0)

Download and install metasploit framwork from here: http://www.metasploit.com/framework/download/

run msfconsole type the following to run exploits:

db_driver sqlite3
db_connect
db_nmap [ip address]
db_autopwn -p -e -t

Unobservable alphabetic password entry

This post has moved HERE

Unobservable pin entry

This post has moved HERE

Grabbing environment variables in gnuplot

This post has moved here

vim magic

So I have a bunch of sequences in fasta format - and I need to rearrange the mo'fos.

To start with the entries look like this:

>F13C5.1        CE19383 WBGene00017422  status:Partially_confirmed      UniProt:O76564  protein_id:AAC64611.1

Run this vim command:

:%s:>\(\S\{4,}\)\t.*UniProt\:\(\S\{6,}\).*$:>\1_CAEEL__\2:g

And now they look like this:

>geneName_OrgID__UniProtAccNo

LaCie Network Space

So.. you've hacked your network space and got root. But your lazy and you can't find up to date instructions on installing ipkg?

wget http://ipkg.nslu2-linux.org/feeds/optware/mssii/cross/stable/mssii-bootstrap_1.2-7_arm.xsh

and then sh mssii-bootstrap_1.2-7_arm.xsh

Now you have ipkg goodness.

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;
}

Floating point precision

I was recently asked a question about the precision of floating point numbers. As most people know the representation used by most computers to store real numbers does not have continuous precision across the number line.

I thought it would be interesting therefore to plot the density of numbers across the number line. This is a nice illustration of the fact that numbers near 0 have higher precision that large numbers. The following plots are for 32bit floats (this was on a 64bit Snow Leopard Mac):




The following code was used to generate the data, which was plotted with Gnuplot:

#include <iostream>
#include <limits>

using namespace std;

int main() {


  //float smallest = 0.00000000000000000000001;
  float smallest = 0.00000001;

  float block_size = 10;
  for(float range_start=-100000;range_start<100000;range_start+=block_size) {
    float range_end = range_start+block_size;

    size_t different_n=0;
    float n;
    for(float n=range_start;n<range_end;) {

      float new_n=n;
      for(float i=1;new_n <= n;i++) {
        new_n = n + (smallest*i);
      }
      n = new_n;
      different_n++;
    }

    cout << range_start+((range_end-range_start)/2) << " " << different_n << endl;
  }
}

mpt2sas raid card from Dell r710 driver installation with Centos 5.3

The REDACTED at Dell don't support Centos 5.3 on r610s raid card (Dell H200). I seriously REDRACTED hate propriety drivers and it's going to be a constant nightmare as we change kernel versions, but whatever. The problem is basically that the drivers are signed for Redhat Enterprise and Centos doesn't like that. So you have to crack open the driver package and remove the signing. A fixed package is available here: i686 x86_64.

I robbed most of the instructions to do this from here:

http://centos.org/modules/newbb/viewtopic.php?topic_id=22284&forum=37

Replace all instance of string with an incrementing number

Perl solution with thanks to rmp. In this case replace all instances of MAGIC in the file small.sam:

cat small.sam | perl -pe 's/MAGIC/$c++/eg'

Pausing/Suspending a running process

I had a couple of analysis jobs running they appeared to be stepping all over each other, with both jobs trashing a lot. I wanted to suspend one of the jobs to allow the other to continue in peace. You can do this using kill:

kill -STOP processid

and to continue it later:

kill -CONT processid

Found here:

http://linuxpoison.blogspot.com/2007/11/temporarily-suspend-process.html

Very simple OpenCL example

This blog post has been moved HERE

.