source: trunk/src/Merge.cc @ 106

Last change on this file since 106 was 106, checked in by cecilia, 18 years ago

handling missing values, if missing in all positions
:wq

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