This file is indexed.

/usr/include/libphylo/likeDistPropEB.h is in rate4site 3.0.0-5.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
// $Id: likeDistProp.h 962 2006-11-07 15:13:34Z privmane $

#ifndef ___LIKE_DIST_PROP_EB
#define ___LIKE_DIST_PROP_EB

#include "definitions.h"
#include "countTableComponent.h"
#include "multipleStochasticProcess.h"
#include "gammaDistribution.h"
#include "logFile.h"
#include <cmath>

class likeDistPropEB {
private:
	multipleStochasticProcess * _msp;
	const gammaDistribution* _pProportionDist;
	const MDOUBLE _maxPairwiseDistance;
	const MDOUBLE _minPairwiseDistance;
	const MDOUBLE _toll;
public:
	const MDOUBLE giveDistance(	const vector< vector<countTableComponentGamProportional> >& ctc,const int nodeID,
								MDOUBLE& resL,const MDOUBLE initialGuess= 0.03) const;
	explicit likeDistPropEB(multipleStochasticProcess * msp,
							const gammaDistribution* pProportionDist,
							const MDOUBLE toll =0.0001,
							const MDOUBLE maxPairwiseDistance = 5.0,
							const MDOUBLE minPairwiseDistance = 0.0000001) : _msp(msp) ,_pProportionDist(pProportionDist), _maxPairwiseDistance(maxPairwiseDistance), _minPairwiseDistance(minPairwiseDistance),_toll(toll){
	}
	likeDistPropEB(const likeDistPropEB & other)
		: _msp(other._msp),_pProportionDist(other._pProportionDist),_maxPairwiseDistance(other._maxPairwiseDistance),_minPairwiseDistance(other._minPairwiseDistance),_toll(other._toll){}
	virtual likeDistPropEB* clone() const {return new likeDistPropEB(*this);}
};



class C_evallikeDistPropEB_d{ // derivative.
public:
  C_evallikeDistPropEB_d(const vector< vector<countTableComponentGamProportional> >& ctc,
				 multipleStochasticProcess* msp,const gammaDistribution* pProportionDist,const int nodeID) : _ctc(ctc), _msp(msp), _pProportionDist(pProportionDist), _nodeID(nodeID) {};
private:
	const vector< vector<countTableComponentGamProportional> >& _ctc;
	multipleStochasticProcess* _msp;
	const gammaDistribution* _pProportionDist;
	const int _nodeID;
public:
	MDOUBLE operator() (MDOUBLE dist) {
		const MDOUBLE epsilonPIJ = 1e-10;
		MDOUBLE sumDL = 0.0;
		for (int gene=0; gene < _msp->getSPVecSize(); ++gene) {
			for (int alph1=0; alph1 <  _ctc[gene][_nodeID].alphabetSize(); ++alph1){
				for (int alph2=0; alph2 <  _ctc[gene][_nodeID].alphabetSize(); ++alph2){
					for(int globalRateCategor = 0;globalRateCategor < _pProportionDist->categories();++globalRateCategor){
						_msp->getSp(gene)->setGlobalRate(_pProportionDist->rates(globalRateCategor));
						MDOUBLE globalRate = _pProportionDist->rates(globalRateCategor);
						for (int localRateCategor = 0; localRateCategor < _msp->getSp(gene)->categories(); ++localRateCategor) {
							MDOUBLE localRate = _msp->getSp(gene)->rates(localRateCategor);
							MDOUBLE pij= _msp->getSp(gene)->Pij_t(alph1,alph2,dist*globalRate*localRate);
							if (pij<epsilonPIJ) {
								pij = epsilonPIJ;
							}
							MDOUBLE dpij = _msp->getSp(gene)->dPij_dt(alph1,alph2,dist*globalRate*localRate);
							
							//sumDL+= _ctc[gene][_nodeID].getCounts(alph1,alph2,globalRateCategor,localRateCategor)*dpij*_pProportionDist->ratesProb(globalRateCategor)*sp->ratesProb(localRateCategor)
							//			*globalRate*localRate/pij;
							sumDL+= _ctc[gene][_nodeID].getCounts(alph1,alph2,globalRateCategor,localRateCategor)*dpij*globalRate*localRate/pij;
						}
					}
				}
			}
		}
		LOG(12,<<"check bl="<<dist<<" gives sumDL "<<sumDL<<endl);
		return -sumDL;
	};
};



class C_evallikeDistPropEB{
private:
	const vector< vector<countTableComponentGamProportional> >& _ctc;
	multipleStochasticProcess* _msp;
	const gammaDistribution* _pProportionDist;
	const int _nodeID;
public:
	C_evallikeDistPropEB(const vector< vector<countTableComponentGamProportional> >& ctc,
					multipleStochasticProcess* msp,const gammaDistribution* pProportionDist,const int nodeID):_ctc(ctc), _msp(msp), _pProportionDist(pProportionDist), _nodeID(nodeID) {};

	MDOUBLE operator() (MDOUBLE dist) {
		const MDOUBLE epsilonPIJ = 1e-10;
		MDOUBLE sumL = 0.0;
		for (int gene=0; gene < _msp->getSPVecSize(); ++gene) {
			for (int alph1=0; alph1 < _ctc[gene][_nodeID].alphabetSize(); ++alph1){
				for (int alph2=0; alph2 <  _ctc[gene][_nodeID].alphabetSize(); ++alph2){
					for(int globalRateCategor = 0;globalRateCategor < _pProportionDist->categories();++globalRateCategor){
						_msp->getSp(gene)->setGlobalRate(_pProportionDist->rates(globalRateCategor));
						MDOUBLE globalRate = _pProportionDist->rates(globalRateCategor);
						for (int localRateCategor = 0; localRateCategor < _msp->getSp(gene)->categories(); ++localRateCategor) {
							MDOUBLE localRate = _msp->getSp(gene)->rates(localRateCategor);
							MDOUBLE pij= _msp->getSp(gene)->Pij_t(alph1,alph2,dist*globalRate*localRate);
							if (pij<epsilonPIJ) {
								pij = epsilonPIJ;
							}
							sumL += _ctc[gene][_nodeID].getCounts(alph1,alph2,globalRateCategor,localRateCategor)*(log(pij)-log(_msp->getSp(gene)->freq(alph2)));//*_pProportionDist->ratesProb(globalRateCategor)*sp->ratesProb(localRateCategor);
						}
					}
				}
			}
		}
		LOG(8,<<"check bl="<<dist<<" gives sumL "<<sumL<<endl);
		return -sumL;
	};
};

#endif