dataset.c

Go to the documentation of this file.

/*
*  Copyright (C) 2007 Vladimir Dergachev
*
*  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
*/

#define _GNU_SOURCE
#define _FILE_OFFSET_BITS 64
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <errno.h>
#include <alloca.h>
#include <math.h>
#include <sys/time.h>
#include <time.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_integration.h>

#include <sys/types.h>
#include <sys/stat.h>
#include <dirent.h>
#include <fcntl.h>

#include "dataset.h"
#include "statistics.h"
#include "cmdline.h"
#include "intervals.h"
#include "hookup.h"
#include "grid.h"
#include "rastermagic.h"
#include "hookup.h"
#include "util.h"
#include "jobs.h"

#include <lal/GeneratePulsarSignal.h>

// typedef double REAL8;
// typedef int INT4;
// typedef float REAL4;

extern FILE *LOG;
extern struct gengetopt_args_info args_info;
extern int first_bin;
extern int nbins;

extern int ntotal_polarizations;
extern SKY_GRID *fine_grid;
extern SKY_GRID *patch_grid;

extern INT64 spindown_start;

extern int fake_injection;
extern double spindown;

char s[20000];

DATASET * datasets=NULL;
int d_free=0;
int d_size=0;

extern int do_CutOff;

int lock_file=-1;

/* this is a busy wait.. as condor does not allow sleep() */
static void spin_wait(int seconds)
{
time_t start_time, end_time;
int i;
time(&start_time);
i=0;
while(1) {
        i++;
        if(i>1000000) {
                time(&end_time);
                if(end_time>start_time+seconds)return;
                i=0;
                }
        }
}

/* we need the structure to be global as a workaround for a bug in condor libraries:
   this way the pointer is 32bit even on 64bit machines */

static void acquire_lock(char *filename)
{
int i;
struct flock fl;

lock_file=open(filename, O_CREAT | O_RDWR, 0666);

memset(&fl, 0, sizeof(fl));
fl.l_type=F_WRLCK;
fl.l_whence=SEEK_SET;
fl.l_start=0;
fl.l_len=1;
errno=0;
i=0;
while((lock_file<0) || (direct_fcntl(lock_file, F_SETLK, &fl)<0)){
        if(lock_file<0) {
                fprintf(stderr, "Could not open file \"%s\" for locking\n", filename);
                }
        if(i>10) {
                fprintf(stderr, "Waiting for lock: %s\n", strerror(errno));
                i=0;
                }
        if(lock_file>=0)close(lock_file);
        condor_safe_sleep(args_info.lock_retry_delay_arg);
        lock_file=open(filename, O_CREAT | O_RDWR, 0666);
        i++;
        }
}

static void release_lock(void)
{
if(lock_file>=0)close(lock_file);
lock_file=-1;
}

typedef struct {
        float e[26];
        int bin;

        double a_plus;
        double a_cross;
        double cos_e;
        double sin_e;

        double segment_start;

        double ref_time;

        double freq;
        double spindown;

        double ra;
        double dec;

        double coherence_time;

        double extra_phase; /* phase accumulated from the start of the run */
        double phase_integration_factor;

        } SIGNAL_PARAMS;

static void compute_signal(double *re, double *im, double *f, double t, SIGNAL_PARAMS *p)
{
float det_vel[3]={NAN, NAN, NAN};
float f_plus, f_cross;
double doppler, omega_t, c_omega_t, s_omega_t;
double hann;

get_detector_vel(round(t), det_vel);

doppler=p->e[0]*det_vel[0]+p->e[1]*det_vel[1]+p->e[2]*det_vel[2];

/*get_AM_response(round(t), p->dec, p->ra, 0.7853982, &f_plus, &f_cross);
fprintf(stderr, "%d f_plus=%f f_cross=%f\n", (int)round(t), f_plus, f_cross);*/
get_AM_response(round(t), p->dec, p->ra, 0.0, &f_plus, &f_cross);
// fprintf(stderr, "%d f_plus=%f f_cross=%f\n", (int)round(t), f_plus, f_cross);

*f=p->freq*(1.0+doppler)+p->spindown*(t-p->ref_time);

//*f=p->freq*(1.0+doppler);

omega_t=2.0*M_PI*(*f-p->bin/p->coherence_time)*(t-p->segment_start)+p->extra_phase;
hann=0.5*(1.0-cos(2.0*M_PI*(t-p->segment_start)/p->coherence_time));

c_omega_t=cos(omega_t);
s_omega_t=sin(omega_t);

/* 
        cos(a)*cos(-b)=0.5*(cos(a-b)+cos(a+b))
        sin(a)*cos(-b)=0.5*(sin(a-b)+sin(a+b))
        cos(a)*sin(-b)=0.5*(sin(a-b)-sin(a+b))
        sin(a)*sin(-b)=0.5*(cos(a-b)-cos(a+b))

        |f-w|<<1
        cos(ft)*exp(-iwt)~= 0.5*(cos((f-w)t)+i*sin((f-w)t))
        sin(ft)*exp(-iwt)~= 0.5*(sin((f-w)t)-i*cos((f-w)t))
*/

// fprintf(stderr, "a_plus=%f a_cross=%f cos_e=%f sin_e=%f\n", p->a_plus, p->a_cross, p->cos_e, p->sin_e);

/* Note: we do not have an extra factor of 0.5 because of normalization convention used, 
  there is an extra factor of 2.0 to compensate for Hann windowing */
*re=hann*(f_plus*(p->a_plus*c_omega_t*p->cos_e-p->a_cross*s_omega_t*p->sin_e)+
        f_cross*(p->a_plus*c_omega_t*p->sin_e+p->a_cross*s_omega_t*p->cos_e));

*im=hann*(f_plus*(p->a_plus*s_omega_t*p->cos_e+p->a_cross*c_omega_t*p->sin_e)+
        f_cross*(p->a_plus*s_omega_t*p->sin_e-p->a_cross*c_omega_t*p->cos_e));

//fprintf(stderr, "response=%g a_plus=%g a_cross=%g \n", response, p->a_plus, p->a_cross);

// fprintf(stderr, "f=%.4f (%.3f) re=%.2f im=%.2f omega_t=%.2f t=%.1f (%f) doppler=%g\n", (*f), (*f)*1800.0-p->bin*1.0, *re, *im, omega_t, t-p->ref_time, t, doppler);

}

static double signal_re(double t, void *params)
{
double re, im, f;
compute_signal(&re, &im, &f, t, params);
return(re);
}

static double signal_im(double t, void *params)
{
double re, im, f;
compute_signal(&re, &im, &f, t, params);
return(im);
}

static double signal_omega(double t, void *params)
{
SIGNAL_PARAMS *p=params;
float det_vel[3];
double doppler, omega;

get_detector_vel(round(t), det_vel);

doppler=p->e[0]*det_vel[0]+p->e[1]*det_vel[1]+p->e[2]*det_vel[2];
omega=2.0*M_PI*p->freq*doppler-p->phase_integration_factor;

return(omega);
}

static double compute_phase(gsl_integration_workspace *giw, int giw_size, double start_time, double end_time)
{
SIGNAL_PARAMS p;
gsl_function F; 
int err;
double result, abserr;
double T=end_time-start_time;

if(!(fabs(T)>0))return 0.0;

F.params=&p;
p.bin=0;

p.ra=args_info.fake_ra_arg;
p.dec=args_info.fake_dec_arg;
p.freq=args_info.fake_freq_arg;
p.spindown=args_info.fake_spindown_arg;
p.ref_time=args_info.fake_ref_time_arg;
p.segment_start=0;
p.coherence_time=0;

p.cos_e=cos(2.0*(args_info.fake_psi_arg-args_info.fake_phi_arg));
p.sin_e=sin(2.0*(args_info.fake_psi_arg-args_info.fake_phi_arg));

precompute_am_constants(p.e, args_info.fake_ra_arg, args_info.fake_dec_arg);

F.function=signal_omega;
p.phase_integration_factor=0.0;

while(1) {
        err=gsl_integration_qag(&F, start_time, end_time,
                1e-6, 1e-5,
                giw_size,
                GSL_INTEG_GAUSS61,
                giw,
                &result,
                &abserr
                );

        if(fabs(result)>10) {
                p.phase_integration_factor+=result/T;
                continue;
                }

        p.extra_phase=result+p.phase_integration_factor*T+
                2*M_PI*(p.freq*T+0.5*p.spindown*T*T+p.spindown*(start_time-p.ref_time)*T);

/*      if(err)fprintf(stderr, "phase %d result=%g (%g) pif=%g abserr=%g %s\n", segment, result, result/(d->gps[segment]-p.ref_time), p.phase_integration_factor, abserr,  gsl_strerror(err)); */
        break;
        }
return(p.extra_phase);
}

static void test_inject_fake_signal(void)
{
SIGNAL_PARAMS p;
double cos_i;
double re, im, f, re2, im2, f2;
int status=0;

p.bin=0;

p.ra=2.56;
p.dec=1.43;
p.freq=500.12;
p.spindown=1.23e-9;
p.ref_time=793154935;
p.segment_start=793161245;
p.coherence_time=1800.0;
p.extra_phase=0;

cos_i=cos(1.23);

p.a_plus=(1+cos_i*cos_i)*0.5;
p.a_cross=cos_i;

p.cos_e=cos(2.0*(1.23));
p.sin_e=sin(2.0*(1.23));
precompute_am_constants(p.e, p.ra, p.dec);

get_detector("LHO");

compute_signal(&re, &im, &f, 793161250.0, &p);
fprintf(stderr, "compute_signal_test1: %g %g %f\n", re, im, f);
fprintf(LOG, "compute_signal_test1: %g %g %f\n", re, im, f);

if(abs(re-9.59669e-06)>1e-11 ||
   abs(im-1.83957e-05)>1e-11 ||
   abs(f-500.101774)>1e-5) status|=1;

p.cos_e=1.0;
p.sin_e=0.0;

compute_signal(&re, &im, &f, 793161250.0, &p);
fprintf(stderr, "compute_signal_test2: %g %g %f\n", re, im, f);
fprintf(LOG, "compute_signal_test2: %g %g %f\n", re, im, f);

p.cos_e=1.0;
p.sin_e=0.0;
p.extra_phase=M_PI*0.5;

compute_signal(&re2, &im2, &f2, 793161250.0, &p);
fprintf(stderr, "compute_signal_test3: %g %g %f\n", re2, im2, f2);
fprintf(LOG, "compute_signal_test3: %g %g %f\n", re2, im2, f2);

if(abs(f2-f)>1e-5 ||
   abs(re-im2)>1e-11 ||
   abs(im+re2)>1e-11) status|=2;

if(status) {
        fprintf(stderr, "compute_signal_test: failed %d\n", status);
        fprintf(LOG, "compute_signal_test: failed %d\n", status);
        exit(-1);
        }
}

static void inject_fake_signal(DATASET *d, int segment, double accumulated_phase)
{
double result, abserr;
size_t neval;
int err;
int window=5, bin;
int i;
double re, im, f, cos_i;
SIGNAL_PARAMS p;
gsl_function F; 

F.params=&p;
p.bin=0;

p.ra=args_info.fake_ra_arg;
p.dec=args_info.fake_dec_arg;
p.freq=args_info.fake_freq_arg;
p.spindown=args_info.fake_spindown_arg;
p.ref_time=args_info.fake_ref_time_arg;
p.segment_start=d->gps[segment];
p.coherence_time=d->coherence_time;

//fprintf(stderr, "frequency=%f spindown=%g ra=%f dec=%f\n", p.freq, p.spindown, p.ra, p.dec);

cos_i=cos(args_info.fake_iota_arg);

p.a_plus=(1+cos_i*cos_i)*0.5;
p.a_cross=cos_i;

p.cos_e=cos(2.0*(args_info.fake_psi_arg-args_info.fake_phi_arg));
p.sin_e=sin(2.0*(args_info.fake_psi_arg-args_info.fake_phi_arg));

// fprintf(stderr, "a_plus=%f a_cross=%f cos_e=%f sin_e=%f\n", p.a_plus, p.a_cross, p.cos_e, p.sin_e);

precompute_am_constants(p.e, args_info.fake_ra_arg, args_info.fake_dec_arg);
compute_signal(&re, &im, &f, d->gps[segment]+(int)(d->coherence_time/2), &p);


bin=round(f*1800.0-d->first_bin);
if((bin+window>nbins) || (bin<window))fprintf(stderr, "Injected signal outside loaded band: bin=%d, segment=%d\n", bin, segment);
if(bin+window>nbins)bin=nbins-window;
if(bin<window)bin=window;

p.extra_phase=accumulated_phase;


for(i=bin-window; i<=bin+window; i++) {
        p.bin=i+d->first_bin;

        F.function=signal_re;
        err=gsl_integration_qng(&F, d->gps[segment], d->gps[segment]+d->coherence_time,
                1, 1e-3,
                &result, &abserr, &neval);
/*      fprintf(stderr, "re %d %d result=%g abserr=%g %d %s\n", segment, i, result, abserr, neval, gsl_strerror(err)); */
        d->re[segment*d->nbins+i]+=args_info.fake_strain_arg*result*16384.0/args_info.strain_norm_factor_arg;


        F.function=signal_im;
        err=gsl_integration_qng(&F, d->gps[segment], d->gps[segment]+d->coherence_time,
                1, 1e-3,
                &result, &abserr, &neval);
/*      fprintf(stderr, "im %d %d result=%g abserr=%g %d %s\n", segment, i, result, abserr, neval, gsl_strerror(err)); */
        d->im[segment*d->nbins+i]+=args_info.fake_strain_arg*result*16384.0/args_info.strain_norm_factor_arg;
        }

}

static void inject_fake_signal2(DATASET *d, int segment)
{
SFTandSignalParams sft_params;
PulsarSignalParams pulsar_params;

sft_params.resTrig=0;
sft_params.Dterms=5;
sft_params.nSamples=5;

}

static void init_dataset(DATASET *d)
{
int i;
memset(d, 0, sizeof(*d));
d->name="unknown";
for(i=0;i<MAX_LOCKS;i++)
        d->lock_file[i]=NULL;
d->validated=0;
d->detector="unknown";
d->gps_start=0;
d->gps_stop=-1;

d->coherence_time=1800.0;
d->nbins=nbins;
d->first_bin=first_bin;

d->segment_list=NULL;
d->veto_segment_list=NULL;
d->no_duplicate_gps=1;

d->size=1000;
d->free=0;
d->gps=do_alloc(d->size, sizeof(*d->gps));
d->dc_factor=1.0;
d->dc_factor_touched=0;
d->dc_factor_blocked=0;
d->re=do_alloc(d->size*d->nbins, sizeof(*d->re));
d->im=do_alloc(d->size*d->nbins, sizeof(*d->im));
d->sft_veto=NULL;
d->veto_level=1e-2; /* this takes care of SFTs that have 1/100 weight... */
d->veto_spike_level=1.7; /* this takes care of SFTs with spikes 50 times median level */

d->weight=1.0;

d->polarizations=NULL;
d->polarizations=do_alloc(ntotal_polarizations, sizeof(*(d->polarizations)));
memset(d->polarizations, 0, (ntotal_polarizations)*sizeof(*(d->polarizations)));
}

static void compute_detector_speed(DATASET *d)
{
int i;
float det_vel_ra, det_vel_dec;
double orbital_axis[3];
float band_axis_ra, band_axis_dec;
fprintf(stderr,"Computing detector speed for dataset %s\n", d->name);
d->detector_velocity=do_alloc(3*d->free, sizeof(*d->detector_velocity));
d->average_detector_velocity[0]=0.0;
d->average_detector_velocity[1]=0.0;
d->average_detector_velocity[2]=0.0;
for(i=0;i<d->free;i++){
        /* middle of the 30min interval */
        get_detector_vel(d->gps[i]+(int)(d->coherence_time*0.5),&(d->detector_velocity[3*i]));
                
        d->average_detector_velocity[0]+=d->detector_velocity[3*i];
        d->average_detector_velocity[1]+=d->detector_velocity[3*i+1];
        d->average_detector_velocity[2]+=d->detector_velocity[3*i+2];
        }
d->average_detector_velocity[0]/=d->free;
d->average_detector_velocity[1]/=d->free;
d->average_detector_velocity[2]/=d->free;
fprintf(LOG,"dataset %s average detector velocity: %g %g %g\n", d->name, 
        d->average_detector_velocity[0],
        d->average_detector_velocity[1],
        d->average_detector_velocity[2]);
det_vel_dec=atan2f(d->average_detector_velocity[2], 
        sqrt(d->average_detector_velocity[0]*d->average_detector_velocity[0]+
                d->average_detector_velocity[1]*d->average_detector_velocity[1]));
det_vel_ra=atan2f(d->average_detector_velocity[1], d->average_detector_velocity[0]);
if(det_vel_ra<0)det_vel_ra+=2.0*M_PI;
fprintf(LOG,"dataset %s average detector velocity RA (degrees) : %f\n", d->name, det_vel_ra*180.0/M_PI);
fprintf(LOG,"dataset %s average detector velocity DEC (degrees): %f\n", d->name, det_vel_dec*180.0/M_PI);
fprintf(stderr,"dataset %s average detector velocity RA (degrees) : %f\n", d->name, det_vel_ra*180.0/M_PI);
fprintf(stderr,"dataset %s average detector velocity DEC (degrees): %f\n", d->name, det_vel_dec*180.0/M_PI);

orbital_axis[0]=0.0;
orbital_axis[1]=-sin(M_PI*23.44/180.0);
orbital_axis[2]=cos(M_PI*23.44/180.0);

/* crossproduct gives the vector perpedicular to both the average doppler shift and
  orbital axis */
d->band_axis[0]=orbital_axis[1]*d->average_detector_velocity[2]-orbital_axis[2]*d->average_detector_velocity[1];
d->band_axis[1]=orbital_axis[2]*d->average_detector_velocity[0]-orbital_axis[0]*d->average_detector_velocity[2];
d->band_axis[2]=orbital_axis[0]*d->average_detector_velocity[1]-orbital_axis[1]*d->average_detector_velocity[0];

/* Normalize */
d->band_axis_norm=sqrt(d->band_axis[0]*d->band_axis[0]+
        d->band_axis[1]*d->band_axis[1]+
        d->band_axis[2]*d->band_axis[2]);
/* replace 0.0 with something more reasonable later */
if(d->band_axis_norm<=0.0){
        d->band_axis[0]=0.0;
        d->band_axis[1]=0.0;
        d->band_axis[2]=1.0;
        } else {
        d->band_axis[0]/=d->band_axis_norm;
        d->band_axis[1]/=d->band_axis_norm;
        d->band_axis[2]/=d->band_axis_norm;
        }
/* 
  Normalize band_axis_norm so it matches the definition of \vec{u}
*/
d->band_axis_norm*=2.0*M_PI/(365.0*24.0*3600.0);

fprintf(LOG, "dataset %s auto band axis norm: %g\n", d->name, d->band_axis_norm);
fprintf(LOG, "dataset %s maximum S contribution from Doppler shifts: %g\n", d->name, d->band_axis_norm*(d->first_bin+d->nbins*0.5)/d->coherence_time);

if(args_info.band_axis_norm_given) {
        d->band_axis_norm=args_info.band_axis_norm_arg;
        }

fprintf(LOG, "dataset %s actual band axis norm: %g\n", d->name, d->band_axis_norm);

d->large_S=6.0/(d->coherence_time*(d->gps[d->free-1]-d->gps[0]+d->coherence_time));

fprintf(LOG, "dataset %s auto large S: %g\n", d->name, d->large_S);

if(args_info.large_S_given){
        d->large_S=args_info.large_S_arg;
        }
fprintf(LOG, "dataset %s large S: %g\n", d->name, d->large_S);
        

fprintf(LOG,"dataset %s auto band axis: %g %g %g\n", d->name, 
        d->band_axis[0],
        d->band_axis[1],
        d->band_axis[2]);
fprintf(stderr,"dataset %s auto band axis: %g %g %g\n", d->name, 
        d->band_axis[0],
        d->band_axis[1],
        d->band_axis[2]);

fprintf(LOG, "dataset %s band_axis: %s\n", d->name, args_info.band_axis_arg);
fprintf(stderr, "dataset %s band_axis: %s\n", d->name, args_info.band_axis_arg);
if(!strcasecmp(args_info.band_axis_arg, "equatorial")){
        d->band_axis[0]=0.0;
        d->band_axis[1]=0.0;
        d->band_axis[2]=1.0;    
        } else
if(!strncasecmp(args_info.band_axis_arg, "explicit", 8)){
        int q;
        q=sscanf(args_info.band_axis_arg+8, "(%lf,%lf,%lf)", 
                &(d->band_axis[0]),
                &(d->band_axis[1]),
                &(d->band_axis[2]));
        if(q!=3){
                fprintf(stderr,"Warning ! In dataset %s not all explicit band axis values were assigned. Format error ?\n", d->name);
                fprintf(LOG,"Warning ! In dataset %s not all explicit band axis values were assigned. Format error ?\n", d->name);
                }
        }
        
fprintf(LOG,"dataset %s actual band axis: %g %g %g\n", d->name, 
        d->band_axis[0],
        d->band_axis[1],
        d->band_axis[2]);
fprintf(stderr,"dataset %s actual band axis: %g %g %g\n", d->name, 
        d->band_axis[0],
        d->band_axis[1],
        d->band_axis[2]);

band_axis_dec=atan2f(d->band_axis[2], 
        sqrt(d->band_axis[0]*d->band_axis[0]+d->band_axis[1]*d->band_axis[1]));
band_axis_ra=atan2f(d->band_axis[1], d->band_axis[0]);

if(band_axis_ra<0)band_axis_ra+=2.0*M_PI;
fprintf(stderr,"dataset %s band axis RA (degrees) : %f\n", d->name, band_axis_ra*180.0/M_PI);
fprintf(stderr,"dataset %s band axis DEC (degrees): %f\n", d->name, band_axis_dec*180.0/M_PI);
fprintf(LOG,"dataset %s band axis RA (degrees) : %f\n", d->name, band_axis_ra*180.0/M_PI);
fprintf(LOG,"dataset %s band axis DEC (degrees): %f\n", d->name, band_axis_dec*180.0/M_PI);
}

/* Note: this messes with tm variable,
  tm[i]=exp(log(10)*TMedians[i])/Modulation[i] */
float FindCutOff(float *tm, int nsegments)
{
int i;
double sum,sum_squared,mean,sigma;
double best_cutoff,smallest_sigma;
double a;
int best_i;
sort_floats(tm, nsegments);
sum=0;
sum_squared=0;
best_i=0;
smallest_sigma=10;
best_cutoff=tm[0];
for(i=0;i<nsegments;i++) {
        a=tm[i]/tm[0];
        sum+=a;
        sum_squared+=a*a;
        mean=sum/(i+1);
        sigma=10*sqrt((sum_squared))/(i+1);
        
        if(sigma<smallest_sigma) {
                smallest_sigma=sigma;
                best_i=i;
                best_cutoff=tm[i];
                }
        }
//fprintf(stderr,"Cutoff: i=%d sigma=%g cutoff=%g (%g,%g,..)\n",best_i, smallest_sigma, best_cutoff,tm[0],tm[1]);
return best_cutoff;
}

void apply_hanning_filter(DATASET *d)
{
int i,k;
float *tmp, *p;
tmp=do_alloc(d->nbins, sizeof(*tmp));
for(i=0;i<d->free;i++) {
        p=&(d->re[i*d->nbins]);

        /* wrap around, just so that we don't have 0 exactly */
        tmp[0]=p[0]-(p[d->nbins-1]+p[1])*0.5;
        tmp[d->nbins-1]=p[d->nbins-1]-(p[d->nbins-2]+p[0])*0.5;

        for(k=1;k<(d->nbins-1);k++) {
                tmp[k]=p[k]-(p[k-1]+p[k+1])*0.5;
                }
        memcpy(p, tmp, nbins*sizeof(*p));       
        }
for(i=0;i<d->free;i++) {
        p=&(d->im[i*d->nbins]);

        /* wrap around, just so that we don't have 0 exactly */
        tmp[0]=p[0]-(p[d->nbins-1]+p[1])*0.5;
        tmp[d->nbins-1]=p[d->nbins-1]-(p[d->nbins-2]+p[0])*0.5;

        for(k=1;k<(d->nbins-1);k++) {
                tmp[k]=p[k]-(p[k-1]+p[k+1])*0.5;
                }
        memcpy(p, tmp, nbins*sizeof(*p));       
        }
free(tmp);
}

typedef struct {
        int index;
        INT64 gps;
        } ELEMENT;

int element_cmp(ELEMENT *e1, ELEMENT *e2)
{
if(e1->gps<e2->gps)return(-1);
if(e1->gps>e2->gps)return(1);
return 0;
}

void sort_dataset(DATASET *d)
{
ELEMENT *p;
float *bin_re, *bin_im;
INT64 *gps;
int i, k;

p=do_alloc(d->free, sizeof(*p));
for(i=0;i<d->free;i++) {
        p[i].index=i;
        p[i].gps=d->gps[i];
        }

qsort(p, d->free, sizeof(*p), element_cmp);

bin_re=do_alloc(d->free*d->nbins, sizeof(*bin_re));
for(i=0;i<d->free;i++) {
        k=p[i].index;
        memcpy(&(bin_re[i*d->nbins]), &(d->re[k*d->nbins]), d->nbins*sizeof(*bin_re));
        }
free(d->re);
d->re=bin_re;

bin_im=do_alloc(d->free*d->nbins, sizeof(*bin_im));
for(i=0;i<d->free;i++) {
        k=p[i].index;
        memcpy(&(bin_im[i*d->nbins]), &(d->im[k*d->nbins]), d->nbins*sizeof(*bin_im));
        }
free(d->im);
d->im=bin_im;

gps=do_alloc(d->free, sizeof(*gps));
for(i=0;i<d->free;i++) {
        k=p[i].index;
        gps[i]=d->gps[k];
        }
free(d->gps);
d->gps=gps;

d->size=d->free;

free(p);
}

void veto_sfts(DATASET *d)
{
int i, j;
float power, *x, *y;
d->sft_veto=do_alloc(d->free, sizeof(*d->sft_veto));
for(i=0;i<d->free;i++) {
        d->sft_veto[i]=0;

        if(d->expTMedians[i]<d->veto_level) {
                d->sft_veto[i]=1;
                continue;
                }

        x=&(d->re[i*d->nbins]);
        y=&(d->im[i*d->nbins]);
        for(j=0;j<d->nbins;j++) {
                power=(*x)*(*x)+(*y)*(*y);
                if(log10(power)>d->TMedians[i]+d->veto_spike_level) {
                        d->sft_veto[i]=2;
                        break;
                        }
                x++;
                y++;
                }
        }
}

typedef struct {
        int point;
        DATASET *d;
        
        } CUTOFF_DATA;

void cutoff_cruncher(int thread_id, CUTOFF_DATA *data)
{
int j, m;
float *tm;
tm=do_alloc(data->d->free, sizeof(*tm));
for(m=0;m<ntotal_polarizations;m++) {
/*      if(patch_grid->band[i]<0){
                d->polarizations[m].patch_CutOff[i]=0.0;
                continue;
                }*/
        if(!do_CutOff) {
                data->d->polarizations[m].patch_CutOff[data->point]=1e30;
                continue;
                }
        for(j=0;j<data->d->free;j++)tm[j]=1.0/(sqrt(data->d->expTMedians[j])*AM_response(j, patch_grid, data->point, data->d->polarizations[m].AM_coeffs));
        data->d->polarizations[m].patch_CutOff[data->point]=FindCutOff(tm, data->d->free);
        }
free(tm);
}

static int validate_dataset(DATASET *d)
{
int i,j, m, k;
float *tm;
CUTOFF_DATA *cd;

if(d->validated)return 1;
fprintf(stderr, "Validating dataset \"%s\"\n", d->name);
fprintf(stderr, "Found %d segments\n", d->free);

fprintf(LOG, "dataset %s nsegments : %d\n", d->name, d->free);
fprintf(LOG, "dataset %s detector : %s\n", d->name, d->detector);

if(d->free<1){
        /* free large chunks of ram that might have been reserved */
        if(d->size>0) {
                free(d->gps);
                free(d->re);
                free(d->im);
                d->gps=NULL;
                d->re=NULL;
                d->im=NULL;
                }
        return 0;
        }

if(!strcmp(d->detector, "unknown")) {
        fprintf(stderr, "Each dataset must specify detector being used\n");
        return 0;
        }

get_detector(d->detector);

sort_dataset(d);

if(fake_injection) {
        gsl_integration_workspace *giw;
        double accumulated_phase, current_phase, start_time;
        int workspace_size=10000;

        fprintf(stderr, "Injecting fake signal.\n");
        fprintf(LOG, "Injecting fake signal.\n");
        giw=gsl_integration_workspace_alloc(workspace_size);

        start_time=args_info.fake_ref_time_arg;
        accumulated_phase=0;
        for(i=0;i<d->free;i++) {
                current_phase=compute_phase(giw, workspace_size, start_time, d->gps[i]);
                inject_fake_signal(d, i, accumulated_phase+current_phase);
                if(fabs(d->gps[i]-start_time)>PHASE_ACCUMULATION_TIMEOUT) {
                        start_time=d->gps[i];
                        accumulated_phase+=current_phase;
                        }
                }
        gsl_integration_workspace_free(giw);
        }

/* compute_power(d); */

/* compute time from the start of the run */
d->hours=do_alloc(d->free, sizeof(*(d->hours)));
d->hours_d=do_alloc(d->free, sizeof(*(d->hours_d)));
for(i=0;i<d->free;i++){
        d->hours[i]=(1.0*(d->gps[i]-d->gps[0]))/3600.0;
        d->hours_d[i]=(1.0*(d->gps[i]-d->gps[0]))/3600.0;
        }
        
/* compute frequency array */
d->frequencies=do_alloc(d->nbins, sizeof(*(d->frequencies)));
d->freq_d=do_alloc(d->nbins, sizeof(*(d->freq_d)));
for(i=0;i<d->nbins;i++){
        d->freq_d[i]=(1.0*(d->first_bin+i))/d->coherence_time;
        d->frequencies[i]=(1.0*(d->first_bin+i))/d->coherence_time;
        }

compute_detector_speed(d);

d->TMedians=do_alloc(d->free, sizeof(*(d->TMedians)));
d->expTMedians=do_alloc(d->free, sizeof(*(d