Commit 7546f0a3 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

adding copyright info to src/*

parent 0b7ca88a
/*
* Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
/*********************************************************************************/
/* */
/* BayesLine fits the LIGO/Virgo power spectra using a model made up */
......@@ -93,7 +112,7 @@ static double logprior(double *lower, double *upper, double *Snf, int ilow, int
{
if(Snf[i]>upper[i] || Snf[i]<lower[i])
{
lgp=-1e60;
lgp=-1e60;
}
}
return (lgp);
......@@ -641,15 +660,15 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
s3 = 0.5;
s4 = 0.5;
for(i=0; i<lines_x->n; i++)
{
lines_x->larray[i] = i;
lines_y->larray[i] = i;
}
// FILE *temp=fopen("start_psd.dat","w");
// for(i=0; i<ncut; i++) fprintf(temp,"%i %lg %lg %lg\n",i,Snf[i],priors->lower[i],priors->upper[i]);
// fclose(temp);
for(i=0; i<lines_x->n; i++)
{
lines_x->larray[i] = i;
lines_y->larray[i] = i;
}
// FILE *temp=fopen("start_psd.dat","w");
// for(i=0; i<ncut; i++) fprintf(temp,"%i %lg %lg %lg\n",i,Snf[i],priors->lower[i],priors->upper[i]);
// fclose(temp);
baseav = 0.0;
......@@ -668,7 +687,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
full_spectrum_spline(Sline, Sbase, freq, data, lines_x);
for(i=0; i< ncut; i++) Sn[i] = Sbase[i]+Sline[i];
for(i=0; i<ncut; i++) Snx[i] = Sn[i];
for(i=0; i<ncut; i++) Snx[i] = Sn[i];
if(!bayesline->constantLogLFlag)
logLx = loglike(power, Sn, ncut);
......@@ -683,14 +702,14 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
{
logPsx = logprior_gaussian(priors->mean, priors->sigma, Snx, 0, ncut);
}
if(logPsx<0)
{
printf("how did PSD come in outside of the prior?\n");
FILE *temp=fopen("failed_psd.dat","w");
for(i=0; i<ncut; i++) fprintf(temp,"%i %lg %lg %lg\n",i,Snx[i],priors->lower[i],priors->upper[i]);
exit(1);
if(logPsx<0)
{
printf("how did PSD come in outside of the prior?\n");
FILE *temp=fopen("failed_psd.dat","w");
for(i=0; i<ncut; i++) fprintf(temp,"%i %lg %lg %lg\n",i,Snx[i],priors->lower[i],priors->upper[i]);
exit(1);
}
// set up proposal for frequency jumps
......@@ -1163,7 +1182,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
full_spectrum_spline(Sline, Sbase, freq, data, lines_y);
for(i=0; i< ncut; i++) Sn[i] = Sbase[i]+Sline[i];
logLy = loglike(power, Sn, ncut);
if(bayesline->constantLogLFlag) logLy = 1.0;
if(bayesline->constantLogLFlag) logLy = 1.0;
}
else logLy = -1e60;
}
......@@ -1172,23 +1191,23 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
{
full_spectrum_single(Sn, Snx, Sbasex, freq, data, lines_x, lines_y, ii, &ilowx, &ihighx, &ilowy, &ihighy);
logLy = logLx + loglike_single(power, Sn, Snx, ilowx, ihighx, ilowy, ihighy);
if(bayesline->constantLogLFlag) logLy = 1.0;
if(bayesline->constantLogLFlag) logLy = 1.0;
}
if(typ == 2) // add new line with label ii
{
full_spectrum_add_or_subtract(Sn, Snx, Sbasex, freq, data, lines_y, ii, &ilowy, &ihighy, 1);
logLy = logLx + loglike_pm(power, Sn, Snx, ilowy, ihighy);
if(bayesline->constantLogLFlag) logLy = 1.0;
if(bayesline->constantLogLFlag) logLy = 1.0;
}
if(typ == 3) // remove line with label ki
{
full_spectrum_add_or_subtract(Sn, Snx, Sbasex, freq, data, lines_x, ki, &ilowx, &ihighx, -1);
logLy = logLx + loglike_pm(power, Sn, Snx, ilowx, ihighx);
if(bayesline->constantLogLFlag) logLy = 1.0;
if(bayesline->constantLogLFlag) logLy = 1.0;
}
// prior on line number e(-ZETA * n). (this is a prior, not a proposal, so opposite sign)
// effectively an SNR cut on lines
......@@ -1215,7 +1234,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
alpha = log(gsl_rng_uniform(r));
if(logH > alpha)
{
{
if(typ == 0) ac0++;
if(typ == 1) ac1++;
if(typ == 4) ac2++;
......@@ -1245,54 +1264,54 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
}//end prior check
//Every 1000 steps update focus region
if(mc%1000 == 0)
{
pmax = -1.0;
for(i=0; i< ncut; i++)
{
x = power[i]/Snx[i];
if(x > pmax)
{
pmax = x;
k = i;
}
}
// define the focus region (only used if focus flag = 1)
fcl = freq[k]-dff;
fch = freq[k]+dff;
Ac = power[k];
if(fcl < freq[0]) fcl = freq[0];
if(fch > freq[ncut-1]) fch = freq[ncut-1];
//if(focus==1)fprintf(stdout,"Focusing on [%f %f] Hz...\n", fcl, fch);
// set up proposal for frequency jumps
xsm =0.0;
baseav = 0.0;
for(i=0; i< ncut; i++)
{
x = power[i]/Snx[i];
if(x < 16.0) x = 1.0;
if(x >= 16.0) x = 100.0;
fprop[i] = x;
xsm += x;
baseav += log(Sbasex[i]);
}
baseav /= (double)(ncut);
pmax = -1.0;
for(i=0; i< ncut; i++)
{
fprop[i] /= xsm;
if(fprop[i] > pmax) pmax = fprop[i];
}
}//end focus update
//Every 1000 steps update focus region
if(mc%1000 == 0)
{
pmax = -1.0;
for(i=0; i< ncut; i++)
{
x = power[i]/Snx[i];
if(x > pmax)
{
pmax = x;
k = i;
}
}
// define the focus region (only used if focus flag = 1)
fcl = freq[k]-dff;
fch = freq[k]+dff;
Ac = power[k];
if(fcl < freq[0]) fcl = freq[0];
if(fch > freq[ncut-1]) fch = freq[ncut-1];
//if(focus==1)fprintf(stdout,"Focusing on [%f %f] Hz...\n", fcl, fch);
// set up proposal for frequency jumps
xsm =0.0;
baseav = 0.0;
for(i=0; i< ncut; i++)
{
x = power[i]/Snx[i];
if(x < 16.0) x = 1.0;
if(x >= 16.0) x = 100.0;
fprop[i] = x;
xsm += x;
baseav += log(Sbasex[i]);
}
baseav /= (double)(ncut);
pmax = -1.0;
for(i=0; i< ncut; i++)
{
fprop[i] /= xsm;
if(fprop[i] > pmax) pmax = fprop[i];
}
}//end focus update
}//End MCMC loop
......@@ -1310,7 +1329,7 @@ void BayesLineLorentzSplineMCMC(struct BayesLineParams *bayesline, double heat,
// return updated spectral model
for(i=0; i< ncut; i++) Snf[i] = Snx[i];
// re-map the array to 0..mx ordering
for(i=0; i< lines_x->n; i++)
{
......@@ -1516,42 +1535,42 @@ void BayesLineInitialize(struct BayesLineParams *bayesline)
/* */
/******************************************************************************/
int k;
double mdn;
double x,z;
/* Make local pointers to BayesLineParams structure members */
dataParams *data = bayesline->data;
splineParams *spline = bayesline->spline;
double *Sna = bayesline->Sna;
double *fa = bayesline->fa;
//initialize spline grid & PSD values
k=0;
for(j=0; j<bayesline->spline->n; j++)
{
x = data->fmin+(data->fhigh-data->fmin)*(double)(j)/(double)(bayesline->spline->n-1);
if (x <= fa[0]) z = Sna[0];
else if(x >= fa[jj-1]) z = Sna[jj-1];
else
{
i=k-10;
mdn = 0.0;
do
{
if(i>=0) mdn = fa[i];
i++;
} while(x > mdn);
k = i;
z = Sna[i];
}
spline->points[j] = x;
spline->data[j] = z;
}
int k;
double mdn;
double x,z;
/* Make local pointers to BayesLineParams structure members */
dataParams *data = bayesline->data;
splineParams *spline = bayesline->spline;
double *Sna = bayesline->Sna;
double *fa = bayesline->fa;
//initialize spline grid & PSD values
k=0;
for(j=0; j<bayesline->spline->n; j++)
{
x = data->fmin+(data->fhigh-data->fmin)*(double)(j)/(double)(bayesline->spline->n-1);
if (x <= fa[0]) z = Sna[0];
else if(x >= fa[jj-1]) z = Sna[jj-1];
else
{
i=k-10;
mdn = 0.0;
do
{
if(i>=0) mdn = fa[i];
i++;
} while(x > mdn);
k = i;
z = Sna[i];
}
spline->points[j] = x;
spline->data[j] = z;
}
for(j=0; j<bayesline->spline->n; j++)
{
......@@ -1622,7 +1641,7 @@ void copy_bayesline_params(struct BayesLineParams *origin, struct BayesLineParam
{
copy->data->tmax = origin->data->tmax;
int i;
int i;
for(i=0; i<origin->data->ncut; i++)
{
copy->spow[i] = origin->spow[i];
......@@ -1638,12 +1657,12 @@ void copy_bayesline_params(struct BayesLineParams *origin, struct BayesLineParam
copy->sfreq[i] = origin->sfreq[i];
copy->Sbase[i] = origin->Sbase[i];
copy->Sline[i] = origin->Sline[i];
}
}
copy_splineParams(origin->spline, copy->spline);
copy_splineParams(origin->spline_x,copy->spline_x);
copy_lorentzianParams(origin->lines_x, copy->lines_x);
}
void print_line_model(FILE *fptr, struct BayesLineParams *bayesline)
......@@ -2408,12 +2427,12 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
// Tukey window parameter. Flat for (1-alpha) of data
//t_rise = 0.4; // Standard LAL setting
/*
t_rise increased because 0.4 was found to be insufficient.
Perhaps there is an additional data conditioning step in
LALInferenceReadData?
*/
t_rise = 1.0;
/*
t_rise increased because 0.4 was found to be insufficient.
Perhaps there is an additional data conditioning step in
LALInferenceReadData?
*/
t_rise = 1.0;
alpha = (2.0*t_rise/Tobs);
df = 1.0/Tobs; // frequency resolution
......@@ -2492,22 +2511,22 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
// re-compute the power spectrum using the cleaned data
tukey(D, alpha, N);
gsl_fft_real_radix2_transform(D, 1, N);
gsl_fft_real_radix2_transform(D, 1, N);
// Form spectral model for whitening data (lines plus a smooth component)
spectrum(D, Sn, specD, sspecD, df, N);
fac = Tobs*Tobs/((double)(N)*(double)(N))/2.;
//copy output of clean (D) into residual for output
residual[0] = 0.0;
residual[1] = 0.0;
for(i=1; i< N/2; i++)
{
residual[2*i] = D[i] * sqrt(fac);
residual[2*i+1] = D[N-i] * sqrt(fac);
}
//copy output of clean (D) into residual for output
residual[0] = 0.0;
residual[1] = 0.0;
for(i=1; i< N/2; i++)
{
residual[2*i] = D[i] * sqrt(fac);
residual[2*i+1] = D[N-i] * sqrt(fac);
}
if(modelprint == 1)
{
......@@ -2580,7 +2599,7 @@ void blstart(double *data, double *residual, int N, double dt, double fmin, int
lineh[j] = (max-1.0)*sspecD[ii]*fac;
Abar = 0.5*(specD[ii+1]/sspecD[ii+1]+specD[ii-1]/sspecD[ii-1]);
lineQ[j] = sqrt((max/Abar-1.0))*linef[j]*Tobs;
j++;
j++;
//printf("%d %e %e %e\n", j, linef[j], linew[j], lineQ[j]);
}
}
......@@ -2713,10 +2732,10 @@ static void Lorentzian(double *Slines, double Tobs, lorentzianParams *lines, int
}
void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double *freqData)
{
//Initialize BayesLine data structures
BayesLineInitialize(bayesline);
{
//Initialize BayesLine data structures
BayesLineInitialize(bayesline);
/*********************************************************************************************
The pspline array holds the locations of the control points (knots) used in the spline
......@@ -2731,7 +2750,7 @@ void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double
lorentzianParams *lines = bayesline->lines_x;
splineParams *spline = bayesline->spline_x;
int i,j;
int i,j;
//number of samples in the data
double Tobs = data->Tobs;
......@@ -2759,34 +2778,34 @@ void BayesLineBurnin(struct BayesLineParams *bayesline, double *timeData, double
bayesline->Sline[i] = 0.0;
}
//interpolate cubic spline points
CubicSplineGSL(spline->n,spline->points,spline->data,N/2-imin,f+imin,logSn+imin);
//initialize work space for spline model
bayesline->data->flow = bayesline->data->fmin;
for(i=0; i< bayesline->data->ncut; i++)
{
j = i + imin - bayesline->data->nmin;
bayesline->power[j] = freqData[2*j]*freqData[2*j]+freqData[2*j+1]*freqData[2*j+1];
bayesline->spow[i] = bayesline->power[j];
bayesline->sfreq[i] = bayesline->freq[j];
}
//form-up line model
for(i=0; i<lines->n; i++)
{
lines->larray[i] = i;
Lorentzian(Sl, Tobs, lines, lines->larray[i], N);
}
//combine spline and line model
for(i=0; i<N/2-imin; i++)
{
bayesline->Sbase[i] = exp(logSn[i+imin]);
bayesline->Sline[i] = Sl[i+imin];
bayesline->Snf[i] = bayesline->Sbase[i] + bayesline->Sline[i];
}
//interpolate cubic spline points
CubicSplineGSL(spline->n,spline->points,spline->data,N/2-imin,f+imin,logSn+imin);
//initialize work space for spline model
bayesline->data->flow = bayesline->data->fmin;
for(i=0; i< bayesline->data->ncut; i++)
{
j = i + imin - bayesline->data->nmin;
bayesline->power[j] = freqData[2*j]*freqData[2*j]+freqData[2*j+1]*freqData[2*j+1];
bayesline->spow[i] = bayesline->power[j];
bayesline->sfreq[i] = bayesline->freq[j];
}
//form-up line model
for(i=0; i<lines->n; i++)
{
lines->larray[i] = i;
Lorentzian(Sl, Tobs, lines, lines->larray[i], N);
}
//combine spline and line model
for(i=0; i<N/2-imin; i++)
{
bayesline->Sbase[i] = exp(logSn[i+imin]);
bayesline->Sline[i] = Sl[i+imin];
bayesline->Snf[i] = bayesline->Sbase[i] + bayesline->Sline[i];
}
//Use initial estimate of PSD to set priors
for(i=0; i<N/2-imin; i++)
......
/*
* Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
......@@ -79,7 +98,7 @@ struct BayesLineParams
dataParams *data;
splineParams *spline;
splineParams *spline_x;
lorentzianParams *lines_x;
lorentzianParams *lines_x;
BayesLinePriors *priors;
double *Snf;
......
/*
Compile: See Makefile in bayeswave directory
* Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg, Margaret Millhouse,Katerina Chatziioannou
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
/*************************** REQUIRED LIBRARIES ***************************/
#include <stdio.h>
......@@ -374,9 +390,9 @@ int main(int argc, char *argv[])
timeData[ifo][2*i] = data->s[ifo][2*(imin+1)];
timeData[ifo][2*i+1] = data->s[ifo][2*(imin+1)+1];
}
fftw_wrapper(timeData[ifo],data->N,-1);
//normalize time series as expected by BayesLine
fftw_wrapper(timeData[ifo],data->N,-1);
//normalize time series as expected by BayesLine
for(i=0; i<N; i++) timeData[ifo][i]/=data->Tobs;
fprintf(stdout,"BayesLine search phase for IFO %s\n", data->ifos[ifo]);
......
/*
* Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
......
/* Compile:
gcc -O2 deglitcher.c -lm -o deglitcher -lgsl
/*
* Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg, James A. Clark, Jonah B. Kanner
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA