238 const ali_t& __ali,
double __cp,
239 bool __veczw,
int __num_threads
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;
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
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
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()));
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;
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";
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);
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); }
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);
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);
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 );