bmalgo.h

ncbi源码

/*
 * ===========================================================================
 * PRODUCTION $Log: bmalgo.h,v $
 * PRODUCTION Revision 1000.0  2004/04/21 15:59:59  gouriano
 * PRODUCTION PRODUCTION: IMPORTED [CATCHUP_003] Dev-tree R1.1
 * PRODUCTION
 * ===========================================================================
 */

/*
Copyright (c) 2003 Anatoliy Kuznetsov.

Permission is hereby granted, free of charge, to any person 
obtaining a copy of this software and associated documentation 
files (the "Software"), to deal in the Software without restriction, 
including without limitation the rights to use, copy, modify, merge, 
publish, distribute, sublicense, and/or sell copies of the Software, 
and to permit persons to whom the Software is furnished to do so, 
subject to the following conditions:

The above copyright notice and this permission notice shall be included 
in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, 
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES 
OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. 
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, 
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR 
OTHER DEALINGS IN THE SOFTWARE.
*/

#ifndef BMALGO__H__INCLUDED__
#define BMALGO__H__INCLUDED__

#include "bm.h"
#include "bmfunc.h"
#include "bmdef.h"

namespace bm
{

/*! \defgroup setalgo Set algorithms 
 *  Set algorithms 
 *
 */

/*! \defgroup distance Distance metrics 
 *  Algorithms to compute binary distance metrics
 *  \ingroup setalgo
 */


/*! 
    \brief    Distance metrics codes defined for vectors A and B
    \ingroup  distance
*/
enum distance_metric
{
    COUNT_AND,       //!< (A & B).count()
    COUNT_XOR,       //!< (A ^ B).count()
    COUNT_OR,        //!< (A | B).count()
    COUNT_SUB_AB,    //!< (A - B).count()
    COUNT_SUB_BA,    //!< (B - A).count()
    COUNT_A,         //!< A.count()
    COUNT_B          //!< B.count()
};

/*! 
    \brief Distance metric descriptor, holds metric code and result.
    \sa distance_operation
*/

struct distance_metric_descriptor
{
     distance_metric   metric;
     bm::id_t          result;
     
     distance_metric_descriptor(distance_metric m)
     : metric(m),
       result(0)
    {}
    distance_metric_descriptor()
    : metric(bm::COUNT_XOR),
      result(0)
    {}
    
    /*! 
        \brief Sets metric result to 0
    */
    void reset()
    {
        result = 0;
    }
};


/*!
    \brief Internal function computes different distance metrics.
    \internal 
    \ingroup  distance
     
*/
inline
void combine_count_operation_with_block(const bm::word_t* blk,
                                        unsigned gap,
                                        const bm::word_t* arg_blk,
                                        int arg_gap,
                                        bm::word_t* temp_blk,
                                        distance_metric_descriptor* dmit,
                                        distance_metric_descriptor* dmit_end)
                                            
{
     gap_word_t* res=0;
     
     gap_word_t* g1 = BMGAP_PTR(blk);
     gap_word_t* g2 = BMGAP_PTR(arg_blk);
     
     if (gap) // first block GAP-type
     {
         if (arg_gap)  // both blocks GAP-type
         {
             gap_word_t tmp_buf[bm::gap_max_buff_len * 3]; // temporary result
             
             for (distance_metric_descriptor* it = dmit;it < dmit_end; ++it)
             {
                 distance_metric_descriptor& dmd = *it;
                 
                 switch (dmd.metric)
                 {
                 case bm::COUNT_AND:
                     res = gap_operation_and(g1, g2, tmp_buf);
                     break;
                 case bm::COUNT_OR:
                     res = gap_operation_or(g1, g2, tmp_buf);
                     break;
                 case bm::COUNT_SUB_AB:
                     res = gap_operation_sub(g1, g2, tmp_buf); 
                     break;
                 case bm::COUNT_SUB_BA:
                     res = gap_operation_sub(g2, g1, tmp_buf); 
                     break;
                 case bm::COUNT_XOR:
                     res = gap_operation_xor(g1, g2, tmp_buf); 
                    break;
                 case bm::COUNT_A:
                    res = g1;
                    break;
                 case bm::COUNT_B:
                    res = g2;
                    break;
                 } // switch
                 
                 if (res)
                     dmd.result += gap_bit_count(res);
                    
             } // for it
             
             return;

         }
         else // first block - GAP, argument - BITSET
         {
             for (distance_metric_descriptor* it = dmit;it < dmit_end; ++it)
             {
                 distance_metric_descriptor& dmd = *it;
                 
                 switch (dmd.metric)
                 {
                 case bm::COUNT_AND:
                     if (arg_blk)
                        dmd.result += gap_bitset_and_count(arg_blk, g1);
                     break;
                 case bm::COUNT_OR:
                     if (!arg_blk)
                        dmd.result += gap_bit_count(g1);
                     else
                        dmd.result += gap_bitset_or_count(arg_blk, g1); 
                     break;
                 case bm::COUNT_SUB_AB:
                     gap_convert_to_bitset((bm::word_t*) temp_blk, 
                                           g1, bm::set_block_size);
                     dmd.result += 
                       bit_operation_sub_count((bm::word_t*)temp_blk, 
                          ((bm::word_t*)temp_blk) + bm::set_block_size,
                           arg_blk);
                 
                     break;
                 case bm::COUNT_SUB_BA:
                     dmd.metric = bm::COUNT_SUB_AB; // recursive call to SUB_AB
                     combine_count_operation_with_block(arg_blk,
                                                        arg_gap,
                                                        blk,
                                                        gap,
                                                        temp_blk,
                                                        it, it+1);
                     dmd.metric = bm::COUNT_SUB_BA; // restore status quo
                     break;
                 case bm::COUNT_XOR:
                     if (!arg_blk)
                        dmd.result += gap_bit_count(g1);
                     else
                        dmd.result += gap_bitset_xor_count(arg_blk, g1);
                     break;
                 case bm::COUNT_A:
                    if (g1)
                        dmd.result += gap_bit_count(g1);
                    break;
                 case bm::COUNT_B:
                    if (arg_blk)
                    {
                        dmd.result += 
                          bit_block_calc_count(arg_blk, 
                                               arg_blk + bm::set_block_size);
                    }
                    break;
                 } // switch
                                     
             } // for it
             
             return;
         
         }
     } 
     else // first block is BITSET-type
     {     
         if (arg_gap) // second argument block is GAP-type
         {
             for (distance_metric_descriptor* it = dmit;it < dmit_end; ++it)
             {
                 distance_metric_descriptor& dmd = *it;
                 
                 switch (dmd.metric)
                 {
                 case bm::COUNT_AND:
                     if (blk) 
                        dmd.result += gap_bitset_and_count(blk, g2);                         
                     break;
                 case bm::COUNT_OR:
                     if (!blk)
                        dmd.result += gap_bit_count(g2);
                     else
                        dmd.result += gap_bitset_or_count(blk, g2);
                     break;
                 case bm::COUNT_SUB_AB:
                     if (blk)
                        dmd.result += gap_bitset_sub_count(blk, g2);
                     break;
                 case bm::COUNT_SUB_BA:
                     dmd.metric = bm::COUNT_SUB_AB; // recursive call to SUB_AB
                     combine_count_operation_with_block(arg_blk,
                                                        arg_gap,
                                                        blk,
                                                        gap,
                                                        temp_blk,
                                                        it, it+1);
                     dmd.metric = bm::COUNT_SUB_BA; // restore status quo
                     break;
                 case bm::COUNT_XOR:
                     if (!blk)
                        dmd.result += gap_bit_count(g2);
                     else
                        dmd.result += gap_bitset_xor_count(blk, g2); 
                    break;
                 case bm::COUNT_A:
                    if (blk)
                    {
                        dmd.result += 
                            bit_block_calc_count(blk, 
                                                 blk + bm::set_block_size);
                    }
                    break;
                 case bm::COUNT_B:
                    if (g2)
                        dmd.result += gap_bit_count(g2);
                    break;
                 } // switch
                                     
             } // for it
             
             return;
         }
     }

     // --------------------------------------------
     //
     // Here we combine two plain bitblocks 

     const bm::word_t* blk_end;
     const bm::word_t* arg_end;

     blk_end = blk + (bm::set_block_size);
     arg_end = arg_blk + (bm::set_block_size);


     for (distance_metric_descriptor* it = dmit; it < dmit_end; ++it)
     {
         distance_metric_descriptor& dmd = *it;

         switch (dmd.metric)
         {
         case bm::COUNT_AND:
             dmd.result += 
                bit_operation_and_count(blk, blk_end, arg_blk);
             break;
         case bm::COUNT_OR:
             dmd.result += 
                bit_operation_or_count(blk, blk_end, arg_blk);
             break;
         case bm::COUNT_SUB_AB:
             dmd.result += 
                bit_operation_sub_count(blk, blk_end, arg_blk);
             break;
         case bm::COUNT_SUB_BA:
             dmd.result += 
                bit_operation_sub_count(arg_blk, arg_end, blk);
             break;
         case bm::COUNT_XOR:
             dmd.result += 
                bit_operation_xor_count(blk, blk_end, arg_blk);
             break;
         case bm::COUNT_A:
            if (blk)
                dmd.result += bit_block_calc_count(blk, blk_end);
            break;
         case bm::COUNT_B:
            if (arg_blk)
                dmd.result += bit_block_calc_count(arg_blk, arg_end);
            break;
         } // switch

     } // for it
     
     

}

/*!
    \brief Distance computing template function.

    Function receives two bitvectors and an array of distance metrics
    (metrics pipeline). Function computes all metrics saves result into
    corresponding pipeline results (distance_metric_descriptor::result)
    An important detail is that function reuses metric descriptors, 
    incrementing received values. It allows you to accumulate results 
    from different calls in the pipeline.
    
    \param bv1      - argument bitvector 1 (A)
    \param bv2      - argument bitvector 2 (B)
    \param dmit     - pointer to first element of metric descriptors array
                      Input-Output parameter, receives metric code as input,
                      computation is added to "result" field
    \param dmit_end - pointer to (last+1) element of metric descriptors array
    \ingroup  distance
    
*/

template<class BV>
void distance_operation(const BV& bv1, 
                        const BV& bv2, 
                        distance_metric_descriptor* dmit,
                        distance_metric_descriptor* dmit_end)
{
    const typename BV::blocks_manager& bman1 = bv1.get_blocks_manager();
    const typename BV::blocks_manager& bman2 = bv2.get_blocks_manager();
    
    bm::word_t* temp_blk = 0;
    
    {
        for (distance_metric_descriptor* it = dmit; it < dmit_end; ++it)
        {
            if (it->metric == bm::COUNT_SUB_AB || 
                it->metric == bm::COUNT_SUB_BA)
            {
                temp_blk = BV::allocate_tempblock();
                break;
            }
        }
    }
    
  
    bm::word_t*** blk_root = bman1.get_rootblock();
    unsigned block_idx = 0;
    unsigned i, j;
    
    const bm::word_t* blk;
    const bm::word_t* arg_blk;
    bool  blk_gap;
    bool  arg_gap;

    BM_SET_MMX_GUARD

    for (i = 0; i < bman1.top_block_size(); ++i)
    {
        bm::word_t** blk_blk = blk_root[i];

        if (blk_blk == 0) // not allocated
        {
            const bm::word_t* const* bvbb = bman2.get_topblock(i);
            if (bvbb == 0) 
            {
                block_idx += bm::set_array_size;
                continue;
            }

            blk = 0;
            blk_gap = false;

            for (j = 0; j < bm::set_array_size; ++j,++block_idx)