source: trunk/src/Merge.cc @ 83

// $Id: Merge.cc 83 2004-05-26 14:42:25Z cecilia $

#include <algorithm>
#include <cmath>
#include <fstream>
#include <map>

#include "stl_utility.h"
#include "Merge.h"

namespace theplu {
namespace cpptools {

  using namespace std;

  Merge::Merge(istream& data, istream& error, istream& annotations, int nof_assays)
    : data_(data), error_(error), nof_assays_(nof_assays)
  {
    // skapa annot-vektor
    while(!annotations.eof()){
      string str;
      getline(annotations,str,'\n');
      annots_.push_back(str);
    }
  }

  gslapi::matrix Merge::calculate_weights(void) const
  {
    gslapi::matrix weights(error_);
    for (unsigned int i=0; i<weights.rows(); i++)
      for (unsigned int j=0; j<weights.columns(); j++)
        weights(i,j) = 1/(1+weights(i,j));
    // weights = w_i=1/(1+u_i)
    return weights;
  }

  void Merge::do_something(void)
  {
    using namespace std;
    gslapi::matrix weights;
    weights = calculate_weights();

    // för varje annot: kolla vilka rader den finns på
    map<string, vector<u_int> > name;
    for (u_int i=0; i<annots_.size(); i++)
      name[annots_[i]].push_back(i);
    
    // för varje annot
    for (map<string,vector<u_int> >::const_iterator k=name.begin(); k!=name.end(); k++)
      merging((*k).second, weights);
  }

  void Merge::merging(vector<u_int> index, gslapi::matrix& weights){

    // merge genes for each assay
    for (u_int j=0; j<nof_assays_; j++){

      // calc. m
      double sum_wx = 0;
      double sum_w = 0;
      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
        sum_wx += weights(*i,j)*data_(*i,j);
        sum_w += weights(*i,j);
      }
      double m = sum_wx/sum_w;
      
      // calc. U
      double sum_u = 0;
      double sum_var = 0;
      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
        sum_u += 1/error_(*i,j);
        double w = weights(*i,j);
        double diff = data_(*i,j)-m;
        sum_var += w*w*diff*diff;
      }
      double U = 1/sum_u+sum_var/(sum_w*sum_w);
      
      // insert new values
      for (vector<u_int>::const_iterator i=index.begin(); i!=index.end(); i++){
        data_(*i,j) = m;
        error_(*i,j) = U;
      }
      
    }
    
  }

}} // of namespace cpptools and namespace theplu