source: trunk/test/alignment_test.cc @ 1659

Last change on this file since 1659 was 1487, checked in by Jari Häkkinen, 13 years ago

Addresses #436. GPL license copy reference should also be updated.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 3.3 KB
Line 
1// $Id: alignment_test.cc 1487 2008-09-10 08:41:36Z jari $
2
3/*
4  Copyright (C) 2005 Jari Häkkinen, Peter Johansson
5  Copyright (C) 2006 Jari Häkkinen
6  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
7
8  This file is part of the yat library, http://dev.thep.lu.se/yat
9
10  The yat library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU General Public License as
12  published by the Free Software Foundation; either version 3 of the
13  License, or (at your option) any later version.
14
15  The yat library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with yat. If not, see <http://www.gnu.org/licenses/>.
22*/
23
24#include "Suite.h"
25
26#include "yat/utility/Alignment.h"
27
28#include "yat/utility/Matrix.h"
29
30#include <gsl/gsl_cdf.h>
31
32#include <cassert>
33#include <cmath>
34#include <fstream>
35#include <iostream>
36#include <sstream>
37#include <string>
38#include <vector>
39
40using namespace theplu::yat;
41
42void read_vector(std::vector<double>& v,std::string& str)
43{
44  std::istringstream s(str);
45  while (!s.eof()) {
46    double d;
47    s >> d;
48    if (!s.fail())
49      v.push_back(d);
50  }
51}
52
53std::istream& operator>>(std::istream& s,
54                         std::vector<std::vector<double> >& m)
55{
56  while (!s.eof()) {
57    std::string row;
58    char ch;
59    while (((ch=s.get())!='\n') && s.good())
60      row+=ch;
61    if (row.size() || s.good()) { // if stream good, read an empty peak set
62      std::vector<double> v;
63      read_vector(v,row);
64      m.push_back(v);
65    }
66  }
67  return s;
68}
69
70double match(const double x, const double y, const double s)
71{ 
72  return 2*gsl_cdf_gaussian_Q(std::abs(x-y),sqrt(2)*s); 
73}
74 
75
76double score(const std::vector<double>& l1,
77             const std::vector<double>& l2,
78             const double sigma, 
79             theplu::yat::utility::Matrix& dot_matrix,
80             std::vector<std::pair<size_t,size_t> >& path)
81{
82  dot_matrix.resize(l1.size(),l2.size());
83  for (size_t i=0; i<l1.size(); i++) 
84    for (size_t j=0; j<l2.size(); j++) {
85      assert(i<dot_matrix.rows());
86      assert(j<dot_matrix.columns());
87      dot_matrix(i,j)=match(l1[i],l2[j],sigma);
88    }
89  return theplu::yat::utility::NeedlemanWunsch(dot_matrix,path,0);
90}
91
92
93int main(int argc, char* argv[])
94{
95  test::Suite suite(argc, argv);
96
97  std::ifstream s(test::filename("data/isoform.peaks").c_str());
98  std::vector<std::vector<double> > peaksets;
99  s >> peaksets;
100
101  for (size_t i=0; i<peaksets.size()-1; i++) 
102    for (size_t j=i+1; j<peaksets.size(); j++) {
103      utility::Matrix dot_m;
104      std::vector<std::pair<size_t,size_t> > path;
105      score(peaksets[i], peaksets[j], 1.0, dot_m, path);
106    }
107
108  std::string a("AGGUUGUCCGUGGUGAGUUCGCA");
109  std::string b("GAGGUUGUCCGUGGUGAGUUCG");
110  utility::Matrix m(a.size(), b.size());
111  for (size_t j=0; j<a.size(); ++j)
112    for (size_t k=0; k<b.size(); ++k)
113      m(j,k) = (a[j]==b[k] ? 1 : -1000);
114  double score=utility::SmithWaterman(m, 100, 100);
115  if (score!=21)
116    suite.add(false);
117
118  // testing ssearch
119  if (utility::ssearch("Hello", "Hll", 0.0, 1.0)!=2){
120    suite.err() << "aligning 'Hello' and 'Hll' gives score " 
121                << utility::ssearch("Hello", "Hll", 0.0, 1.0)
122                << " expected " << 2 << std::endl;
123    suite.add(false);
124  }
125  if (utility::ssearch("Hello", "Peter said you can't say 'allo", 1, 1)!=3)
126    suite.add(false);
127
128  return suite.return_value();
129}
Note: See TracBrowser for help on using the repository browser.