File: Sampler.h

package info (click to toggle)
bitseq 0.7.5+dfsg-5
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, sid
  • size: 676 kB
  • sloc: cpp: 7,043; python: 562; makefile: 150; sh: 52
file content (99 lines) | stat: -rw-r--r-- 3,056 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
#ifndef SAMPLER_H
#define SAMPLER_H

#include<vector>
#include<fstream>
#include "boost/random/mersenne_twister.hpp"
#include "boost/random/gamma_distribution.hpp"
#include "boost/random/uniform_01.hpp"

using namespace std;

#include "GibbsParameters.h"
#include "TagAlignments.h"

// compute statistics
//#define DoSTATS
//#define DoDebug

typedef pair<double,double> pairD;

class Sampler{
   protected:
   long m, samplesN, samplesLogged, samplesTotal, samplesOut, Nmap, Nunmap;
   const distributionParameters *beta,*dir;
   const TagAlignments *alignments;
   const vector<double> *isoformLengths;
   boost::random::mt11213b rng_mt;
   boost::random::gamma_distribution<double> gammaDistribution;
   typedef boost::random::gamma_distribution<double>::param_type gDP;
   // Need by children:
   boost::random::uniform_01<double> uniformDistribution;
   
   bool doLog,save;
   string saveType;
   ofstream *outFile;
   double saveNorm,logRate;
#ifdef DoSTATS   
   long long nT,nZ,nTa;
   double tT,tZ,tTa;
#endif

   vector<long> C;
   double sumC0;
   vector<double> theta;
   vector<double> thetaActLog;
   vector<pairD> thetaSum;
   vector<pairD> thetaSqSum;
   pairD sumNorm;

   // Sample theta.
   void sampleTheta();
   // Compute tau.
   void getTau(vector <double> &tau, double norm);
   // Append current expression samples into file opened for saving samples.
   void appendFile();
   // Update sums of theta and theta^2.
   void updateSums();

   public:

   Sampler();
   virtual ~Sampler();
   // Initialize sampler, set seed and use it to generate new seed.
   void init(long m, long samplesTotal, long samplesOut, long Nunmap,
             const TagAlignments *alignments,
             const distributionParameters &betaPar, 
             const distributionParameters &dirPar, 
             long &seed);
   // Reset sampler's stats before new iteration
   void resetSampler(long samplesTotal);
   // Return mean C[0].
   long getAverageC0();
   // Get vector of mean theta expression. Has "two columns" first is calculated
   // from all samples, the second only from the "thinned" samples.
   void getAverage(vector<pairD> &av);
   // Get mean for transcript i.
   pairD getAverage(long i);
   // Get within variance for transcript i.
   pairD getWithinVariance(long i);
   // Get sum of theta^2, sum of theta, and their norm for transcript i.
   void getThetaSums(long i, double *thSqSum, double *thSum, double *sumN);
   // Return norms for theta sums.
   pairD getSumNorms() const { return sumNorm; }
   // Set sampler into state where samples are saved into the outFile.
   void saveSamples(ofstream *outFile, const vector<double> *isoformLengths,
                    const string &saveType, double norm = 0);
   // Stop saving samples into the file.
   void noSave();
   // Get theta act logged values.
   const vector<double>& getThetaActLog(){return thetaActLog;}

   // Produce new McMc samples.
   virtual void sample();
   // If necessary ("thinned sample") sample theta; update sums.
   virtual void update();
};


#endif