22 #include <boost/format.hpp>
40 extern "C" void glassofast_(
int *,
g_fp_t *,
g_fp_t *,
g_fp_t *,
int *,
int *,
int *,
g_fp_t *,
g_fp_t *,
g_fp_t *,
g_fp_t *,
int *,
int *);
44 void spotrf_(
const char* UPLO,
const int* N,
float* A,
const int* LDA,
int* INFO);
47 void sgetrf_(
int* M,
int* N,
float* A,
int* LDA,
int* IPIV,
int* INFO);
50 void sgetri_(
int* N,
float* A,
int* lda,
int* IPIV,
float* WORK,
int* LWORK,
int* INFO);
61 template<
typename _Tp>
66 ct_vector(uint16_t __alilen = 0) :
std::vector<_Tp>(__alilen*__alilen), alilen(__alilen) {}
71 return this->_M_impl._M_start[ __i*alilen + __j ];
74 inline const _Tp&
operator()(uint16_t __i, uint16_t __j)
const
76 return this->_M_impl._M_start[ __i*alilen + __j ];
80 template<
typename _Tp>
86 af_vector(uint16_t __alilen, uint8_t __q, _Tp __v) :
std::vector<_Tp>(__alilen*__q, __v), alilen(__alilen), q(__q) {}
90 return this->_M_impl._M_start[ __i*q + __ai ];
92 inline const _Tp&
operator()(uint16_t __i, uint8_t __ai)
const
94 return this->_M_impl._M_start[ + __i*q + __ai ];
98 class pf_vector :
public std::vector<predictor::pairfreq_t> {
105 pf_vector() :
std::vector<_Tp>(), alilen(0), q(0), d3(0), d2(0), d1(0) {}
106 pf_vector(uint16_t __alilen, uint8_t __q, _Tp __v) :
std::vector<_Tp>(__alilen*__alilen*__q*__q, __v), alilen(__alilen), q(__q), d3(q), d2(q*d3), d1(alilen*d2) {}
108 inline _Tp&
operator()(uint16_t __i, uint16_t __j, uint8_t __ai, uint8_t __aj)
110 return this->_M_impl._M_start[ (__i)*d1 + (__j)*d2 + (__ai)*d3 + (__aj) ];
112 inline const _Tp&
operator()(uint16_t __i, uint16_t __j, uint8_t __ai, uint8_t __aj)
const
114 return this->_M_impl._M_start[ (__i)*d1 + (__j)*d2 + (__ai)*d3 + (__aj) ];
118 template<
typename _Tp>
126 cov_vector(uint16_t __alilen, uint8_t __q) :
std::vector<_Tp>(__alilen*__q*__alilen*__q),
127 alilen(__alilen), q(__q), d1(alilen*q) {}
128 cov_vector(uint16_t __alilen, uint8_t __q, _Tp __v) :
std::vector<_Tp>(__alilen*__q*__alilen*__q, __v),
129 alilen(__alilen), q(__q), d1(alilen*q) {}
134 return this->_M_impl._M_start[ __row*d1 + __col ];
136 inline const _Tp&
operator()(
size_t __row,
size_t __col)
const
138 return this->_M_impl._M_start[ __row*d1 + __col ];
141 inline _Tp&
operator()(uint16_t __i, uint16_t __j, uint8_t __ai, uint8_t __aj)
143 return this->_M_impl._M_start[ (__i*q+__ai)*d1 + __j*q+__aj ];
145 inline const _Tp&
operator()(uint16_t __i, uint16_t __j, uint8_t __ai, uint8_t __aj)
const
147 return this->_M_impl._M_start[ (__i*q+__ai)*d1 + __j*q+__aj ];
154 virtual double operator()(
const uint16_t i,
const uint16_t j)
const = 0;
158 template<
typename _Tp,
typename _Tri =
size_t,
typename _Tci =
size_t>
164 d2matrix(_Tri __rows, _Tci __cols) :
std::vector<_Tp>(__rows*__cols), rows(__rows), cols(__cols)
166 memset(this->_M_impl._M_start, 0, rows*cols*
sizeof(_Tp));
170 return this->_M_impl._M_start[ __r*cols + __c ];
174 return this->_M_impl._M_start[ __r*cols + __c ];
180 const uint16_t __alilen,
const string &__key,
const _rawscore_calc_t &__fo );
182 static vector<contact_t>
183 _apc(
const ct_vector<fp_t>& __raw_ctscore,
const vector<double>& __apc_bg,
const double __apc_mean,
const uint16_t __mincontsep,
bool __filt );
185 static inline vector<contact_t>
188 static inline uint32_t
192 static inline __m128i
_mm_setzero_si128(){ __m128i m;
long long *mp = (
long long *)&m; mp[1] = mp[0] = 0;
return m; }
199 static void _brk(sigval_t __sigval)
201 *
static_cast<int*
>(__sigval.sival_ptr) = 1;
209 _glasso_timer(
int *__sival_ptr, time_t __tv_sec,
bool __dbg =
false) : dbg(__dbg)
212 sev.sigev_notify = SIGEV_THREAD;
213 sev.sigev_notify_function = _brk;
214 sev.sigev_value.sival_ptr = __sival_ptr;
215 sev.sigev_notify_attributes = NULL;
216 if (timer_create(CLOCK_PROCESS_CPUTIME_ID, &sev, &_timerid) == -1)
218 if (dbg) cerr <<
"timer_create(CLOCK_PROCESS_CPUTIME_ID, ...): " << strerror(errno) <<
" at " << __FILE__ <<
":" << __LINE__ <<
"\n";
220 if (timer_create(CLOCK_MONOTONIC, &sev, &_timerid) == -1)
throw std::runtime_error(bo::str(bo::format(
"%s at %s:%d") % strerror(errno) % __FILE__ % __LINE__ ));
222 struct itimerspec its;
223 its.it_value.tv_sec = __tv_sec;
224 its.it_value.tv_nsec = 0;
225 its.it_interval.tv_sec = 0;
226 its.it_interval.tv_nsec = 0;
227 if (timer_settime(_timerid, 0, &its, NULL) == -1)
throw std::runtime_error(bo::str(bo::format(
"%s at %s:%d") % strerror(errno) % __FILE__ % __LINE__ ));
232 if (timer_delete(_timerid) == -1)
throw std::runtime_error(strerror(errno));
236 void predictor::get_seq_weights(
238 const ali_t& __ali,
double __cp,
239 bool __veczw,
int __num_threads
242 if (__ali.seqcnt < 2)
throw alilen_error(bo::str(bo::format(
"alignment size (%d) < 2") % __ali.seqcnt));
244 if (__ali.alilenpad > 4080)
throw alilen_error(bo::str(bo::format(
"alignment (%d) longer than 4080") % __ali.alilenpad));
245 if (__cp < 0 || __cp > 1)
throw range_error(bo::str(bo::format(
"clustering threshold is out of range [0-1] %g") % __cp ));
246 if (__num_threads < 0)
throw range_error(bo::str(bo::format(
"number of threads is out of range [0-) %d") % __num_threads ));
249 int num_threads = __num_threads ? __num_threads : omp_get_max_threads();
250 if(dbg) cerr <<
"will use " << num_threads <<
" OMP threads\n";
253 if(dbg) cerr <<
"no OMP thread support\n";
260 time_t t0 = time(NULL);
262 uint16_t mismthresh = __ali.alilen - ceil( __ali.alilen * __cp );
263 int matchthresh = __ali.alilenpad - mismthresh;
266 simcnt_t simcnt(__ali.seqcnt, num_threads);
280 #pragma omp parallel for num_threads(num_threads), schedule(dynamic)
282 for(uint32_t i=0; i<__ali.seqcnt-1; ++i)
286 && omp_get_thread_num() == 0
289 { time_t t1 = time(NULL);
if( t1 != t0 ) { t0 = t1; fprintf(stderr,
"%s%d/%d",
"[1G[K", i, __ali.seqcnt); } }
291 const uint8_t *a0 = &__ali(i,0);
293 for(uint32_t j=i+1; j<__ali.seqcnt; ++j)
295 const uint8_t *a = a0, *b = &__ali(j,0);
296 uint16_t matc128 = 0;
298 for(uint16_t k = 0, k_end = __ali.alilenpad; k < k_end; ++k)
299 matc128 += ( a[k] == b[k] );
301 if( matc128 >= matchthresh )
306 int thread_num = omp_get_thread_num();
308 ++simcnt(i,thread_num);
309 ++simcnt(j,thread_num);
319 uint32_t wchunk = cache_holds_nseq / 2;
321 if(dbg) cerr <<
"SSE2 veczweight, wchunk = "<< wchunk <<
"\n";
324 __m128i m1 = _mm_cmpeq_epi8( m0, m0 );
326 uint32_t chunk_e = (__ali.seqcnt-2)/wchunk + 1;
328 #pragma omp parallel for num_threads(num_threads), schedule(dynamic)
330 for(uint32_t chunki=0; chunki<chunk_e; ++chunki)
336 && omp_get_thread_num() == 0
339 { time_t t1 = time(NULL);
if( t1 != t0 ) { t0 = t1; fprintf(stderr,
"%s%d,%d/%d",
"[1G[K", chunki, chunki, chunk_e); } }
341 uint32_t i_b = chunki*wchunk;
342 uint32_t i_e = min( i_b+wchunk, __ali.seqcnt-1 );
343 uint32_t j_e = min( i_b+wchunk, __ali.seqcnt );
344 for(uint32_t i=i_b; i<i_e; ++i)
346 __m128i *a0 = (__m128i*)&__ali(i,0);
352 #define _count_matches_sse2_helper(__omp_get_thread_num) {\
353 __m128i *a = a0, *b = (__m128i*)&__ali(j,0);\
354 __m128i matc128 = m1;\
356 for(uint16_t k = 0, k_end = __ali.alilen16; k < k_end; ++k) \
358 __m128i m = _mm_cmpeq_epi8( a[k], b[k] );\
359 matc128 = _mm_add_epi8( m, matc128 );\
362 matc128 = _mm_sad_epu8( m1, matc128 );\
363 matc128 = _mm_add_epi32( matc128,\
364 _mm_srli_si128( matc128, 8 ) );\
366 if( _mm_cvtsi128_si32(matc128) >= matchthresh )\
368 int thread_num = (__omp_get_thread_num);\
369 ++simcnt(i,thread_num);\
370 ++simcnt(j,thread_num);\
373 #define _count_matches_sse2() _count_matches_sse2_helper(omp_get_thread_num())
375 #define _count_matches_sse2() _count_matches_sse2_helper(0)
377 for(uint32_t j = i+1; j<j_e; ++j)
378 _count_matches_sse2();
383 for(uint32_t chunkj=chunki+1; chunkj<chunk_e; ++chunkj)
387 && omp_get_thread_num() == 0
390 { time_t t1 = time(NULL);
if( t1 != t0 ) { t0 = t1; fprintf(stderr,
"%s%d,%d/%d",
"[1G[K", chunki, chunkj, chunk_e); } }
392 uint32_t i_b = chunki*wchunk;
393 uint32_t i_e = min( i_b+wchunk, __ali.seqcnt-1 );
394 uint32_t j_b = chunkj*wchunk;
395 uint32_t j_e = min( j_b+wchunk, __ali.seqcnt );
396 for(uint32_t i=i_b; i<i_e; ++i)
398 __m128i *a0 = (__m128i*)&__ali(i,0);
400 for(uint32_t j=j_b; j<j_e; ++j)
401 _count_matches_sse2();
404 #undef _count_matches_sse2
405 #undef _count_matches_sse2_helper
408 #else // should never happen unless bug
409 else throw std::runtime_error(bo::str(bo::format(
"ooops, there is a bug in %s around %s:%d") %
PACKAGE % __FILE__ % __LINE__ ));
415 for(uint32_t i = 0; i<__ali.seqcnt; ++i)
419 uint32_t sci = simcnt(i,thread_num);
422 for(
int thread_num = 0; thread_num < num_threads; ++thread_num)
423 sci += simcnt(i,thread_num);
425 __wtot += (__aliw[i] = 1.0 / (1 + sci));
428 if(dbg) fprintf(stderr,
"total weight (variation) of alignment = %16.15g\n", __wtot );
435 double __dens,
double __gapth, uint16_t __mincontsep,
436 double __pseudocnt,
double __pscnt_weight,
bool __estimate_ivcov,
double __shrink_lambda,
437 bool __cov20,
bool __apply_gapth,
double __rho,
438 int __num_threads, time_t __icme_timeout,
time_res_t *__timing
441 if (__ali.seqcnt < 2)
throw alilen_error(bo::str(bo::format(
"alignment size (%d) < 2") % __ali.seqcnt));
442 if (__ali.seqcnt != __aliw.size())
throw alilen_error(bo::str(bo::format(
"alignment size (%d) != alignment weight vector size (%d)") % __ali.seqcnt % __aliw.size()));
444 if (__ali.alilen < 4)
throw alilen_error(bo::str(bo::format(
"alignment (%d) shorter than 4") % __ali.alilen));
445 if (__ali.alilen < __mincontsep+1)
throw alilen_error(bo::str(bo::format(
"alignment (%d) shorter than __mincontsep (%d) + 1") % __ali.alilen % __mincontsep));
448 if(__dens < 0 || __dens > 1)
throw range_error(bo::str(bo::format(
"contact density is out of range [0-1] %g") % __dens ));
449 if(__gapth < 0 || __gapth > 1)
throw range_error(bo::str(bo::format(
"gap threshold is out of range [0-1] %g") % __gapth ));
450 if(__mincontsep < 1)
throw range_error(bo::str(bo::format(
"minimum contact separation is out of range [1-) %d") % __mincontsep ));
451 if(__num_threads < 0)
throw range_error(bo::str(bo::format(
"number of threads is out of range [0-) %d") % __num_threads ));
452 if(__pseudocnt < 0)
throw range_error(bo::str(bo::format(
"pseudo count is out of range [0-) %g") % __pseudocnt ));
453 if(__pscnt_weight < 0 || __pscnt_weight > 1)
throw range_error(bo::str(bo::format(
"pseudo count weight is out of range [0-1] %g") % __pscnt_weight ));
454 if(__shrink_lambda < 0 || __shrink_lambda > 1)
throw range_error(bo::str(bo::format(
"shrinkage lambda is out of range [0-1] %g") % __shrink_lambda ));
457 int num_threads = __num_threads ? __num_threads : omp_get_max_threads();
458 if(dbg) cerr <<
"will use " << num_threads <<
" OMP threads\n";
461 if(dbg) cerr <<
"no OMP thread support\n";
464 if(__timing) (*__timing)[
"num_threads"] = num_threads;
468 const double ps_wtot = __pseudocnt * 21 + __wtot;
471 freq_t ps_aliw[__ali.seqcnt]; memcpy(ps_aliw, __aliw.data(), __ali.seqcnt*
sizeof(
freq_t));
472 cblas_dscal(__ali.seqcnt, 1/ps_wtot, ps_aliw, 1);
474 for(uint32_t k = 0; k < __ali.seqcnt; ++k)
476 for(uint16_t j = 0; j < __ali.alilen; ++j)
478 uint8_t aa = __ali(k,j);
479 if(aa < 21) aafreq(j,aa) += ps_aliw[k];
483 vector<uint16_t> gap_cols;
485 for(uint16_t i = 0; i < __ali.alilen; ++i)
487 if( aafreq(i,
aamap_gapidx) > __gapth ) gap_cols.push_back(i);
490 if(dbg) cerr <<
"calculated column aa frequencies, gap cols = " << gap_cols.size() <<
"\n";
493 pf_vector pairfreq( __ali.alilen, 21, __pseudocnt/21.0 / ps_wtot );
496 gettimeofday(&t_before, NULL);
498 for(uint32_t k = 0; k<__ali.seqcnt; ++k)
500 freq_t ps_aliwk = ps_aliw[k];
501 const uint8_t *alik = &__ali(k,0);
504 #pragma omp parallel for num_threads(num_threads), schedule(static)
506 for(uint16_t i = 0; i < (__ali.alilen-2)/2+1; ++i)
508 #define _fill_pfi(__i) {\
509 pairfreq_t *pfi = &pairfreq(__i,0,0,0);\
510 for(uint16_t j = __i+1; j < __ali.alilen; ++j)\
512 uint8_t ai = alik[__i];\
513 uint8_t aj = alik[j];\
514 if(ai < 21 && aj < 21)\
515 pfi[ j*pairfreq.d2 + ai*pairfreq.d3 + aj ] += ps_aliwk;\
520 uint16_t ri = __ali.alilen-2-i;
527 for(uint16_t i = 0; i < __ali.alilen-1; ++i)
528 for(uint16_t j = i+1; j < __ali.alilen; ++j)
530 pairfreq_t *pfij = &pairfreq(i,j,0,0), *pfji = &pairfreq(j,i,0,0);
531 for(uint8_t ai = 0; ai < 21; ++ai)
532 for(uint8_t aj = 0; aj < 21; ++aj)
533 pfji[ aj*pairfreq.
d3 + ai ] = pfij[ ai*pairfreq.
d3 + aj ];
536 for(uint16_t i = 0; i < __ali.alilen; ++i)
539 memset(pairfreq.data() + i*pairfreq.
d1 + i*pairfreq.
d2, 0, 21*21*
sizeof(
pairfreq_t));
540 for(uint8_t aa = 0; aa < 21; ++aa)
541 pairfreq(i,i,aa,aa) = aafreq(i,aa);
544 { timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_before, &t_diff);
545 if(__timing) (*__timing)[
"pairfreq"] = t_diff.tv_sec+t_diff.tv_usec/1e6;
if(dbg) fprintf(stderr,
"calculated pair frequency table in %g secs\n", t_diff.tv_sec+t_diff.tv_usec/1e6 ); }
548 map<string, ct_vector<fp_t> > raw_ctscore;
549 map<string, vector<double> > apc_bg;
550 map<string, double> apc_mean;
556 _rawscore_MI_true_t(
const af_vector<freq_t> &__Pi_true,
const pf_vector &__Pij_true ) : _Pi_true(__Pi_true), _Pij_true(__Pij_true) {}
557 virtual double operator()(
const uint16_t i,
const uint16_t j)
const
560 for(uint8_t ai = 0; ai < 21; ++ai)
561 for(uint8_t aj = 0; aj < 21; ++aj)
562 if(_Pij_true(i,j,ai,aj) > 0)
563 raw_ij += _Pij_true(i,j,ai,aj) * log( _Pij_true(i,j,ai,aj) / _Pi_true(i,ai) / _Pi_true(j,aj) );
566 } _rawscore_MI_true(aafreq, pairfreq);
567 _raw_score_matrix(raw_ctscore, apc_bg, apc_mean, __ali.alilen,
"MI", _rawscore_MI_true );
568 if(dbg) fprintf(stderr,
"collected apc_mean[MI] = %16.15g\n", apc_mean[
"MI"]);
573 if( __pscnt_weight != 0.0 )
575 for(
size_t i = 0; i < pairfreq.size(); ++i) pairfreq[i] = (1-__pscnt_weight)*pairfreq[i] + __pscnt_weight/21/21;
576 for(
size_t i = 0; i < aafreq.size(); ++i) aafreq[i] = (1-__pscnt_weight)*aafreq[i] + __pscnt_weight/21;
578 for(uint16_t i = 0; i < __ali.alilen; ++i)
579 for(uint8_t ai = 0; ai < 21; ++ai)
580 for(uint8_t aj = 0; aj < 21; ++aj)
581 pairfreq(i,i,ai,aj) = ai==aj ? (1-__pscnt_weight)*Pij_true(i,i,ai,aj) + __pscnt_weight/21 : 0;
590 for(uint16_t i = 0; i < __ali.alilen; ++i)
591 for(uint16_t j = 0; j < __ali.alilen; ++j)
592 for(uint8_t ai = 0; ai < 21; ++ai)
594 if(j==0) aafs += aafreq(i,ai);
595 for(uint8_t aj = 0; aj < 21; ++aj)
596 pfs += pairfreq(i,j,ai,aj);
598 fprintf(stderr,
"aa freq sum (cell) = %.15g, pairfreq sum (cell) = %.15g\n", aafs/__ali.alilen, pfs/__ali.alilen/__ali.alilen);
601 #ifdef RETURN_AFTER_PAIRFREQ
602 cerr <<
"returning because of RETURN_AFTER_PAIRFREQ " << __FILE__ <<
":" << __LINE__ <<
"\n";
603 return map<string, vector<contact_t> >();
611 __apply_gapth ? __ali.alilen - gap_cols.size() : __ali.alilen,
616 uint16_t i = 0, covi = 0;
617 vector<uint16_t>::const_iterator gc_bi = ( __apply_gapth ? gap_cols.begin() : gap_cols.end() ), gc_ei = gap_cols.end();
618 for(; i < __ali.alilen; ++i)
620 if(gc_bi!=gc_ei && *gc_bi == i){ ++gc_bi;
continue; }
622 uint16_t j = 0, covj = 0;
623 vector<uint16_t>::const_iterator gc_bj = ( __apply_gapth ? gap_cols.begin() : gap_cols.end() ), gc_ej = gap_cols.end();
624 for(; j < __ali.alilen; ++j)
626 if(gc_bj!=gc_ej && *gc_bj == j){ ++gc_bj;
continue; }
628 for(uint8_t ai = 0; ai < covm.q; ++ai)
629 for(uint8_t aj = 0; aj < covm.q; ++aj)
630 covm(covi,covj,ai,aj) = pairfreq(i,j,ai,aj) - aafreq(i,ai) * aafreq(j,aj);
635 {
pf_vector _empty; pairfreq.swap(_empty); }
639 for(uint16_t i = 0; i < covm.alilen; ++i)
640 for(uint8_t ai = 0; ai < covm.q; ++ai)
641 for(uint8_t aj = 0; aj < covm.q; ++aj)
642 if(ai != aj) covm(i,i,ai,aj) = 0;
643 if(dbg) cerr<<
"formed covariance matrix ("<<covm.
alilen<<
"/"<<__ali.alilen<<
","<<gap_cols.size()<<
")\n";
648 if(__shrink_lambda>0)
650 gettimeofday(&t_before, NULL);
653 for(
size_t i = 0; i < covm.d1; ++i) dmean += covm(i, i);
655 const double& lambda = __shrink_lambda;
659 if(dbg) cerr <<
"cov matrix Cholesky...";
663 vector<float> tempmat = covm;
665 spotrf_(
"L",&N,tempmat.data(),&N,&cholres);
666 if(cholres < 0)
throw std::runtime_error(bo::str(bo::format(
"illegal value for argument %d") % -cholres ));
668 if(cholres==0)
break;
670 if(dbg) cerr <<
" not positive definite";
672 for(uint16_t i = 0; i < covm.alilen; ++i)
673 for(uint16_t j = 0; j < covm.alilen; ++j)
674 for(uint8_t ai = 0; ai < covm.q; ++ai)
675 for(uint8_t aj = 0; aj < covm.q; ++aj)
677 covm(i,j,ai,aj) *= (1.0 - lambda);
679 covm(i,j,ai,aj) = lambda*dmean + (1.0-lambda)*covm(i,j,ai,aj);
683 double dmean_dev = 0;
684 for(
size_t i = 0; i < covm.d1; ++i) dmean_dev += pow(dmean - covm(i,i), 2);
685 dmean_dev /= covm.d1; dmean_dev = sqrt(dmean_dev);
686 fprintf(stderr,
", dmean dev: %g\n", dmean_dev );
689 { timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_before, &t_diff);
690 if(__timing) (*__timing)[
"shrink"] = t_diff.tv_sec+t_diff.tv_usec/1e6;
if(dbg) fprintf(stderr,
"\nshrank covariance matrix in %g secs\n", t_diff.tv_sec+t_diff.tv_usec/1e6 ); }
693 #ifdef LIBFREEC_CHKPT
695 FILE *f = fopen(
"covm.dat",
"w");
697 if( fwrite(&covd,
sizeof(covd), 1, f) != 1 ) { cerr <<
"error writing 'covm.dat': " << strerror(errno) <<
"\n"; fclose(f); unlink(
"covm.dat");
break; }
698 if( fwrite(&galilen,
sizeof(galilen), 1, f) != 1 ) { cerr <<
"error writing 'covm.dat': " << strerror(errno) <<
"\n"; fclose(f); unlink(
"covm.dat");
break; }
699 if( fwrite( covm.data(),
sizeof(*covm.data()), covm.size(), f ) != covm.size() )
701 cerr <<
"error writing 'covm.dat': " << strerror(errno) <<
"\n";
702 fclose(f); unlink(
"covm.dat");
705 cerr <<
"wrote 'covm.dat'\n";
716 gettimeofday(&t_before, NULL);
719 fp_t density = -1, prev_dens = -1;
722 if( rho < 0 ) rho = std::max(0.001, 1.0 / __wtot);
723 else if(density >= 0)
725 if( density == 0 ) rhofac = 0.5;
728 if( rhofac != 0 && prev_dens > 0 && density > 0 )
729 rhofac = pow( rhofac, log( __dens / density ) / log( density / prev_dens ) );
731 rhofac = __dens > density ? 0.9 : 1.1;
737 if(dbg) fprintf(stderr,
"will try rho %g, %s solution\n", rho, (g_ia ?
"approximate" :
"exact"));
739 if(rho <= 0 || rho >= 1)
throw std::runtime_error(bo::str(bo::format(
"regularization parameter rho is out of expected range (0-1): %g") % rho ));
743 for(uint16_t i = 0; i < rho_m.
alilen; ++i)
744 for(uint16_t j = 0; j < rho_m.
alilen; ++j)
745 for(uint8_t ai = 0; ai < rho_m.
q; ++ai)
746 for(uint8_t aj = 0; aj < rho_m.
q; ++aj)
749 (i == j && ai != aj )
751 rho_m(i,j,ai,aj) = 1e9;
753 #ifdef LIBFREEC_CHKPT
755 FILE *f = fopen(
"rho_m.dat",
"w");
757 if( fwrite( rho_m.data(),
sizeof(*rho_m.data()), rho_m.size(), f ) != rho_m.size() )
759 cerr <<
"error writing 'rho_m.dat': " << strerror(errno) <<
"\n";
760 fclose(f); unlink(
"rho_m.dat");
763 cerr <<
"wrote 'rho_m.dat'\n";
771 int g_n = covm.d1, g_is = 0, g_msg = dbg, g_maxit = 1e5, g_jerr, g_brk = 0;
774 vector<g_fp_t> ww( g_n*g_n );
776 timeval t_before;
if(dbg) gettimeofday(&t_before, NULL);
784 vector<float> g_Wd(g_n), g_WXj(g_n);
785 glassofast_(&g_n, covm.data(), rho_m.data(), &g_thr, &g_maxit, &g_msg, &g_is, wwi.data(), ww.data(), g_Wd.data(), g_WXj.data(), &g_jerr, &g_brk);
790 throw icme_timeout_error(bo::str(bo::format(
"inverse covariance matrix estimation exceeded time limit (%d sec)") % __icme_timeout ));
break;
792 throw std::runtime_error(
"an error occurred in glassofast_");
798 if(dbg){ timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_before, &t_diff);
799 fprintf(stderr,
"%d iterations of glassofast took %g secs\n", g_niter, t_diff.tv_sec+t_diff.tv_usec/1e6 ); }
802 if( __dens == 0 )
break;
807 for(
size_t i = 0; i < wwi.
d1; ++i)
808 for(
size_t j = i+1; j < wwi.
d1; ++j)
809 if( wwi(i, j) != 0 ) ++nz;
811 density = (
fp_t)nz / ( wwi.
d1 * (wwi.
d1-1) / 2 );
813 if(dbg) fprintf(stderr,
"density = %7g, target sp = %g, |density-target sp|/target sp = %g\n", density, __dens, fabs(density - __dens)/__dens);
815 while( __dens == 0 || fabs(__dens - density)/__dens > 0.01 );
823 int N = wwi.
d1, INFO;
826 timeval t_before;
if(dbg) gettimeofday(&t_before, NULL);
828 sgetrf_(&N,&N,wwi.data(),&N,IPIV,&INFO);
829 if(INFO)
throw std::runtime_error(bo::str(bo::format(
"LU factorization error %d") % INFO ));
832 if(dbg){ timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_before, &t_diff);
833 fprintf(stderr,
"LU factorization took %g secs, ", t_diff.tv_sec+t_diff.tv_usec/1e6 ); }
835 int LWORK = -1; vector<float> WORK(1);
836 sgetri_(&N,wwi.data(),&N,IPIV,WORK.data(),&LWORK,&INFO);
837 LWORK = N*N < WORK[0] ? N*N : WORK[0];
840 sgetri_(&N,wwi.data(),&N,IPIV,WORK.data(),&LWORK,&INFO);
841 if(INFO)
throw std::runtime_error(bo::str(bo::format(
"matrix inversion error %d") % INFO ));
843 if(dbg){ timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_before, &t_diff);
844 fprintf(stderr,
"inverted matrix (incl LUf) in %g secs\n", t_diff.tv_sec+t_diff.tv_usec/1e6 ); }
847 if(__timing){ timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_before, &t_diff);
848 (*__timing)[
"inv"] = t_diff.tv_sec+t_diff.tv_usec/1e6; }
854 for(
size_t i = 0; i < wwi.
d1; ++i)
855 for(
size_t j = i+1; j < wwi.
d1; ++j)
856 if( wwi(i, j) != 0 ){ ++nz; checksum += wwi(i, j); }
858 fp_t density = (
fp_t)nz / wwi.
d1 / (wwi.
d1-1) * 2;
859 fprintf(stderr,
"density of inverse covariance matrix = %g (cksum %.7g)\n", density, checksum);
867 uint16_t i = 0, covi = 0;
868 vector<uint16_t>::const_iterator gc_bi = ( __apply_gapth ? gap_cols.begin() : gap_cols.end() ), gc_ei = gap_cols.end();
869 for(; i < __ali.alilen; ++i)
871 if(gc_bi!=gc_ei && *gc_bi == i){ ++gc_bi;
continue; }
873 uint16_t j = 0, covj = 0;
874 vector<uint16_t>::const_iterator gc_bj = ( __apply_gapth ? gap_cols.begin() : gap_cols.end() ), gc_ej = gap_cols.end();
875 for(; j < __ali.alilen; ++j)
877 if(gc_bj!=gc_ej && *gc_bj == j){ ++gc_bj;
continue; }
879 for(uint8_t ai = 0; ai < wwiwg.
q; ++ai)
880 for(uint8_t aj = 0; aj < wwiwg.
q; ++aj)
881 wwiwg(i,j,ai,aj) = wwi(covi,covj,ai,aj);
886 if(dbg) cerr<<
"went back to gapped ("<<wwiwg.
alilen<<
") wwi matrix\n";
899 virtual double operator()(
const uint16_t i,
const uint16_t j)
const
902 for(uint8_t ai = 0; ai < 20; ++ai)
903 for(uint8_t aj = 0; aj < 20; ++aj)
904 raw_ij += fabs(_wwiwg(i,j,ai,aj));
907 } _rawscore_l1norm(wwiwg);
908 _raw_score_matrix(raw_ctscore, apc_bg, apc_mean, __ali.alilen,
"l1norm", _rawscore_l1norm );
909 if(dbg) fprintf(stderr,
"collected apc_mean[l1norm] = %16.15g\n", apc_mean[
"l1norm"]);
919 virtual double operator()(
const uint16_t i,
const uint16_t j)
const
921 vector<float> W_mf( 21*21, 0.0 );
922 vector<double> rmean_W_mf( 21, 0.0 );
923 vector<double> cmean_W_mf( 21, 0.0 );
924 double mean2_W_mf = 0;
925 for(uint8_t ai = 0; ai < 20; ++ai)
926 for(uint8_t aj = 0; aj < 20; ++aj)
928 const double cell = ( W_mf[ ai*21+aj ] = -_wwiwg(i,j,ai,aj) );
930 rmean_W_mf[ai] += cell/21;
931 cmean_W_mf[aj] += cell/21;
933 mean2_W_mf /= W_mf.size();
936 double sum_sum_W_mf_zsg2 = 0.0;
937 for(uint8_t ai = 0; ai < 21; ++ai)
938 for(uint8_t aj = 0; aj < 21; ++aj)
940 double W_mf_zsg = (double)W_mf[ ai*21+aj ] - rmean_W_mf[ai] - cmean_W_mf[aj] + mean2_W_mf;
941 sum_sum_W_mf_zsg2 += W_mf_zsg*W_mf_zsg;
945 double raw_ij = sqrt(sum_sum_W_mf_zsg2);
949 } _rawscore_fro(wwiwg);
950 _raw_score_matrix(raw_ctscore, apc_bg, apc_mean, __ali.alilen,
"fro", _rawscore_fro );
951 if(dbg) fprintf(stderr,
"collected apc_mean[fro] = %16.15g\n", apc_mean[
"fro"]);
956 map<string, vector<contact_t> > res;
958 res[
"l1norm"] =
_apc( raw_ctscore[
"l1norm"], apc_bg[
"l1norm"], apc_mean[
"l1norm"], __mincontsep,
true );
960 res[
"MI"] =
_raw_as_is( raw_ctscore[
"MI"], __mincontsep );
962 res[
"fro"] =
_apc( raw_ctscore[
"fro"], apc_bg[
"fro"], apc_mean[
"fro"], __mincontsep,
false );
971 predictor::run(
const ali_t& __ali,
double __clustpc,
972 double __density,
double __gapth, uint16_t __mincontsep,
973 double __pseudocnt,
double __pscnt_weight,
bool __estimate_ivcov,
double __shrink_lambda,
974 bool __cov20,
bool __apply_gapth,
double __rho,
975 bool __veczw,
int __num_threads, time_t __icme_timeout,
time_res_t *__timing
978 timeval t_top; gettimeofday(&t_top, NULL);
979 timeval t_before; gettimeofday(&t_before, NULL);
982 get_seq_weights(aliw, wtot, __ali, __clustpc, __veczw, __num_threads);
984 { timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_before, &t_diff);
985 if(__timing) (*__timing)[
"seqw"] = t_diff.tv_sec+t_diff.tv_usec/1e6;
986 if(dbg) fprintf(stderr,
"\nseq weight loop for %d seqs took %g secs\n", __ali.seqcnt, t_diff.tv_sec+t_diff.tv_usec/1e6 ); }
989 __density, __gapth, __mincontsep,
990 __pseudocnt, __pscnt_weight, __estimate_ivcov, __shrink_lambda,
991 __cov20, __apply_gapth, __rho,
992 __num_threads, __icme_timeout, __timing);
994 { timeval t_now, t_diff; gettimeofday(&t_now, NULL); timersub(&t_now, &t_top, &t_diff);
995 if(__timing) (*__timing)[
"all"] = t_diff.tv_sec+t_diff.tv_usec/1e6;
if(dbg) fprintf(stderr,
"all done in %g secs\n", t_diff.tv_sec+t_diff.tv_usec/1e6); }
1003 if( __al.size() != alilen )
throw alilen_error(bo::str(bo::format(
"alignment length mismatch, expected %d, got %d") % alilen % __al.size() ));
1006 uint8_t *alptr =
reinterpret_cast<uint8_t*
>(data()) + (seqcnt-1) * alilenpad;
1007 memcpy(alptr, __al.data(), __al.size());
1009 if( capacity() == size() )
1010 reserve( size() + 1024*alilen16 );
1024 const uint16_t __alilen,
const string &__key,
const _rawscore_calc_t &__fo )
1027 vector<double>& apc_bg = __apc_bg[__key] = vector<double>( __alilen );
1028 double& mean = __apc_mean[__key] = 0;
1030 for(uint16_t i = 0; i < __alilen; ++i)
1031 for(uint16_t j = i+1; j < __alilen; ++j)
1033 double raw_ij = __fo(i, j);
1035 raw_cts[ i*__alilen + j ] = raw_cts[ j*__alilen + i ] = raw_ij;
1038 double bg_ij = raw_ij / (__alilen-1);
1042 mean /= __alilen * (__alilen-1) / 2;
1045 vector<contact_t>
_apc(
const ct_vector<fp_t>& __raw_ctscore,
const vector<double>& __apc_bg,
const double __apc_mean,
const uint16_t __mincontsep,
bool __filt )
1047 const uint16_t alilen = __apc_bg.size();
1050 vector<contact_t> res;
1051 res.reserve( alilen*(alilen-1)/2*0.03 );
1053 for(
int i = 0, e = alilen-__mincontsep; i < e; ++i)
1054 for(uint16_t j = i+__mincontsep; j < alilen; ++j)
1055 if(!__filt || __raw_ctscore(i,j) > 0 )
1058 __raw_ctscore(i,j) - __apc_bg[i] * __apc_bg[j] / __apc_mean
1066 const uint16_t alilen = __raw_ctscore.
alilen;
1068 vector<contact_t> res;
1069 res.reserve( alilen*(alilen-1)/2 );
1071 for(
int i = 0, e = alilen-__mincontsep; i < e; ++i)
1072 for(uint16_t j = i+__mincontsep; j < alilen; ++j)
1081 static inline uint32_t
1084 long l1_dcache_size = sysconf(_SC_LEVEL1_DCACHE_SIZE);
1086 long l2_cache_size = sysconf(_SC_LEVEL2_CACHE_SIZE);
1088 long l3_cache_size = sysconf(_SC_LEVEL3_CACHE_SIZE);
1091 uint32_t wchunk = 0;
1093 if(l1_dcache_size > 0) wchunk = 0.5*l1_dcache_size / __seqsize;
1094 if(wchunk < 4 && l2_cache_size > 0) wchunk = 0.5*l2_cache_size / __seqsize;
1095 if(wchunk < 4 && l3_cache_size > 0) wchunk = 0.5*l3_cache_size / __seqsize;
1096 if(wchunk < 4){ wchunk = 0.5*1048576 / __seqsize;
1097 if(wchunk < 4) wchunk = 4;