raredisplay.cpp

   1 /*
   2  *  raredisplay.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 11/18/08.
   6  *  Copyright 2008 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "raredisplay.h"
  11
  12 /***********************************************************************/
  13
  14 void RareDisplay::init(string label){
  15         try {
  16                 this->label = label;
  17         }
  18         catch(exception& e) {
  19                 m->errorOut(e, "RareDisplay", "init");
  20                 exit(1);
  21         }
  22 }
  23
  24 /***********************************************************************/
  25
  26 void RareDisplay::update(SAbundVector* rank){
  27         try {
  28                 int newNSeqs = rank->getNumSeqs();
  29                 vector<double> data = estimate->getValues(rank);
  30
  31                 if(nIters != 1){
  32
  33                         double oldS = var[index] * ( nIters - 2 );
  34                         double delta = data[0] - results[index];
  35                         double newMean = results[index] + delta / nIters;
  36                         double newS = oldS + delta * ( data[0] - newMean );
  37                         double newVar = newS / ( nIters - 1 );
  38
  39                         seqs[index] = newNSeqs;
  40                         results[index] = newMean;
  41                         var[index] = newVar;
  42
  43                         index++;
  44                 }
  45                 else{
  46                         seqs.push_back(newNSeqs);
  47                         results.push_back(data[0]);
  48                         var.push_back(0.0);
  49                 }
  50         }
  51         catch(exception& e) {
  52                 m->errorOut(e, "RareDisplay", "update");
  53                 exit(1);
  54         }
  55 }
  56
  57 /***********************************************************************/
  58 void RareDisplay::update(vector<SharedRAbundVector*> shared, int numSeqs, int numGroupComb) {
  59         try {
  60                 vector<double> data = estimate->getValues(shared);
  61                 double newNSeqs = data[0];
  62
  63                 if(nIters != 1){
  64
  65                         double oldS = var[index] * ( nIters - 2 );
  66                         double delta = data[0] - results[index];
  67                         double newMean = results[index] + delta / nIters;
  68                         double newS = oldS + delta * ( data[0] - newMean );
  69                         double newVar = newS / ( nIters - 1 );
  70                         seqs[index] = newNSeqs;
  71                         results[index] = newMean;
  72                         var[index] = newVar;
  73
  74                         index++;
  75                 }
  76                 else{
  77
  78                         seqs.push_back(newNSeqs);
  79                         results.push_back(data[0]);
  80                         var.push_back(0.0);
  81
  82                 }
  83         }
  84         catch(exception& e) {
  85                 m->errorOut(e, "RareDisplay", "update");
  86                 exit(1);
  87         }
  88 }
  89
  90 /***********************************************************************/
  91
  92 void RareDisplay::reset(){
  93         try {
  94                 nIters++;
  95                 index = 0;
  96         }
  97         catch(exception& e) {
  98                 m->errorOut(e, "RareDisplay", "reset");
  99                 exit(1);
 100         }
 101 }
 102
 103 /***********************************************************************/
 104
 105 void RareDisplay::close(){
 106         try {
 107
 108                 output->initFile(label);
 109
 110                 for (int i = 0; i < seqs.size(); i++) {
 111
 112                         vector<double> data(3,0);
 113                         double variance = var[i];
 114
 115                         data[0] = results[i];
 116
 117                         double ci = 1.96 * pow(variance, 0.5);
 118                         data[1] = data[0] - ci;
 119                         data[2] = data[0] + ci;
 120
 121                         output->output(seqs[i], data);
 122                 }
 123
 124                 nIters = 1;
 125                 index = 0;
 126
 127                 seqs.clear();
 128                 results.clear();
 129                 var.clear();
 130
 131                 output->resetFile();
 132         }
 133         catch(exception& e) {
 134                 m->errorOut(e, "RareDisplay", "close");
 135                 exit(1);
 136         }
 137 }
 138 /***********************************************************************/
 139
 140 void RareDisplay::inputTempFiles(string filename){
 141         try {
 142                 ifstream in;
 143                 m->openInputFile(filename, in);
 144
 145                 int thisIters;
 146                 in >> thisIters; m->gobble(in);
 147
 148                 for (int i = 0; i < seqs.size(); i++) {
 149                         double tempresult, tempvar;
 150                         in >> tempresult >> tempvar; m->gobble(in);
 151
 152                         //find weighted result
 153                         results[i] = ((nIters * results[i]) + (thisIters * tempresult)) / (float)(nIters + thisIters);
 154
 155                         var[i] = ((nIters * var[i]) + (thisIters * tempvar)) / (float)(nIters + thisIters);
 156                 }
 157
 158                 in.close();
 159         }
 160         catch(exception& e) {
 161                 m->errorOut(e, "RareDisplay", "inputTempFiles");
 162                 exit(1);
 163         }
 164 }
 165
 166 /***********************************************************************/
 167
 168 void RareDisplay::outputTempFiles(string filename){
 169         try {
 170                 ofstream out;
 171                 m->openOutputFile(filename, out);
 172
 173                 out << nIters << endl;
 174
 175                 for (int i = 0; i < seqs.size(); i++) {
 176                         out << results[i] << '\t' << var[i] << endl;
 177                 }
 178
 179                 out.close();
 180         }
 181         catch(exception& e) {
 182                 m->errorOut(e, "RareDisplay", "outputTempFiles");
 183                 exit(1);
 184         }
 185 }
 186
 187
 188 /***********************************************************************/
 189