1 /*
2 Copyright 2009-2012 Andreas Biegert, Christof Angermueller
3
4 This file is part of the CS-BLAST package.
5
6 The CS-BLAST package is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 The CS-BLAST package is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 #ifndef CS_PROFILE_INL_H_
21 #define CS_PROFILE_INL_H_
22
23 #include "profile.h"
24 #include <cstring>
25
26 namespace cs {
27
28 template <class Abc>
Profile()29 Profile<Abc>::Profile() : nn(0), v(NULL) {}
30
31 template <class Abc>
Profile(size_t n)32 Profile<Abc>::Profile(size_t n) : nn(n), v(n>0 ? new double*[n] : NULL) {
33 size_t i,nel=n*Abc::kSizeAny;
34 if (v) v[0] = nel>0 ? new double[nel] : NULL;
35 for (i=1;i<n;i++) v[i] = v[i-1] + Abc::kSizeAny;
36 }
37
38 template <class Abc>
Profile(size_t n,const double & a)39 Profile<Abc>::Profile(size_t n, const double &a)
40 : nn(n), v(n>0 ? new double*[n] : NULL) {
41 size_t i,j,nel=n*Abc::kSizeAny;
42 if (v) v[0] = nel>0 ? new double[nel] : NULL;
43 for (i=1; i< n; i++) v[i] = v[i-1] + Abc::kSizeAny;
44 for (i=0; i< n; i++) for (j=0; j<Abc::kSizeAny; j++) v[i][j] = a;
45 }
46
47 template <class Abc>
Profile(size_t n,const double * a)48 Profile<Abc>::Profile(size_t n, const double *a)
49 : nn(n), v(n>0 ? new double*[n] : NULL) {
50 size_t i,nel=n*Abc::kSizeAny;
51 if (v) v[0] = nel>0 ? new double[nel] : NULL;
52 for (i=1; i< n; i++) v[i] = v[i-1] + Abc::kSizeAny;
53 memcpy(v[0], a, nel * sizeof(double));
54 }
55
56 template <class Abc>
Profile(const Profile & rhs)57 Profile<Abc>::Profile(const Profile &rhs)
58 : nn(rhs.nn), v(nn>0 ? new double*[nn] : NULL) {
59 size_t i,nel=nn*Abc::kSizeAny;
60 if (v) v[0] = nel>0 ? new double[nel] : NULL;
61 for (i=1; i< nn; i++) v[i] = v[i-1] + Abc::kSizeAny;
62 memcpy(v[0], rhs[0], nel * sizeof(double));
63 }
64
65 template <class Abc>
Profile(const CountProfile<Abc> & cp)66 Profile<Abc>::Profile(const CountProfile<Abc>& cp)
67 : nn(cp.length()), v(nn>0 ? new double*[nn] : NULL) {
68 size_t i,j,nel=nn*Abc::kSizeAny;
69 if (v) v[0] = nel>0 ? new double[nel] : NULL;
70 for (i=1; i< nn; i++) v[i] = v[i-1] + Abc::kSizeAny;
71 for (i=0; i< nn; i++)
72 for (j=0; j<Abc::kSizeAny; j++)
73 v[i][j] = cp.neff[i] > 0 ? cp.counts[i][j] / cp.neff[i] : 0.0;
74 }
75
76
77 template <class Abc>
78 Profile<Abc> & Profile<Abc>::operator=(const Profile<Abc> &rhs) {
79 if (this != &rhs) {
80 size_t i,nel;
81 nel = rhs.nn*Abc::kSizeAny;
82 if (nn != rhs.nn) {
83 if (v != NULL) {
84 delete[] v[0];
85 delete[] v;
86 }
87 nn=rhs.nn;
88 v = nn>0 ? new double*[nn] : NULL;
89 if (v) v[0] = nel>0 ? new double[nel] : NULL;
90 for (i=1; i< nn; i++) v[i] = v[i-1] + Abc::kSizeAny;
91 }
92 memcpy(v[0], rhs[0], nel * sizeof(double));
93 }
94 return *this;
95 }
96
97 template <class Abc>
98 inline double* Profile<Abc>::operator[](const size_t i) {
99 #ifdef CHECKBOUNDS
100 if (i<0 || i>=nn) {
101 throw Exception("Profile subscript out of bounds");
102 }
103 #endif
104 return v[i];
105 }
106
107 template <class Abc>
108 inline const double* Profile<Abc>::operator[](const size_t i) const {
109 #ifdef CHECKBOUNDS
110 if (i<0 || i>=nn) {
111 throw Exception("Profile subscript out of bounds");
112 }
113 #endif
114 return v[i];
115 }
116
117 template <class Abc>
Resize(size_t newn)118 void Profile<Abc>::Resize(size_t newn) {
119 size_t i,nel;
120 if (newn != nn) {
121 if (v != NULL) {
122 delete[] v[0];
123 delete[] v;
124 }
125 nn = newn;
126 v = nn>0 ? new double*[nn] : NULL;
127 nel = nn*Abc::kSizeAny;
128 if (v) v[0] = nel>0 ? new double[nel] : NULL;
129 for (i=1; i< nn; i++) v[i] = v[i-1] + Abc::kSizeAny;
130 }
131 }
132
133 template <class Abc>
Assign(size_t newn,const double & a)134 void Profile<Abc>::Assign(size_t newn, const double& a) {
135 size_t i,j,nel;
136 if (newn != nn) {
137 if (v != NULL) {
138 delete[] v[0];
139 delete[] v;
140 }
141 nn = newn;
142 v = nn>0 ? new double*[nn] : NULL;
143 nel = nn*Abc::kSizeAny;
144 if (v) v[0] = nel>0 ? new double[nel] : NULL;
145 for (i=1; i< nn; i++) v[i] = v[i-1] + *Abc::kSizeAny;
146 }
147 for (i=0; i< nn; i++) for (j=0; j<*Abc::kSizeAny; j++) v[i][j] = a;
148 }
149
150 template<class Abc>
Insert(size_t idx,const Profile<Abc> & other)151 void Profile<Abc>::Insert(size_t idx, const Profile<Abc>& other) {
152 size_t n = MIN(other.length(), length() - idx);
153 memcpy(v[idx], other[0], n * Abc::kSizeAny * sizeof(double));
154 }
155
156 template <class Abc>
~Profile()157 Profile<Abc>::~Profile() {
158 if (v != NULL) {
159 delete[] v[0];
160 delete[] v;
161 }
162 }
163
164 // Assigns given constant value or default to all entries in matrix
165 template <class Abc>
Assign(Profile<Abc> & p,double val)166 inline void Assign(Profile<Abc>& p, double val) {
167 for (size_t i = 0; i < p.length(); ++i)
168 for (size_t j = 0; j < Abc::kSizeAny; ++j)
169 p[i][j] = val;
170 }
171
172 // Normalizes all profile columns to fixed value. Iff 'incl_any' is true,
173 // normalization also includes values for ANY letter
174 template <class Abc>
Normalize(Profile<Abc> & p,double val,bool incl_any)175 inline void Normalize(Profile<Abc>& p, double val, bool incl_any) {
176 const size_t abc_size = incl_any ? Abc::kSizeAny : Abc::kSize;
177 for (size_t i = 0; i < p.length(); ++i) {
178 double sum = 0.0;
179 for (size_t a = 0; a < abc_size; ++a) sum += p[i][a];
180 if (fabs(val - sum) > kNormalize && sum != 0.0) {
181 double fac = val / sum;
182 for (size_t a = 0; a < abc_size; ++a) p[i][a] *= fac;
183 }
184 }
185 }
186
187 // Normalizes all profile columns to corresponding value in vector 'norm'.
188 // Iff 'incl_any' is true, normalization also includes values for ANY letter.
189 template <class Abc>
Normalize(Profile<Abc> & p,const Vector<double> & norm,bool incl_any)190 inline void Normalize(Profile<Abc>& p, const Vector<double>& norm, bool incl_any) {
191 const size_t abc_size = incl_any ? Abc::kSizeAny : Abc::kSize;
192 for (size_t i = 0; i < p.length(); ++i) {
193 double sum = 0.0;
194 for (size_t a = 0; a < abc_size; ++a) sum += p[i][a];
195 if (fabs(norm[i] - sum) > kNormalize && sum != 0.0) {
196 double fac = norm[i] / sum;
197 for (size_t a = 0; a < abc_size; ++a) p[i][a] *= fac;
198 }
199 }
200 }
201
202 // Prints profile in human-readable format for debugging.
203 template<class Abc>
204 std::ostream& operator<< (std::ostream& out, const Profile<Abc>& p) {
205 out << "Profile" << std::endl;
206 for (size_t a = 0; a < Abc::kSizeAny; ++a)
207 out << "\t" << Abc::kIntToChar[a];
208 out << std::endl;
209 for (size_t i = 0; i < p.length(); ++i) {
210 out << i+1;
211 for (size_t a = 0; a < Abc::kSizeAny; ++a)
212 out << strprintf("\t%6.4f", p[i][a]);
213 out << std::endl;
214 }
215 return out;
216 }
217
218 // Calculates the entropy of the given profile using logarithm base 2.
219 template <class Abc>
Entropy(const Profile<Abc> & p)220 inline double Entropy(const Profile<Abc>& p) {
221 double rv = 0.0;
222 for (size_t i = 0; i < p.length(); ++i) {
223 for (size_t a = 0; a < Abc::kSize; ++a)
224 if (p[i][a] > FLT_MIN) rv -= p[i][a] * log2(p[i][a]);
225 }
226 return rv;
227 }
228
229 // Calculates the Neff in the given profile.
230 template <class Abc>
Neff(const Profile<Abc> & p)231 inline double Neff(const Profile<Abc>& p) {
232 return p.length() > 0 ? pow(2, Entropy(p) / p.length()) : 0;
233 }
234
235
236 } // namespace cs
237
238 #endif // CS_PROFILE_INL_H_
239