source: branches/lib/GDE/MUSCLE/src/fastdistkbit.cpp

#include "muscle.h"
#include "distfunc.h"
#include "seqvect.h"
#include <math.h>

#define MIN(x, y)       ((x) < (y) ? (x) : (y))

static void SetKmerBitVector(const Seq &s, byte Bits[])
        {
        const unsigned uLength = s.Length();
        const unsigned k = 3;   // kmer length
        unsigned i = 0;
        unsigned c = 0;
        unsigned h = 0;
        for (unsigned j = 0; j < k - 1; ++j)
                {
                unsigned x = CharToLetterEx(s[i++]);
                if (x <= AX_Y)
                        c = c*20 + x;
                else
                        {
                        c = 0;
                        h = j + 1;
                        }
                }
        for ( ; i < uLength; ++i)
                {
                unsigned x = CharToLetterEx(s[i++]);
                if (x <= AX_Y)
                        c = (c*20 + x)%8000;
                else
                        {
                        c = 0;
                        h = i + k;
                        }
                if (i >= h)
                        {
                        unsigned ByteOffset = c/8;
                        unsigned BitOffset = c%8;
                        Bits[ByteOffset] |= (1 << BitOffset);
                        }
                }
        }

static unsigned CommonBitCount(const byte Bits1[], const byte Bits2[])
        {
        const byte * const p1end = Bits1 + 1000;
        const byte *p2 = Bits2;

        unsigned uCount = 0;
        for (const byte *p1 = Bits1; p1 != p1end; ++p1)
                {
        // Here is a cute trick for efficiently counting the
        // bits common between two bytes by combining them into
        // a single word.
                unsigned b = *p1 | (*p2 << 8);
                while (b != 0)
                        {
                        if (b & 0x101)
                                ++uCount;
                        b >>= 1;
                        }
                ++p2;
                }
        return uCount;
        }

void DistKbit20_3(const SeqVect &v, DistFunc &DF)
        {
        const unsigned uSeqCount = v.Length();
        DF.SetCount(uSeqCount);

// There are 20^3 = 8,000 distinct kmers in the 20-letter alphabet.
// For each sequence, we create a bit vector of length 8,000, i.e.
// 1,000 bytes, having one bit per kmer. The bit is set to 1 if the
// kmer is present in the sequence.
        const unsigned uBytes = uSeqCount*1000;
        byte *BitVector = new byte[uBytes];
        memset(BitVector, 0, uBytes);

        SetProgressDesc("K-bit distance matrix");
        for (unsigned uSeqIndex = 0; uSeqIndex < uSeqCount; ++uSeqIndex)
                SetKmerBitVector(*v[uSeqIndex], BitVector + uSeqIndex*1000);

        unsigned uDone = 0;
        const unsigned uTotal = (uSeqCount*(uSeqCount - 1))/2;
        for (unsigned uSeqIndex1 = 0; uSeqIndex1 < uSeqCount; ++uSeqIndex1)
                {
                const byte *Bits1 = BitVector + uSeqIndex1*1000;
                const unsigned uLength1 = v[uSeqIndex1]->Length();
                for (unsigned uSeqIndex2 = 0; uSeqIndex2 < uSeqIndex1; ++uSeqIndex2)
                        {
                        const byte *Bits2 = BitVector + uSeqIndex2*1000;
                        const unsigned uLength2 = v[uSeqIndex2]->Length();
                        const float fCount = (float) CommonBitCount(Bits1, Bits2);

                // Distance measure = K / min(L1, L2)
                // K is number of distinct kmers that are found in both sequences
                        const float fDist = fCount / MIN(uLength1, uLength2);
                        DF.SetDist(uSeqIndex1, uSeqIndex2, fDist);
                        if (uDone%10000 == 0)
                                Progress(uDone, uTotal);
                        ++uDone;
                        }
                }
        ProgressStepsDone();

        delete[] BitVector;
        }