Commit 30492a87 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

paving the way for relaxed polarization constraint

parent 270766d3
......@@ -1546,19 +1546,19 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
xsm =0.0;
pmax = -1.0;
for(i=0; i< ncut; i++)
{
{
x = power[i]/Sbase[i];
if(x < 10.0) x = 1.0;
if(x >= 10.0) x = 100.0;
if(x > pmax)
{
pmax = x;
k = i;
}
if(x < 10.0) x = 1.0;
if(x >= 10.0) x = 100.0;
fprop[i] = x;
xsm += x;
}
// define the focus region (only used if focus flag = 1)
fcl = freq[k]-dff;
fch = freq[k]+dff;
......
......@@ -1070,6 +1070,7 @@ void print_help_message(void)
fprintf(stdout," --noSignalAmplitudePrior use same SNR-dependent prior for signal & glitch model\n");
fprintf(stdout," --noAmplitudeProposal don't draw from SNR-dependent amplitude prior\n");
fprintf(stdout," --varyExtrinsicAmplitude update wavelet amplitudes in extrinsic update (FIXED?)\n");
fprintf(stdout," --ellipticalPolarization assume elliptical polarization for signal model\n");
fprintf(stdout," --noClusterProposal disable clustering & TF density proposal \n");
fprintf(stdout," --clusterWeight fractional weight for TF to proximity proposal (0.5)\n");
fprintf(stdout," --ampPriorPeak SNR where amplitude prior peaks (5)\n");
......
......@@ -56,7 +56,6 @@
#define CLIGHT 2.99792458e8 // Speed of light (m/s)
#define NE 6 //Number of extrinsic parameters
// #define NW 5 //Number of intrinsic parameters
//Tuning parameters for wavelet proximity proposal
//TODO: These should not be global!
......@@ -75,6 +74,7 @@ struct Wavelet
int size;
int smax;
int *index;
int dimension;
//clustering prior
int Ncluster;
......@@ -90,6 +90,7 @@ struct Wavelet
struct Model
{
int size;
int Npol;
double *Sn;
double **Snf;
......@@ -101,12 +102,13 @@ struct Model
double logLnorm;
double *detLogL;
double *extParams;
double **h;
double **response;
int chirpletFlag;
int NW;
struct Wavelet *signal;
struct Wavelet **signal;
struct Wavelet **glitch;
struct Network *projection;
......@@ -269,6 +271,7 @@ struct Data
int signalFlag;
int noiseFlag;
int polarizationFlag;
int fixIntrinsicFlag;
int fixExtrinsicFlag;
int geocenterPSDFlag;
......@@ -338,7 +341,7 @@ struct Data
double (*glitch_amplitude_prior) (double *, double *, double, double);
double (*intrinsic_likelihood)(int, int, int, struct Model*, struct Model*, struct Prior*, struct Chain*, struct Data*);
double (*extrinsic_likelihood)(struct Network*, double *, double **, double *, double *, double **, struct Data*, double, double);
double (*extrinsic_likelihood)(struct Network*, double *, double **, double *, struct Wavelet **, double **, struct Data*, double, double);
char runName[100];
char outputDirectory[1000];
......
......@@ -2324,7 +2324,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
/*********************************/
logP = model->logL;
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) logP += model->glitch[ifo]->logp;
if(data->signalFlag) logP += model->signal->logp;
if(data->signalFlag) logP += model->signal[0]->logp;
/*********************************/
/* */
......@@ -2336,11 +2336,11 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
extrinsic_fisher_update(data, model);
for(j=0; j<4; j++) detC -= log(model->fisher->evalue[j])/2.;
//printf(" MAP signal model dimension: %i wavelets\n",model->signal->size);
for(i=1; i<=model->signal->size; i++)
//printf(" MAP signal model dimension: %i wavelets\n",model->signal[0]->size);
for(i=1; i<=model->signal[0]->size; i++)
{
ii = model->signal->index[i];
intrinsic_fisher_update(model->signal->intParams[ii], sigmas, model->SnGeo, 1.0, data->Tobs, data->NW);
ii = model->signal[0]->index[i];
intrinsic_fisher_update(model->signal[0]->intParams[ii], sigmas, model->SnGeo, 1.0, data->Tobs, data->NW,1 );
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
}
......@@ -2353,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, data->NW);
intrinsic_fisher_update(model->glitch[ifo]->intParams[ii], sigmas, model->Snf[ifo], 1.0, data->Tobs, data->NW, 1);
//for(j=0; j<NW; j++)printf(" %i: %g\n",j,sigmas[j]);
for(j=0; j<NW; j++) detC += log(sigmas[j]);
}
......@@ -2381,7 +2381,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
logV += log(2.); //ecc
//intrinsic
dim+=NW*model->signal->size;
dim+=NW*model->signal[0]->size;
logVsignal = 0;
......@@ -2394,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,data->NW);
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];
......@@ -2403,7 +2403,7 @@ void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *
logVsignal+=log(maxAmp-minAmp); //amplitude
logVsignal+=log(prior->range[4][1]-prior->range[4][0]); //phase
logV += model->signal->size*logVsignal;
logV += model->signal[0]->size*logVsignal;
}
if(data->glitchFlag)
{
......@@ -2421,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,data->NW);
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];
......
......@@ -51,7 +51,7 @@ static void Lorentzian(double *Slines, double Tobs, lorentzianParams *lines, int
// add or remove the old line
f2=f*f;
f4=f2*f2;
amplitude = A*f4/(f2*Q*Q);
amplitude = A*f4;//(f2*Q*Q);
for(i=istart; i<istop; i++)
{
freq = (double)i/Tobs;
......@@ -60,7 +60,8 @@ static void Lorentzian(double *Slines, double Tobs, lorentzianParams *lines, int
z = 1.0;
if(x > deltafmax) z = exp(-(x-deltafmax)/deltafmax);
Slines[i] += z*amplitude/(fsq*(fsq-f2)*(fsq-f2));
//Slines[i] += z*amplitude/(fsq*(fsq-f2)*(fsq-f2));
Slines[i] += z*amplitude/(f2*fsq+Q*Q*(fsq-f2)*(fsq-f2));
}
}
......@@ -909,6 +910,11 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " polarization constraint is ");
if(data->polarizationFlag) fprintf(fptr, "ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " background prior is ...... ");
if(data->backgroundPriorFlag && data->glitchFlag) fprintf(fptr, "ENABLED");
else fprintf(fptr, "DISABLED");
......@@ -952,11 +958,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, " chirplet option is ...... ");
if(data->chirpletFlag) fprintf(fptr,"ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, "\n");
}
......@@ -977,7 +983,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], data->NW);
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) print_glitch_model(chain->intGlitchChainFile[ifo][ic], model_x->glitch[ifo]);
//print PSD parameters
if(data->bayesLineFlag)
......@@ -1022,7 +1028,7 @@ void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Mode
int NI = data->NI;
fprintf(fptr,"%d %f ", chain->mc, model->logL-data->logLc+model->logLnorm);
fprintf(fptr,"%d ", model->signal->size);
fprintf(fptr,"%d ", model->signal[0]->size);
for(ifo=0; ifo<NI; ifo++) fprintf(fptr,"%d ", model->glitch[ifo]->size);
for(ifo=0; ifo<NI; ifo++) fprintf(fptr, "%f ",model->Sn[ifo]);
// fprintf(fptr, "%f %f ",logZ, sqrt(varZ));
......@@ -1035,21 +1041,20 @@ 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);
fprintf(fptr,"%i ",model->signal[0]->size);
for(k=0; k<NE; k++) fprintf(fptr,"%.12g ", model->extParams[k]);
for(j=1;j<=model->signal->size; j++) // print out first few sine-gaussians
for(j=1;j<=model->signal[0]->size; j++) // print out first few sine-gaussians
{
i = model->signal->index[j];
i = model->signal[0]->index[j];
//print all parameters
for(k=0; k< NW; k++) fprintf(fptr,"%.12g ", model->signal->intParams[i][k]);
for(k=0; k<model->signal[0]->dimension; k++) fprintf(fptr,"%.12g ", model->signal[0]->intParams[i][k]);
}
fprintf(fptr,"\n");
}
void print_glitch_model(FILE *fptr, struct Wavelet *glitch, int NW)
void print_glitch_model(FILE *fptr, struct Wavelet *glitch)
{
int i,j,k;
fprintf(fptr,"%i ",glitch->size);
......@@ -1058,7 +1063,7 @@ void print_glitch_model(FILE *fptr, struct Wavelet *glitch, int NW)
i = glitch->index[j];
//print all parameters
for(k=0; k< NW; k++) fprintf(fptr,"%.12g ", glitch->intParams[i][k]);
for(k=0; k<glitch->dimension; k++) fprintf(fptr,"%.12g ", glitch->intParams[i][k]);
}
fprintf(fptr,"\n");
}
......@@ -1088,7 +1093,7 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
//print intrinsic parameter stats
fprintf(stdout, " DIM: ");
fprintf(stdout, "DGW=%d ", model[chain->index[0]]->signal->size);
fprintf(stdout, "DGW=%d ", model[chain->index[0]]->signal[0]->size);
for(ifo=0; ifo<NI; ifo++) fprintf(stdout, "D%i=%d ", ifo, model[chain->index[0]]->glitch[ifo]->size);
fprintf(stdout, "m0=%i m1=%i m2=%i m3=%i ", chain->mod0, chain->mod1, chain->mod2, chain->mod3);
fprintf(stdout, "\n");
......@@ -1165,7 +1170,7 @@ void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf,
fclose(waveout);
}
void print_colored_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_colored_time_domain_waveforms(char filename[], double *h, int N, double Tobs, int imin, int imax, double tmin, double tmax)
{
/*
TODO: invFFT comes out time-reversed.
......@@ -1186,9 +1191,7 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub
{
if(i>imin && i<imax)
{
x = 1.0;//sqrt(Snf[i]*eta);
ht[2*i] /= x;
ht[2*i+1] /= -x;
ht[2*i+1] *= -1.0;
}
else
{
......@@ -1500,6 +1503,10 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
}
else data->backgroundPriorFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--ellipticalPolarization");
if(ppt) data->polarizationFlag = 1;
else data->polarizationFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--fixExtrinsicParams");
if(ppt) data->fixExtrinsicFlag = 1;
else data->fixExtrinsicFlag = 0;
......@@ -1797,7 +1804,6 @@ 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;
......@@ -1815,7 +1821,7 @@ void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **para
for(d=1; d<=glitch[ifo]->size; d++)
{
glitch[ifo]->index[d]=d;
for(i=0; i<NW; i++)
for(i=0; i<data->NW; i++)
fscanf(paramsfile[ifo],"%lg ", &glitch[ifo]->intParams[d][i]);
//Append new wavelet to glitch[ifo]->templates array
......@@ -1833,9 +1839,9 @@ 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;
//HACKED SIGNAL PARAMETERS
struct Wavelet *signal = model->signal[0];
//Zero out signal->templates array (holds the linear combination of wavelets)
for(i=0; i<data->N; i++) signal->templates[i] = 0.0;
......@@ -1853,15 +1859,18 @@ void parse_signal_parameters(struct Data *data, struct Model *model, FILE *param
for(d=1; d<=signal->size; d++)
{
signal->index[d]=d;
for(i=0; i<NW; i++) fscanf(paramsfile,"%lg ", &signal->intParams[d][i]);
for(i=0; i<data->NW; i++) fscanf(paramsfile,"%lg ", &signal->intParams[d][i]);
//Append new wavelet to signal->templates array
model->wavelet(signal->templates, signal->intParams[d], data->N, 1, data->Tobs);
}
//Combine template waveforms into polarization array
combinePolarizations(data, model->signal, model->h, model->extParams, model->Npol);
//Project geocenter waveform onto network
waveformProject(data, model->projection, model->extParams, model->response, signal->templates, data->fmin, data->fmax);
waveformProject(data, model->projection, model->extParams, model->response, model->h, data->fmin, data->fmax);
//Get + and x polarizations of hrec
compute_reconstructed_plus_and_cross(data, model->projection, model->extParams, hrecPlus, hrecCross, signal->templates, data->fmin, data->fmax);
......@@ -1877,7 +1886,8 @@ void dump_time_domain_waveform(struct Data *data, struct Model *model, struct Pr
int ifo;
char filename[128];
struct Wavelet *signal = model->signal;
//TODO: HACKED signal pointer
struct Wavelet *signal = model->signal[0];
//Zero out signal->templates array (holds the linear combination of wavelets)
for(i=0; i<data->N; i++) signal->templates[i] = 0.0;
......@@ -1893,8 +1903,11 @@ void dump_time_domain_waveform(struct Data *data, struct Model *model, struct Pr
}
//Combine template waveforms into polarization array
combinePolarizations(data, model->signal, model->h, model->extParams, model->Npol);
//Project geocenter waveform onto network
waveformProject(data, model->projection, model->extParams, model->response, signal->templates, data->fmin, data->fmax);
waveformProject(data, model->projection, model->extParams, model->response, model->h, data->fmin, data->fmax);
for(ifo=0; ifo<data->NI; ifo++)
{
......
......@@ -40,13 +40,13 @@ 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, int NW);
void print_glitch_model(FILE *fptr, struct Wavelet *glitch);
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);
void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax);
void print_frequency_domain_data(char filename[], double *h, int N, double Tobs, int imin, int imax);
void print_colored_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_colored_time_domain_waveforms(char filename[], double *h, int N, double Tobs, int imin, int imax, double tmin, double tmax);
void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *prior, ProcessParamsTable *commandLine);
void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **paramsfile, double **grec);
......
This diff is collapsed.
......@@ -18,13 +18,12 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
double Sum_Extreme(double **a, double **b, double **Sn, int n, double *delt, double *pshift, double Tobs, int NI, int imin, int imax, double t0);
double EvaluateSearchLogLikelihood (int typ, int ii, int det, struct Model *mx, struct Model *my, struct Prior *prior, struct Chain *chain, struct Data *data);
double EvaluateConstantLogLikelihood (int typ, int ii, int det, UNUSED struct Model *mx, struct Model *my, struct Prior *prior, UNUSED struct Chain *chain, UNUSED struct Data *data);
double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx, struct Model *my, struct Prior *prior, struct Chain *chain, struct Data *data);
double EvaluateExtrinsicSearchLogLikelihood (struct Network *projection, double *params, double **invSnf, double *Sn, double *geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED double *Sn, UNUSED double *geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, double *geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicSearchLogLikelihood (struct Network *projection, double *params, double **invSnf, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
double EvaluateExtrinsicConstantLogLikelihood (UNUSED struct Network *projection, double *params, UNUSED double **invSnf, UNUSED double *Sn, UNUSED struct Wavelet **geo, UNUSED double **g, UNUSED struct Data *data, UNUSED double fmin, UNUSED double fmax);
double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, double *params, double **invSnf, double *Sn, struct Wavelet **geo, double **g, struct Data *data, double fmin, double fmax);
......
This diff is collapsed.
......@@ -42,12 +42,12 @@ void Shf_Geocenter(struct Data *data, struct Model *model, double *SnGeo, double
double AntennaPattern;
double ecc = params[3];
int nmax = model->signal->size;
int index[model->signal->size];
int nmax = model->signal[0]->size;
int index[model->signal[0]->size];
for(n=0; n<nmax; n++)
{
i = (int)(model->signal->intParams[n+1][1]*data->Tobs);
i = (int)(model->signal[0]->intParams[n+1][1]*data->Tobs);
index[n] = i;
SnGeo[i] = 0.0;
}
......@@ -366,9 +366,8 @@ 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;
int NW = data->NW;
/* Model dimension parameters */
// number of sineGaussians we print to file
prior->omax = omax;
......@@ -382,7 +381,7 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
double b = 2.9;
//TODO Nwavelet unstable at i=0?
for(i=1; i<200; i++) prior->Nwavelet[i] = log(((double)i*4.0*sqrt(3.)/(LAL_PI*pow(b,2.)))/(3.0+pow((double)i/b,4.)));
prior->Nwavelet[1]=prior->Nwavelet[0];
prior->Nwavelet[0]=prior->Nwavelet[1];
prior->smin = 0;
prior->smax = data->Dmax;
......@@ -440,8 +439,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;
double betamin = -1.0;
double betamax = 1.0;
/*
if(data->amplitudePriorFlag)
......@@ -457,12 +456,13 @@ void initialize_priors(struct Data *data, struct Prior *prior, int omax)
prior->range[1][1] = fmax;
prior->range[2][0] = Qmin;
prior->range[2][1] = Qmax;
prior->range[3][0] = Amin;
prior->range[3][1] = Amax;
prior->range[3][0] = Amin;
prior->range[3][1] = Amax;
prior->range[4][0] = 0.0;
prior->range[4][1] = LAL_TWOPI;
if (NW==6) {
if(data->chirpletFlag)
{
prior->range[5][0] = betamin;
prior->range[5][1] = betamax;
}
......@@ -620,8 +620,8 @@ 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;
double betamin = -1.0;
double betamax = 1.0;
//restrict prior range for wavelet time to 1 second around trigger
......@@ -647,6 +647,11 @@ void reset_priors(struct Data *data, struct Prior *prior)
prior->range[3][1] = Amax;
prior->range[4][0] = 0.0;
prior->range[4][1] = LAL_TWOPI;
if(data->chirpletFlag)
{
prior->range[5][0] = betamin;
prior->range[5][1] = betamax;
}
prior->TFV = (tmax-tmin)*(fmax-fmin)*(Qmax-Qmin);
prior->logTFV = log(prior->TFV);
......@@ -721,7 +726,13 @@ 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;
int Npol = 2; // (+) and (x) polarization
int NW = data->NW; //number of intrinsic parameters / frame (5 for wavelet)
// if(data->polarizationFlag) Npol = 1;
// else Npol = 2;
//gaussian noise model
model->Sn = dvector(0,NI-1);
model->SnS = dmatrix(0,NI-1,0,halfN-1);
......@@ -743,13 +754,15 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
model->SnGeo[i] = 1./model->SnGeo[i];
}
//extrinsic parameters
model->extParams = dvector(0,NE);
//full instrument response model (projected geocenter + glitches)
model->response = dmatrix(0,NI-1,0,N-1);
//incident signal model
model->h = dmatrix(0,Npol-1,0,N-1);
//likelhood is a pointer to a function which returns logL
model->logL = 0.0;
......@@ -765,13 +778,25 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//allocate memory for signal wavelet structure
model->signal = malloc(sizeof(struct Wavelet));
/*
New for O3: signal structure is a 2-element vector for (+) and (x) polarizations
*/
model->Npol = Npol;
//Setup wavelet structure holding signal model
if(data->signalFlag)
initialize_wavelet(model->signal,N, NW,prior->smax, 1);
else
initialize_wavelet(model->signal,N, NW,prior->smax, 0);
model->signal = malloc(Npol * sizeof(struct Wavelet *));
for(i=0; i<Npol; i++)
{
model->signal[i] = malloc(sizeof(struct Wavelet));
//assign wavelet dimension for signal model
model->signal[i]->dimension = NW; //wavelets
//Setup wavelet structure holding signal model
if(data->signalFlag)
initialize_wavelet(model->signal[i],N, prior->smax, 1);
else
initialize_wavelet(model->signal[i],N, prior->smax, 0);
}
//allocate memory for structure holding projection coefficients
model->projection = malloc(sizeof(struct Network));
......@@ -781,13 +806,18 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//Allocate and initialize glitch structures
model->glitch = malloc(NI * sizeof(struct Wavelet *));
for(i=0; i<NI; i++) model->glitch[i] = malloc(sizeof(struct Wavelet));
for(i=0; i<NI; i++)
{
//assign wavelet dimension for glitch model
model->glitch[i]->dimension = NW; //wavelets
if(data->glitchFlag)
initialize_wavelet(model->glitch[i],N, NW,prior->smax, 1);
initialize_wavelet(model->glitch[i],N, prior->smax, 1);
else
initialize_wavelet(model->glitch[i],N, NW,prior->smax, 0);
initialize_wavelet(model->glitch[i],N, prior->smax, 0);
}
//allocate memory for structure holding extrinsic fisher matrix
......@@ -799,16 +829,17 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
initialize_fisher(model->intfisher, NW);
//select wavelet basis function
if (data->chirpletFlag) {
model->wavelet = ChirpletFourier;
model->wavelet_bandwidth = ChirpletBandwidth;
}
else{
model->wavelet = SineGaussianFourier;
model->wavelet_bandwidth = SineGaussianBandwidth;
}