source: trunk/lib/utility/Merge.cc @ 301

Last change on this file since 301 was 301, checked in by Peter, 18 years ago

modified includes in tests

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.9 KB
Line 
1// $Id: Merge.cc 301 2005-04-30 13:39:27Z peter $
2
3
4#include <c++_tools/utility/Merge.h>
5
6#include <c++_tools/statistics/Averager.h>
7#include <c++_tools/statistics/AveragerWeighted.h>
8#include <c++_tools/utility/stl_utility.h>
9
10#include <algorithm>
11#include <cmath>
12#include <fstream>
13#include <map>
14
15namespace theplu {
16namespace utility {
17
18  using namespace std;
19
20  Merge::Merge(gslapi::matrix& data, gslapi::matrix& error, 
21               vector<string>& annotations)
22    : data_(data), error_(error), annots_(annotations)
23  {
24   
25    if(data_.columns() == error_.columns()){
26      nof_assays_ = data_.columns();
27    }else{
28      cerr << "not the same number of assays in data and error matrix";
29    }
30   
31    if(data_.rows() != error_.rows() || data_.rows() != annots_.size()){
32      cerr << "wrong number of gene annotations (!= nof genes)";
33    }
34  }
35
36  Merge::Merge(istream& data, istream& error, istream& annotations)
37    : data_(data), error_(error)
38  {
39    // skapa annot-vektor
40    while(!annotations.eof()){
41      string str;
42      getline(annotations,str,'\n');
43      annots_.push_back(str);
44    }
45   
46    if(data_.columns() == error_.columns()){
47      nof_assays_ = data_.columns();
48    }else{
49      cerr << "not the same number of assays in data and error matrix";
50    }
51   
52    if(data_.rows() != error_.rows() || data_.rows() != annots_.size()){
53      cerr << "wrong number of gene annotations (!= nof genes)";
54    }
55  }
56 
57  void Merge::calculate_weights(gslapi::matrix& weights) {
58   
59    for (unsigned int i=0; i<weights.rows(); i++)
60      for (unsigned int j=0; j<weights.columns(); j++){
61        // if error < 0 (missing value) weight=0
62        if(weights(i,j) >= 0){
63          weights(i,j) = 1/(1+weights(i,j));
64        }else{
65          weights(i,j) = 0;
66        }
67      }
68  }
69
70  void Merge::do_something(void){
71
72    using namespace std;
73    gslapi::matrix weights(error_);
74    calculate_weights(weights);
75
76    // för varje annot: kolla vilka rader den finns på
77    map<string, vector<u_int> > name;
78    for (u_int i=0; i<annots_.size(); i++){
79      name[annots_[i]].push_back(i);
80    }
81   
82    // för varje annot, merge
83    for (map<string,vector<u_int> >::const_iterator k=name.begin(); 
84         k!=name.end(); k++){
85      if((*k).second.size()>1){
86        merging((*k).second, weights);
87      }
88    }
89  }
90
91  void Merge::merging(vector<u_int> index, gslapi::matrix& weights){
92    statistics::AveragerWeighted w_av;
93    statistics::Averager av;
94    // merge genes for each assay
95    for (u_int j=0; j<nof_assays_; j++){
96      // calc. m
97      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
98        w_av.add(data_(*i,j),weights(*i,j));
99        if(error_(*i,j)>0){
100          av.add(1/error_(*i,j));
101        }
102      }
103      double m = w_av.mean();
104      // calc. U
105      double U = -1;
106      if(av.n()>0){ // if av not empty, else missing
107        U = 1/av.mean()+w_av.squared_error();
108      }
109      // insert new values
110      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
111        data_(*i,j) = m;
112        error_(*i,j) = U;
113      }
114      w_av.reset();
115      av.reset();
116    }
117  }
118
119}} // of namespace utility and namespace theplu
Note: See TracBrowser for help on using the repository browser.