@@ -1,5 +1,5 @@

 Package: seqbias

-Version: 0.99.6

+Version: 0.99.7

 Date: 25-12-2010

 Title: Estimation of per-position bias in high-throughput sequencing data

 Description: This package implements a model of per-position sequencing bias in

@@ -13,7 +13,7 @@ YAML_CPP_OBJ = aliascontent.o conversion.o emitter.o emitterstate.o \

                stream.o tag.o

 SEQBIAS_OBJ = sequencing_bias.o kmers.o miscmath.o common.o \

-			  table.o superfasthash.o logger.o seqbias.o \

+			  table.o logger.o seqbias.o \

 GNULIB_OBJ = asprintf.o printf-args.o printf-parse.o vasnprintf.o vasprintf.o

@@ -5,6 +5,8 @@

 #include <algorithm>

 #include <cstdlib>

 #include "samtools/faidx.h"

+#include "samtools/sam.h"

+#include "samtools/bam.h"

 typedef long          pos;     /* genomic position */

@@ -24,6 +26,12 @@ void rev( T* xs, int n ) {

     while( i < j ) std::swap(xs[i++],xs[j--]);

 }

+

+/* defined in bam_aux.c */

+extern "C" {

+void bam_init_header_hash(bam_header_t *header);

+}

+

 /* Change the behavior of the faidx_fetch_seq function to be more useful. If

  * coordinates are outside the actual sequence, write N's, rather than adjusting

  * the start,end. */

@@ -3,7 +3,6 @@

 #include "common.hpp"

 #include "miscmath.hpp"

 #include "logger.h"

-#include "asprintf.h"

 #include <string>

 #include <cstring>

@@ -24,8 +23,6 @@ kmer_matrix::kmer_matrix( const YAML::Node& node )

     for( i = 0; i < n*m; i++ ) {

         node_A[i] >> A[i];

     }

-

-    stored_row = new double[ m ];

 }

 kmer_matrix::kmer_matrix( size_t n, size_t k )

@@ -34,7 +31,6 @@ kmer_matrix::kmer_matrix( size_t n, size_t k )

     m = 1<<(2*k); /* == 4^k */

     A = new double[ n * m ];

-    stored_row = new double[ m ];

 }

 kmer_matrix::kmer_matrix( const kmer_matrix& M )

@@ -45,10 +41,6 @@ kmer_matrix::kmer_matrix( const kmer_matrix& M )

     m = M.m;

     A = new double[ n * m ];

     memcpy( (void*)A, (void*)M.A, n * m * sizeof(double) );

-

-

-    stored_row = new double[ m ];

-    memcpy( (void*)stored_row, (void*)M.stored_row, m * sizeof(double) );

 }

 size_t kmer_matrix::getn() const { return n; }

@@ -66,13 +58,9 @@ void kmer_matrix::operator=( const kmer_matrix& M )

         delete[] A;

         A = new double[ n * m ];

-

-        delete[] stored_row;

-        stored_row = new double[ m ];

     }

     memcpy( (void*)A, (void*)M.A, n * m * sizeof(double) );

-    memcpy( (void*)stored_row, (void*)M.stored_row, m * sizeof(double) );

 }

@@ -81,11 +69,11 @@ void kmer_matrix::to_yaml( YAML::Emitter& out ) const

     out << YAML::BeginMap;

     out << YAML::Key   << "k";

-    out << YAML::Value << (unsigned int)k;

+    out << YAML::Value << k;

     out << YAML::Key   << "n";

-    out << YAML::Value << (unsigned int)n;

+    out << YAML::Value << n;

     out << YAML::Key   << "m";

-    out << YAML::Value << (unsigned int)m;

+    out << YAML::Value << m;

     out << YAML::Key   << "A";

     out << YAML::Flow;

     out << YAML::Value;

@@ -103,7 +91,6 @@ void kmer_matrix::to_yaml( YAML::Emitter& out ) const

 kmer_matrix::~kmer_matrix()

 {

     delete[] A;

-    delete[] stored_row;

 }

@@ -119,7 +106,7 @@ void kmer_matrix::setall( double x )

     for( i = 0; i < n * m; i++ ) A[i] = x;

 }

-void kmer_matrix::setrow( size_t i, double x )

+void kmer_matrix::setrowall( size_t i, double x )

 {

     size_t j;

     for( j = 0; j < m; j++ ) {

@@ -127,6 +114,16 @@ void kmer_matrix::setrow( size_t i, double x )

     }

 }

+void kmer_matrix::getrow( size_t i, double* xs )

+{

+    memcpy( (void*)xs, (void*)(A + i * m), m * sizeof(double) );

+}

+

+void kmer_matrix::setrow( size_t i, double* xs )

+{

+    memcpy( (void*)(A + i * m ), (void*)xs, m * sizeof(double) );

+}

+

 void kmer_matrix::dist_normalize()

 {

@@ -155,21 +152,30 @@ void kmer_matrix::dist_conditionalize( int effective_k )

 }

-void kmer_matrix::dist_conditionalize_row( size_t i, int effective_k  )

+void kmer_matrix::dist_conditionalize_row( size_t i, size_t j, int effective_k  )

 {

     if( effective_k <= 0 ) effective_k = m;

-

-    kmer L;

-    kmer L_max = 1 << (2*(effective_k-1));

+    kmer L, R;

+    kmer L_max = 1 << (2*j);

+    kmer R_max = 1 << (2*(effective_k-j-1));

     kmer K;

-    kmer nt;

     double z;

+    kmer nt;

     for( L = 0; L < L_max; L++ ) {

-        K = L<<2;

-        z = 0.0;

-        for( nt = 0; nt < 4; nt++ ) z += A[ i * m + (nt | K) ];

-        for( nt = 0; nt < 4; nt++ ) A[ i * m + (nt | K) ] /= z;

+        for( R = 0; R < R_max; R++ ) { 

+            z = 0.0;

+

+            for( nt = 0; nt < 4; nt++ ) {

+                K = (L<<(2*(effective_k-j))) | (nt<<(2*(effective_k-j-1))) | R;

+                z += A[ i*m + K ];

+            }

+

+            for( nt = 0; nt < 4; nt++ ) {

+                K = (L<<(2*(effective_k-j))) | (nt<<(2*(effective_k-j-1))) | R;

+                A[ i*m + K ] /= z;

+            }

+        }

     }

 }

@@ -186,24 +192,11 @@ void kmer_matrix::log_transform_row( size_t i, int effective_k )

 }

-void kmer_matrix::store_row( size_t i )

-{

-    memcpy( (void*)stored_row, (void*)(A + i * m), m * sizeof(double) );

-    stored_row_index = i;

-}

-

-

-void kmer_matrix::restore_stored_row()

-{

-    memcpy( (void*)(A + stored_row_index * m), (void*)stored_row, m * sizeof(double) );

-}

-

-

 const size_t sequence::kmer_max_k = 4*sizeof(kmer);

-sequence::sequence( const char* s, int c )

-    : c(c), xs(NULL), n(0)

+sequence::sequence( const char* s, int c, double w )

+    : c(c), w(w), xs(NULL), n(0)

 {

     n = strlen(s);

     if( n > 0 ) {

@@ -223,6 +216,7 @@ sequence::sequence( const sequence& s )

     : xs(NULL), n(0)

 {

     c = s.c;

+    w = s.w;

     n = s.n;

     if( n > 0 ) {

         xs = new kmer[ n/kmer_max_k + 1 ];

@@ -291,6 +285,10 @@ motif::motif( const YAML::Node& node )

     }

     P = new kmer_matrix( node["P"] );

+

+    R = new bool[n*n];

+    memset( R, 0, n*n*sizeof(bool) );

+    compute_reachability();

 }

 motif::motif( size_t n, size_t k, int c )

@@ -301,6 +299,9 @@ motif::motif( size_t n, size_t k, int c )

     parents = new bool[n*n];

     memset( parents, 0, n*n*sizeof(bool) );

+

+    R = new bool[n*n];

+    memset( R, 0, n*n*sizeof(bool) );

 }

@@ -313,6 +314,9 @@ motif::motif( const motif& M )

     parents = new bool[n*n];

     memcpy( parents, M.parents, n*n*sizeof(bool) );

+

+    R = new bool[n*n];

+    memset( R, 0, n*n*sizeof(bool) );

 }

@@ -320,6 +324,7 @@ motif::motif( const motif& M )

 motif::~motif()

 {

     delete[] parents;

+    delete[] R;

     delete P;

 }

@@ -329,10 +334,10 @@ void motif::to_yaml( YAML::Emitter& out ) const

     out << YAML::BeginMap;

     out << YAML::Key   << "n";

-    out << YAML::Value << (unsigned int)n;

+    out << YAML::Value << n;

     out << YAML::Key   << "k";

-    out << YAML::Value << (unsigned int)k;

+    out << YAML::Value << k;

     out << YAML::Key   << "c";

     out << YAML::Value << c;

@@ -385,7 +390,7 @@ size_t motif::num_parents( size_t i ) const

 {

     size_t j;

     size_t M = 0;

-    for( j = 0; j <= i; j++ ) {

+    for( j = 0; j < n; j++ ) {

         if( parents[i*n+j] ) M++;

     }

@@ -402,19 +407,29 @@ void motif::set_edge( size_t i, size_t j, bool x )

     parents[j*n+i] = x;

 }

-

-void motif::store_row( size_t i )

+bool motif::reachable( size_t i, size_t j )

 {

-    P->store_row( i );

+    return R[j*n+i];

 }

-

-void motif::restore_stored_row()

+void motif::compute_reachability()

 {

-    P->restore_stored_row();

+    /* find the transitive closure of the parents matrix via flyod-warshall */

+

+    memcpy( R, parents, n*n*sizeof(bool) );

+

+    size_t k, i, j;

+    for( k = 0; k < n; k++ ) {

+        for( i = 0; i < n; i++ ) {

+            for( j = 0; j < n; j++ ) {

+                R[j*n+i] = R[j*n+i] || (R[k*n+i] && R[j*n+k]);

+            }

+        }

+    }

 }

+

 char* motif::print_model_graph( int offset )

 {

     std::string graph_str;

@@ -439,7 +454,8 @@ char* motif::print_model_graph( int offset )

     for( j = 0; j < (int)n; j++ ) {

         if( !parents[j*n+j] ) continue;

-        for( i = 0; i < j; i++ ) {

+        for( i = 0; i < (int)n; i++ ) {

+            if( i == j ) continue;

             if( parents[j*n+i] ) {

                 r = asprintf( &tmp, "n%d -> n%d;\n", i, j );

                 graph_str += tmp;

@@ -455,50 +471,13 @@ char* motif::print_model_graph( int offset )

 }

-void motif::add_all_edges( const std::deque<sequence*>* data )

-{

-    size_t i, j;

-    size_t n_parents;

-    size_t m;

-    kmer K;

-    std::deque<sequence*>::const_iterator seq;

-

-    for( j = 0; j < n; j++ ) {

-        for( i = k-1 <= j ? j-(k-1) : 0; i <= j; i++ ) {

-            set_edge( i, j, true );

-        }

-        P->setrow( j, 0.0 );

-        n_parents = num_parents(j);

-        m = 1 << (2*n_parents);

-

-        for( K = 0; K < m; K++ ) {

-            (*P)( j, K ) = pseudocount;

-        }

-

-        for( seq = data->begin(); seq != data->end(); seq++ ) {

-            if( (*seq)->c == c && (*seq)->get( parents + j*n, n, K ) ) {

-                (*P)( j, K )++;

-            }

-        }

-

-        P->dist_normalize_row( j );

-        P->dist_conditionalize_row( j, n_parents );

-        P->log_transform_row( j, n_parents );

-    }

-}

-

-

 /* make an edge i --> j

  * That is, condition j on i. */

 void motif::add_edge( size_t i, size_t j, const std::deque<sequence*>* data )

 {

-    if( i > j ) {

-        failf( "Invalid motif edge (%zu, %zu)\n", i, j );

-    }

-

     set_edge( i, j, true );

-    P->setrow( j, 0.0 );

+    P->setrowall( j, 0.0 );

     size_t n_parents = num_parents(j);

     size_t m = 1 << (2*n_parents);

     kmer K;

@@ -510,26 +489,30 @@ void motif::add_edge( size_t i, size_t j, const std::deque<sequence*>* data )

     std::deque<sequence*>::const_iterator seq;

     for( seq = data->begin(); seq != data->end(); seq++ ) {

         if( (*seq)->c == c && (*seq)->get( parents + j*n, n, K ) ) {

-            (*P)( j, K )++;

+            (*P)( j, K ) += (*seq)->w;

         }

     }

+    size_t n_pred_parents = 0;

+    size_t u;

+    for( u = 0; u < j; u++ ) {

+        if( has_edge( u, j ) ) n_pred_parents++;

+    }

+

     P->dist_normalize_row( j );

-    P->dist_conditionalize_row( j, n_parents );

+    P->dist_conditionalize_row( j, n_pred_parents, n_parents );

     P->log_transform_row( j, n_parents );

+

+    compute_reachability();

 }

 void motif::remove_edge( size_t i, size_t j, const std::deque<sequence*>* data )

 {

-    if( i > j ) {

-        failf( "Invalid motif edge (%zu, %zu)\n", i, j );

-    }

-

     set_edge( i, j, false );

-    P->setrow( j, 0.0 );

+    P->setrowall( j, 0.0 );

     size_t n_parents = num_parents(j);

     size_t m = 1 << (2*n_parents);

     kmer K;

@@ -545,38 +528,22 @@ void motif::remove_edge( size_t i, size_t j, const std::deque<sequence*>* data )

         }

     }

+    size_t n_pred_parents = 0;

+    size_t u;

+    for( u = 0; u < j; u++ ) {

+        if( has_edge( u, j ) ) n_pred_parents++;

+    }

+

     P->dist_normalize_row( j );

-    P->dist_conditionalize_row( j, n_parents );

+    P->dist_conditionalize_row( j, n_pred_parents, n_parents );

     P->log_transform_row( j, n_parents );

-}

+    compute_reachability();

+}

-double motif_log_likelihood( const motif& M0, const motif& M1,

-                             const std::deque<sequence*>* training_seqs )

-{

-    double L, L0, L1;

-

-    L = 0.0;

-    

-    std::deque<sequence*>::const_iterator i;

-    for( i = training_seqs->begin(); i != training_seqs->end(); i++ ) {

-        L0 = M0.eval( **i );

-        L1 = M1.eval( **i );

-        if( (*i)->c == 0 ) {

-            L += L0 - logaddexp( L0, L1 );

-        }

-        else if( (*i)->c == 1 ) {

-            L += L1 - logaddexp( L0, L1 );

-        }

-    }

-

-    return L;

-}

-

-

 /* Let x_i be the sequence of training example i, c be the class, and M be the

  * model. Then given l0 and l1, where,

@@ -590,15 +557,19 @@ double motif_log_likelihood( const motif& M0, const motif& M1,

  *                      = sum_i log( P( x_i | c_i, M ) / ( P( x_i | c = 0, M ) + P( x_i | c = 1, M ) )

  *

  */

-double conditional_likelihood( size_t n, const double* l0, const double* l1, const int* c )

+double conditional_likelihood( size_t n, const double* l0, const double* l1, const int* c,

+                               double prior )

 {

     double l;

     size_t i;

+    double p = log(prior);

+    double q = log(1-prior);

+

     l = 0.0;

     for( i = 0; i < n; i++ ) {

-        l += (c[i] == 0 ? l0[i] : l1[i])

-                    - logaddexp( l0[i], l1[i] );

+        l += (c[i] == 0 ? (q+l0[i]) : (p+l1[i]))

+                    - logaddexp( q+l0[i], p+l1[i] );

     }

     return l;

@@ -615,7 +586,7 @@ void motif::update_likelihood_column( double* L, size_t n, size_t m, size_t j,

     for( i = 0; i < n; i++ )  {

         if( (*training_seqs)[i]->get( parents + j * m, m, K ) ) {

-            L[ i * m + j ] = (*P)( j, K );

+            L[ i * m + j ] = (*training_seqs)[i]->w * (*P)( j, K );

         }

     }

 }

@@ -628,24 +599,24 @@ void motif::update_likelihood_column( double* L, size_t n, size_t m, size_t j,

 /* Akaike Information Criterion */

 double aic( double L, double n_obs, double n_params, double c = 1.0 )

 {

-    return L - c*n_params*n_obs / (n_obs - n_params - 1.0);

+    return L - c*n_params;

 }

-/* Bayesian (Schwarz) Information Criterion */

-double bic( double L, double n_obs, double n_params, double c = 1.0 )

+/* Akaike Information Criterion with second order small-sample bias correction.

+ * (From Hurvich and Tsai (1989)) */

+double aicc( double L, double n_obs, double n_params, double c = 1.0 )

 {

-    return 2.0*L - c*n_params*log(n_obs);

+    return L - c*n_params

+             - (n_params + 1) * (n_params + 2) / (n_obs - n_params - 2);

 }

-

-/* Hannan-Quinn Information Criterion */

-double hqic( double L, double n_obs, double n_params, double c = 1.0 )

+/* Bayesian (Schwarz) Information Criterion */

+double bic( double L, double n_obs, double n_params, double c = 1.0 )

 {

-    return L - c*n_params*log(log(n_obs));

+    return 2.0*L - c*n_params*log(n_obs);

 }

-

 void train_motifs( motif& M0, motif& M1,

                    const std::deque<sequence*>* training_seqs,

                    size_t max_dep_dist, double complexity_penalty )

@@ -658,11 +629,11 @@ void train_motifs( motif& M0, motif& M1,

         failf( "Motif models of mismatching size. (%zu != %zu)\n", M0.n, M1.n );

     }

-    double (*compute_ic)( double, double, double, double ) = aic;

+    double (*compute_ic)( double, double, double, double ) = aicc;

-    size_t i, j;

-    size_t i_start;

+    int i, j;

+    int i_start, i_end;

     const size_t n = training_seqs->size();

     const size_t m = M0.n;

@@ -680,23 +651,51 @@ void train_motifs( motif& M0, motif& M1,

     double* l1 = new double[ n ]; memset( (void*)l1, 0, n * sizeof(double) );

+    /* vectors for saving and restoring parameters, to avoid recomputing

+     * frequencies */

+    double* M0_row_j = new double[ M0.P->getm() ];

+    double* M1_row_j = new double[ M0.P->getm() ];

+    double* M0_row_i = new double[ M0.P->getm() ];

+    double* M1_row_i = new double[ M0.P->getm() ];

+

+

     /* 0-1 vector giving labeling each sequence as foreground or background */

     int* cs = new int[ n ];

-    for( i = 0; i < n; i++ ) {

+    for( i = 0; i < (int)n; i++ ) {

         cs[i] = (*training_seqs)[i]->c;

     }

+    /* set prior probability of an example being foreground */

+    double ws[2];

+    for( i = 0; i < (int)n; i++ ) {

+        if( cs[i] <= 1 ) ws[cs[i]] += (*training_seqs)[i]->w;

+    }

+    double prior = ws[1] / (ws[0] + ws[1]);

+

     /* backup likelihood matrix columns, to restore state after trying a new edge */

-    double* b0 = new double[ n ];

-    double* b1 = new double[ n ];

+    double* b0_j = new double[ n ];

+    double* b1_j = new double[ n ];

+    double* b0_i = new double[ n ];

+    double* b1_i = new double[ n ];

     /* keeping track of the optimal edge */

-    double ic, ic_curr, ic_best;

-    size_t j_best, i_best;

+    double ic, ic_curr;

+

+    double ic_forw_best;

+    int j_forw_best, i_forw_best;

+    double ic_back_best;

+    int j_back_best, i_back_best;

+

+    double ic_rev_best;

+    int j_rev_best, i_rev_best;

+

+

+    /* for cycle detection */

+    bool has_ij_path;

     /* parameters to compute information criterion */

     double n_obs    = training_seqs->size();

@@ -706,51 +705,81 @@ void train_motifs( motif& M0, motif& M1,

     /* log conditional likelihood */

     double l; 

-

     /* baseline ic */

-    l = conditional_likelihood( n, l0, l1, cs );

+    l = conditional_likelihood( n, l0, l1, cs, prior );

     ic_curr = compute_ic( l, n_obs, n_params, complexity_penalty );

+    /* for pretty output */

+    size_t col;

+    const char* col_base = "\n%34s";

+    const size_t col_max = 30;

+

+

     size_t round_num = 0;

     while( true ) {

         round_num++;

         log_printf( LOG_MSG, "round %4d (ic = %0.4e) ", round_num, ic_curr );

+        col = 0;

+

+        ic_forw_best = ic_back_best = ic_rev_best = -HUGE_VAL;

+        i_forw_best = i_back_best = i_rev_best = 0;

+        j_forw_best = j_back_best = j_rev_best = 0;

-        ic_best = -HUGE_VAL;

-        j_best = i_best = 0;

-        for( j = 0; j < M0.n; j++ ) {

+        /* phase 1: try all possible edge additions */

+        for( j = M0.n-1; j >= 0; j-- ) {

             if( M0.has_edge( j, j ) ) {

-                if( max_dep_dist == 0 || j <= max_dep_dist ) i_start = 0;

-                else i_start = i - max_dep_dist;

+                if( max_dep_dist == 0 || j <= (int)max_dep_dist ) i_start = 0;

+                else i_start = j - max_dep_dist;

+

+                if( max_dep_dist == 0 ) i_end = M0.n - 1;

+                else i_end = std::min( M0.n-1, j + max_dep_dist );

             }

-            else i_start = j;

+            else i_start = i_end = j;

+

-            for( i = i_start; i <= j; i++ ) {

+            for( i = i_start; i <= i_end; i++ ) {

+                /* skip existing edges  */

                 if( M0.has_edge( i, j ) ) {

                     continue;

                 }

+                /* skip edges that would introduce cycles */

+                if( M0.reachable( j, i ) ) {

+                    continue;

+                }

+

+                /* skip edges that would exceed the parent limit */

                 if( M0.num_parents(j) >= M0.k ) {

                     continue;

                 }

-                log_puts( LOG_MSG, "." );

+                /* skip edges that are equivalent to one already tried */

+                if( i > j && M0.num_parents(j) == 1 && M0.num_parents(i) == 1 ) {

+                    continue;

+                }

-                /* keep track of the old parameters to avoid retraining */

-                M0.store_row(j);

-                M1.store_row(j);

+

+                log_puts( LOG_MSG, "+" );

+                if( ++col > col_max ) {

+                    col = 0;

+                    log_printf( LOG_MSG, col_base, "" );

+                }

+                /* keep track of the old parameters to avoid retraining */

+                M0.P->getrow( j, M0_row_j );

+                M1.P->getrow( j, M1_row_j );

+

                 /* keep track of the old likelihoods to avoid reevaluating */

-                colcpy( b0, L0, j, n, m  );

-                colcpy( b1, L1, j, n, m );

+                colcpy( b0_j, L0, j, n, m  );

+                colcpy( b1_j, L1, j, n, m );

                 /* add edge */

                 M0.add_edge( i, j, training_seqs );

@@ -762,254 +791,313 @@ void train_motifs( motif& M0, motif& M1,

                 /* update training example likelihoods */

-                vecsub( l0, b0, n );

+                vecsub( l0, b0_j, n );

                 vecaddcol( l0, L0, n, m, j );

-                vecsub( l1, b1, n );

+                vecsub( l1, b1_j, n );

                 vecaddcol( l1, L1, n, m, j );

-                l        = conditional_likelihood( n, l0, l1, cs );

+                l        = conditional_likelihood( n, l0, l1, cs, prior );

                 n_params = M0.num_params() + M1.num_params();

                 ic       = compute_ic( l, n_obs, n_params, complexity_penalty );

-                if( ic > ic_best ) {

-                    ic_best = ic;

-                    i_best = i;

-                    j_best = j;

+

+                if( ic >= ic_forw_best ) {

+                    ic_forw_best = ic;

+                    i_forw_best = i;

+                    j_forw_best = j;

                 }

+                /* TODO: delete me */

+                //log_printf( LOG_MSG, "\n(%zu,%zu) : %0.4e\n", i, j, ic );

                 /* remove edge */

                 M0.set_edge( i, j, false );

                 M1.set_edge( i, j, false );

                 /* restore previous parameters */

-                M0.restore_stored_row();

-                M1.restore_stored_row();

+                M0.P->setrow( j, M0_row_j );

+                M1.P->setrow( j, M1_row_j );

                 /* restore previous likelihoods */

                 vecsubcol( l0, L0, n, m, j );

-                vecadd( l0, b0, n );

+                vecadd( l0, b0_j, n );

                 vecsubcol( l1, L1, n, m, j );

-                vecadd( l1, b1, n );

+                vecadd( l1, b1_j, n );

-                matsetcol( L0, b0, n, m, j );

-                matsetcol( L1, b1, n, m, j );

+                matsetcol( L0, b0_j, n, m, j );

+                matsetcol( L1, b1_j, n, m, j );

             }

         }

-        log_puts( LOG_MSG, "\n" );

-

-        if( ic_best <= ic_curr || ic_best == -HUGE_VAL ) break;

-

-        ic_curr = ic_best;

-

-        M0.add_edge( i_best, j_best, training_seqs );

-        M1.add_edge( i_best, j_best, training_seqs );

-

-        vecsubcol( l0, L0, n, m, j_best );

-        vecsubcol( l1, L1, n, m, j_best );

-

-        M0.update_likelihood_column( L0, n, m, j_best, training_seqs );

-        M1.update_likelihood_column( L1, n, m, j_best, training_seqs );

-

-        vecaddcol( l0, L0, n, m, j_best );

-        vecaddcol( l1, L1, n, m, j_best );

-    }

+        /* phase 2: try all possible edge removals */

+        for( j = 0; j < (int)M0.n; j++ ) {

+            for( i = 0; i < (int)M0.n; i++ ) {

+                if( !M0.has_edge( i, j ) ) continue;

+                if( i == j && M0.num_parents(j) > 1 ) continue;

+                log_puts( LOG_MSG, "-" );

+                if( ++col > col_max ) {

+                    col = 0;

+                    log_printf( LOG_MSG, col_base, "" );

+                }

-    delete[] cs;

-    delete[] L0;

-    delete[] L1;

-    delete[] l0;

-    delete[] l1;

-    delete[] b0;

-    delete[] b1;

-    

-    log_unindent();

-}

-

-

-

-

-void train_motifs_backwards( motif& M0, motif& M1,

-                             const std::deque<sequence*>* training_seqs,

-                             size_t max_dep_dist, double complexity_penalty )

-{

-    log_puts( LOG_MSG, "training motifs (backwards) ...\n" );

-    log_indent();

-

-

-    if( M0.n != M1.n ) {

-        failf( "Motif models of mismatching size. (%zu != %zu)\n", M0.n, M1.n );

-    }

-

-    double (*compute_ic)( double, double, double, double ) = aic;

-

-    size_t i, j;

-    size_t i_last;

-

-    const size_t n = training_seqs->size();

-    const size_t m = M0.n;

-

-    /* likelihood matrices: 

-     * Lk_ij gives the likelihood of training example i on node j in model k.

-     */

-    double* L0 = new double[ n * m ]; memset( (void*)L0, 0, n * m * sizeof(double) );

-    double* L1 = new double[ n * m ]; memset( (void*)L1, 0, n * m * sizeof(double) );

-

-    /* likelihood vectors:

-     * summed columns of L0, L1, maintained to minimize redundant computation.

-     * */

-    double* l0 = new double[ n ]; memset( (void*)l0, 0, n * sizeof(double) );

-    double* l1 = new double[ n ]; memset( (void*)l1, 0, n * sizeof(double) );

+                /* keep track of the old parameters to avoid retraining */

+                M0.P->getrow( j, M0_row_j );

+                M1.P->getrow( j, M1_row_j );

+                /* keep track of the old likelihoods to avoid reevaluating */

+                colcpy( b0_j, L0, j, n, m  );

+                colcpy( b1_j, L1, j, n, m );

-    /* 0-1 vectors giving labeling each sequence as foreground or background */

-    int* cs = new int[ training_seqs->size() ];

+                /* remove edge */

+                M0.remove_edge( i, j, training_seqs );

+                M1.remove_edge( i, j, training_seqs );

-    for( i = 0; i < training_seqs->size(); i++ ) {

-        cs[i] = (*training_seqs)[i]->c;

-    }

+                /* evaluate likelihoods for that column */

+                M0.update_likelihood_column( L0, n, m, j, training_seqs );

+                M1.update_likelihood_column( L1, n, m, j, training_seqs );

+                /* update training example likelihoods */

+                vecsub( l0, b0_j, n );

+                vecaddcol( l0, L0, n, m, j );

-    /* backup likelihood matrix columns, to restore state after trying a new edge */

-    double* b0 = new double[ n ];

-    double* b1 = new double[ n ];

-    

+                vecsub( l1, b1_j, n );

+                vecaddcol( l1, L1, n, m, j );

+                /* evaluate likelihood / ic */

+                l        = conditional_likelihood( n, l0, l1, cs, prior );

+                n_params = M0.num_params() + M1.num_params();

+                ic       = compute_ic( l, n_obs, n_params, complexity_penalty );

-    /* keeping track of the optimal edge */

-    double ic, ic_curr, ic_best;

-    size_t j_best, i_best;

+                if( ic > ic_back_best ) {

+                    ic_back_best = ic;

+                    i_back_best = i;

+                    j_back_best = j;

+                }

+                /* TODO: delete me */

+                //log_printf( LOG_MSG, "\n(%zu,%zu) : %0.4e\n", i, j, ic );

-    /* start with all edges */

-    M0.add_all_edges( training_seqs );

-    M1.add_all_edges( training_seqs );

+                /* replace edge */

+                M0.set_edge( i, j, true );

+                M1.set_edge( i, j, true );

+                /* restore previous parameters */

+                M0.P->setrow( j, M0_row_j );

+                M1.P->setrow( j, M1_row_j );

-    /* parameters to compute information criterion */

-    double n_obs    = training_seqs->size();

-    double n_params = M0.num_params() + M1.num_params();