source: trunk/test/test_alignment.cc @ 261

Last change on this file since 261 was 261, checked in by Peter, 17 years ago

fixed bug and added test for alignment

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 1.9 KB
Line 
1// $Id: test_alignment.cc 261 2005-03-09 17:06:19Z peter $
2
3#include "Alignment.h"
4
5#include "matrix.h"
6
7#include <gsl/gsl_cdf.h>
8
9#include <cmath>
10#include <fstream>
11#include <iostream>
12#include <sstream>
13#include <string>
14#include <vector>
15
16void read_vector(std::vector<double>& v,std::string& str)
17{
18  std::istringstream s(str);
19  while (!s.eof()) {
20    double d;
21    s >> d;
22    if (!s.fail())
23      v.push_back(d);
24  }
25}
26
27std::istream& operator>>(std::istream& s,
28                         std::vector<std::vector<double> >& m)
29{
30  while (!s.eof()) {
31    std::string row;
32    char ch;
33    while (((ch=s.get())!='\n') && s.good())
34      row+=ch;
35    if (row.size() || s.good()) { // if stream good, read an empty peak set
36      std::vector<double> v;
37      read_vector(v,row);
38      m.push_back(v);
39    }
40  }
41  return s;
42}
43
44double match(const double x, const double y, const double s)
45{ 
46  return 2*gsl_cdf_gaussian_Q(fabs(x-y),sqrt(2)*s); 
47}
48 
49
50double score(const std::vector<double>& l1,
51             const std::vector<double>& l2,
52             const double sigma, 
53             theplu::gslapi::matrix& dot_matrix,
54             std::vector<std::pair<size_t,size_t> >& path)
55{
56  dot_matrix  = theplu::gslapi::matrix(l1.size(),l2.size());
57  for (size_t i=0; i<l1.size(); i++) 
58    for (size_t j=0; j<l2.size(); j++) {
59      assert(i<dot_matrix.rows());
60      assert(j<dot_matrix.columns());
61      dot_matrix(i,j)=match(l1[i],l2[j],sigma);
62    }
63  return theplu::cpptools::alignment::NeedlemanWunsch(dot_matrix,path,0);
64}
65
66
67int main()
68{
69  bool ok=true;
70
71  std::ifstream s(std::string("data/isoform.peaks").c_str());
72  std::vector<std::vector<double> > peaksets;
73  s >> peaksets;
74
75  for (size_t i=0; i<peaksets.size()-1; i++) 
76    for (size_t j=i+1; j<peaksets.size(); j++) {
77      theplu::gslapi::matrix dot_m;
78      std::vector<std::pair<size_t,size_t> > path;
79      score(peaksets[i], peaksets[j], 1.0, dot_m, path);
80    }
81
82
83  if (ok) {
84    std::cout << "test_alignment: SUCCESS\n";
85    return 0;
86  }
87  std::cerr << "test_alignment: FAILED\n";
88  return -1;
89
90}
Note: See TracBrowser for help on using the repository browser.