test/alignment2.cc

Code
Comments
Other
Rev Date Author Line
2847 22 Sep 12 peter 1 // $Id$
2847 22 Sep 12 peter 2
2847 22 Sep 12 peter 3 /*
2847 22 Sep 12 peter 4   Copyright (C) 2012 Peter Johansson
2847 22 Sep 12 peter 5
2847 22 Sep 12 peter 6   This file is part of the yat library, http://dev.thep.lu.se/yat
2847 22 Sep 12 peter 7
2847 22 Sep 12 peter 8   The yat library is free software; you can redistribute it and/or
2847 22 Sep 12 peter 9   modify it under the terms of the GNU General Public License as
2847 22 Sep 12 peter 10   published by the Free Software Foundation; either version 3 of the
2847 22 Sep 12 peter 11   License, or (at your option) any later version.
2847 22 Sep 12 peter 12
2847 22 Sep 12 peter 13   The yat library is distributed in the hope that it will be useful,
2847 22 Sep 12 peter 14   but WITHOUT ANY WARRANTY; without even the implied warranty of
2847 22 Sep 12 peter 15   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
2847 22 Sep 12 peter 16   General Public License for more details.
2847 22 Sep 12 peter 17
2847 22 Sep 12 peter 18   You should have received a copy of the GNU General Public License
2847 22 Sep 12 peter 19   along with yat. If not, see <http://www.gnu.org/licenses/>.
2847 22 Sep 12 peter 20 */
2847 22 Sep 12 peter 21
2881 18 Nov 12 peter 22 #include <config.h>
2881 18 Nov 12 peter 23
2847 22 Sep 12 peter 24 #include "Suite.h"
2847 22 Sep 12 peter 25
2847 22 Sep 12 peter 26 #include "yat/utility/Aligner.h"
2847 22 Sep 12 peter 27 #include "yat/utility/Alignment.h"
2847 22 Sep 12 peter 28
2847 22 Sep 12 peter 29 #include "yat/utility/Matrix.h"
2847 22 Sep 12 peter 30 #include "yat/utility/utility.h"
2847 22 Sep 12 peter 31
2847 22 Sep 12 peter 32 #include <gsl/gsl_cdf.h>
2847 22 Sep 12 peter 33
2847 22 Sep 12 peter 34 #include <cassert>
2847 22 Sep 12 peter 35 #include <cmath>
2847 22 Sep 12 peter 36 #include <fstream>
2847 22 Sep 12 peter 37 #include <string>
2847 22 Sep 12 peter 38 #include <utility>
2847 22 Sep 12 peter 39 #include <vector>
2847 22 Sep 12 peter 40
2847 22 Sep 12 peter 41 using namespace theplu::yat;
2847 22 Sep 12 peter 42
2847 22 Sep 12 peter 43 int main(int argc, char* argv[])
2847 22 Sep 12 peter 44 {
2847 22 Sep 12 peter 45   test::Suite suite(argc, argv);
2847 22 Sep 12 peter 46
2847 22 Sep 12 peter 47   std::ifstream is(test::filename("data/isoform.peaks").c_str());
2847 22 Sep 12 peter 48   std::vector<std::vector<double> > data;
2847 22 Sep 12 peter 49   utility::load(is, data, '\0', '\n', true, false);
2847 22 Sep 12 peter 50   is.close();
2847 22 Sep 12 peter 51
2847 22 Sep 12 peter 52   assert(data.size()==4);
2847 22 Sep 12 peter 53
2847 22 Sep 12 peter 54   const std::vector<double>& x(data[2]);
2847 22 Sep 12 peter 55   const std::vector<double>& y(data[3]);
2847 22 Sep 12 peter 56
2847 22 Sep 12 peter 57   // create score matrix
2847 22 Sep 12 peter 58   suite.out() << "calculate score matrix\n";
2847 22 Sep 12 peter 59   utility::Matrix score(x.size(), y.size());
2847 22 Sep 12 peter 60   double sqrt2 = std::sqrt(2);
2847 22 Sep 12 peter 61   for (size_t i=0; i<x.size(); ++i)
2847 22 Sep 12 peter 62     for (size_t j=0; j<y.size(); ++j) {
2847 22 Sep 12 peter 63       score(i,j) = 2*gsl_cdf_gaussian_Q(std::abs(x[i]-y[j]), sqrt2);
2847 22 Sep 12 peter 64       if (score(i,j) > 0.001)
2847 22 Sep 12 peter 65         suite.out() << i << " " << j << " "
2847 22 Sep 12 peter 66                     << x[i] << " " << y[j] << " "
2847 22 Sep 12 peter 67                     << score(i, j) << "\n";
2847 22 Sep 12 peter 68     }
2847 22 Sep 12 peter 69
2847 22 Sep 12 peter 70   std::vector<std::pair<size_t, size_t> > path;
2847 22 Sep 12 peter 71   double nw = utility::NeedlemanWunsch(score, path, 0);
2847 22 Sep 12 peter 72
2847 22 Sep 12 peter 73   suite.out() << "NW score: " << nw << "\n";
2847 22 Sep 12 peter 74
2847 22 Sep 12 peter 75   if (!suite.equal_fix(nw, 2.50342, 1e-5)) {
2847 22 Sep 12 peter 76     suite.add(false);
2847 22 Sep 12 peter 77     suite.err() << "unexpected NW score\n";
2847 22 Sep 12 peter 78   }
2847 22 Sep 12 peter 79
2847 22 Sep 12 peter 80   suite.out() << "path:\n";
2847 22 Sep 12 peter 81   for (size_t i=0; i<path.size(); ++i)
2847 22 Sep 12 peter 82     suite.out() << path[i].first << " " << path[i].second << "\n";
2847 22 Sep 12 peter 83
2847 22 Sep 12 peter 84   std::vector<std::pair<size_t, size_t> > must_be;
2847 22 Sep 12 peter 85   must_be.push_back(std::make_pair(0, 1));
2847 22 Sep 12 peter 86   must_be.push_back(std::make_pair(3, 5));
2847 22 Sep 12 peter 87   must_be.push_back(std::make_pair(4, 10));
2847 22 Sep 12 peter 88   must_be.push_back(std::make_pair(5, 12));
2847 22 Sep 12 peter 89
2847 22 Sep 12 peter 90   // check that these entries are aligned, i.e., that they are in
2847 22 Sep 12 peter 91   // returned path
2847 22 Sep 12 peter 92   for (size_t i=0; i<must_be.size(); ++i)
2847 22 Sep 12 peter 93     if (std::find(path.begin(), path.end(), must_be[i])==path.end()) {
2847 22 Sep 12 peter 94       suite.add(false);
2847 22 Sep 12 peter 95       suite.err() << "error: did not find entry "
2847 22 Sep 12 peter 96                   << must_be[i].first << " " << must_be[i].second
2847 22 Sep 12 peter 97                   << " in path\n";
2847 22 Sep 12 peter 98     }
2847 22 Sep 12 peter 99
2847 22 Sep 12 peter 100   return suite.return_value();
2847 22 Sep 12 peter 101 }

Generated on Sat Oct 14 02:20:36 2023 (UTC) by svndigest 0.11.2pre (dev build 1672)