/usr/include/libphylo/likeDistPropEB.h is in rate4site 3.0.0-2.
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
|