BitString.h

Go to the documentation of this file.

/*
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or
 * (at your option) any later version.
 *
 * Written (W) 2009 Soeren Sonnenburg
 * Copyright (C) 2009 Fraunhofer Institute FIRST and Max Planck Society
 */
#ifndef __BITSTRING_H__
#define __BITSTRING_H__

#include <shogun/features/Alphabet.h>
#include <shogun/lib/common.h>
#include <shogun/lib/io.h>
#include <shogun/lib/MemoryMappedFile.h>
#include <shogun/lib/Mathematics.h>

namespace shogun
{
class CBitString : public CSGObject
{
    public:
        CBitString(EAlphabet alpha, int32_t width=1) : CSGObject(), string(NULL), length(0)
        {
            alphabet=new CAlphabet(alpha);
            int32_t nbits=alphabet->get_num_bits();
            word_len = width*nbits;

            mask=0;
            for (int32_t j=0; j<word_len; j++)
                mask=(mask<<1) | (uint64_t) 1;
            mask<<=sizeof(uint64_t)*8-word_len;

            single_mask=0;
            for (int32_t j=0; j<nbits; j++)
                mask=(mask<<1) | (uint64_t) 1;
        }

        ~CBitString()
        {
            cleanup();
            SG_UNREF(alphabet);
        }

        void cleanup()
        {
            delete[] string;
            string=NULL;
            length=0;
        }

        void obtain_from_char(char* str, uint64_t len)
        {
            cleanup();
            uint64_t stream_len=len/sizeof(uint64_t)+1;
            string=new uint64_t[stream_len];
            length=len;

            uint64_t w=0;
            int32_t nbits=alphabet->get_num_bits();
            uint64_t j=0;
            int32_t nfit=8*sizeof(w)/nbits;
            for (uint64_t i=0; i<len; i++)
            {
                w= (w << nbits) | alphabet->remap_to_bin((uint8_t) str[j]);

                if (i % nfit == nfit-1)
                {
                    string[j]=w;
                    j++;
                }
            }

            if (j<stream_len)
            {
                string[j]=w;
                j++;
            }

            ASSERT(j==stream_len);
        }

        void load_fasta_file(const char* fname, bool ignore_invalid=false)
        {
            cleanup();

            uint64_t len=0;
            uint64_t offs=0;

            CMemoryMappedFile<char> f(fname);

            uint64_t id_len=0;
            char* id=f.get_line(id_len, offs);

            if (!id_len || id[0]!='>')
                SG_SERROR("No fasta hunks (lines starting with '>') found\n");

            if (offs==f.get_size())
                SG_SERROR("Empty file?\n");

            char* fasta=NULL;
            char* s=NULL;
            int32_t spanned_lines=0;
            int32_t fasta_len=0;

            while (true)
            {
                s=f.get_line(len, offs);

                if (!fasta)
                    fasta=s;

                if (!s || len==0)
                    SG_SERROR("Error reading fasta entry in line %d len=%ld", spanned_lines+1, len);

                if (s[0]=='>')
                    SG_SERROR("Multiple fasta hunks (lines starting with '>') are not supported!\n");

                if (offs==f.get_size())
                {
                    offs-=len+1; // seek to beginning
                    fasta_len+=len;

                    uint64_t w=0;
                    int32_t nbits=alphabet->get_num_bits();
                    uint64_t nfit=8*sizeof(w)/nbits;

                    len = fasta_len-spanned_lines;
                    uint64_t stream_len=len/(nfit)+1;
                    string=new uint64_t[stream_len];
                    length=len;

                    int32_t idx=0;
                    int32_t k=0;

                    for (int32_t j=0; j<fasta_len; j++, k++)
                    {
                        if (fasta[j]=='\n')
                        {
                            k--;
                            continue;
                        }

                        w= (w << nbits) | alphabet->remap_to_bin((uint8_t) fasta[j]);

                        if (k % nfit == nfit-1)
                        {
                            string[idx]=w;
                            w=0;
                            idx++;
                        }
                    }

                    if (idx<stream_len)
                        string[idx]=w<<(nbits*(nfit - k%nfit));

                    break;
                }

                spanned_lines++;
                fasta_len+=len+1; // including '\n'
            }
        }

        void set_string(uint64_t* str, uint64_t len)
        {
            cleanup();
            string=str;
            length=len;
        }

        void create(uint64_t len)
        {
            cleanup();
            uint64_t stream_len=len/sizeof(uint64_t)+1;
            string=new uint64_t[stream_len];
            CMath::fill_vector(string, (int32_t) stream_len, (uint64_t) 0);
            length=len;
        }

        /*
        inline uint64_t condense(uint64_t bits, uint64_t mask) const
            {
                uint64_t res = 0 ;
                uint64_t m=mask ;
                uint64_t b=bits ;
                char cnt=0 ;
                
                char ar[256][256] ;

                for (int i=0; i<8; i++)
                {
                    //fprintf(stdout, "%i %lx %lx %lx %i\n", i, res, m, b, (int)cnt) ;
                    if (m&1)
                        res = res>>8 | ar[b&255][m&255] ;
                    //else
                    //  cnt++ ;
                    m=m>>8 ;
                    b=b>>8 ;
                }
                res=res>>cnt ;
                //fprintf(stdout, "%lx %lx %lx\n", res, m, b) ;
                
                //res = res & bits & mask ;
                
                return res ;
            }
        */

        inline uint64_t operator[](uint64_t index) const
        {
            ASSERT(index<length);

            uint64_t bitindex=alphabet->get_num_bits()*index;
            int32_t ws=8*sizeof(uint64_t);
            uint64_t i=bitindex/ws;
            int32_t j=bitindex % ws;
            int32_t missing=word_len-(ws-j);

            //SG_SPRINT("i=%lld j=%d ws=%d word_len=%d missing=%d left=%llx shift=%d\n", i, j, ws, word_len, missing, ( string[i] << j ) & mask, ws-word_len);
            uint64_t res= ((string[i] << j) & mask ) >> (ws-word_len);

            if (missing>0)
                res|= ( string[i+1] >> (ws-missing) );

            //res = condense(res, 1<<31|1<<24|1<<10|1<<8|1<<4|1<<2|1<<1|1) ;

            return res;
        }

        inline void set_binary_word(uint16_t word, uint64_t index)
        {
            ASSERT(index<length);

            uint64_t bitindex=alphabet->get_num_bits()*index;
            int32_t ws=8*sizeof(uint64_t);
            uint64_t i=bitindex/ws;
            int32_t j=bitindex % ws;
            int32_t missing=word_len-(ws-j);

            uint64_t sl = j-word_len;
            uint64_t ml;
            uint64_t wl;
            uint64_t mr;
            uint64_t wr;

            if (sl>0)
            {
                ml=mask<<sl;
                wl=word<<sl ;
            }
            else
            {
                sl=-sl;
                ml=mask>>sl;
                wl=word>>sl;
            }

            string[i] = ( string[i] & (~ml) ) | ( wl & ml);

            if (missing>0)
            {
                mr = mask<<missing;
                wr = word<<missing;
                string[i+1] = ( string[i+1] & (~mr) ) | ( wr & mr);
            }

        }

        inline uint64_t get_length() const { return length-word_len/alphabet->get_num_bits()+1; }

        inline virtual const char* get_name() const { return "BitString"; }

    private:
        CAlphabet* alphabet;
        uint64_t* string;
        uint64_t length;
        int32_t word_len;
        uint64_t mask;
        uint64_t single_mask;
};
}
#endif //__BITSTRING_H__