source: trunk/src/Merge.cc @ 83

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

added test output

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.0 KB
Line 
1// $Id: Merge.cc 83 2004-05-26 14:42:25Z 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 annot
49    for (map<string,vector<u_int> >::const_iterator k=name.begin(); k!=name.end(); k++)
50      merging((*k).second, weights);
51  }
52
53  void Merge::merging(vector<u_int> index, gslapi::matrix& weights){
54
55    // merge genes for each assay
56    for (u_int j=0; j<nof_assays_; j++){
57
58      // calc. m
59      double sum_wx = 0;
60      double sum_w = 0;
61      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
62        sum_wx += weights(*i,j)*data_(*i,j);
63        sum_w += weights(*i,j);
64      }
65      double m = sum_wx/sum_w;
66     
67      // calc. U
68      double sum_u = 0;
69      double sum_var = 0;
70      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
71        sum_u += 1/error_(*i,j);
72        double w = weights(*i,j);
73        double diff = data_(*i,j)-m;
74        sum_var += w*w*diff*diff;
75      }
76      double U = 1/sum_u+sum_var/(sum_w*sum_w);
77     
78      // insert new values
79      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
80        data_(*i,j) = m;
81        error_(*i,j) = U;
82      }
83     
84    }
85   
86  }
87
88}} // of namespace cpptools and namespace theplu
Note: See TracBrowser for help on using the repository browser.