Commit 9b0b7e45 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

fixed memory problem in TF proposal

relaxed dependency on offsource PSD (fingers crossed)


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@617 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 0b0784a0
......@@ -1533,23 +1533,9 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
else
logLx = 1.0;
// logPsx = logprior(priors->invsigma, priors->mean, Snx, 0, ncut);
// if(priorFlag)logPsx = logprior(priors->sigma, priors->mean, Snx, 0, ncut);
if(priorFlag==1)
{
logPsx = logprior(priors->lower, priors->upper, Snx, 0, ncut);
// printf("starting position for PSD: %g\n",logPsx);
// if(logPsx < -1e50)
// {
// FILE *temp = fopen("temp.dat","w");
// for(i=0; i<ncut; i++)
// {
// fprintf(temp,"%lg %lg %lg %lg\n",(double)i/4.,priors->lower[i],Snx[i],priors->upper[i]);
// }
// fclose(temp);
// system_pause();
// }
// logPsx = logprior_gaussian_model(priors->mean, priors->sigma, Snx, spline_x->points, spline_x->n, lines_x->f, lines_x->n, data);
}
if(priorFlag==2)
{
......@@ -2533,8 +2519,8 @@ void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double
}
else
{
psd[i] = bayesline->priors->SAmax;
splinePSD[i] = bayesline->priors->SAmax;
psd[i] = 1.0;
splinePSD[i] = 1.0;
}
invpsd[i] = 1./psd[i];
}
......
......@@ -323,9 +323,16 @@ int main(int argc, char *argv[])
bayesline[0][ifo]->power[i] = (asd[2*j]*asd[2*j]+asd[2*j+1]*asd[2*j+1]);
}
//Quit lying to BayesLine about the data. Poor BayesLine.
for(i=0; i<(imax-imin); i++)
{
j = i+imin;
bayesline[0][ifo]->power[i] = (data->s[ifo][2*j]*data->s[ifo][2*j]+data->s[ifo][2*j+1]*data->s[ifo][2*j+1]);
}
fprintf(stdout,"BayesLine search phase for IFO %i\n", ifo);
BayesLineSearch(bayesline[0][ifo]);
//TODO: Does BayesLineRJMCMC use asd?
fprintf(stdout,"BayesLine characterization phase for IFO %i\n", ifo);
BayesLineRJMCMC(bayesline[0][ifo], asd, model[0]->Snf[ifo], model[0]->invSnf[ifo], model[0]->SnS[ifo], N, 100000, 1.0, 2);
......@@ -713,6 +720,8 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing signal model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZsignal, &data->varZsignal);
//shut off signal model for checkpointing
......@@ -747,6 +756,8 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing glitch model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
/*if(data->glitchOnlyFlag==1 && strcmp(data->channels[0],"AdLIGO"))
......@@ -938,6 +949,7 @@ void print_help_message(void)
fprintf(stdout," --clusterBeta cluster prior exp(-beta K) (4)\n");
fprintf(stdout," --clusterGamma occam penalty exp(gamma*J) (0)\n");
fprintf(stdout," --updateGeocenterPSD geocenter PSD depends on extrinsic parameters\n");
fprintf(stdout," --waveletPrior use empirical distribution on number of wavelets (O1)\n");
fprintf(stdout," --backgroundPrior name of 2-column bkg frequency distribution file\n");
fprintf(stdout," --noOrientationProposal disable MCMC proposal for psi/ecc\n");
fprintf(stdout," --uniformAmplitudePrior don't use SNR-dependent amplitude prior\n");
......
......@@ -1141,13 +1141,13 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
for(ifo=0; ifo<NI; ifo++)
{
printf("Setting up time-frequency proposal for IFO%i...\n",ifo);
TFprop_setup(data, model[0], prior->range, tf, ifo);
TFprop_setup(data, model[chain->index[0]], prior->range, tf, ifo);
if(data->gnuplotFlag) TFprop_plot(prior->range, tf, data->Tobs/(double)data->N, ifo);
}
if(data->signalFlag)
{
printf("Setting up time-frequency proposal for signal model...\n");
TFprop_signal(data, prior->range, tf, model[0]->projection);
TFprop_signal(data, prior->range, tf, model[chain->index[0]]->projection);
}
}
......@@ -1491,8 +1491,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(data->signalFlag || data->glitchFlag)
{
free_TF_proposal(data, tf);
free(tf);
}
free(tf);
*logEvidence = logZ;
*varLogEvidence = varZ;
......@@ -2169,9 +2169,12 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
/*
if(model_x->size != model_y->size)
{
printf("%i: %i->%i logLy=%g logLx=%g py=%g px=%g qy=%g qx=%g, dx=%i,dy=%i,num=%g, den=%g, logH=%g\n",typ,model_x->size,model_y->size,model_y->logL,model_x->logL,wave_y->logp,wave_x->logp,logqy,logqx,wave_x->size,wave_y->size,wave_y->logp+ logqx,- wave_x->logp - logqy, logH);
}
*/
printf("%i: %i->%i logLy=%g logLx=%g py=%g px=%g qy=%g qx=%g, dx=%i,dy=%i,num=%g, den=%g, logH=%g\n",typ,model_x->size,model_y->size,model_y->logL,model_x->logL,wave_y->logp,wave_x->logp,logqy,logqx,wave_x->size,wave_y->size,wave_y->logp+ logqx,- wave_x->logp - logqy, logH);
printf(" %lg\n",prior->Nwavelet[wave_x->size]);
printf(" %lg\n",prior->Nwavelet[wave_y->size]);
printf("dif=%lg\n",prior->Nwavelet[wave_y->size]-prior->Nwavelet[wave_x->size]);
}
*/
//if(chain->index[ic]==0 && data->stochasticFlag)printf("logLy=%g logLx=%g Ay=%g Ax=%g detC=%g logH=%g\n",model_y->logL-data->logLc,model_x->logL-data->logLc,model_x->background->logamp,model_y->background->logamp, model_y->background->detCij[data->N/4],logH);
......
......@@ -896,8 +896,8 @@ void initialize_bayesline(struct BayesLineParams **bayesline, struct Data *data,
//bayesline[ifo]->priors->invsigma[i] = 1./(psd[ifo][i+(int)(Tobs*fmin)]*1.0);
bayesline[ifo]->priors->sigma[i] = psd[ifo][i+imin];
bayesline[ifo]->priors->mean[i] = psd[ifo][i+imin];
bayesline[ifo]->priors->lower[i] = psd[ifo][i+imin]/10.;
bayesline[ifo]->priors->upper[i] = psd[ifo][i+imin]*10.;
bayesline[ifo]->priors->lower[i] = psd[ifo][i+imin]/100.;
bayesline[ifo]->priors->upper[i] = psd[ifo][i+imin]*2.;
}
......@@ -1212,7 +1212,7 @@ void initialize_TF_proposal(struct Data *data, struct Prior *prior, struct TimeF
tf->N = data->N;
//4s grid for Q
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0]);
tf->nQ = (int)(prior->range[2][1] - prior->range[2][0]);
//4 Hz grid for frequency
tf->nf = (int)(prior->range[1][1] - prior->range[1][0])/4;
......
......@@ -3902,4 +3902,4 @@ int acor(double *mean, double *sigma, double *tau, double X[], int L, int maxlag
free(C);
return 0;
}
\ No newline at end of file
}
......@@ -1083,7 +1083,7 @@ void TFprop_draw(double *params, double **range, struct TimeFrequencyMap *tf, in
t = tmin + (tmax-tmin)*ran2(seed);
i = (int)((f-fmin)/df);
j = (int)((t)/dt);
j = (int)((t-tmin)/dt);
k = (int)((Q-Qmin)/dQ);
alpha = tf->max[ifo][k]*ran2(seed);
......@@ -1092,7 +1092,7 @@ void TFprop_draw(double *params, double **range, struct TimeFrequencyMap *tf, in
params[0] = t;
params[1] = f;
params[2] = Q;
params[2] = Q;
}
double TFprop_density(double *params, double **range, struct TimeFrequencyMap *tf, int ifo)
......@@ -1111,7 +1111,7 @@ double TFprop_density(double *params, double **range, struct TimeFrequencyMap *t
double dQ = (Qmax-Qmin)/(double)tf->nQ;
int i = (int)((params[0])/dt);
int i = (int)((params[0]-tmin)/dt);
int j = (int)((params[1]-fmin)/df);
int k = (int)((params[2]-Qmin)/dQ);
......@@ -1258,7 +1258,9 @@ void TFprop_setup(struct Data *data, struct Model *model, double **range, struct
free_dvector(AF,0,n-1);
free_dvector(h, 0,n-1);
free_dvector(d, 0,n-1);
free_dvector(params,0,4);
}
void TFprop_plot(double **range, struct TimeFrequencyMap *tf, double dt, int ifo)
......
......@@ -18,7 +18,7 @@
LIBS = Frame metaio gsl fftw3 fftw3f gslcblas m
#CCFLAGS = -g -Wall -O3 -std=gnu99
#CCFLAGS = -g
#DEBUGGING OPTIONS
CCFLAGS = -g -Wall -Wunused -Wextra -O3
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment