source: trunk/src/Merge.cc @ 82

Last change on this file since 82 was 82, checked in by cecilia, 19 years ago

annot-map

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.2 KB
Line 
1// $Id: Merge.cc 82 2004-05-26 14:16:03Z 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
11namespace theplu {
12namespace cpptools {
13
14  using namespace std;
15
16  Merge::Merge(istream& data, istream& error, istream& annotations, int nof_assays)
17    : data_(data), error_(error), nof_assays_(nof_assays)
18  {
19    // skapa annot-vektor
20    while(!annotations.eof()){
21      string str;
22      getline(annotations,str,'\n');
23      annots_.push_back(str);
24    }
25  }
26
27  gslapi::matrix Merge::calculate_weights(void) const
28  {
29    gslapi::matrix weights(error_);
30    for (unsigned int i=0; i<weights.rows(); i++)
31      for (unsigned int j=0; j<weights.columns(); j++)
32        weights(i,j) = 1/(1+weights(i,j));
33    // weights = w_i=1/(1+u_i)
34    return weights;
35  }
36
37  void Merge::do_something(void)
38  {
39    using namespace std;
40    gslapi::matrix weights;
41    weights = calculate_weights();
42
43    // för varje annot: kolla vilka rader den finns på
44    map<string, vector<u_int> > name;
45    for (u_int i=0; i<annots_.size(); i++)
46      name[annots_[i]].push_back(i);
47   
48    // för varje namn:
49   
50    for (map<string,vector<u_int> >::const_iterator k=name.begin(); k!=name.end(); k++)
51      merging((*k).second, weights);
52  }
53
54  void Merge::merging(vector<u_int> index, gslapi::matrix& weights){
55    // nya rader i error_ och data_
56   
57    vector<double> new_data;
58    vector<double> new_error;
59    // merge genes for each assay
60    for (u_int j=0; j<nof_assays_; j++){
61
62      // calc. m
63      double sum_wx = 0;
64      double sum_w = 0;
65      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
66        sum_wx += weights(*i,j)*data_(*i,j);
67        sum_w += weights(*i,j);
68      }
69      double m = sum_wx/sum_w;
70      new_data.push_back(m);
71     
72      // calc. U
73      double sum_u = 0;
74      double sum_var = 0;
75      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
76        sum_u += 1/error_(*i,j);
77        double w = weights(*i,j);
78        double diff = data_(*i,j)-m;
79        sum_var += w*w*diff*diff;
80      }
81      double U = 1/sum_u+sum_var/(sum_w*sum_w);
82      new_error.push_back(U);
83     
84      // insert new values
85      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
86        data_(*i,j) = m;
87        error_(*i,j) = U;
88      }
89     
90    }
91   
92  }
93
94}} // of namespace cpptools and namespace theplu
Note: See TracBrowser for help on using the repository browser.