source: trunk/src/Merge.cc @ 81

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

merge object + test program + data files

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 2.7 KB
Line 
1// $Id: Merge.cc 81 2004-05-26 13:30:01Z 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 (string?)
20    while(annotations){
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    // för varje annot: kolla vilka rader den finns på
43    map<string, vector<u_int> > name;
44    cout << annots_.size() << " size" << endl;
45    for (u_int i=0; i<annots_.size(); i++) {
46      map<string,vector<u_int> >::iterator p;
47      p = name.find(annots_[i]);
48      vector<u_int> index;
49      if(p != name.end()){
50        index = (*p).second;
51      }
52      index.push_back(i); 
53      for (vector<u_int>::const_iterator j=index.begin(); j!=index.end(); j++){
54        cout << *j << "* ";
55      }
56      cout << endl;
57      name.insert(pair<string, vector<u_int> >(annots_[i], index));
58    }
59   
60    // för varje namn:
61   
62    for (map<string,vector<u_int> >::const_iterator k=name.begin(); k!=name.end(); k++) {
63//      merging((*k).second, weights);
64      vector<u_int> t = (*k).second;
65      for (vector<u_int>::const_iterator i=t.begin(); i!=t.end(); i++){
66        cout << *i << " ";
67      }
68      cout << endl;
69    }
70  }
71
72  void Merge::merging(vector<u_int> index, gslapi::matrix& weights){
73    // nya rader i error_ och data_
74   
75    vector<double> new_data;
76    vector<double> new_error;
77    // merge genes for each assay
78    for (u_int j=0; j<nof_assays_; j++){
79
80      // calc. m
81      double sum_wx = 0;
82      double sum_w = 0;
83      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
84        sum_wx += weights(*i,j)*data_(*i,j);
85        sum_w += weights(*i,j);
86      }
87      double m = sum_wx/sum_w;
88      new_data.push_back(m);
89     
90      // calc. U
91      double sum_u = 0;
92      double sum_var = 0;
93      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
94        sum_u += 1/error_(*i,j);
95        double w = weights(*i,j);
96        double diff = data_(*i,j)-m;
97        sum_var += w*w*diff*diff;
98      }
99      double U = 1/sum_u+sum_var/(sum_w*sum_w);
100      new_error.push_back(U);
101     
102      // insert new values
103      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
104        data_(*i,j) = m;
105        error_(*i,j) = U;
106      }
107     
108    }
109   
110  }
111
112}} // of namespace cpptools and namespace theplu
Note: See TracBrowser for help on using the repository browser.