Alphabet.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) 2006-2009 Soeren Sonnenburg
 * Copyright (C) 2006-2009 Fraunhofer Institute FIRST and Max-Planck-Society
 */

#ifndef _CALPHABET__H__
#define _CALPHABET__H__

#include "base/SGObject.h"
#include "lib/Mathematics.h"
#include "lib/common.h"

namespace shogun
{
enum EAlphabet
{
    DNA=0,

    RAWDNA=1,

    RNA=2,

    PROTEIN=3,

    // BINARY just 0 and 1
    BINARY=4,

    ALPHANUM=5,

    CUBE=6,

    RAWBYTE=7,

    IUPAC_NUCLEIC_ACID=8,

    IUPAC_AMINO_ACID=9,

    NONE=10,

    DIGIT=11,

    DIGIT2=12,

    RAWDIGIT=13,

    RAWDIGIT2=14,

    UNKNOWN=15,

    SNP=16,

    RAWSNP=17
};


class CAlphabet : public CSGObject
{
    public:

        //CAlphabet();


        CAlphabet(char* alpha, int32_t len);

        CAlphabet(EAlphabet alpha);

        CAlphabet(CAlphabet* alpha);
        virtual ~CAlphabet();

        bool set_alphabet(EAlphabet alpha);

        inline EAlphabet get_alphabet() const
        {
            return alphabet;
        }

        inline int32_t get_num_symbols() const
        {
            return num_symbols;
        }

        inline int32_t get_num_bits() const
        {
            return num_bits;
        }

        inline uint8_t remap_to_bin(uint8_t c)
        {
            return maptable_to_bin[c];
        }

        inline uint8_t remap_to_char(uint8_t c)
        {
            return maptable_to_char[c];
        }

        void clear_histogram();

        template <class T>
        void add_string_to_histogram(T* p, int64_t len)
        {
            for (int64_t i=0; i<len; i++)
                add_byte_to_histogram((uint8_t) (p[i]));
        }

        inline void add_byte_to_histogram(uint8_t p)
        {
            histogram[p]++;
        }

        void print_histogram();

        inline void get_hist(int64_t** h, int32_t* len)
        {
            int32_t hist_size=(1 << (sizeof(uint8_t)*8));
            ASSERT(h && len);
            *h=(int64_t*) malloc(sizeof(int64_t)*hist_size);
            ASSERT(*h);
            *len=hist_size;
            ASSERT(*len);
            memcpy(*h, &histogram[0], sizeof(int64_t)*hist_size);
        }

        inline const int64_t* get_histogram()
        {
            return &histogram[0];
        }

        bool check_alphabet(bool print_error=true);

        inline bool is_valid(uint8_t c)
        {
            return valid_chars[c];
        }

        bool check_alphabet_size(bool print_error=true);

        int32_t get_num_symbols_in_histogram();

        int32_t get_max_value_in_histogram();

        int32_t get_num_bits_in_histogram();

        static const char* get_alphabet_name(EAlphabet alphabet);


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

        template <class ST>
        static void translate_from_single_order(ST* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val)
        {
            int32_t i,j;
            ST value=0;

            for (i=sequence_length-1; i>= p_order-1; i--) //convert interval of size T
            {
                value=0;
                for (j=i; j>=i-p_order+1; j--)
                    value= (value >> max_val) | (obs[j] << (max_val * (p_order-1)));

                obs[i]= (ST) value;
            }

            for (i=p_order-2;i>=0;i--)
            {
                if (i>=sequence_length)
                    continue;

                value=0;
                for (j=i; j>=i-p_order+1; j--)
                {
                    value= (value >> max_val);
                    if (j>=0 && j<sequence_length)
                        value|=obs[j] << (max_val * (p_order-1));
                }
                obs[i]=value;
            }

            // TODO we should get rid of this loop!
            if (start>0)
            {
                for (i=start; i<sequence_length; i++)
                    obs[i-start]=obs[i];
            }
        }

        template <class ST>
        static void translate_from_single_order_reversed(ST* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val)
        {
            int32_t i,j;
            ST value=0;

            for (i=sequence_length-1; i>= p_order-1; i--) //convert interval of size T
            {
                value=0;
                for (j=i; j>=i-p_order+1; j--)
                    value= (value << max_val) | obs[j];

                obs[i]= (ST) value;
            }

            for (i=p_order-2;i>=0;i--)
            {
                if (i>=sequence_length)
                    continue;

                value=0;
                for (j=i; j>=i-p_order+1; j--)
                {
                    value= (value << max_val);
                    if (j>=0 && j<sequence_length)
                        value|=obs[j];
                }
                obs[i]=value;
            }

            // TODO we should get rid of this loop!
            if (start>0)
            {
                for (i=start; i<sequence_length; i++)
                    obs[i-start]=obs[i];
            }
        }

        template <class ST>
        static void translate_from_single_order(ST* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
        {
            ASSERT(gap>=0);

            const int32_t start_gap=(p_order-gap)/2;
            const int32_t end_gap=start_gap+gap;

            int32_t i,j;
            ST value=0;

            // almost all positions
            for (i=sequence_length-1; i>=p_order-1; i--) //convert interval of size T
            {
                value=0;
                for (j=i; j>=i-p_order+1; j--)
                {
                    if (i-j<start_gap)
                    {
                        value= (value >> max_val) | (obs[j] << (max_val * (p_order-1-gap)));
                    }
                    else if (i-j>=end_gap)
                    {
                        value= (value >> max_val) | (obs[j] << (max_val * (p_order-1-gap)));
                    }
                }
                obs[i]= (ST) value;
            }

            // the remaining `order` positions
            for (i=p_order-2;i>=0;i--)
            {
                if (i>=sequence_length)
                    continue;

                value=0;
                for (j=i; j>=i-p_order+1; j--)
                {
                    if (i-j<start_gap)
                    {
                        value= (value >> max_val);
                        if (j>=0 && j<sequence_length)
                            value|=obs[j] << (max_val * (p_order-1-gap));
                    }
                    else if (i-j>=end_gap)
                    {
                        value= (value >> max_val);
                        if (j>=0 && j<sequence_length)
                            value|=obs[j] << (max_val * (p_order-1-gap));
                    }
                }
                obs[i]=value;
            }

            // TODO we should get rid of this loop!
            if (start>0)
            {
                for (i=start; i<sequence_length; i++)
                    obs[i-start]=obs[i];
            }
        }

        template <class ST>
        static void translate_from_single_order_reversed(ST* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
        {
            ASSERT(gap>=0);

            const int32_t start_gap=(p_order-gap)/2;
            const int32_t end_gap=start_gap+gap;

            int32_t i,j;
            ST value=0;

            // almost all positions
            for (i=sequence_length-1; i>=p_order-1; i--) //convert interval of size T
            {
                value=0;
                for (j=i; j>=i-p_order+1; j--)
                {
                    if (i-j<start_gap)
                        value= (value << max_val) | obs[j];
                    else if (i-j>=end_gap)
                        value= (value << max_val) | obs[j];
                }
                obs[i]= (ST) value;
            }

            // the remaining `order` positions
            for (i=p_order-2;i>=0;i--)
            {
                if (i>=sequence_length)
                    continue;

                value=0;
                for (j=i; j>=i-p_order+1; j--)
                {
                    if (i-j<start_gap)
                    {
                        value= value << max_val;
                        if (j>=0 && j<sequence_length)
                            value|=obs[j];
                    }
                    else if (i-j>=end_gap)
                    {
                        value= value << max_val;
                        if (j>=0 && j<sequence_length)
                            value|=obs[j];
                    }
                }
                obs[i]=value;
            }

            // TODO we should get rid of this loop!
            if (start>0)
            {
                for (i=start; i<sequence_length; i++)
                    obs[i-start]=obs[i];
            }
        }


    protected:
        void init_map_table();

        void copy_histogram(CAlphabet* src);



#ifdef HAVE_BOOST_SERIALIZATION
    protected:

        friend class ::boost::serialization::access;

        template<class Archive>

            void serialize(Archive & ar, const unsigned int archive_version)
            {

                SG_DEBUG("archiving CAlphabet\n");

                ar & ::boost::serialization::base_object<CSGObject>(*this);

                ar & num_bits;

                SG_DEBUG("done with CAlphabet\n");

            }


#endif //HAVE_BOOST_SERIALIZATION


    public:
        static const uint8_t B_A;
        static const uint8_t B_C;
        static const uint8_t B_G;
        static const uint8_t B_T;
        static const uint8_t B_0;
        static const uint8_t MAPTABLE_UNDEF;
        static const char* alphabet_names[18];

    protected:
        EAlphabet alphabet;
        int32_t num_symbols;
        int32_t num_bits;
        bool valid_chars[1 << (sizeof(uint8_t)*8)];
        uint8_t maptable_to_bin[1 << (sizeof(uint8_t)*8)];
        uint8_t maptable_to_char[1 << (sizeof(uint8_t)*8)];
        int64_t histogram[1 << (sizeof(uint8_t)*8)];
};


#ifndef DOXYGEN_SHOULD_SKIP_THIS
template<> inline void CAlphabet::translate_from_single_order(float32_t* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
{
}

template<> inline void CAlphabet::translate_from_single_order(float64_t* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
{
}

template<> inline void CAlphabet::translate_from_single_order(floatmax_t* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
{
}

template<> inline void CAlphabet::translate_from_single_order_reversed(float32_t* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
{
}

template<> inline void CAlphabet::translate_from_single_order_reversed(float64_t* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
{
}

template<> inline void CAlphabet::translate_from_single_order_reversed(floatmax_t* obs, int32_t sequence_length, int32_t start, int32_t p_order, int32_t max_val, int32_t gap)
{
}
#endif

}

#ifdef HAVE_BOOST_SERIALIZATION

namespace boost {
namespace serialization {
   template<class Archive>
    inline void save_construct_data(Archive & ar, shogun::CAlphabet* t, const unsigned int archive_version)
    {
      std::cout << "archiving construction data CAlphabet" << std::endl;
      ar << t->alphabet;
      std::cout << "done archiving construction data CAlphabet" << std::endl;
    }

    template<class Archive>
    inline void load_construct_data(Archive & ar, shogun::CAlphabet * t, const unsigned int archive_version)
    {
      //SG_DEBUG("archiving construction data CAlphabet\n");
      shogun::EAlphabet a;
      ar >> a;
      new(t)shogun::CAlphabet(a);
      //SG_DEBUG("done construction data archiving CAlphabet\n");
    }
} // serialization
} // namespace boost
#endif //HAVE_BOOST_SERIALIZATION



#endif