Commit 1fa769d6 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

Polarization branch recovers elliptic polarization behavior (at least for GW150914...)

parent 1134ed2b
......@@ -231,7 +231,7 @@ int main(int argc, char *argv[])
//give chain buffer in case we need more
int NC=chain->NC;
chain->NC+=10;
allocate_chain(chain,NI);
allocate_chain(chain,NI,data->Npol);
for(ic=0; ic<chain->NC; ic++) chain->index[ic]=ic;
/*
......@@ -990,7 +990,7 @@ int main(int argc, char *argv[])
free(model[ic]);
}
free(model);
free_chain(chain,NI);
free_chain(chain,NI,data->Npol);
return 0;
......
......@@ -191,7 +191,7 @@ struct Chain
FILE *runFile;
FILE **intChainFile;
FILE **intWaveChainFile;
FILE ***intWaveChainFile;
FILE ***intGlitchChainFile;
FILE ***lineChainFile;
FILE ***splineChainFile;
......@@ -246,6 +246,7 @@ struct Data
{
int N;
int NI;
int Npol;
int tsize;
int fsize;
int imin;
......
......@@ -380,23 +380,29 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
initialize_model(model, data, prior, psd, &chain->seed);
for(ifo=0; ifo<NI; ifo++)model->Sn[ifo]=1.0;
FILE **signalParams = malloc(data->Npol*sizeof(FILE *));
//File containing signal-model waveform parameters
sprintf(filename,"%s/%ssignal_params.dat.0",injpath,injname);
FILE *signalParams = fopen(filename,"r");
for(n=0; n<Nsample+1; n++) fgets(burnrows,100000,signalParams);
sprintf(filename,"%sbayeswave.run",data->runName);
fprintf(stdout,"%s\n",filename);
runFile = fopen(filename,"a");
fprintf(runFile, " ========== Signal Injection ========= \n");
fprintf(runFile,"%s\n",burnrows);
fclose(runFile);
rewind(signalParams);
for(n=0; n<Nsample; n++) fgets(burnrows,100000,signalParams);
for(n=0; n<data->Npol; n++)
{
sprintf(filename,"%s/%ssignal_params_h%i.dat.0",injpath,injname,n);
signalParams[n] = fopen(filename,"r");
for(n=0; n<Nsample+1; n++) fgets(burnrows,100000,signalParams[n]);
fprintf(runFile,"%s\n",burnrows);
rewind(signalParams[n]);
for(n=0; n<Nsample; n++) fgets(burnrows,100000,signalParams[n]);
}
fclose(runFile);
parse_signal_parameters(data, model, signalParams, hinj, ginj, ginj);
double SNRnet;
......@@ -425,7 +431,8 @@ void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, stru
//Add model to data
for(ifo=0; ifo<NI; ifo++) for (i=0; i<data->N; i++) data->s[ifo][i] += hinj[ifo][i];
fclose(signalParams);
for(n=0; n<data->Npol; n++)fclose(signalParams[n]);
free(signalParams);
free_model(model, data, prior);
}
......@@ -969,7 +976,7 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic)
{
int ifo;
int n;
int NI = data->NI;
struct Model *model_x = model[chain->index[ic]];
......@@ -980,18 +987,18 @@ void print_chain_files(struct Data *data, struct Chain *chain, struct Model **mo
print_model(chain->intChainFile[ic], data, chain, model_x);
//print signal model parameters
if(data->signalFlag) print_signal_model(chain->intWaveChainFile[ic], model_x);
if(data->signalFlag) for(n=0; n<model_x->Npol; n++) print_signal_model(chain->intWaveChainFile[n][ic], model_x, n);
//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(n=0; n<NI; n++) print_glitch_model(chain->intGlitchChainFile[n][ic], model_x->glitch[n]);
//print PSD parameters
if(data->bayesLineFlag)
{
for(ifo=0; ifo<NI; ifo++)
for(n=0; n<NI; n++)
{
print_line_model(chain->lineChainFile[ic][ifo],psd_x[ifo]);
print_spline_model(chain->splineChainFile[ic][ifo],psd_x[ifo]);
print_line_model(chain->lineChainFile[ic][n],psd_x[n]);
print_spline_model(chain->splineChainFile[ic][n],psd_x[n]);
}
}
......@@ -1000,24 +1007,25 @@ void print_chain_files(struct Data *data, struct Chain *chain, struct Model **mo
void flush_chain_files(struct Data *data, struct Chain *chain, int ic)
{
int ifo;
int n;
int NP = data->Npol;
int NI = data->NI;
fflush(chain->intChainFile[ic]);
//print signal model parameters
if(data->signalFlag) fflush(chain->intWaveChainFile[ic]);
if(data->signalFlag) for(n=0; n<NP; n++) fflush(chain->intWaveChainFile[n][ic]);
//print glitch model parameters
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) fflush(chain->intGlitchChainFile[ifo][ic]);
if(data->glitchFlag) for(n=0; n<NI; n++) fflush(chain->intGlitchChainFile[n][ic]);
//print PSD parameters
if(data->bayesLineFlag)
{
for(ifo=0; ifo<NI; ifo++)
for(n=0; n<NI; n++)
{
fflush(chain->lineChainFile[ic][ifo]);
fflush(chain->splineChainFile[ic][ifo]);
fflush(chain->lineChainFile[ic][n]);
fflush(chain->splineChainFile[ic][n]);
}
}
}
......@@ -1038,18 +1046,18 @@ void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Mode
fprintf(fptr,"\n");
}
void print_signal_model(FILE *fptr, struct Model *model)
void print_signal_model(FILE *fptr, struct Model *model, int n)
{
int i,j,k;
fprintf(fptr,"%i ",model->signal[0]->size);
fprintf(fptr,"%i ",model->signal[n]->size);
for(k=0; k<NE; k++) fprintf(fptr,"%.12g ", model->extParams[k]);
for(j=1;j<=model->signal[0]->size; j++) // print out first few sine-gaussians
for(j=1;j<=model->signal[n]->size; j++) // print out first few sine-gaussians
{
i = model->signal[0]->index[j];
i = model->signal[n]->index[j];
//print all parameters
for(k=0; k<model->signal[0]->dimension; k++) fprintf(fptr,"%.12g ", model->signal[0]->intParams[i][k]);
for(k=0; k<model->signal[n]->dimension; k++) fprintf(fptr,"%.12g ", model->signal[n]->intParams[i][k]);
}
fprintf(fptr,"\n");
}
......@@ -1834,36 +1842,40 @@ void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **para
}
}
void parse_signal_parameters(struct Data *data, struct Model *model, FILE *paramsfile, double **hrec, double **hrecPlus, double **hrecCross)
void parse_signal_parameters(struct Data *data, struct Model *model, FILE **paramsfile, double **hrec, double **hrecPlus, double **hrecCross)
{
int ifo;
int d;
int i;
int n;
//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;
//Get the number of wavelet basis functions for the current sample
fscanf(paramsfile, "%i", &signal->size);
//First in the wavechain file are the extrinsic parameters
for(i=0; i<NE; i++) fscanf(paramsfile,"%lg", &model->extParams[i]);
//Compute network projection coefficients
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++)
for(n=0; n<data->Npol; n++)
{
signal->index[d]=d;
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);
struct Wavelet *signal = model->signal[n];
//Zero out signal->templates array (holds the linear combination of wavelets)
for(i=0; i<data->N; i++) signal->templates[i] = 0.0;
//Get the number of wavelet basis functions for the current sample
fscanf(paramsfile[n], "%i", &signal->size);
//First in the wavechain file are the extrinsic parameters
for(i=0; i<NE; i++) fscanf(paramsfile[n],"%lg", &model->extParams[i]);
//Compute network projection coefficients
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++)
{
signal->index[d]=d;
for(i=0; i<data->NW; i++) fscanf(paramsfile[n],"%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
......@@ -1873,7 +1885,7 @@ void parse_signal_parameters(struct Data *data, struct Model *model, FILE *param
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);
//compute_reconstructed_plus_and_cross(data, model->projection, model->extParams, hrecPlus, hrecCross, signal->templates, data->fmin, data->fmax);
//fill hrec with reconstructed waveform
for(ifo=0; ifo<data->NI; ifo++) for(i=0; i<data->N; i++) {hrec[ifo][i] = model->response[ifo][i];}
......
......@@ -39,7 +39,7 @@ void print_chain_status(struct Data *data, struct Chain *chain, struct Model **m
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_signal_model(FILE *fptr, struct Model *model, int n);
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);
......@@ -50,7 +50,7 @@ void print_colored_time_domain_waveforms(char filename[], double *h, int N, doub
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);
void parse_signal_parameters(struct Data *data, struct Model *model, FILE *paramsfile, double **hrec, double **hrecPlus, double **hrecCross);
void parse_signal_parameters(struct Data *data, struct Model *model, FILE **paramsfile, double **hrec, double **hrecPlus, double **hrecCross);
void parse_bayesline_parameters(struct Data *data, struct Model *model, FILE **splinechain, FILE **linechain, double **psd);
void dump_time_domain_waveform(struct Data *data, struct Model *model, struct Prior *prior, char root[]);
......
......@@ -257,8 +257,11 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
if(data->signalFlag)
{
//intrinsic
sprintf(filename,"chains/%s_params.dat.%i",modelname,ic);
chain->intWaveChainFile[ic] = fopen(filename,"w");
for(n=0; n<data->Npol; n++)
{
sprintf(filename,"chains/%s_params_h%i.dat.%i",modelname,n,ic);
chain->intWaveChainFile[n][ic] = fopen(filename,"w");
}
}
//Parameters for reproducing PSD model
......@@ -281,7 +284,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
static void resume_sampler(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline)
{
//TODO: set internal mcmc flags right based on chain->mc
int i;
int i,n;
int ic,ifo;
int NI = data->NI;
......@@ -463,16 +466,21 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
// Signal model
if(data->signalFlag)
{
sprintf(filename,"checkpoint/params.dat.%i",ic);
if( (fptr = fopen(filename,"r")) == NULL)
fptr2 = malloc(data->Npol*sizeof(FILE *));
for(n=0; n<data->Npol; n++)
{
fprintf(stderr,"Error: Could not checkpoint signal model\n");
fprintf(stderr," Parameter file %s does not exist\n",filename);
fprintf(stderr," Exiting to system (1)\n");
exit(1);
sprintf(filename,"checkpoint/params_h%i.dat.%i",n,ic);
if( (fptr2[n] = fopen(filename,"r")) == NULL)
{
fprintf(stderr,"Error: Could not checkpoint signal model\n");
fprintf(stderr," Parameter file %s does not exist\n",filename);
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
parse_signal_parameters(data, model[ic], fptr2, h, h, h);
for(n=0; n<data->Npol; n++) fclose(fptr2[n]);
free(fptr2);
}
parse_signal_parameters(data, model[ic], fptr, h, h, h);
fclose(fptr);
}
// PSD model
......@@ -619,8 +627,11 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
if(data->signalFlag)
{
//intrinsic
sprintf(filename,"chains/%s_params.dat.%i",modelname,ic);
chain->intWaveChainFile[ic] = fopen(filename,"a");
for(n=0; n<data->Npol; n++)
{
sprintf(filename,"chains/%s_params_h%i.dat.%i",modelname,n,ic);
chain->intWaveChainFile[n][ic] = fopen(filename,"a");
}
}
//Parameters for reproducing PSD model
......@@ -658,7 +669,7 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
static void save_sampler(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, char modelname[])
{
int i;
int i,n;
int ic,ifo;
int NI = data->NI;
......@@ -778,10 +789,13 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
// Signal model
if(data->signalFlag)
{
sprintf(filename,"checkpoint/params.dat.%i",ic);
fptr = fopen(filename,"w");
print_signal_model(fptr, model[ic]);
fclose(fptr);
for(n=0; n<data->Npol; n++)
{
sprintf(filename,"checkpoint/params_h%i.dat.%i",n,ic);
fptr = fopen(filename,"w");
print_signal_model(fptr, model[ic],n);
fclose(fptr);
}
}
// PSD model
......@@ -859,7 +873,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
int M = chain->count;
int NC = chain->NC;
int i, ic, ifo;
int i, ic, ifo, n;
double logZ = 0.0;
double varZ = 0.0;
......@@ -1385,8 +1399,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if(data->glitchFlag || data->signalFlag)
{
fclose(chain->intChainFile[ic]);
if(data->signalFlag)
fclose(chain->intWaveChainFile[ic]);
if(data->signalFlag) for(n=0; n<data->Npol; n++)
fclose(chain->intWaveChainFile[n][ic]);
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++)
fclose(chain->intGlitchChainFile[ifo][ic]);
}
......@@ -1602,6 +1616,22 @@ static void add_wavelet(struct Wavelet *wave_x, struct Wavelet *wave_y, int *ii)
*ii = index;
}
void constrain_hplus_hcross(struct Wavelet **wave, int i)
{
/*
give hx special treatment (wavelets at same TFQ as h+, but independenta A,phi
*/
wave[1]->intParams[i][0] = wave[0]->intParams[i][0]; //t0
wave[1]->intParams[i][1] = wave[0]->intParams[i][1]; //f0
wave[1]->intParams[i][2] = wave[0]->intParams[i][2]; //Q
// printf("A+ = %g, Ax = %g\n",wave[0]->intParams[i][3],wave[1]->intParams[i][3]); //A
// printf("Phi+ = %g, Phix = %g\n",wave[0]->intParams[i][4],wave[1]->intParams[i][4]); //A
if(wave[0]->dimension > 5)
{
wave[1]->intParams[i][5] = wave[0]->intParams[i][5]; //beta (~fdot for chirplets)
}
}
/**************************************************************/
void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Model **model, struct Chain *chain, struct TimeFrequencyMap *tf, long *seed, int ic)
......@@ -1883,6 +1913,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
paramsy = wave_y[n]->intParams;
draw_wavelet_params(paramsy[ii], range, seed, wave_y[n]->dimension);
}
if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
/* TF density for proposal */
if(data->clusterProposalFlag )
......@@ -1898,7 +1929,8 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
else draw_time_frequency(det, ii, wave_x[n], wave_y[n], range, seed, TFtoProx, tf, &k);
}
if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
if(ic==0)
{
if(k==0)
......@@ -1987,6 +2019,13 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
else for(n=nmin; n<nmax; n++) logqy += -log(prior->range[3][1]-prior->range[3][0]);
//check priors
for(n=nmin; n<nmax; n++)
{
paramsy = wave_y[n]->intParams;
test += checkrange(paramsy[ii],prior->range, wave_y[n]->dimension);
}
}//end birth move
}//end test condition
}//end trans-dimensional update
......@@ -2016,11 +2055,9 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
i = (int)(ran2(seed)*(double)(wave_x[nmin]->size))+1; // pick a term to update
ii = wave_x[nmin]->index[i];
alpha = ran2(seed);
// replace existing wavelet (birth/death move)
//TODO: Birth/Death move is currently disabled
if(alpha > 1.0)//chain->intPropRate/chain->temperature[ic])
if(ran2(seed) > 1.0)//chain->intPropRate/chain->temperature[ic])
{
// Draw all wavelet parameters from prior
......@@ -2100,6 +2137,14 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
intrinsic_halfcycle_proposal(paramsx[ii],paramsy[ii],seed);
}
if(det==-1 && model_y->Npol > 1)
{
//give hx the same phase shift as h+
wave_y[1]->intParams[ii][4] = wave_x[1]->intParams[ii][4] + (wave_y[0]->intParams[ii][4]-wave_x[0]->intParams[ii][4]);
constrain_hplus_hcross(wave_y, ii);
}
}
/*
......@@ -2132,6 +2177,8 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], 1.0, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
}
}
if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
}
/*
Draw from prior:
......@@ -2159,8 +2206,15 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
}
}
if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
}
}//end perturb step
//check priors
for(n=nmin; n<nmax; n++)
{
paramsy = wave_y[n]->intParams;
test += checkrange(paramsy[ii],prior->range, wave_y[n]->dimension);
}
}//end model/dimension check
else test = 1;
}//end fixed-dimension update
......
......@@ -251,7 +251,7 @@ void initialize_chain(struct Chain *chain, int flag)
chain->burnFlag = flag; //Flag for burn-in (1 = yes)
}
void free_chain(struct Chain *chain, int NI)
void free_chain(struct Chain *chain, int NI, int NP)
{
free(chain->ptacc);
free(chain->ptprop);
......@@ -265,19 +265,20 @@ void free_chain(struct Chain *chain, int NI)
free(chain->intChainFile);
int n;
for(n=0; n<NP; n++) free(chain->intWaveChainFile[n]);
free(chain->intWaveChainFile);
int ifo;
for(ifo=0; ifo<NI; ifo++) free(chain->intGlitchChainFile[ifo]);
for(n=0; n<NI; n++) free(chain->intGlitchChainFile[n]);
free(chain->intGlitchChainFile);
//logL chain for thermodynamic integration
int i;
for(i=0; i<chain->NC; i++) free(chain->logLchain[i]);
for(n=0; n<chain->NC; n++) free(chain->logLchain[n]);
free(chain->logLchain);
}
void allocate_chain(struct Chain *chain, int NI)
void allocate_chain(struct Chain *chain, int NI, int NP)
{
chain->ptacc = malloc(chain->NC*sizeof(int));
chain->ptprop = malloc(chain->NC*sizeof(int));
......@@ -291,7 +292,10 @@ void allocate_chain(struct Chain *chain, int NI)
chain->intChainFile = malloc(chain->NC*sizeof(FILE*));
chain->intWaveChainFile = malloc(chain->NC*sizeof(FILE*));
int n;
chain->intWaveChainFile = malloc(NP*sizeof(FILE**));
for(n=0; n<NP; n++) chain->intWaveChainFile[n] = malloc(chain->NC*sizeof(FILE*));
int ifo;
chain->intGlitchChainFile = malloc(NI*sizeof(FILE**));
......@@ -307,7 +311,7 @@ void allocate_chain(struct Chain *chain, int NI)
}
//logL chain for thermodynamic integration
int i,n;
int i;
chain->nPoints = 3*chain->count/4/chain->cycle;
chain->logLchain = malloc(chain->NC*sizeof(double *));
for(i=0; i<chain->NC; i++)
......@@ -340,9 +344,9 @@ void resize_model(struct Data *data, struct Chain *chain, struct Prior *prior, s
}
//choose new number of chains
free_chain(chain, data->NI);
free_chain(chain, data->NI, data->Npol);
chain->NC = NC;
allocate_chain(chain, data->NI);
allocate_chain(chain, data->NI, data->Npol);
for(ic=0; ic<chain->NC; ic++) chain->index[ic]=ic;
*model = realloc(*model, chain->NC * sizeof(struct Model*) );
......@@ -680,6 +684,10 @@ void initialize_data(struct Data *data, double **s, int N, int tsize, double Tob
data->r = dmatrix(0,NI-1,0,N-1); // residual
data->rh = dmatrix(0,NI-1,0,N-1);
// number of signal polarizations (0==h+, 1==hx, ...)
data->Npol = 2;
if(data->polarizationFlag) data->Npol = 1;
// initialize data vectors
for(i=0; i<N; i++)
......@@ -726,7 +734,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 Npol = 2; // (+) and (x) polarization
int Npol = data->Npol;
int NW = data->NW; //number of intrinsic parameters / frame (5 for wavelet)
......@@ -758,7 +766,7 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
model->response = dmatrix(0,NI-1,0,N-1);
//incident signal model
model->h = dmatrix(0,Npol-1,0,N-1);
model->h = dmatrix(0,Npol,0,N-1);
//likelhood is a pointer to a function which returns logL
model->logL = 0.0;
......@@ -1020,12 +1028,9 @@ void copy_int_model(struct Model *origin, struct Model *copy, int N, int NI, int
if(det==-1)
{
for(ip=0; ip<origin->Npol; ip++)
{
copy_wavelet(origin->signal[ip], copy->signal[ip], N);
for(i=0; i<N; i++) copy->h[ip][i] = origin->h[ip][i];
}
for(ifo=0; ifo<NI; ifo++) for(i=0; i<N; i++) copy->response[ifo][i] = origin->response[ifo][i];
for(ip=0; ip<origin->Npol; ip++) copy_wavelet(origin->signal[ip], copy->signal[ip], N);
for(ip=0; ip<=origin->Npol; ip++) for(i=0; i<N; i++) copy->h[ip][i] = origin->h[ip][i];
for(ifo=0; ifo<NI; ifo++) for(i=0; i<N; i++) copy->response[ifo][i] = origin->response[ifo][i];
}
else
{
......@@ -1091,7 +1096,7 @@ void free_model(struct Model *model, struct Data *data, struct Prior *prior)
free_dvector(model->detLogL,0,