Commit 4d93d2a4 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

found memory errors


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@596 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent a751e0ed
......@@ -260,13 +260,6 @@ int main(int argc, char *argv[])
int omax = 10;
if(data->Dmax < omax) omax = data->Dmax;
// maximum number of sineGaussians
prior->gmin = ivector(0,data->NI-1);
prior->gmax = ivector(0,data->NI-1);
/* Wavelet parameters */
prior->range = dmatrix(0,4,0,1);
initialize_priors(data, prior, omax);
data->glitch_amplitude_proposal = draw_glitch_amplitude;
......@@ -358,31 +351,9 @@ int main(int argc, char *argv[])
fprintf(stdout,"BayesLine search phase for IFO %i\n", ifo);
BayesLineSearch(bayesline[0][ifo], asd, fmin, fmax, deltaT, Tobs);
// FILE *temp=fopen("temp.dat","w");
// //at frequencies below fmin output SAmax (killing lower frequencies in any future inner products)
// for(i=(int)(floor(bayesline[0][ifo]->data->fmin * bayesline[0][ifo]->data->Tobs)); i<(int)(floor(bayesline[0][ifo]->data->fmax * bayesline[0][ifo]->data->Tobs)); i++)
// {
// fprintf(temp,"%lg %lg %lg %lg\n",(double)i/data->Tobs, psd[ifo][i], asd[2*i]*asd[2*i] + asd[2*i+1]*asd[2*i+1],bayesline[0][ifo]->Snf[i-imin]);
// }
// fclose(temp);
// system_pause();
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);
// FILE *temp=fopen("temp3.dat","w");
// //at frequencies below fmin output SAmax (killing lower frequencies in any future inner products)
// for(i=(int)(floor(bayesline[0][ifo]->data->fmin * bayesline[0][ifo]->data->Tobs)); i<(int)(floor(bayesline[0][ifo]->data->fmax * bayesline[0][ifo]->data->Tobs)); i++)
// {
// fprintf(temp,"%lg %lg %lg %lg\n",(double)i/data->Tobs, asd[2*i]*asd[2*i] + asd[2*i+1]*asd[2*i+1],bayesline[0][ifo]->Snf[i-imin], model[0]->Snf[ifo][i]);
// }
// fclose(temp);
// system_pause();
//Quit lying to BayesLine about the data. Poor BayesLine.
for(i=0; i<(imax-imin); i++)
{
......
......@@ -434,8 +434,8 @@ void BayesWaveBurstInject(ProcessParamsTable *commandLine, struct Data *data, st
void getChannels(ProcessParamsTable *commandLine, char **channel)
{
ProcessParamsTable *ppt = NULL;
char ifo[10];
char arg[20];
char ifo[4];
char arg[128];
ProcessParamsTable *argument=commandLine;
......
......@@ -1235,7 +1235,12 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
if(params[5] > Tobs) params[5] -= Tobs;
// check that we haven't mapped out of range
if(extrinsic_checkrange(params)) return -1.0e60;
if(extrinsic_checkrange(params))
{
free_dmatrix(h,0,NI-1,0,N-1);
free_dmatrix(r,0,NI-1,0,N-1);
return -1.0e60;
}
//Now update the full model with current extrinsic parameters
computeProjectionCoeffs(data, projection, params);
......
......@@ -166,7 +166,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
current geocenter waveform and the detector response function.
*/
computeProjectionCoeffs(data, model[ic]->projection, model[ic]->extParams);
Shf_Geocenter(data, model[ic]->projection, model[ic]->Snf, model[ic]->SnGeo, model[ic]->extParams);
Shf_Geocenter_inv(data, model[ic]->projection, model[ic]->invSnf, model[ic]->SnGeo, model[ic]->extParams);
for(j=1; j<=signal->size; j++)
{
......@@ -1557,7 +1557,7 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
BayesLineRJMCMC(bl_x[ifo], data->r[ifo], model_x->Snf[ifo], model_x->invSnf[ifo], model_x->SnS[ifo], N, chain->cycle, chain->beta, priorFlag);
}
Shf_Geocenter(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
Shf_Geocenter_inv(data, model_x->projection, model_x->invSnf, model_x->SnGeo, model_x->extParams);
//recompute likelihoods & priors of current chain
model_x->logLnorm = 0.0;
......@@ -1681,7 +1681,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
int *larrayy;
double **paramsy;
if(data->signalFlag) Shf_Geocenter(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
if(data->signalFlag) Shf_Geocenter_inv(data, model_x->projection, model_x->invSnf, model_x->SnGeo, model_x->extParams);
acc=1;
if(data->signalFlag) copy_int_model(model_x, model_y, N, NI,-1);
......@@ -2322,7 +2322,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(data->amplitudePriorFlag)
{
//model_x->projection has be updated by data->extrinsic_likelihood, so it stores F+ and Fx for params_y.
Shf_Geocenter(data, model_x->projection, model_x->Snf, SnGeox, paramsx);
Shf_Geocenter_inv(data, model_x->projection, model_x->invSnf, SnGeox, paramsx);
for(i=1; i<ienddim; i++)
{
logpx += (data->signal_amplitude_prior(smodel->intParams[smodel->index[i]],SnGeox, 1.0, data->Tobs, prior->sSNRpeak));
......@@ -2398,7 +2398,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(data->amplitudePriorFlag)
{
//model_x->projection has be updated by data->extrinsic_likelihood, so it stores F+ and Fx for params_y.
Shf_Geocenter(data, model_x->projection, model_x->Snf, SnGeoy, paramsy);
Shf_Geocenter_inv(data, model_x->projection, model_x->invSnf, SnGeoy, paramsy);
for(i=1; i<ienddim; i++) logpy += (data->signal_amplitude_prior(intParams[i],SnGeoy, 1.0,data->Tobs, prior->sSNRpeak));
}
//logpy += dim*log(paramsy[5])/2.;
......@@ -2557,7 +2557,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
computeProjectionCoeffs(data, model_x->projection, model_x->extParams);
waveformProject(data, model_x->projection, model_x->extParams, model_x->response, model_x->signal->templates, data->fmin, data->fmax);
Shf_Geocenter(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
Shf_Geocenter_inv(data, model_x->projection, model_x->invSnf, model_x->SnGeo, model_x->extParams);
free(glitch);
......
......@@ -48,6 +48,30 @@ void Shf_Geocenter(struct Data *data, struct Network *projection, double **Snf,
for(i=imin; i<imax; i++) SnGeo[i] = 1./SnGeo[i];
}
void Shf_Geocenter_inv(struct Data *data, struct Network *projection, double **invSnf, double *SnGeo, double *params)
{
int ifo,i, halfN=data->N/2;
int NI = data->NI;
int imin = data->imin;
int imax = data->imax;
double AntennaPattern;
double ecc = params[3];
for(i=0; i<halfN; i++) SnGeo[i] = 0.0;
for(ifo=0; ifo<NI; ifo++)
{
AntennaPattern = projection->Fplus[ifo]*projection->Fplus[ifo] + ecc*ecc*projection->Fcross[ifo]*projection->Fcross[ifo];
for(i=imin; i<imax; i++) SnGeo[i] += AntennaPattern*invSnf[ifo][i];
}
for(i=0; i<imin; i++)SnGeo[i]=1.0;
for(i=imax; i<halfN; i++)SnGeo[i]=1.0;
for(i=imin; i<imax; i++) SnGeo[i] = 1./SnGeo[i];
}
void OverlapReductionFunction(struct Data *data)
{
......@@ -649,6 +673,11 @@ void initialize_data(struct Data *data, double **s, int N, int tsize, double Tob
data->varZnoise = 0.0;
data->varZsignal = 0.0;
data->varZfull = 0.0;
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 0;
}
void initialize_model(struct Model *model, struct Data *data, struct Prior *prior, double **psd, long *seed)
......
......@@ -5,6 +5,8 @@
/* ********************************************************************************** */
void Shf_Geocenter(struct Data *data, struct Network *projection, double **Snf, double *SnGeo, double *params);
void Shf_Geocenter_inv(struct Data *data, struct Network *projection, double **invSnf, double *SnGeo, double *params);
void OverlapReductionFunction(struct Data *data);
void ComputeNoiseCorrelationMatrix(struct Data *data, double **Snf, double *Sn, struct Background *background);
......
......@@ -408,15 +408,18 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
void extrinsic_fisher_update(struct Data *data, struct Model *model)
{
int i;
int i,j;
double *params = model->extParams;
double *geo = model->signal->templates;
double *Sn = model->Sn;
double **glitch = NULL;
glitch = malloc(data->NI*sizeof(double*));
for(i=0; i<data->NI; i++) glitch[i]=model->glitch[i]->templates;
double **glitch = malloc(data->NI*sizeof(double*));
for(i=0; i<data->NI; i++)
{
glitch[i]=malloc(data->N*sizeof(double));
for(j=0; j<data->N; j++)glitch[i][j] = model->glitch[i]->templates[j];
}
struct FisherMatrix *fisher = model->fisher;
......@@ -425,6 +428,7 @@ void extrinsic_fisher_update(struct Data *data, struct Model *model)
matrix_eigenstuff(fisher->matrix, fisher->evector, fisher->evalue, fisher->N);
for(i=0; i<data->NI; i++) free(glitch[i]);
free(glitch);
}
......@@ -438,6 +442,8 @@ void fisher_matrix_proposal(struct FisherMatrix *fisher, double *params, long *s
double Amps[N];
double jump[N];
for(i=0; i<N; i++) Amps[i]=jump[i]=0.0;
// not used since we jump along eigendirections
double sqD = 1./sqrt((double)(N));
......
......@@ -229,7 +229,6 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
int ifo, iend=N/2;
double df,freq;
double toff;
double Fplus, Fcross;
// common extrinsic parameters
......@@ -245,7 +244,8 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
double dim,dre;
double trigArgument;
double toff = 0.0;
double trigArgument=0.0;
for(ifo = 0; ifo<NI; ifo++)
{
......
......@@ -83,4 +83,4 @@ BayesWavePost: BayesWavePost.c $(OBJS)
$(CC) $(CCFLAGS) $(LDFLAGS) -Wl,-rpath -Wl,$(LAL_LIBS) -o bayeswave_post BayesWavePost.c $(OBJS) $(LIBS:%=-l%)
clean:
rm *.o bayeswave bayeswave_post
rm -f *.o bayeswave bayeswave_post
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