Commit d139ec6f authored by Meg Millhouse's avatar Meg Millhouse
Browse files

merge chirplets

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@825 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 655d2c65
......@@ -107,6 +107,13 @@ int main(int argc, char *argv[])
struct Prior *prior = malloc(sizeof(struct Prior));
parse_command_line(data, chain, prior, runState->commandLine);
if (data->chirpletFlag) {
data->NW = 6;
}
else{
data->NW = 5;
}
// Move to the directory specified by the user, or stay in ./
chdir(data->outputDirectory);
......@@ -231,7 +238,7 @@ int main(int argc, char *argv[])
Setup PRIOR structure
*/
// int Dmax = data->Dmax;
// data->Dmax = 100;
// data->Dmax = 200;
int omax = 10;
if(data->Dmax < omax) omax = data->Dmax;
......
......@@ -56,7 +56,7 @@
#define CLIGHT 2.99792458e8 // Speed of light (m/s)
#define NE 6 //Number of extrinsic parameters
#define NW 5 //Number of intrinsic parameters
// #define NW 5 //Number of intrinsic parameters
//Tuning parameters for wavelet proximity proposal
//TODO: These should not be global!
......@@ -102,6 +102,9 @@ struct Model
double *detLogL;
double *extParams;
double **response;
int chirpletFlag;
int NW;
struct Wavelet *signal;
struct Wavelet **glitch;
......@@ -292,7 +295,11 @@ struct Data
int nukeFlag; //tell BWP to remove chains/ and checkpoint/ directories
int srate;
int srate;
int chirpletFlag;
int NW;
double dt;
double df;
double fmin;
......
......@@ -2303,7 +2303,8 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
{
int i,j,ii, ifo;
int NI = data->NI;
int NW = data->NW;
double logP;
double detC=0;
double logV=0;
......@@ -2339,7 +2340,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
for(i=1; i<=model->signal->size; i++)
{
ii = model->signal->index[i];
intrinsic_fisher_update(model->signal->intParams[ii], sigmas, model->SnGeo, 1.0, data->Tobs);
intrinsic_fisher_update(model->signal->intParams[ii], sigmas, model->SnGeo, 1.0, data->Tobs, data->NW);
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
}
......@@ -2352,7 +2353,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
{
ii = model->glitch[ifo]->index[i];
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,model->glitch[ifo]->intParams[ii][j]);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], 1.0, data->Tobs);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], 1.0, data->Tobs, data->NW);
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,sigmas[j]);
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
......@@ -2393,7 +2394,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
maxAmp=-1.e60;
for(i=0; i<100; i++)
{
draw_wavelet_params(params, prior->range, &chain->seed);
draw_wavelet_params(params, prior->range, &chain->seed,data->NW);
data->signal_amplitude_proposal(params, model->SnGeo, 1.0, &chain->seed, data->Tobs, prior->range,prior->sSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
......@@ -2420,7 +2421,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
maxAmp=-1.e60;
for(i=0; i<100; i++)
{
draw_wavelet_params(params, prior->range, &chain->seed);
draw_wavelet_params(params, prior->range, &chain->seed,data->NW);
data->glitch_amplitude_proposal(params, model->Snf[ifo], 1.0, &chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
if(params[3]<minAmp)minAmp=params[3];
if(params[3]>maxAmp)maxAmp=params[3];
......
......@@ -951,6 +951,11 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
if(data->gnuplotFlag) fprintf(fptr,"ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " chirplet option is ...... ");
if(data->chirpletFlag) fprintf(fptr,"ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, "\n");
}
......@@ -971,7 +976,7 @@ void print_chain_files(struct Data *data, struct Chain *chain, struct Model **mo
if(data->signalFlag) print_signal_model(chain->intWaveChainFile[ic], model_x);
//print glitch model parameters
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) print_glitch_model(chain->intGlitchChainFile[ifo][ic], model_x->glitch[ifo]);
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) print_glitch_model(chain->intGlitchChainFile[ifo][ic], model_x->glitch[ifo], data->NW);
//print PSD parameters
if(data->bayesLineFlag)
......@@ -1029,6 +1034,7 @@ void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Mode
void print_signal_model(FILE *fptr, struct Model *model)
{
int i,j,k;
int NW = model->NW;
fprintf(fptr,"%i ",model->signal->size);
for(k=0; k<NE; k++) fprintf(fptr,"%.12g ", model->extParams[k]);
......@@ -1037,12 +1043,12 @@ void print_signal_model(FILE *fptr, struct Model *model)
i = model->signal->index[j];
//print all parameters
for(k=0; k<=4; k++) fprintf(fptr,"%.12g ", model->signal->intParams[i][k]);
for(k=0; k< NW; k++) fprintf(fptr,"%.12g ", model->signal->intParams[i][k]);
}
fprintf(fptr,"\n");
}
void print_glitch_model(FILE *fptr, struct Wavelet *glitch)
void print_glitch_model(FILE *fptr, struct Wavelet *glitch, int NW)
{
int i,j,k;
fprintf(fptr,"%i ",glitch->size);
......@@ -1051,7 +1057,7 @@ void print_glitch_model(FILE *fptr, struct Wavelet *glitch)
i = glitch->index[j];
//print all parameters
for(k=0; k<=4; k++) fprintf(fptr,"%.12g ", glitch->intParams[i][k]);
for(k=0; k< NW; k++) fprintf(fptr,"%.12g ", glitch->intParams[i][k]);
}
fprintf(fptr,"\n");
}
......@@ -1620,6 +1626,11 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
ppt = LALInferenceGetProcParamVal(commandLine, "--nuke");
if(ppt) data->nukeFlag = 1;
else data->nukeFlag = 0;
/* chirplet basis */
ppt = LALInferenceGetProcParamVal(commandLine, "--chirplets");
if(ppt) data->chirpletFlag = 1;
else data->chirpletFlag = 0;
/*** Get GMST for the trigger time ***/
......@@ -1785,6 +1796,7 @@ void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **para
int ifo;
int d;
int i;
int NW = data->NW;
struct Wavelet **glitch = model->glitch;
......@@ -1820,7 +1832,8 @@ void parse_signal_parameters(struct Data *data, struct Model *model, FILE *param
int ifo;
int d;
int i;
int NW = data->NW;
struct Wavelet *signal = model->signal;
//Zero out signal->templates array (holds the linear combination of wavelets)
......
......@@ -40,7 +40,7 @@ void flush_chain_files(struct Data *data, struct Chain *chain, int ic);
void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Model *model);
void print_signal_model(FILE *fptr, struct Model *model);
void print_glitch_model(FILE *fptr, struct Wavelet *glitch);
void print_glitch_model(FILE *fptr, struct Wavelet *glitch, int NW);
void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
......
......@@ -354,6 +354,7 @@ double EvaluateSearchLogLikelihood(int typ, int ii, int det, struct Model *mx, s
int NI = data->NI;
int imin = data->imin;
int imax = data->imax;
int NW = data->NW;
double Tobs = data->Tobs;
double **r = data->r;
......@@ -413,7 +414,7 @@ double EvaluateSearchLogLikelihood(int typ, int ii, int det, struct Model *mx, s
// Rejection sample on prior range
if(wave_y->size>0 )
{
if(checkrange(paramsy,prior->range) )
if(checkrange(paramsy,prior->range,NW) )
{
free(gx);
free(gy);
......@@ -616,7 +617,7 @@ double EvaluateSearchLogLikelihood(int typ, int ii, int det, struct Model *mx, s
// check that we haven't mapped out of range
if(wave_y->size>0)
{
if(checkrange(paramsy,prior->range))
if(checkrange(paramsy,prior->range,NW))
{
free(gx);
free(gy);
......@@ -661,7 +662,7 @@ double EvaluateConstantLogLikelihood(int typ, int ii, int det, UNUSED struct Mod
*/
if(typ!=3 && wave->size > 0)
{
if( checkrange(wave->intParams[ii],prior->range)) return -1.0e60;
if( checkrange(wave->intParams[ii],prior->range,data->NW)) return -1.0e60;
}
return 0.0;
......@@ -683,6 +684,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
int imax = data->imax;
int jmin = imin;
int jmax = imax;
int NW = data->NW;
double **r = data->r;
double **s = data->s;
......@@ -740,7 +742,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
*/
if( (data->signalFlag || data->glitchFlag) && typ!=3 && wave_y->size>0 )
{
if(checkrange(paramsy,prior->range))
if(checkrange(paramsy,prior->range,NW))
{
free(gx);
free(gy);
......@@ -888,7 +890,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
wBirth->size++;
// shift time of signal
for(i=0; i<5; i++) wBirth->intParams[jj][i] = wDeath->intParams[ii][i];
for(i=0; i<NW; i++) wBirth->intParams[jj][i] = wDeath->intParams[ii][i];
if(det==-1)
{
wBirth->intParams[jj][0] += my->projection->deltaT[ifo];
......
......@@ -148,7 +148,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
model[ic]->size++;
for(j=1; j<=glitch[ifo]->size; j++)
{
draw_wavelet_params(glitch[ifo]->intParams[j], prior->range, &chain->seed);
draw_wavelet_params(glitch[ifo]->intParams[j], prior->range, &chain->seed, data->NW);
if(data->amplitudePriorFlag) data->glitch_amplitude_proposal(glitch[ifo]->intParams[j], model[ic]->Snf[ifo], 1.0, &chain->seed, data->Tobs, prior->range, prior->gSNRpeak);
model[ic]->wavelet(glitch[ifo]->templates, glitch[ifo]->intParams[j], data->N, 1, Tobs);
......@@ -175,7 +175,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
for(j=1; j<=signal->size; j++)
{
draw_wavelet_params(signal->intParams[j], prior->range, &chain->seed);
draw_wavelet_params(signal->intParams[j], prior->range, &chain->seed, data->NW);
if(data->amplitudePriorFlag) data->signal_amplitude_proposal(signal->intParams[j], model[ic]->SnGeo, 1.0, &chain->seed, data->Tobs, prior->range, prior->sSNRpeak);
model[ic]->wavelet(signal->templates, signal->intParams[j], data->N, 1, Tobs);
}
......@@ -760,7 +760,7 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
{
sprintf(filename,"checkpoint/params_ifo%i.dat.%i",ifo,ic);
fptr = fopen(filename,"w");
print_glitch_model(fptr, glitch[ifo]);
print_glitch_model(fptr, glitch[ifo],data->NW);
fclose(fptr);
}
}
......@@ -995,8 +995,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
//Set MAP model to initial state
copy_psd_model(model[chain->index[0]], model_MAP, N, NI);
copy_ext_model(model[chain->index[0]], model_MAP, N, NI);
copy_int_model(model[chain->index[0]], model_MAP, N, NI,-1);
for(ifo=0; ifo<NI; ifo++) copy_int_model(model[chain->index[0]], model_MAP, N, NI,ifo);
copy_int_model(model[chain->index[0]], model_MAP, N, NI, data->NW, -1);
for(ifo=0; ifo<NI; ifo++) copy_int_model(model[chain->index[0]], model_MAP, N, NI,data->NW, ifo);
FILE *script = NULL;
FILE *psdscript = NULL;
......@@ -1269,8 +1269,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
logPostMap = logPost;
copy_psd_model(model[chain->index[0]], model_MAP, N, NI);
copy_ext_model(model[chain->index[0]], model_MAP, N, NI);
copy_int_model(model[chain->index[0]], model_MAP, N, NI,-1);
for(ifo=0; ifo<NI; ifo++) copy_int_model(model[chain->index[0]], model_MAP, N, NI,ifo);
copy_int_model(model[chain->index[0]], model_MAP, N, NI, data->NW, -1);
for(ifo=0; ifo<NI; ifo++) copy_int_model(model[chain->index[0]], model_MAP, N, NI, data->NW, ifo);
// LaplaceApproximation(data, model_MAP, chain, prior, &logZ);
}
}
......@@ -1558,6 +1558,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
int clusterFlag;
int uniformFlag;
int densityFlag;
int NW = data->NW;
//ratio of TF proposals to proximity proposals
double TFtoProx = data->TFtoProx;
......@@ -1603,7 +1604,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
struct Model *model_y = malloc(sizeof(struct Model));
initialize_model(model_y, data, prior, model_x->Snf, seed);
copy_ext_model(model_x, model_y, N, NI);
copy_int_model(model_x, model_y, N, NI,-1);
copy_int_model(model_x, model_y, N, NI, data->NW, -1);
copy_psd_model(model_x, model_y, N, NI);
model_y->size=model_x->size;
......@@ -1617,8 +1618,8 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(data->signalFlag) Shf_Geocenter_full(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
acc=1;
if(data->signalFlag) copy_int_model(model_x, model_y, N, NI,-1);
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) copy_int_model(model_x, model_y, N, NI, ifo);
if(data->signalFlag) copy_int_model(model_x, model_y, N, NI, data->NW, -1);
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) copy_int_model(model_x, model_y, N, NI, data->NW, ifo);
for(mc=1; mc<=M; mc++)
{
......@@ -1630,7 +1631,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
phaseFlag=0;
//acc gets set to 1 if model_y gets accepted (no need to copy at next iteration)
if(acc==0) copy_int_model(model_x, model_y, N, data->NI, det);
if(acc==0) copy_int_model(model_x, model_y, N, data->NI, data->NW, det);
if(data->psdFitFlag) for(ifo=0; ifo<NI; ifo++) Sny[ifo] = Snx[ifo];
else for(ifo=0; ifo<NI; ifo++) Sny[ifo] = 1.0;
......@@ -1724,7 +1725,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
k++;
larrayy[k] = larrayx[j];
jj = larrayx[j];
for(kk=0; kk <= 4; kk++) paramsy[jj][kk] = paramsx[jj][kk];
for(kk=0; kk < NW; kk++) paramsy[jj][kk] = paramsx[jj][kk];
}
if(j == i) ii = larrayx[j]; // take note of who's demise is being proposed
}
......@@ -1779,7 +1780,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
larrayy[wave_y->size] = ii;
// Draw all wavelet parameters from prior
draw_wavelet_params(paramsy[ii], range, seed);
draw_wavelet_params(paramsy[ii], range, seed, data->NW);
// TF density for proposal
if(data->clusterProposalFlag )
......@@ -1891,7 +1892,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
// Draw all wavelet parameters from prior
draw_wavelet_params(paramsy[ii], range, seed);
draw_wavelet_params(paramsy[ii], range, seed, data->NW);
// TF density for proposal
if(data->clusterProposalFlag)
......@@ -1966,15 +1967,15 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
//TODO: Try new fisher proposal
if(det==-1) intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, seed, Tobs, &test);
else intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, Tobs, &test);
if(det==-1) intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, seed, Tobs, &test, data->NW);
else intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, Tobs, &test, data->NW);
//if(det==-1) intrinsic_fisher_proposal_2(model_x,range, paramsx[ii], paramsy[ii], model_x->SnGeo, 1.0, &logqx, &logqy, seed, Tobs, &test);
//else intrinsic_fisher_proposal_2(model_x,range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, Tobs, &test);
}
else
{
draw_wavelet_params(paramsy[ii], range, seed);
draw_wavelet_params(paramsy[ii], range, seed, data->NW);
if(ran2(seed)<0.5 && data->amplitudeProposalFlag)
{
......@@ -2108,7 +2109,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
acc = 1;
copy_int_model(model_y, model_x, N, data->NI, det);
copy_int_model(model_y, model_x, N, data->NI, data->NW, det);
if(data->stochasticFlag) copy_background(model_y->background, model_x->background, data->NI, data->N/2);
......@@ -2154,6 +2155,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//Unpack structures and use convenient (and familiar) names
int NI = data->NI;
int NW = data->NW;
/* CHAIN */
int M = 3*chain->cycle;
......@@ -2173,7 +2175,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
dim = smodel->size;
ienddim = dim+1;
iendN = data->N/2;
intParams = dmatrix(1,dim,0,4);
intParams = dmatrix(1,dim,0,NW-1);
//Set extrinsic parameter fisher matrix (numerical derivatives of response to geocenter waveform)
......@@ -2237,7 +2239,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
logJ=0.;
sky=0;
test=0;
for(i=1;i<ienddim;i++) for(j=0; j<5; j++) intParams[i][j] = smodel->intParams[smodel->index[i]][j];
for(i=1;i<ienddim;i++) for(j=0; j<NW; j++) intParams[i][j] = smodel->intParams[smodel->index[i]][j];
for(i=0; i<iendN; i++) SnGeoy[i]=SnGeox[i];
//Initialize parameter vectors & likelihood for extrinsic parameter MCMC
......@@ -2281,7 +2283,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(i=1; i<ienddim; i++)
{
intParams[i][3] = smodel->intParams[smodel->index[i]][3]*paramsy[5];//amplitude
test += checkrange(intParams[i],prior->range);
test += checkrange(intParams[i],prior->range,NW);
}
if(test==0)
......@@ -2336,7 +2338,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
free(glitch);
free_dvector(paramsy,0,NE);
free_dvector(SnGeoy,0,iendN);
free_dmatrix(intParams,1,dim,0,4);
free_dmatrix(intParams,1,dim,0,NW-1);
}
/* ********************************************************************************** */
......
......@@ -366,7 +366,9 @@ void resize_model(struct Data *data, struct Chain *chain, struct Prior *prior, s
void initialize_priors(struct Data *data, struct Prior *prior, int omax)
{
int i,j,n;
int NW = data->NW;
/* Model dimension parameters */
// number of sineGaussians we print to file
prior->omax = omax;
......@@ -399,7 +401,7 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
/* Wavelet parameters */
prior->range = dmatrix(0,4,0,1);
prior->range = dmatrix(0,NW-1,0,1);
double Amin=0.0, Amax=1.0;
//Set Amplitude priors based on the data channel being used
......@@ -438,6 +440,8 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
double tmax = data->Tobs;//-tmin;
double fmin = data->fmin;
double fmax = data->fmax;
double betamin = -1.0;
double betamax = 1.0;
/*
if(data->amplitudePriorFlag)
......@@ -457,6 +461,11 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
prior->range[3][1] = Amax;
prior->range[4][0] = 0.0;
prior->range[4][1] = LAL_TWOPI;
if (NW==6) {
prior->range[5][0] = betamin;
prior->range[5][1] = betamax;
}
prior->TFV = (tmax-tmin)*(fmax-fmin)*(Qmax-Qmin);
prior->logTFV = log(prior->TFV);
......@@ -611,7 +620,10 @@ void reset_priors(struct Data *data, struct Prior *prior)
double tmax = data->Tobs;//-tmin;
double fmin = data->fmin;
double fmax = data->fmax;
double betamin = -1.0;
double betamax = 1.0;
//restrict prior range for wavelet time to 1 second around trigger
if(data->runPhase==1)
{
......@@ -709,6 +721,7 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
int N = data->N;
int halfN = data->N/2;
int NI = data->NI;
int NW = data->NW;
//gaussian noise model
model->Sn = dvector(0,NI-1);
model->SnS = dmatrix(0,NI-1,0,halfN-1);
......@@ -756,9 +769,9 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//Setup wavelet structure holding signal model
if(data->signalFlag)
initialize_wavelet(model->signal,N, prior->smax, 1);
initialize_wavelet(model->signal,N, NW,prior->smax, 1);
else
initialize_wavelet(model->signal,N, prior->smax, 0);
initialize_wavelet(model->signal,N, NW,prior->smax, 0);
//allocate memory for structure holding projection coefficients
model->projection = malloc(sizeof(struct Network));
......@@ -772,9 +785,9 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
for(i=0; i<NI; i++)
{
if(data->glitchFlag)
initialize_wavelet(model->glitch[i],N, prior->smax, 1);
initialize_wavelet(model->glitch[i],N, NW,prior->smax, 1);
else
initialize_wavelet(model->glitch[i],N, prior->smax, 0);
initialize_wavelet(model->glitch[i],N, NW,prior->smax, 0);
}
//allocate memory for structure holding extrinsic fisher matrix
......@@ -786,8 +799,16 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
initialize_fisher(model->intfisher, NW);
//select wavelet basis function
model->wavelet = SineGaussianFourier;
model->wavelet_bandwidth = SineGaussianBandwidth;
if (data->chirpletFlag) {
model->wavelet = ChirpletFourier;
model->wavelet_bandwidth = ChirpletBandwidth;
}
else{
model->wavelet = SineGaussianFourier;
model->wavelet_bandwidth = SineGaussianBandwidth;
}
//Stochastic background
if(data->stochasticFlag)
......@@ -795,6 +816,14 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
model->background = malloc(sizeof(struct Background));
initialize_background(model->background, NI, halfN);
}
// Chirplet param check
if (data->chirpletFlag) {
model->NW = 6;
}
else{
model->NW = 5;
}
}
......@@ -931,10 +960,12 @@ void reset_likelihood(struct Data *data)
}
}
void reset_params(struct Wavelet *wave_x, struct Wavelet *wave_y, double **range)
// Unused
/*
void reset_params(struct Wavelet *wave_x, struct Wavelet *wave_y, double **range, int NW)
{
int i,j,k;
int sx = wave_x->size;
int smax = wave_x->smax;
......@@ -959,9 +990,10 @@ void reset_params(struct Wavelet *wave_x, struct Wavelet *wave_y, double **range
{
larrayy[j] = larrayx[j];
i = larrayx[j];
for(k=0; k <= 4; k++) paramsy[i][k] = paramsx[i][k];
for(k=0; k < NW; k++) paramsy[i][k] = paramsx[i][k];
}
}
*/
void initialize_network(struct Network *projection, int N, int NI)
{
......@@ -987,7 +1019,8 @@ void free_network(struct Network *projection, int NI)
free(projection->dtimes);
}
void copy_int_model(struct Model *origin, struct Model *copy, int N, int NI, int det)
void copy_int_model(struct Model *origin, struct Model *copy, int N, int NI, int NW, int det)
{
int i,ifo;
......@@ -997,12 +1030,12 @@ void copy_int_model(struct Model *origin, struct Model *copy, int N, int NI, int
if(det==-1)
{
copy_wavelet(origin->signal, copy->signal, N);
copy_wavelet(origin->signal, copy->signal, N, NW);
for(ifo=0; ifo<NI; ifo++) for(i=0; i<N; i++) copy->response[ifo][i] = origin->response[ifo][i];
}
else
{
copy_wavelet(origin->glitch[det],copy->glitch[det], N);
copy_wavelet(origin->glitch[det],copy->glitch[det], N, NW);
}
for(ifo=0; ifo<NI; ifo++)
......@@ -1052,6 +1085,7 @@ void free_model(struct Model *model, struct Data *data, struct Prior *prior)
int N = data->N;
int NI = data->NI;
int halfN = data->N/2;
int NW = data->NW;