00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #ifndef OST_SEQ_ALG_DISTANCE_MAP_HH
00025 #define OST_SEQ_ALG_DISTANCE_MAP_HH
00026
00027 #include <numeric>
00028 #include <cmath>
00029 #include <algorithm>
00030
00031 #include <ost/tri_matrix.hh>
00032 #include <ost/integrity_error.hh>
00033 #include <ost/seq/alignment_handle.hh>
00034 #include <ost/seq/alg/module_config.hh>
00035
00036 namespace ost { namespace seq { namespace alg {
00037
00039
00040 namespace {
00041
00042 template <typename Pair>
00043 struct sum_first {
00044 typename Pair::first_type operator()(const typename Pair::first_type& x,
00045 const Pair& y) const
00046 {
00047 return x + y.first;
00048 }
00049 };
00050
00051 template <typename Pair>
00052 struct CompareFirst {
00053 bool operator()(const Pair& x, const Pair& y) const {
00054 return (x.first < y.first);
00055 }
00056 };
00057
00058 }
00060
00063 struct DLLEXPORT_OST_SEQ_ALG Distances {
00064
00065 typedef std::vector<std::pair<Real,int> > FloatArray;
00066
00067 uint GetDataSize() const { return values_.size(); }
00068
00069 void Add(Real dist, int index) {
00070 values_.push_back(std::make_pair(dist,index));
00071 }
00072
00073 Real GetAverage() const {
00074 if (values_.empty()) throw IntegrityError("List of distances empty!");
00075 Real sum = std::accumulate(values_.begin(), values_.end(), 0.0,
00076 sum_first<FloatArray::value_type>());
00077 return sum/values_.size();
00078 }
00079
00080 std::pair<Real,int> GetMin() const {
00081 if (values_.empty()) throw IntegrityError("List of distances empty!");
00082 return *std::min_element(values_.begin(), values_.end(),
00083 CompareFirst<FloatArray::value_type>());
00084 }
00085
00086 std::pair<Real,int> GetMax() const {
00087 if (values_.empty()) throw IntegrityError("List of distances empty!");
00088 return *std::max_element(values_.begin(), values_.end(),
00089 CompareFirst<FloatArray::value_type>());
00090 }
00091
00092 std::pair<Real,int> GetDataElement(uint index) const {
00093 if (values_.size() != 0) {
00094 if (values_.size() > index) {
00095 return values_[index];
00096 } else {
00097 std::stringstream error_message;
00098 error_message << "size of distances data vector: " << values_.size()
00099 << " is smaller than requested: " << index+1;
00100 throw IntegrityError(error_message.str());
00101 }
00102 } else {
00103 throw IntegrityError("List of distances empty!");
00104 }
00105 }
00106
00107 Real GetStdDev() const {
00108 if (values_.empty()) throw IntegrityError("List of distances empty!");
00109 const Real avg = this->GetAverage();
00110 return this->do_get_stddev(avg);
00111 }
00112
00113 Real GetWeightedStdDev(Real sigma) const {
00114 if (values_.empty()) throw IntegrityError("List of distances empty!");
00115 const Real avg = this->GetAverage();
00116 const Real weight = std::exp(avg/(-2*sigma));
00117 return this->do_get_stddev(avg) * weight;
00118 }
00119
00120 Real GetNormStdDev() const {
00121 if (values_.empty()) throw IntegrityError("List of distances empty!");
00122 const Real avg = this->GetAverage();
00123 return this->do_get_stddev(avg) / avg;
00124 }
00125
00126 bool operator==(const Distances& rhs) const { return values_ == rhs.values_; }
00127 bool operator!=(const Distances& rhs) const { return values_ != rhs.values_; }
00128
00129 private:
00130
00131 Real do_get_stddev(Real avg) const {
00132 Real var = 0.0;
00133 for (FloatArray::const_iterator i = values_.begin(),
00134 e = values_.end(); i != e; ++i) {
00135 const Real diff = i->first - avg;
00136 var += diff*diff;
00137 }
00138 return std::sqrt(var/values_.size());
00139 }
00140
00141 FloatArray values_;
00142 };
00143
00147 class DistanceMap;
00148 typedef boost::shared_ptr<DistanceMap> DistanceMapPtr;
00149
00150 class DLLEXPORT_OST_SEQ_ALG DistanceMap {
00151 public:
00152 DistanceMap(uint nresidues, uint num_structures)
00153 : dist_(nresidues), num_structures_(num_structures) { }
00154
00155 void AddDistance(uint i, uint j, Real distance, int index) {
00156 dist_(i ,j).Add(distance,index);
00157 }
00158
00159 const Distances& GetDistances(uint i, uint j) const {
00160 return dist_(i, j);
00161 }
00162
00163 uint GetSize() const { return dist_.GetSize(); }
00164
00165 uint GetNumStructures() const { return num_structures_; }
00166
00167 private:
00168 TriMatrix<Distances> dist_;
00169 uint num_structures_;
00170 };
00171
00172
00182 DistanceMapPtr DLLEXPORT_OST_SEQ_ALG
00183 CreateDistanceMap(const seq::AlignmentHandle& alignment);
00184
00185 }}}
00186
00187 #endif