#include"Prune.h"
/*****************************************************************/
// Polar Method (Box, Muller and Marsaglia) for two independent
// normal distributions
/*****************************************************************/
twodouble Prune :: normal (double mu1, double var1, double mu2,
double var2)
{
double v1, v2, s;
do {
v1 = 2 * drand48() - 1;
v2 = 2 * drand48() - 1;
s = v1*v1 + v2*v2;
}
while (s >= 1 || s < 1e-30);
s = sqrt((-2*log(s))/s);
v1 *= s;
v2 *= s;
twodouble x;
x.x1=v1*sqrt(var1)+mu1;
x.x2=v2*sqrt(var2)+mu2;
return (x);
}
/*****************************************************************/
// Taylor Series expansion to find area under Normal distribution
// P(z\le x)=0.5 + {{1}\over{\sqrt{2\pi}}} \sum_{k=0}^\infty
// {{(-1)^k x^{2*k+1}} \over {(2k+1) 2^k k!}}
// z is standard normal.
/*****************************************************************/
double Prune :: pnorm(double x) // Gaussian CDF
{
if(x<-6.5) return 0;
if(x>6.5) return 1;
double factK=1;
double sum=0.0;
double term=1.0;
int k=0;
while(fabs(term)>exp(-23)){
term=(1.0/sqrt(2.0*M_PI))*pow(-1,k)*pow(x,2*k+1)/
((2*k+1)*pow(2,k)*factK);
sum+=term;
k++;
factK*=k;
}
sum+=0.5;
if(sum<1e-10) sum=0;
return sum;
}
/*****************************************************************/
// Returns pdf of (a+x)/(b+y), where x and y are independent
// standard normal distributed.
/*****************************************************************/
double Prune :: ratio(double a, double b, double t)
{
double q,tf;
q=(b+a*t)/sqrt(1+t*t);
tf=(exp(-0.5*(a*a+b*b))/(M_PI*(1+t*t)))*(1+(q/dnorm(q))*
(pnorm(q)-pnorm(0)));
return tf;
}
/*****************************************************************/
double Prune :: dnorm(double x) // Gaussian PDF
{
return ((1.0/sqrt(2.0*M_PI))*exp(-x*x/2));
}
/*****************************************************************/
// The actual values of pruning are generated here
/*****************************************************************/
void Prune :: setZeta(char * basefilename)
{
int i,k;
struct timeval tv;
struct timezone tz;
gettimeofday(&tv,&tz);
int SEED=tv.tv_usec;
srand48(SEED);
char * filename=new char[35];
sprintf(filename,"%s.cb",basefilename);
ifstream cbfile (filename); // Control, BP
sprintf(filename,"%s.cd",basefilename);
ifstream cdfile (filename); // Control, DL
sprintf(filename,"%s.sb",basefilename);
ifstream sbfile (filename); // Stress, BP
sprintf(filename,"%s.sd",basefilename);
ifstream sdfile (filename); // Stress, DL
if(!cbfile || !cdfile || !sbfile || !sdfile){
cerr << "\nOne or More of the Stats files not found.\n";
exit(1);
}
// ************* CONVENTION *************/
// xyzp; x \in {c,s} - control and stress
// y \in {b,d} - B.P. or D.L.
// z \in {a,b} - apical or basal
// p \in {m,v} - mean or variance
// *************************************/
double * cdam; double * cdav; double * cdbm; double * cdbv;
double * cbam; double * cbav; double * cbbm; double * cbbv;
double * sdam; double * sdav; double * sdbm; double * sdbv;
double * sbam; double * sbav; double * sbbm; double * sbbv;
cdam=new double[antra]; cdav=new double[antra];
cbam=new double[antra]; cbav=new double[antra];
cdbm=new double[bntra]; cdbv=new double[bntra];
cbbm=new double[bntra]; cbbv=new double[bntra];
sdam=new double[antra]; sdav=new double[antra];
sbam=new double[antra]; sbav=new double[antra];
sdbm=new double[bntra]; sdbv=new double[bntra];
sbbm=new double[bntra]; sbbv=new double[bntra];
for (i=0; i<antra; i++){ // Apical means
cdfile >> cdam[i]; cbfile >> cbam[i];
sdfile >> sdam[i]; sbfile >> sbam[i];
}
for (i=0; i<antra; i++){ // Apical variances
cdfile >> cdav[i]; cbfile >> cbav[i];
sdfile >> sdav[i]; sbfile >> sbav[i];
}
for (i=0; i<bntra; i++){ // Basal means
cdfile >> cdbm[i]; cbfile >> cbbm[i];
sdfile >> sdbm[i]; sbfile >> sbbm[i];
}
for (i=0; i<bntra; i++){ // Basal variances
cdfile >> cdbv[i]; cbfile >> cbbv[i];
sdfile >> sdbv[i]; sbfile >> sbbv[i];
}
cdfile.close(); cbfile.close();
sdfile.close(); sbfile.close();
double a,b,max;
double bpsum;
double dlcnt, bpcnt;
dlsum=bpsum=0.0;
dlcnt=bpcnt=0.0;
abpzeta=new double[sdata.na];
bbpzeta=new double[sdata.nb];
adlzeta=new double[sdata.na];
bdlzeta=new double[sdata.nb];
// Find maximum ratio on Apical side for D.L. and B.P. //
// Generate Zeta on Apical Side for DL and BP //
for(i=0; i<sdata.na; i++){
if(i>(antra-1)) k=antra-1;
else k=i;
max=maxratio(sdam[k], sdav[k], cdam[k], cdav[k]);
adlzeta[i]=generateZeta(sdam[k], sdav[k], cdam[k], cdav[k], max);
adlzeta[i] *= sdata.apicalcount[i];
max=maxratio(sbam[k], sbav[k], cbam[k], cbav[k]);
abpzeta[i]=generateZeta(sbam[k], sbav[k], cbam[k], cbav[k], max);
abpzeta[i] *= BPStats.apicalcount[i];
dlsum += adlzeta[i] ;
bpsum += abpzeta[i];
dlcnt += sdata.apicalcount[i];
bpcnt += BPStats.apicalcount[i];
}
// Find maximum ratio on Basal side for DL and BP //
// Generate Zeta on Basal Side for DL and BP //
for(i=0; i<sdata.nb; i++){
if(i>(bntra-1)) k=bntra-1;
else k=i;
max=maxratio(sdbm[k], sdbv[k], cdbm[k], cdbv[k]);
bdlzeta[i]=generateZeta(sdbm[k], sdbv[k], cdbm[k], cdbv[k], max);
bdlzeta[i] *= sdata.basalcount[i];
max=maxratio(sbbm[k], sbbv[k], cbbm[k], cbbv[k]);
bbpzeta[i]=generateZeta(sbbm[k], sbbv[k], cbbm[k], cbbv[k], max);
bbpzeta[i] *= BPStats.basalcount[i];
dlsum += bdlzeta[i] ;
bpsum += bbpzeta[i];
dlcnt += sdata.basalcount[i];
bpcnt += BPStats.basalcount[i];
}
/*
cerr << "\nReduction in DL: " << dlsum ;
cerr << "\nTotal DL: " << dlcnt ;
cerr << "\nReduction in BP: " << bpsum ;
cerr << "\nTotal BP: " << bpcnt ;
cerr << endl ;
*/
delete(cdam); delete(cdav);
delete(cbam); delete(cbav);
delete(cdbm); delete(cdbv);
delete(cbbm); delete(cbbv);
delete(sdam); delete(sdav);
delete(sbam); delete(sbav);
delete(sdbm); delete(sdbv);
delete(sbbm); delete(sbbv);
}
/*****************************************************************/
double Prune :: generateZeta(double sm, double sv, double cm,
double cv, double maxratio)
{
if(!maxratio) return drand48();
double vara=1,varb=1;
double a,b;
if (sv) a=sm/sv; // Avoid division by zero
else if (!sm) return 0; // Variance is zero and mean is zero
else {a=sm; vara=0;} // Variance is zero and mean is nonzero
if (cv) b=cm/cv;
else if (!cm) return 0;
else {b=cm; varb=0;}
twodouble zeta;
double t;
double dzeta;
double min;
if(sm>cm) min=1.0;
else min=sm/cm;
do{
if(!vara && !varb) {
return 1-sm/cm;
}
else if(!vara){
zeta=normal(sm,0,b,1);
t=zeta.x1/zeta.x2;
dzeta=t/cv;
}
else if(!varb){
zeta=normal(a,1,cm,0);
t=zeta.x1/zeta.x2;
dzeta=t*sv;
}
else {
zeta=normal(a,1,b,1); // This is ORIG.
//zeta=normal(sm,sv,cm,cv); // TEST
t=zeta.x1/zeta.x2;
dzeta=(t*sv)/(cv);
}
}
while(ratio(a,b,t)<(0.1*maxratio) || dzeta > 1.0 || dzeta < 0.0);
//while(dzeta > 1.0 || dzeta < (min*0.88));
//while(ratio(a,b,t)<(0.65*maxratio) || dzeta > 1.0 || dzeta < (min*0.85));
return (1-dzeta);
}
/*****************************************************************/
double Prune :: maxratio (double sm, double sv, double cm, double cv)
{
double a,b;
if (sv) a=sm/sv; // Avoid division by zero
else if (!sm) return 0; // Variance is zero and mean is zero
else a=sm; // Variance is zero and mean is nonzero
if (cv) b=cm/cv;
else if (!cm) return 0;
else b=cm;
double max=0.0;
double t;
for(t=-3; t<5; t+=0.01)
if(max<ratio(a,b,t))
max=ratio(a,b,t);
return max;
}
/*****************************************************************/
void Prune :: setProbabilities()
{
setDLProbs();
setBPProbs();
}
/*****************************************************************/
void Prune :: setBPProbs()
{
int i;
double P;
// Apical
for(i=0; i<sdata.na; i++){
if(abpzeta[i]==0 || adlzeta[i]==0)
abpP[i]=0;
else{
if(abpprune[i]==0) abpprune[i]=1;
if(adlprune[i]==0) adlprune[i]=RESLN;
P=(adlzeta[i]*abpprune[i])/(abpzeta[i]*adlprune[i]);
abpP[i]=pow(adlP[i],P);
}
}
// Basal
for(i=0; i<sdata.nb; i++){
if(bbpzeta[i]==0 || bdlzeta[i]==0)
bbpP[i]=0;
else{
if(bbpprune[i]==0) bbpprune[i]=1;
if(bdlprune[i]==0) bdlprune[i]=RESLN;
P=(bdlzeta[i]*bbpprune[i])/(bbpzeta[i]*bdlprune[i]);
bbpP[i]=pow(bdlP[i],P);
}
}
}
/*****************************************************************/
void Prune :: setDLProbs()
{
/* Set Probabilities for DL Pruning by just setting the sum of *
* probabilities to 1 by normalizing (zeta-prune) values. */
double dlSum=0.0;
maxdlP=0;
int i;
for(i=0; i<sdata.na; i++){
adlP[i]=adlzeta[i]-adlprune[i];
dlSum+=adlP[i];
}
for(i=0; i<sdata.nb; i++){
bdlP[i]=bdlzeta[i]-bdlprune[i];
dlSum+=bdlP[i];
}
for(i=0; i<sdata.na; i++){
adlP[i]/=dlSum;
if(maxdlP<adlP[i]) maxdlP=adlP[i];
//cerr << "ADLP: " << i << " " << adlP[i] << endl;
}
for(i=0; i<sdata.nb; i++){
bdlP[i]/=dlSum;
if(maxdlP<bdlP[i]) maxdlP=bdlP[i];
//cerr << "BDLP: " << i << " " << bdlP[i] << endl;
}
maxdlP += 0.1*maxdlP;
}
/*****************************************************************/
