test/alignment.cc

Code
Comments
Other
Rev Date Author Line
261 09 Mar 05 peter 1 // $Id$
261 09 Mar 05 peter 2
675 10 Oct 06 jari 3 /*
2121 13 Dec 09 peter 4   Copyright (C) 2005 Jari Häkkinen, Peter Johansson
2121 13 Dec 09 peter 5   Copyright (C) 2006 Jari Häkkinen
2121 13 Dec 09 peter 6   Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
4359 23 Aug 23 peter 7   Copyright (C) 2012 Peter Johansson
261 09 Mar 05 peter 8
1437 25 Aug 08 peter 9   This file is part of the yat library, http://dev.thep.lu.se/yat
261 09 Mar 05 peter 10
675 10 Oct 06 jari 11   The yat library is free software; you can redistribute it and/or
675 10 Oct 06 jari 12   modify it under the terms of the GNU General Public License as
1486 09 Sep 08 jari 13   published by the Free Software Foundation; either version 3 of the
675 10 Oct 06 jari 14   License, or (at your option) any later version.
675 10 Oct 06 jari 15
675 10 Oct 06 jari 16   The yat library is distributed in the hope that it will be useful,
675 10 Oct 06 jari 17   but WITHOUT ANY WARRANTY; without even the implied warranty of
675 10 Oct 06 jari 18   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
675 10 Oct 06 jari 19   General Public License for more details.
675 10 Oct 06 jari 20
675 10 Oct 06 jari 21   You should have received a copy of the GNU General Public License
1487 10 Sep 08 jari 22   along with yat. If not, see <http://www.gnu.org/licenses/>.
675 10 Oct 06 jari 23 */
675 10 Oct 06 jari 24
2881 18 Nov 12 peter 25 #include <config.h>
2881 18 Nov 12 peter 26
1231 14 Mar 08 peter 27 #include "Suite.h"
1231 14 Mar 08 peter 28
2815 28 Aug 12 peter 29 #include "yat/utility/Aligner.h"
675 10 Oct 06 jari 30 #include "yat/utility/Alignment.h"
675 10 Oct 06 jari 31
1121 22 Feb 08 peter 32 #include "yat/utility/Matrix.h"
675 10 Oct 06 jari 33
261 09 Mar 05 peter 34 #include <gsl/gsl_cdf.h>
261 09 Mar 05 peter 35
442 15 Dec 05 jari 36 #include <cassert>
261 09 Mar 05 peter 37 #include <cmath>
261 09 Mar 05 peter 38 #include <fstream>
261 09 Mar 05 peter 39 #include <iostream>
261 09 Mar 05 peter 40 #include <sstream>
261 09 Mar 05 peter 41 #include <string>
261 09 Mar 05 peter 42 #include <vector>
261 09 Mar 05 peter 43
680 11 Oct 06 jari 44 using namespace theplu::yat;
301 30 Apr 05 peter 45
261 09 Mar 05 peter 46 void read_vector(std::vector<double>& v,std::string& str)
261 09 Mar 05 peter 47 {
261 09 Mar 05 peter 48   std::istringstream s(str);
261 09 Mar 05 peter 49   while (!s.eof()) {
261 09 Mar 05 peter 50     double d;
261 09 Mar 05 peter 51     s >> d;
261 09 Mar 05 peter 52     if (!s.fail())
261 09 Mar 05 peter 53       v.push_back(d);
261 09 Mar 05 peter 54   }
261 09 Mar 05 peter 55 }
261 09 Mar 05 peter 56
261 09 Mar 05 peter 57 std::istream& operator>>(std::istream& s,
261 09 Mar 05 peter 58                          std::vector<std::vector<double> >& m)
261 09 Mar 05 peter 59 {
261 09 Mar 05 peter 60   while (!s.eof()) {
261 09 Mar 05 peter 61     std::string row;
261 09 Mar 05 peter 62     char ch;
261 09 Mar 05 peter 63     while (((ch=s.get())!='\n') && s.good())
261 09 Mar 05 peter 64       row+=ch;
261 09 Mar 05 peter 65     if (row.size() || s.good()) {  // if stream good, read an empty peak set
261 09 Mar 05 peter 66       std::vector<double> v;
261 09 Mar 05 peter 67       read_vector(v,row);
261 09 Mar 05 peter 68       m.push_back(v);
261 09 Mar 05 peter 69     }
261 09 Mar 05 peter 70   }
261 09 Mar 05 peter 71   return s;
261 09 Mar 05 peter 72 }
261 09 Mar 05 peter 73
261 09 Mar 05 peter 74 double match(const double x, const double y, const double s)
4200 19 Aug 22 peter 75 {
4200 19 Aug 22 peter 76   return 2*gsl_cdf_gaussian_Q(std::abs(x-y),sqrt(2)*s);
261 09 Mar 05 peter 77 }
261 09 Mar 05 peter 78
4200 19 Aug 22 peter 79
261 09 Mar 05 peter 80 double score(const std::vector<double>& l1,
261 09 Mar 05 peter 81              const std::vector<double>& l2,
4200 19 Aug 22 peter 82              const double sigma,
1121 22 Feb 08 peter 83              theplu::yat::utility::Matrix& dot_matrix,
261 09 Mar 05 peter 84              std::vector<std::pair<size_t,size_t> >& path)
261 09 Mar 05 peter 85 {
1098 18 Feb 08 jari 86   dot_matrix.resize(l1.size(),l2.size());
4200 19 Aug 22 peter 87   for (size_t i=0; i<l1.size(); i++)
261 09 Mar 05 peter 88     for (size_t j=0; j<l2.size(); j++) {
261 09 Mar 05 peter 89       assert(i<dot_matrix.rows());
261 09 Mar 05 peter 90       assert(j<dot_matrix.columns());
261 09 Mar 05 peter 91       dot_matrix(i,j)=match(l1[i],l2[j],sigma);
261 09 Mar 05 peter 92     }
680 11 Oct 06 jari 93   return theplu::yat::utility::NeedlemanWunsch(dot_matrix,path,0);
261 09 Mar 05 peter 94 }
261 09 Mar 05 peter 95
261 09 Mar 05 peter 96
1231 14 Mar 08 peter 97 int main(int argc, char* argv[])
261 09 Mar 05 peter 98 {
1231 14 Mar 08 peter 99   test::Suite suite(argc, argv);
261 09 Mar 05 peter 100
1251 03 Apr 08 peter 101   std::ifstream s(test::filename("data/isoform.peaks").c_str());
261 09 Mar 05 peter 102   std::vector<std::vector<double> > peaksets;
261 09 Mar 05 peter 103   s >> peaksets;
261 09 Mar 05 peter 104
2853 23 Sep 12 peter 105   utility::Matrix nm_score(peaksets.size(), peaksets.size(), 0);
2853 23 Sep 12 peter 106   utility::Matrix correct(nm_score);
4200 19 Aug 22 peter 107   for (size_t i=0; i<peaksets.size()-1; i++)
261 09 Mar 05 peter 108     for (size_t j=i+1; j<peaksets.size(); j++) {
1121 22 Feb 08 peter 109       utility::Matrix dot_m;
261 09 Mar 05 peter 110       std::vector<std::pair<size_t,size_t> > path;
2853 23 Sep 12 peter 111       nm_score(i,j) = score(peaksets[i], peaksets[j], 1.0, dot_m, path);
2853 23 Sep 12 peter 112       if (i==2 && j==3) {
2853 23 Sep 12 peter 113         for (size_t r=0; r<dot_m.rows(); ++r)
2853 23 Sep 12 peter 114           for (size_t c=0; c<dot_m.columns(); ++c)
2853 23 Sep 12 peter 115             if (dot_m(r,c) > 1e-5)
2853 23 Sep 12 peter 116               std::cerr << r << " " << c << " " << dot_m(r,c) << "\n";
2853 23 Sep 12 peter 117       }
261 09 Mar 05 peter 118     }
2853 23 Sep 12 peter 119   correct(0,1) = 4.947983;
2853 23 Sep 12 peter 120   correct(0,2) = 3.496081;
2853 23 Sep 12 peter 121   correct(0,3) = 4.960331;
2853 23 Sep 12 peter 122   correct(1,2) = 2.486868;
2853 23 Sep 12 peter 123   correct(1,3) = 11.407001;
2853 23 Sep 12 peter 124   correct(2,3) = 2.503420;
261 09 Mar 05 peter 125
2853 23 Sep 12 peter 126   suite.add(suite.equal_range_fix(nm_score.begin(), nm_score.end(),
2853 23 Sep 12 peter 127                                   correct.begin(), 1e-6));
2853 23 Sep 12 peter 128
875 19 Sep 07 peter 129   std::string a("AGGUUGUCCGUGGUGAGUUCGCA");
875 19 Sep 07 peter 130   std::string b("GAGGUUGUCCGUGGUGAGUUCG");
1121 22 Feb 08 peter 131   utility::Matrix m(a.size(), b.size());
875 19 Sep 07 peter 132   for (size_t j=0; j<a.size(); ++j)
875 19 Sep 07 peter 133     for (size_t k=0; k<b.size(); ++k)
875 19 Sep 07 peter 134       m(j,k) = (a[j]==b[k] ? 1 : -1000);
2815 28 Aug 12 peter 135   double score = utility::Aligner(100, 100).smith_waterman(m);
875 19 Sep 07 peter 136   if (score!=21)
1232 15 Mar 08 peter 137     suite.add(false);
875 19 Sep 07 peter 138
867 13 Sep 07 peter 139   // testing ssearch
869 14 Sep 07 peter 140   if (utility::ssearch("Hello", "Hll", 0.0, 1.0)!=2){
4200 19 Aug 22 peter 141     suite.err() << "aligning 'Hello' and 'Hll' gives score "
1231 14 Mar 08 peter 142                 << utility::ssearch("Hello", "Hll", 0.0, 1.0)
1231 14 Mar 08 peter 143                 << " expected " << 2 << std::endl;
1232 15 Mar 08 peter 144     suite.add(false);
869 14 Sep 07 peter 145   }
869 14 Sep 07 peter 146   if (utility::ssearch("Hello", "Peter said you can't say 'allo", 1, 1)!=3)
1232 15 Mar 08 peter 147     suite.add(false);
870 14 Sep 07 peter 148
1231 14 Mar 08 peter 149   return suite.return_value();
261 09 Mar 05 peter 150 }

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