Commit eb6625f3 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

memory fix, optimize integral in ext. likelihood


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@597 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 4d93d2a4
......@@ -1659,7 +1659,7 @@ void parse_signal_parameters(struct Data *data, struct Model *model, FILE *param
for(i=0; i<NE; i++) fscanf(paramsfile,"%lg", &model->extParams[i]);
//Compute network projection coefficients
computeProjectionCoeffs(data, model->projection, model->extParams);
computeProjectionCoeffs(data, model->projection, model->extParams, data->fmin, data->fmax);
//Get each wavelets parameters and fill up signal structure
for(d=1; d<=signal->size; d++)
......@@ -1695,7 +1695,7 @@ void dump_time_domain_waveform(struct Data *data, struct Model *model, struct Pr
for(i=0; i<data->N; i++) signal->templates[i] = 0.0;
//Compute network projection coefficients
computeProjectionCoeffs(data, model->projection, model->extParams);
computeProjectionCoeffs(data, model->projection, model->extParams, data->fmin, data->fmax);
//Get each wavelets parameters and fill up signal structure
for(d=1; d<=signal->size; d++)
......
......@@ -279,7 +279,7 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
model_x->wavelet(model_x->signal->templates, model_x->signal->intParams[i], N, 1, data->Tobs);
}
computeProjectionCoeffs(data, model_x->projection, model_x->extParams);
computeProjectionCoeffs(data, model_x->projection, model_x->extParams, data->fmin, data->fmax);
waveformProject(data, model_x->projection, model_x->extParams, model_x->response, model_x->signal->templates, data->fmin, data->fmax);
}
......@@ -1142,8 +1142,10 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
NI = data->NI;
N = data->N;
imin = data->imin;
imax = data->imax;
//imin = data->imin;
//imax = data->imax;
imin = (int)(fmin*data->Tobs);
imax = (int)(fmax*data->Tobs);
d = data->s;
h = dmatrix(0,NI-1,0,N-1); //template
......@@ -1154,7 +1156,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
for(i=0; i<NI; i++) for(n=0; n<N; n++) h[i][n] = 0.0;
//compute instrument response h[] to geocenter waveform geo[]
computeProjectionCoeffs(data, projection, params);
computeProjectionCoeffs(data, projection, params, fmin, fmax);
waveformProject(data, projection, params, h, geo, fmin, fmax);
......@@ -1205,7 +1207,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
for(i=0; i<NI; i++) for(n=0; n<N; n++) h[i][n] = 0.0;
//compute instrument response h[] to geocenter waveform geo[]
computeProjectionCoeffs(data, projection, params);
computeProjectionCoeffs(data, projection, params, data->fmin, data->fmax);
waveformProject(data, projection, params, h, geo, data->fmin, data->fmax);
//Form up residual
......@@ -1243,8 +1245,8 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
}
//Now update the full model with current extrinsic parameters
computeProjectionCoeffs(data, projection, params);
waveformProject(data, projection, params, h, geo, fmin, fmax);
computeProjectionCoeffs(data, projection, params, data->fmin, data->fmax);
waveformProject(data, projection, params, h, geo, data->fmin, data->fmax);
//Form up residual
for(ifo=0; ifo<NI; ifo++)
{
......
......@@ -165,7 +165,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
-waveformProject uses the stored coeffs and does the convolution between the
current geocenter waveform and the detector response function.
*/
computeProjectionCoeffs(data, model[ic]->projection, model[ic]->extParams);
computeProjectionCoeffs(data, model[ic]->projection, model[ic]->extParams, data->fmin, data->fmax);
Shf_Geocenter_inv(data, model[ic]->projection, model[ic]->invSnf, model[ic]->SnGeo, model[ic]->extParams);
for(j=1; j<=signal->size; j++)
......@@ -2294,6 +2294,8 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(i=0; i<NE; i++) paramsy[i] = paramsx[i];
//Find minimum and maximum frequencies of signal model and only compute logL over that range
// double fmin = data->fmin;
// double fmax = data->fmax;
double fmin = data->fmax;
double fmax = data->fmin;
double fi,ff;
......@@ -2303,13 +2305,12 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(fi<fmin) fmin=fi;
if(ff>fmax) fmax=ff;
}
//Make sure frequency is in range
if(fmin < data->fmin)
fmin = data->fmin;
fmin = data->fmin;
if(fmax > data->fmax)
fmax = data->fmax;
}
fmax = data->fmax;
logLx = data->extrinsic_likelihood(model_x->projection, paramsx, model_x->invSnf, Sn, geo, glitch, data, fmin, fmax);
......@@ -2555,7 +2556,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
model_x->extParams[5] = 1.0;//amplitude
computeProjectionCoeffs(data, model_x->projection, model_x->extParams);
computeProjectionCoeffs(data, model_x->projection, model_x->extParams, data->fmin, data->fmax);
waveformProject(data, model_x->projection, model_x->extParams, model_x->response, model_x->signal->templates, data->fmin, data->fmax);
Shf_Geocenter_inv(data, model_x->projection, model_x->invSnf, model_x->SnGeo, model_x->extParams);
......
......@@ -950,19 +950,26 @@ void reset_params(struct Wavelet *wave_x, struct Wavelet *wave_y, double **range
void initialize_network(struct Network *projection, int N, int NI)
{
projection->expPhase = dmatrix(0,NI-1,0,N-1);
projection->deltaT = dvector(0,NI-1);
projection->Fplus = dvector(0,NI-1);
projection->Fcross = dvector(0,NI-1);
projection->dtimes = dvector(0,NI-1);
projection->expPhase = malloc(NI*sizeof(double *));
projection->deltaT = malloc(NI*sizeof(double));
projection->Fplus = malloc(NI*sizeof(double));
projection->Fcross = malloc(NI*sizeof(double));
projection->dtimes = malloc(NI*sizeof(double));
int i;
for(i=0; i<NI; i++) projection->expPhase[i] = malloc(N*sizeof(double));
}
void free_network(struct Network *projection, int N, int NI)
{
free_dmatrix(projection->expPhase,0,NI-1,0,N-1);
free_dvector(projection->deltaT,0,NI-1);
free_dvector(projection->Fplus,0,NI-1);
free_dvector(projection->Fcross,0,NI-1);
int i;
for(i=0; i<NI; i++) free(projection->expPhase[i]);
free(projection->expPhase);
free(projection->deltaT);
free(projection->Fplus);
free(projection->Fcross);
free(projection->dtimes);
}
void copy_int_model(struct Model *origin, struct Model *copy, int N, int NI, int det)
......@@ -1012,16 +1019,17 @@ void copy_psd_model(struct Model *origin, struct Model *copy, int N, int NI)
void copy_ext_model(struct Model *origin, struct Model *copy, int N, int NI)
{
int i,ifo;
for(i=0; i<NE; i++) copy->extParams[i] = origin->extParams[i];
int i,ifo;
for(i=0; i<NE; i++) copy->extParams[i] = origin->extParams[i];
for(ifo=0; ifo<NI; ifo++)
{
copy->projection->Fplus[ifo] = origin->projection->Fplus[ifo];
copy->projection->Fcross[ifo] = origin->projection->Fcross[ifo];
copy->projection->deltaT[ifo] = origin->projection->deltaT[ifo];
for(i=0; i<N; i++) copy->projection->expPhase[ifo][i] = origin->projection->expPhase[ifo][i];
}
for(ifo=0; ifo<NI; ifo++)
{
copy->projection->Fplus[ifo] = origin->projection->Fplus[ifo];
copy->projection->Fcross[ifo] = origin->projection->Fcross[ifo];
copy->projection->deltaT[ifo] = origin->projection->deltaT[ifo];
copy->projection->dtimes[ifo] = origin->projection->dtimes[ifo];
for(i=0; i<N; i++) copy->projection->expPhase[ifo][i] = origin->projection->expPhase[ifo][i];
}
}
void free_model(struct Model *model, struct Data *data, struct Prior *prior)
......
......@@ -270,7 +270,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
epsilon = dvector(0,NE-1);
// Compute SNR of signal model
computeProjectionCoeffs(data, projectionP, params);
computeProjectionCoeffs(data, projectionP, params, data->fmin, data->fmax);
waveformProject(data, projectionP, params, hP, geo, data->fmin, data->fmax);
//store data
......@@ -319,8 +319,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
paramsP[i] += epsilon[i];
paramsM[i] -= epsilon[i];
computeProjectionCoeffs(data, projectionP, paramsP);
computeProjectionCoeffs(data, projectionM, paramsM);
computeProjectionCoeffs(data, projectionP, paramsP, data->fmin, data->fmax);
computeProjectionCoeffs(data, projectionM, paramsM, data->fmin, data->fmax);
logLP = EvaluateExtrinsicMarkovianLogLikelihood(projectionP, paramsP, invSnf, Sn, geo, glitch, data, data->fmin, data->fmax);
logLM = EvaluateExtrinsicMarkovianLogLikelihood(projectionM, paramsM, invSnf, Sn, geo, glitch, data, data->fmin, data->fmax);
......@@ -344,8 +344,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
paramsP[i] += epsilon[i];
paramsM[i] -= epsilon[i];
computeProjectionCoeffs(data, projectionP, paramsP);
computeProjectionCoeffs(data, projectionM, paramsM);
computeProjectionCoeffs(data, projectionP, paramsP, data->fmin, data->fmax);
computeProjectionCoeffs(data, projectionM, paramsM, data->fmin, data->fmax);
waveformProject(data, projectionP, paramsP, hP, geo, data->fmin, data->fmax);
waveformProject(data, projectionM, paramsM, hM, geo, data->fmin, data->fmax);
......@@ -362,8 +362,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
*/
if(extrinsic_checkrange(paramsP) && extrinsic_checkrange(paramsM))
{
computeProjectionCoeffs(data, projectionP, paramsP);
computeProjectionCoeffs(data, projectionM, paramsM);
computeProjectionCoeffs(data, projectionP, paramsP, data->fmin, data->fmax);
computeProjectionCoeffs(data, projectionM, paramsM, data->fmin, data->fmax);
}
waveformProject(data, projectionP, paramsP, hPP, geo, data->fmin, data->fmax);
......
......@@ -217,7 +217,7 @@ void SineGaussianFisher(double *params, struct FisherMatrix *fisher)
/* ********************************************************************************** */
void computeProjectionCoeffs(struct Data *data, struct Network *projection, double *extParams)
void computeProjectionCoeffs(struct Data *data, struct Network *projection, double *extParams, double fmin, double fmax)
{
int i,re,im;
......@@ -226,7 +226,12 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
int NI = data->NI;
double Tobs = data->Tobs;
// network time-delays and antennae patterns
int ifo, iend=N/2;
// int ifo, iend=N/2;
int ifo;
int istart = (int)(fmin*Tobs);
int iend = (int)(fmax*Tobs);
if(istart<0) istart = 0;
if(iend>N/2) iend = N/2;
double df,freq;
double Fplus, Fcross;
......@@ -263,7 +268,8 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
projection->deltaT[ifo] = toff;
//calculate frequency domain waveforms:
freq = 0.0; // initialize frequency
// freq = 0.0; // initialize frequency
freq = fmin; // initialize frequency
df = 1./Tobs;
//Use recursion relationship for cos(Phase) & sin(Phase)
......@@ -279,7 +285,8 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
dre = sin(dre);
dre = -2.0*dre*dre;
for(i=0; i<iend; i++)
//for(i=0; i<iend; i++)
for(i=istart; i<iend; i++)
{
re = 2*i;
im = re+1;
......@@ -332,7 +339,7 @@ void waveformProject(struct Data *data, struct Network *projection, double *extP
even = 2*i;
odd = even + 1;
if((freq<fmin) || (freq>fmax)) // frequency below LIGO band or above termination point
if((freq<fmin) || (freq>fmax)) // frequency out of band for data or signal
{
response[id][even] = 0.0;
response[id][odd] = 0.0;
......
......@@ -15,6 +15,6 @@ void SineGaussianFisher(double *params, struct FisherMatrix *fisher);
/* */
/* ********************************************************************************** */
void computeProjectionCoeffs(struct Data *data, struct Network *projection, double *extParams);
void computeProjectionCoeffs(struct Data *data, struct Network *projection, double *extParams, double fmin, double fmax);
void waveformProject(struct Data *data, struct Network *projection, double *extParams, double **response, double *geo, double fmin, double fmax);
void compute_reconstructed_plus_and_cross(struct Data *data, struct Network *projection, double *extParams, double **hrecPlus, double **hrecCross, double *geo, double fmin, double fmax);
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