diff --git a/src/main_test.cpp b/src/main_test.cpp
deleted file mode 100644
index 5a9e265..0000000
--- a/src/main_test.cpp
+++ /dev/null
@@ -1,404 +0,0 @@
-#include <iostream>
-#include <string>
-#include <vector>
-#include <utility>
-#include <cmath>
-#include <unordered_map>
-
-#include <Matrix2D.hpp>
-#include <Matrix3D.hpp>
-#include <Matrix4D.hpp>
-#include <ConsensusSequenceModelComputer.hpp>
-#include <EMConsensusSequence.hpp>
-#include <SequenceModelComputer.hpp>
-#include <EMSequence.hpp>
-#include <dna_utility.hpp>
-#include <Statistics.hpp>
-#include <sorting_utility.hpp>
-
-#include <boost/math/distributions.hpp>
-
-#include <KmersStatistics.hpp>
-
-template<class T>
-std::ostream& operator << (std::ostream& stream, const std::vector<T>& v)
-{   for(const auto& x : v)
-    {   stream << x << " " ; }
-    return stream ;
-}
-
-template<class T, class U>
-std::ostream& operator << (std::ostream& stream, const std::unordered_map<T,U>& m)
-{   for(const auto& bucket : m)
-    {   stream << bucket.first << " " << bucket.second << std::endl ; }
-    return stream ;
-}
-
-
-Matrix3D<double> compute_consensus_sequences(const Matrix2D<int>& seq,
-                                             const Matrix4D<double>& prob,
-                                             size_t class_k)
-{
-    size_t n_seq   = seq.get_nrow() ;
-    size_t l_seq   = seq.get_ncol() ;
-    // size_t n_class = prob.get_dim()[1] ;
-    size_t n_shift = prob.get_dim()[2] ;
-    size_t n_flip  = prob.get_dim()[3] ;
-    bool flip      = n_flip == 2 ;
-    size_t l_model = l_seq - n_shift + 1 ;
-
-    Matrix3D<double> cons_seq(n_seq, l_model, 4, 0.) ;
-
-    // aggregate
-    for(size_t i=0; i<n_seq; i++)
-    {   std::cerr << "i " << i << std::endl ;
-        for(size_t s=0; s<n_shift; s++)
-        {   //std::cerr << "s " << s << std::endl ;
-            double tot = 0. ;
-            // aggregate
-            for(size_t j=0; j<l_model; j++)
-            {   // std::cerr << "j " << j << std::endl ;
-                // std::cerr << "seq(" << i << "," << s+j << ")" << std::endl ;
-                int base    = seq(i,s+j) ;
-                // N should participate to A,C,G,T equivalently -> same as doing nothing
-                if(base == dna::char_to_int('N'))
-                {   continue ; }
-                int base_rv = 4 - base - 1 ;
-                // --------- forward ---------
-                // std::cerr << "cons_seq(" << i << "," << j << "," << base << ") += prob(" << i << "," << class_k << "," << 0 << std::endl ;
-                cons_seq(i,j,base)    += prob(i,class_k,s,0) ;
-                tot += prob(i,class_k,s,0) ;
-
-                // --------- reverse ---------
-                if(flip)
-                {   // std::cerr << "cons_seq(" << i << "," << j << "," << base_rv << ") += prob(" << i << "," << class_k << "," << 1 << std::endl ;
-                    cons_seq(i,j,base_rv) += prob(i,class_k,s,1) ;
-                    tot += prob(i,class_k,s,1) ;
-                }
-            }
-        }
-    }
-
-    // normalize
-    for(size_t i=0; i<n_seq; i++)
-    {   for(size_t j=0; j<l_model; j++)
-        {   double tot = 0. ;
-            for(size_t k=0; k<4; k++)
-            {   tot += cons_seq(i,j,k) ; }
-            for(size_t k=0; k<4; k++)
-            {   cons_seq(i,j,k) /= tot ; }
-        }
-    }
-
-    return cons_seq ;
-}
-
-Matrix2D<int> char_to_int(const Matrix2D<char>& seq)
-{
-    size_t n_row = seq.get_nrow() ;
-    size_t n_col = seq.get_ncol() ;
-
-    Matrix2D<int> iseq(n_row, n_col, 0.) ;
-
-    for(size_t i=0; i<n_row; i++)
-    {   for(size_t j=0; j<n_col; j++)
-        {   iseq(i,j) = dna::char_to_int(seq(i,j)) ; }
-    }
-    return iseq ;
-}
-
-Matrix3D<double> sequence_to_consensus(const Matrix2D<int>& seq)
-{
-    size_t n_row  = seq.get_nrow() ;
-    size_t n_col  = seq.get_ncol() ;
-    size_t n_dim3 = 4 ; // A, C, G, T
-
-    int base_n = dna::char_to_int('N') ;
-
-    Matrix3D<double> cseq(n_row, n_col, n_dim3, 1e-100) ;
-
-
-    double prob = 1. - 3*(1e-100) ; // p(base) - p(!base)
-
-    for(size_t i=0; i<n_row; i++)
-    {   for(size_t j=0; j<n_col; j++)
-        {   int base = seq(i,j) ;
-            // N
-            if(base == base_n)
-            {   cseq(i,j,0) = 0.25 ;
-                cseq(i,j,1) = 0.25 ;
-                cseq(i,j,2) = 0.25 ;
-                cseq(i,j,3) = 0.25 ;
-            }
-            // A C G T
-            else
-            {   cseq(i,j, base) =  prob ; }
-        }
-    }
-    return cseq ;
-}
-
-double score_consensus_subseq(const Matrix3D<double>& cons_seq,
-                                                      size_t row,
-                                                      size_t start,
-                                                      const Matrix2D<double>& model)
-{   size_t ncol = cons_seq.get_dim()[1] ;
-    // size_t dim1 = cons_seq.get_dim()[0] ;
-    // size_t dim2 = cons_seq.get_dim()[1] ;
-    size_t dim3 = cons_seq.get_dim()[2] ;
-
-    if(start > ncol - model.get_nrow())
-    {   char msg[4096] ;
-        sprintf(msg, "Error! given start (%zu) is too high. Max value is %zu",
-                start, ncol - model.get_nrow()) ;
-        throw std::invalid_argument(msg) ;
-    }
-    else if(model.get_nrow() > ncol)
-    {   char msg[4096] ;
-        sprintf(msg, "Error! given model is longer than sequences (%zu / %zu)",
-                model.get_nrow(), ncol) ;
-        throw std::invalid_argument(msg) ;
-    }
-    else if(model.get_ncol() != 4)
-    {   char msg[4096] ;
-        sprintf(msg, "Error! given model 2nd dimension is not 4 (%zu)",
-                model.get_ncol()) ;
-        throw std::invalid_argument(msg) ;
-    }
-    else if(dim3 != 4)
-    {   char msg[4096] ;
-        sprintf(msg, "Error! given data 3rd dimension is not 4 (%zu)",
-                dim3) ;
-        throw std::invalid_argument(msg) ;
-    }
-
-    /*
-
-    double ll = 0 ;
-    for(size_t j_seq=from, j_mod=0; j_seq<to; j_seq++, j_mod++)
-    {   for(size_t k=0; k<dim3; k++)
-        {   std::cerr << ll
-                      << " *= "
-                      << "("
-                      << model(j_mod,k)
-                      << " * "
-                      << cons_seq(row,j_seq,k)
-                      << ")"
-                      << std::endl ;
-            ll *= model(j_mod,k) * cons_seq(row,j_seq,k) ;
-        }
-    }*/
-
-    size_t from = start ;
-    size_t to   = from + model.get_nrow() ; // will score [from, to)
-
-    double ll = 0. ;
-
-    for(size_t j_seq=start, j_mod=0; j_seq<to; j_seq++, j_mod++)
-    {   double sum = 0. ;
-        for(size_t k=0; k<dim3; k++)
-        {   std::cerr << sum
-                      << "+= ("
-                      << cons_seq(row, j_seq, k)
-                      << " * "
-                      << model(j_mod, k)
-                      << ")"
-                      << std::endl ;
-            sum += (cons_seq(row, j_seq, k) * model(j_mod, k)) ;
-        }
-        std::cerr << sum << std::endl ;
-        std::cerr << "---------------" << std::endl ;
-        ll += log(sum) ;
-    }
-
-    return ll ;
-}
-
-
-Matrix2D<double> test_consensus(const std::string& file,
-                                size_t n_class,
-                                size_t n_iter,
-                                size_t n_shift,
-                                bool flip,
-                                bool bckg_class,
-                                const std::string& seed,
-                                size_t n_thread)
-{   size_t n_flip = flip + 1 ;
-
-    Matrix2D<int> seq(file) ;
-    // Matrix2D<int> seq(char_to_int(Matrix2D<char>(file))) ;
-    Matrix3D<double> cseq = sequence_to_consensus(seq) ;
-
-    // partition
-    EMConsensusSequence em(cseq,
-                           n_class,
-                           n_iter,
-                           n_shift,
-                           flip,
-                           bckg_class,
-                           seed,
-                           n_thread) ;
-    em.classify() ;
-
-    // get post prob
-    Matrix4D<double> prob = em.get_post_prob() ;
-
-    // compute models
-    ConsensusSequenceModelComputer mc(std::move(cseq),
-                                      prob,
-                                      bckg_class,
-                                      n_thread) ;
-
-    Matrix2D<double> model = mc.get_model() ;
-
-    // compute the class prob
-    size_t n_row   = prob.get_dim()[0] ;
-    std::vector<double> class_prob(n_class, 0.) ;
-    double p_tot = 0. ;
-    for(size_t i=0; i<n_row; i++)
-    {   for(size_t j=0; j<n_class; j++)
-        {   for(size_t k=0; k<n_shift; k++)
-            {   for(size_t l=0; l<n_flip; l++)
-                {   class_prob[j] += prob(i,j,k,l) ;
-                    p_tot         += prob(i,j,k,l) ;
-                }
-            }
-        }
-    }
-    for(auto& prob : class_prob)
-    {   prob /= p_tot ; }
-    // create a matrix containing the class prob in the 1st
-    // column and the model in the remaining columns
-    Matrix2D<double> model_final(model.get_nrow(),
-                                 model.get_ncol() + 1) ;
-    size_t i_class = 0 ;
-    for(size_t i=0; i<model_final.get_nrow(); i++)
-    {   if((i != 0) and (i % 4 == 0))
-        {   i_class++ ; }
-        model_final(i,0) = class_prob[i_class] ;
-    }
-    // fill the remaining with the model parameters
-    for(size_t i=0; i<model.get_nrow(); i++)
-    {   for(size_t j=0; j<model.get_ncol(); j++)
-        {   model_final(i,j+1) = model(i,j) ; }
-    }
-    return model_final ;
-}
-
-Matrix2D<double> test_sequence(const std::string& file,
-                               size_t n_class,
-                               size_t n_iter,
-                               size_t n_shift,
-                               bool flip,
-                               bool bckg_class,
-                               const std::string& seed,
-                               size_t n_thread)
-{   size_t n_flip = flip + 1 ;
-
-    Matrix2D<int> seq(file) ;
-    // Matrix2D<int> seq(char_to_int(Matrix2D<char>(file))) ;
-
-    // partition
-    EMSequence em(seq,
-                  n_class,
-                  n_iter,
-                  n_shift,
-                  flip,
-                  bckg_class,
-                  seed,
-                  n_thread) ;
-    em.classify() ;
-
-    // get post prob
-    Matrix4D<double> prob = em.get_post_prob() ;
-
-    // compute models
-    SequenceModelComputer mc(std::move(seq),
-                             prob,
-                             bckg_class,
-                             n_thread) ;
-
-    Matrix2D<double> model = mc.get_model() ;
-
-    // compute the class prob
-    size_t n_row   = prob.get_dim()[0] ;
-    std::vector<double> class_prob(n_class, 0.) ;
-    double p_tot = 0. ;
-    for(size_t i=0; i<n_row; i++)
-    {   for(size_t j=0; j<n_class; j++)
-        {   for(size_t k=0; k<n_shift; k++)
-            {   for(size_t l=0; l<n_flip; l++)
-                {   class_prob[j] += prob(i,j,k,l) ;
-                    p_tot         += prob(i,j,k,l) ;
-                }
-            }
-        }
-    }
-    for(auto& prob : class_prob)
-    {   prob /= p_tot ; }
-    // create a matrix containing the class prob in the 1st
-    // column and the model in the remaining columns
-    Matrix2D<double> model_final(model.get_nrow(),
-                                 model.get_ncol() + 1) ;
-    size_t i_class = 0 ;
-    for(size_t i=0; i<model_final.get_nrow(); i++)
-    {   if((i != 0) and (i % 4 == 0))
-        {   i_class++ ; }
-        model_final(i,0) = class_prob[i_class] ;
-    }
-    // fill the remaining with the model parameters
-    for(size_t i=0; i<model.get_nrow(); i++)
-    {   for(size_t j=0; j<model.get_ncol(); j++)
-        {   model_final(i,j+1) = model(i,j) ; }
-    }
-    return model_final ;
-}
-
-
-
-
-int main()
-{
-    std::string file = "/local/groux/scATAC-seq/data/"
-                       "10xgenomics_PBMC_5k_motifs/sp1_motifs_10e-7_sequences.mat" ;
-    sequence_to_consensus(Matrix2D<int>(file)).save("sp1_motifs_10e-7_consensussequences.mat3D") ;
-
-    /*
-    size_t n_class   = 2 ;
-    size_t n_iter    = 50 ;
-    size_t n_shift   = 1 ;
-    bool flip        = true ;
-    bool bckg_class  = false ;
-    std::string seed = "" ;
-    size_t n_thread  = 1 ;
-
-    Matrix2D<double> model_final = test_consensus(file,
-                                                   n_class,
-                                                   n_iter,
-                                                   n_shift,
-                                                   flip,
-                                                   bckg_class,
-                                                   seed,
-                                                   n_thread) ;
-    */
-
-    /*
-    Matrix2D<double> model_final = test_sequence(file,
-                                                 n_class,
-                                                 n_iter,
-                                                 n_shift,
-                                                 flip,
-                                                 bckg_class,
-                                                 seed,
-                                                 n_thread) ;
-    */
-
-
-    /*
-    Matrix2D<double> model_final = test_aggregate(file) ;
-    */
-
-    // std::cout << model_final << std::endl ;
-
-    return EXIT_SUCCESS ;
-}