00001 #ifndef OST_MOL_ALG_ADJACENCY_BITMAP_HH
00002 #define OST_MOL_ALG_ADJACENCY_BITMAP_HH
00003
00004 #include <cstring>
00005 #include <boost/shared_ptr.hpp>
00006 #include <ost/img/image_handle.hh>
00007 #include <ost/tri_matrix.hh>
00008 #include <ost/seq/sequence_handle.hh>
00009 #include <ost/seq/alignment_handle.hh>
00010 #include <ost/stdint.hh>
00011 #include "similarity_matrix.hh"
00012 #include "contact_overlap.hh"
00013 #include "module_config.hh"
00014
00015 namespace ost { namespace mol { namespace alg {
00016
00017
00018 struct OverlapResult {
00019 OverlapResult(uint16_t nom, uint16_t denom, uint16_t def):
00020 nominator(nom), denominator(denom), defined(def) {}
00021 OverlapResult(): nominator(0), denominator(0), defined(0) {}
00022 uint16_t nominator;
00023 uint16_t denominator;
00024 uint16_t defined;
00025 };
00026
00027
00028 class AdjacencyBitmap {
00029 public:
00030 explicit AdjacencyBitmap(size_t n_vertices):
00031 num_vertices_(n_vertices),
00032 storage_size_(StorageSize(n_vertices))
00033 {
00034 data_ = static_cast<uint64_t*>(malloc(this->num_bytes()));
00035 end_ = data_+storage_size_*num_vertices_;
00036 memset(data_, 0, this->num_bytes());
00037 }
00038
00039 AdjacencyBitmap(const DistanceMatrix& dmat_a,
00040 const DistanceMatrix& dmat_b,
00041 Real threshold, Real radius,
00042 int start, int end):
00043 num_vertices_(end-start),
00044 storage_size_(StorageSize(end-start)),
00045 data_(static_cast<uint64_t*>(malloc(this->num_bytes()))),
00046 end_(data_+storage_size_*num_vertices_) {
00047
00048 for (int i = start; i < end; ++i) {
00049 this->set(i-start, i-start, false, false);
00050 for (int j = i+1; j < end; ++j) {
00051 Real da = dmat_a.Get(i, j);
00052 Real db = dmat_b.Get(i, j);
00053 bool defined = da>=0 && db>=0 && (db < radius || da < radius);
00054 bool agree = false;
00055 if (defined) {
00056 agree = std::abs(da-db)<threshold;
00057 }
00058 this->set(i-start, j-start, agree, defined);
00059 this->set(j-start, i-start, agree, defined);
00060 }
00061 }
00062 }
00063
00064 AdjacencyBitmap(const SimilarityMatrix& sim, Real threshold):
00065 num_vertices_(sim.GetSize()),
00066 storage_size_(StorageSize(sim.GetSize())),
00067 data_(static_cast<uint64_t*>(malloc(this->num_bytes()))),
00068 end_(data_+storage_size_*num_vertices_) {
00069 memset(data_, 0, this->num_bytes());
00070 for (int i = 0; i < sim.GetSize(); ++i) {
00071 this->set(i, i, false, false);
00072 for (int j = i+1; j < sim.GetSize(); ++j) {
00073 this->set(i, j, sim.Get(i, j)>threshold, sim.Get(i, j)>=0.0);
00074 this->set(j, i, sim.Get(i, j)>threshold, sim.Get(i, j)>=0.0);
00075 }
00076 }
00077 }
00078
00079 AdjacencyBitmap(const AdjacencyBitmap& rhs):
00080 num_vertices_(rhs.num_vertices_),
00081 storage_size_(rhs.storage_size_),
00082 data_(static_cast<uint64_t*>(malloc(this->num_bytes()))),
00083 end_(data_+storage_size_*num_vertices_) {
00084 memcpy(data_, rhs.data_, this->num_bytes());
00085 }
00086
00087 ~AdjacencyBitmap() {
00088 if (data_)
00089 free(data_);
00090 }
00091
00092
00093
00094 OverlapResult Overlap(size_t vert_i, size_t vert_j) const;
00095
00096 void set(size_t i, size_t j, bool val, bool defined=true) {
00097 size_t byte_index = i*storage_size_ + j/(4*sizeof(uint64_t));
00098 assert(byte_index < this->num_bytes());
00099 size_t bit_index = (2*j) % BITS;
00100
00101
00102 uint64_t one = 1;
00103 uint64_t two = 2;
00104 assert(bit_index < sizeof(uint64_t)*8);
00105
00106 if (val) {
00107 data_[byte_index] |= (one << bit_index);
00108 } else {
00109 data_[byte_index] &= ~(one << bit_index);
00110 }
00111
00112 if (defined) {
00113 data_[byte_index] |= (two << bit_index);
00114 } else {
00115 data_[byte_index] &= ~(two << bit_index);
00116 }
00117 }
00118
00119 bool get(size_t i, size_t j) const {
00120 uint64_t one = 1;
00121 size_t byte_index = i*storage_size_ + j/(4*sizeof(uint64_t));
00122 assert(byte_index < storage_size_*num_vertices_);
00123 size_t bit_index = (2*j) % BITS;
00124 assert(bit_index < sizeof(uint64_t)*8);
00125 return data_[byte_index] & (one << bit_index);
00126 }
00127
00128 bool defined(size_t i, size_t j) const {
00129 uint64_t two = 2;
00130 size_t byte_index = i*storage_size_ + j/(4*sizeof(uint64_t));
00131 assert(byte_index < storage_size_*num_vertices_);
00132 size_t bit_index = (2*j) % BITS;
00133 assert(bit_index < sizeof(uint64_t)*8);
00134 return data_[byte_index] & (two << bit_index);
00135 }
00136
00137 size_t size() const {
00138 return num_vertices_;
00139 }
00140
00141 size_t num_bytes() const {
00142 return storage_size_*num_vertices_*sizeof(uint64_t);
00143 }
00144 void dump() const;
00145 private:
00146
00147 AdjacencyBitmap();
00148 AdjacencyBitmap& operator=(const AdjacencyBitmap&);
00149
00150 const static size_t BITS=sizeof(uint64_t)*8;
00151 size_t num_vertices_;
00152 size_t storage_size_;
00153 uint64_t* data_;
00154 uint64_t* end_;
00155
00156 static const uint8_t NUMBER_OF_BITS_SET[];
00157 static const uint8_t NUMBER_OF_EDGES_SET[];
00158 static size_t StorageSize(size_t verts) {
00159 return verts/(4*sizeof(uint64_t)) + (verts % (sizeof(uint64_t)*4) ? 1 : 0);
00160 }
00161 };
00162
00163
00164 }}}
00165 #endif
00166