Select Git revision
FT_example2.m
libFStatNomad.cpp 116.75 KiB
/*
* Copyright (C) 2012 Miroslav Shaltev.
*
* May include code snippets by:
* Holger Pletsch, Karl Wette, Reinhard Prix,
* Badri Krishnan, Alicia Sintes, Bernd Machenschalk.
*
* 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 __STDC_CONSTANT_MACROS
#include "libFStatNomad.h"
static FstatCandidate empty_FstatCandidate;
static ConfigVariables empty_ConfigVariables;
using namespace std;
REAL8 FStatNomad::DoPredictFStat(INT4 fstatsearchtype, LALSegList *segList) {
AntennaPatternMatrix Mmunu;
init_SP(fstatsearchtype,segList);
SkyPosition skypos;
MultiAMCoeffs *multiAMcoef = NULL;
UINT4 X, i, numSFTs;
LALStatus lstatus = blank_status;
/* normalize skyposition: correctly map into [0,2pi]x[-pi/2,pi/2] */
skypos.longitude = SignalPoint.alpha;
skypos.latitude = SignalPoint.delta;
skypos.system = COORDINATESYSTEM_EQUATORIAL;
LALNormalizeSkyPosition ( &lstatus, &skypos, &skypos);
REAL8 rho2 = 0;
for (i = 0; i < Fstat_in_vec->length; i++ ) {
MultiSFTVector* multiSFTs = XLALLoadMultiSFTs( catalogSeq.data + i, freqmin, freqmax );
MultiNoiseWeights *multiNoiseWeights = NULL;
MultiPSDVector *psd = NULL;
LALNormalizeMultiSFTVect ( &lstatus, &psd, multiSFTs, blocksRngMed() );
LALComputeMultiNoiseWeights ( &lstatus, &multiNoiseWeights, psd, blocksRngMed(), 0 );
LALDestroyMultiPSDVector ( &lstatus, &psd );
MultiDetectorStateSeries *multiDetStates = NULL;
LALGetMultiDetectorStates( &lstatus, &multiDetStates, multiSFTs, edat);
LALGetMultiAMCoeffs ( &lstatus, &multiAMcoef, multiDetStates, skypos );
/* noise-weighting of Antenna-patterns and compute A,B,C */
if ( XLALWeightMultiAMCoeffs ( multiAMcoef, multiNoiseWeights ) != XLAL_SUCCESS ) {
LogPrintf (LOG_NORMAL, "XLALWeightMultiAMCoeffs() failed with error = %d\n\n", xlalErrno );
ABORT ( &lstatus, FSTATFCNOMAD_EXLAL, FSTATFCNOMAD_MSGEXLAL );
}
/* OK: we only need the antenna-pattern matrix M_mu_nu */
Mmunu = multiAMcoef->Mmunu;
REAL8 aPlus = 0.5 * SignalPoint.h0 * ( 1.0 + SQ( SignalPoint.cosi) );
REAL8 aCross = SignalPoint.h0 * SignalPoint.cosi;
REAL8 al1, al2, al3;
REAL8 Ap2 = SQ(aPlus);
REAL8 Ac2 = SQ(aCross);
REAL8 cos2psi2 = SQ( cos(2*SignalPoint.psi) );
REAL8 sin2psi2 = SQ( sin(2*SignalPoint.psi) );
al1 = Ap2 * cos2psi2 + Ac2 * sin2psi2; /* A1^2 + A3^2 */
al2 = Ap2 * sin2psi2 + Ac2 * cos2psi2; /* A2^2 + A4^2 */
al3 = ( Ap2 - Ac2 ) * sin(2.0*SignalPoint.psi) * cos(2.0*SignalPoint.psi); /* A1 A2 + A3 A4 */
// /* SNR^2 */
rho2 += Mmunu.Sinv_Tsft * (Mmunu.Ad * al1 + Mmunu.Bd * al2 + 2.0 * Mmunu.Cd * al3 );
}
if (fstatsearchtype == SCSEARCH) {
INT4 nseg = Fstat_in_vec->length;
rho2 = rho2 / (REAL8)nseg;
printf("predicted average (N=%d) rho2: %.16f\n",nseg,rho2);
printf("predicted average (N=%d) 2F: %.16f\n",nseg, 4.0 + rho2);
} else if (fstatsearchtype == FCSEARCH) {
printf("predicted fully coherent rho2: %.16f\n",rho2);
printf("predicted fully coherent 2F: %.16f\n", 4.0 + rho2);
}
// // if (uvar.printFstat)
return rho2;
}
INT4 FStatNomad::write_trials_info(INT4 fstatsearchtype,INT4 cnt)
{
FILE *fpOut;
if ((fpOut = fopen(TrialsFile(fstatsearchtype,cnt), "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", TrialsFile(fstatsearchtype,cnt));
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
std::vector<UINT8> ltrials = trials();
for (INT4 i = 0; i < ltrials.size(); i++) {
fprintf(fpOut,"%d\t%d\n",i,ltrials.at(i));
}
fclose(fpOut);
}
INT4 FStatNomad::write_nomad_output(vector<ExPoint> epv,INT4 fstatsearchtype,INT4 stage) {
UINT4 numpoints;
FILE *fpOut;
INT4 thispoint;
UINT4 startpoint = 0;
numpoints = totalnumsteps;
if ((fpOut = fopen(ToplistFile(fstatsearchtype,stage), "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", ToplistFile(fstatsearchtype,stage));
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
/* get the log string */
fprintf( fpOut, "%%%% User Input:\n");
fprintf( fpOut, "%%%%-------------------------------------------\n");
fprintf( fpOut, "%%%% cmdline: %s\n", LogLine() );
/* add code version ID (only useful for git-derived versions) */
fprintf ( fpOut, "%%%% version: %s\n", VCSInfoString );
tCell best_cell = Cell();
INT4 n = 0;
fprintf( fpOut,"%%%% lower_bound:");
if (FixedDirection.at(PID_ALPHA)) {
fprintf( fpOut, " %.16e",FixedPoint.alpha);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_DELTA)) {
fprintf( fpOut, " %.16e",FixedPoint.delta);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_FREQ)) {
fprintf( fpOut, " %.16e",FixedPoint.freq);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_F1DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f1dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_F2DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f2dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_F3DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f3dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (BinarySearch()) {
if (FixedDirection.at(PID_ASINI)) {
fprintf( fpOut, " %.16e",FixedPoint.asini);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_ARGP)) {
fprintf( fpOut, " %.16e",FixedPoint.argp);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_ECC)) {
fprintf( fpOut, " %.16e",FixedPoint.ecc);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_TPSSB)) {
fprintf( fpOut, " %.16e",XLALGPSGetREAL8(&FixedPoint.TpSSB));
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
if (FixedDirection.at(PID_PERIOD)) {
fprintf( fpOut, " %.16e",FixedPoint.period);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).min);
n++;
}
}
fprintf( fpOut,"\n");
n = 0;
fprintf( fpOut,"%%%% start_point:");
if (FixedDirection.at(PID_ALPHA)) {
fprintf( fpOut, " %.16e",FixedPoint.alpha);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_DELTA)) {
fprintf( fpOut, " %.16e",FixedPoint.delta);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_FREQ)) {
fprintf( fpOut, " %.16e",FixedPoint.freq);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_F1DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f1dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_F2DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f2dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_F3DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f3dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (BinarySearch()) {
if (FixedDirection.at(PID_ASINI)) {
fprintf( fpOut, " %.16e",FixedPoint.asini);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_ARGP)) {
fprintf( fpOut, " %.16e",FixedPoint.argp);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_ECC)) {
fprintf( fpOut, " %.16e",FixedPoint.ecc);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_TPSSB)) {
fprintf( fpOut, " %.16e",XLALGPSGetREAL8(&FixedPoint.TpSSB));
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
if (FixedDirection.at(PID_PERIOD)) {
fprintf( fpOut, " %.16e",FixedPoint.period);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).startpoint);
n++;
}
}
fprintf( fpOut,"\n");
n = 0;
fprintf( fpOut,"%%%% upper_bound:");
if (FixedDirection.at(PID_ALPHA)) {
fprintf( fpOut, " %.16e",FixedPoint.alpha);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_DELTA)) {
fprintf( fpOut, " %.16e",FixedPoint.delta);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_FREQ)) {
fprintf( fpOut, " %.16e",FixedPoint.freq);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_F1DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f1dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_F2DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f2dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_F3DOT)) {
fprintf( fpOut, " %.16e",FixedPoint.f3dot);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (BinarySearch()) {
if (FixedDirection.at(PID_ASINI)) {
fprintf( fpOut, " %.16e",FixedPoint.asini);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_ARGP)) {
fprintf( fpOut, " %.16e",FixedPoint.argp);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_ECC)) {
fprintf( fpOut, " %.16e",FixedPoint.ecc);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_TPSSB)) {
fprintf( fpOut, " %.16e",XLALGPSGetREAL8(&FixedPoint.TpSSB));
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
if (FixedDirection.at(PID_PERIOD)) {
fprintf( fpOut, " %.16e",FixedPoint.period);
}
else {
fprintf( fpOut, " %.16e", directions.at(n).max);
n++;
}
}
fprintf( fpOut,"\n");
if (BinarySearch()) {
fprintf( fpOut, "%%%% nomad_lower: %.16f %.16f %.16f %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",best_cell.lower_bound[0],best_cell.lower_bound[1],best_cell.lower_bound[2],best_cell.lower_bound[3],best_cell.lower_bound[4],best_cell.lower_bound[5],best_cell.lower_bound[6],best_cell.lower_bound[7],best_cell.lower_bound[8],best_cell.lower_bound[9],best_cell.lower_bound[10]);
fprintf( fpOut, "%%%% nomad_start: %.16f %.16f %.16f %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n",best_cell.start_point[0],best_cell.start_point[1],best_cell.start_point[2],best_cell.start_point[3],best_cell.start_point[4],best_cell.start_point[5],best_cell.start_point[6],best_cell.start_point[7],best_cell.start_point[8],best_cell.start_point[9],best_cell.start_point[10]);
fprintf( fpOut, "%%%% nomad_upper: %.16f %.16f %.16f %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e\n\n",best_cell.upper_bound[0],best_cell.upper_bound[1],best_cell.upper_bound[2],best_cell.upper_bound[3],best_cell.upper_bound[4],best_cell.upper_bound[5],best_cell.upper_bound[6],best_cell.upper_bound[7],best_cell.upper_bound[8],best_cell.upper_bound[9],best_cell.upper_bound[10]);
}
else {
fprintf( fpOut, "%%%% nomad_lower: %.16f %.16f %.16f %.16e %.16e %.16e\n",best_cell.lower_bound[0],best_cell.lower_bound[1],best_cell.lower_bound[2],best_cell.lower_bound[3],best_cell.lower_bound[4],best_cell.lower_bound[5]);
fprintf( fpOut, "%%%% nomad_start: %.16f %.16f %.16f %.16e %.16e %.16e\n",best_cell.start_point[0],best_cell.start_point[1],best_cell.start_point[2],best_cell.start_point[3],best_cell.start_point[4],best_cell.start_point[5]);
fprintf( fpOut, "%%%% nomad_upper: %.16f %.16f %.16f %.16e %.16e %.16e\n\n",best_cell.upper_bound[0],best_cell.upper_bound[1],best_cell.upper_bound[2],best_cell.upper_bound[3],best_cell.upper_bound[4],best_cell.upper_bound[5]);
}
for (thispoint=startpoint; thispoint<epv.size(); thispoint++) {
PulsarDopplerParams XLAL_INIT_DECL(cPoint);
cPoint.Alpha = epv.at(thispoint).alpha;
cPoint.Delta = epv.at(thispoint).delta;
cPoint.fkdot[0] = epv.at(thispoint).freq;
cPoint.fkdot[1] = epv.at(thispoint).f1dot;
cPoint.fkdot[2] = epv.at(thispoint).f2dot;
cPoint.fkdot[3] = epv.at(thispoint).f3dot;
if (BinarySearch())
{
cPoint.tp.gpsSeconds = epv.at(thispoint).TpSSB.gpsSeconds;
cPoint.tp.gpsNanoSeconds = epv.at(thispoint).TpSSB.gpsNanoSeconds;
cPoint.argp = epv.at(thispoint).argp;
cPoint.asini = epv.at(thispoint).asini;
cPoint.ecc = epv.at(thispoint).ecc;
cPoint.period = epv.at(thispoint).period;
}
cPoint.refTime = refTimeGPS;
REAL8 Fi[7];
for (INT4 i = 0; i < 7; i++) {
Fi[i] = 0;
}
FstatQuantities Fstat_what;
if (fstatsearchtype == SCSEARCH) {
for (UINT4 k = 0; k < nStacks; k++) {
Fstat_what = FSTATQ_2F;
XLALComputeFstat(&tFstat,Fstat_in_vec->data[k],&cPoint,0,1,Fstat_what);
Fi[0] += tFstat->twoF[0];
Fstat_what = FSTATQ_2F_PER_DET;
XLALComputeFstat(&tFstat,Fstat_in_vec->data[k],&cPoint,0,1,Fstat_what);
for (INT4 i = 0; i < Ndet(); i++) {
Fi[i+1] += tFstat->twoFPerDet[i][0];
}
}
for (INT4 i = 0; i <= Ndet(); i++) {
Fi[i] /= (REAL8)nStacks;
}
}
else
{
Fstat_what = FSTATQ_2F;
XLALComputeFstat(&tFstat,Fstat_in_vec->data[0],&cPoint,0,1,Fstat_what);
Fi[0] = tFstat->twoF[0];
Fstat_what = FSTATQ_2F_PER_DET;
XLALComputeFstat(&tFstat,Fstat_in_vec->data[0],&cPoint,0,1,Fstat_what);
for (INT4 i = 0; i < Ndet(); i++) {
Fi[i+1] = tFstat->twoFPerDet[i][0];
}
}
LogPrintf(LOG_NORMAL,"Write point %d of %d: ",thispoint+1,epv.size());
printf("2F: %.6g",Fi[0]);
for (INT4 i = 1; i <= Ndet(); i++) {
printf(" 2F[%d]: %.6g",i,Fi[i]);
}
printf("\n");
if (BinarySearch()) {
fprintf ( fpOut, "%.16f %.16f %.16f %.16e %.16e %.16e %.16f %.16f %.16f %.16f %.16e",epv.at(thispoint).freq, epv.at(thispoint).alpha, epv.at(thispoint).delta, epv.at(thispoint).f1dot, epv.at(thispoint).f2dot, epv.at(thispoint).f3dot, epv.at(thispoint).period, epv.at(thispoint).asini, XLALGPSGetREAL8(&(epv.at(thispoint).TpSSB)), epv.at(thispoint).argp, epv.at(thispoint).ecc);
}
else {
fprintf ( fpOut, "%.16f %.16f %.16f %.16e %.16e %.16e 1",epv.at(thispoint).freq,epv.at(thispoint).alpha,epv.at(thispoint).delta,epv.at(thispoint).f1dot,epv.at(thispoint).f2dot,epv.at(thispoint).f3dot);
}
for (INT4 i = 0; i <= Ndet(); i++) {
fprintf( fpOut, " %.9g",Fi[i]);
}
fprintf( fpOut, "\n");
}
fclose (fpOut);
}
INT4 FStatNomad::write_targeted_output(ExPoint *p) {
FILE *fpOut;
if ((fpOut = fopen(fnameout(), "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", fnameout());
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
/* get the log string */
fprintf( fpOut, "%%%% User Input:\n");
fprintf( fpOut, "%%%%-------------------------------------------\n");
fprintf( fpOut, "%%%% cmdline: %s\n", LogLine() );
/* add code version ID (only useful for git-derived versions) */
fprintf ( fpOut, "%%%% version: %s\n", VCSInfoString );
if (BinarySearch()) {
fprintf ( fpOut, "%.16f %.16f %.16f %.16e %.16e %.16e %.16f %.16f %.16f %.16f %.16e %.16f\n",p->freq, p->alpha, p->delta, p->f1dot, p->f2dot, p->f3dot, p->period, p->asini, XLALGPSGetREAL8(&p->TpSSB), p->argp, p->ecc, p->stat);
}
else {
fprintf ( fpOut, "%.16g %.16g %.16g %.16e %.16e %.16e %.16g\n",p->freq, p->alpha, p->delta, p->f1dot, p->f2dot, p->f3dot, p->stat);
}
fclose (fpOut);
return 0;
}
INT4 FStatNomad::write_output(INT4 fstatsearchtype,REAL8 freq,REAL8 alpha,REAL8 delta,REAL8 f1dot,REAL8 f2dot,REAL8 f3dot,REAL8 twoF) {
FILE *fpOut;
if ((fpOut = fopen(ToplistFile(fstatsearchtype,0), "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", ToplistFile(fstatsearchtype,0));
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
/* get the log string */
fprintf( fpOut, "%%%% User Input:\n");
fprintf( fpOut, "%%%%-------------------------------------------\n");
fprintf( fpOut, "%%%% cmdline: %s\n", LogLine() );
/* add code version ID (only useful for git-derived versions) */
fprintf ( fpOut, "%%%% version: %s\n", VCSInfoString );
fprintf ( fpOut, "%.16g %.16g %.16g %.16e %.16e %.16e %.16g\n",freq, alpha, delta, f1dot, f2dot, f3dot, fabs(twoF));
fclose (fpOut);
return 0;
}
void FStatNomad::update_to_final_numsteps()
{
if (LoopMeshCoarseningExponent()) {
numsteps = (FinalMaxMeshCoarseningExponent() - FinalMinMeshCoarseningExponent())/FinalStepMeshCoarseningExponent();
}
else {
numsteps = 1;
}
totalnumsteps = numsteps;
if (LTMADS() > 0) {
totalnumsteps = 2*totalnumsteps;
}
if (ReSearch()) {
totalnumsteps = 2*totalnumsteps;
}
printf("totalnumsteps: %d, numsteps: %d\n",totalnumsteps,numsteps);
}
void FStatNomad::ComputeSteps(INT4 fstatsearchtype) {
INT4 multinumsteps = 0;
if ((fstatsearchtype == SCSEARCH) || (fstatsearchtype == FCSEARCH)) {
if (LoopMeshCoarseningExponent()) {
numsteps = (MaxMeshCoarseningExponent() - MinMeshCoarseningExponent())/StepMeshCoarseningExponent();
}
else {
numsteps = 1;
}
totalnumsteps = numsteps;
if (LTMADS() > 0) {
totalnumsteps = 2*totalnumsteps;
}
if (ReSearch()) {
totalnumsteps = 2*totalnumsteps;
}
}
else if (fstatsearchtype == CHSEARCHSC) {
numsteps = (MaxMeshCoarseningExponent() - MinMeshCoarseningExponent())/StepMeshCoarseningExponent();
// numsteps = 1;
totalnumsteps = numsteps;
if (LTMADS() > 0) {
totalnumsteps = 2*totalnumsteps;
}
if (ReSearch()) {
totalnumsteps = 2*totalnumsteps;
}
}
else if (fstatsearchtype == CHSEARCHFC) {
numsteps = (FinalMaxMeshCoarseningExponent() - FinalMinMeshCoarseningExponent())/FinalStepMeshCoarseningExponent();
// numsteps = 1;
totalnumsteps = numsteps;
if (FinalLTMADS() > 0) {
totalnumsteps = 2*totalnumsteps;
}
if (FinalReSearch()) {
totalnumsteps = 2*totalnumsteps;
}
}
printf("totalnumsteps: %d, multinumsteps: %d, numsteps: %d\n",totalnumsteps,multinumsteps,numsteps);
countstep = 0;
trials_reserve(totalnumsteps);
}
void FStatNomad::init_SP(INT4 fstatsearchtype, LALSegList *segmentList) {
const char *fn = __func__;
/* temp loop variables: generally k loops over segments and j over SFTs in a stack */
UINT4 j;
UINT4 k;
/* copy user specified spin variables at reftime */
/* the reference time value in spinRange_refTime will be set in SetUpSFTs() */
SetSpindownRange(&(usefulParams.spinRange_refTime));
/* General GPS times */
XLAL_INIT_DECL(refTimeGPS);
XLAL_INIT_DECL(tStartGPS);
multiSFTs = NULL;
SFTtype *firstSFT;
/* F-statistic computation related stuff */
UINT4 binsFstat1, binsFstatSearch=0;
/* temporary storage for spinrange vector */
/* for 1st stage: read sfts, calculate detector states */
LogPrintf( LOG_NORMAL,"Reading input data ... ");
LALStatus lstatus = blank_status;
LAL_CALL( SetUpSFTs( &lstatus, &Fstat_in_vec, &usefulParams, segmentList), &status);
LogPrintfVerbatim ( LOG_NORMAL, " done.\n");
/* some useful params computed by SetUpSFTs */
var_tStack = usefulParams.tStack;
tObs = usefulParams.tObs;
if (fstatsearchtype == SCSEARCH) {
nStacks = segmentList->length;
}
else {
nStacks = 1;
}
tStartGPS = usefulParams.tStartGPS;
refTimeGPS = usefulParams.spinRange_refTime.refTime;
LogPrintf(LOG_DEBUG, "%% --- GPS reference time = %.4f , GPS data mid time = %.4f\n",
XLALGPSGetREAL8(&refTimeGPS), XLALGPSGetREAL8(&tMidGPS) );
/* print debug info about stacks */
LogPrintf(LOG_DEBUG, "%% --- Setup, N = %d, T = %.0fs, Tobs = %.0fs\n",
nStacks, var_tStack, tObs);
lstatus = blank_status;
INITSTATUS(&lstatus);
ATTATCHSTATUSPTR (&lstatus);
objstatus = &lstatus;
comp_exp2F();
comp_sigma2F();
Ndet(usefulParams.detectorIDs->length);
LogPrintf(LOG_NORMAL,"number of detectors: %d\n",Ndet());
LogPrintf(LOG_DEBUG,"Nseg: %d\texp2F: %f\tsigma2F: %f\n",INT4(Nseg()),exp2F(),sigma2F());
return 0;
} /* SEARCH_SP */
/* ################################### SEARCH_SP ############################### */
void FStatNomad::DoTargetedSearch(INT4 fstatsearchtype) {
REAL8 this_point_avg2F = 0.;
ExPoint p;
SetPoint(&p);
if (fstatsearchtype == FCSEARCH) {
get_FC2F(&p);
PrintPoint(Nseg(),&p,"Targeted search:",BinarySearch());
p.stat = -p.stat;
if (outputLoudest()) {
vector<ExPoint> epv;
epv.reserve(1);
ExPoint tExPoint;
SetPoint(&tExPoint);
epv.push_back(tExPoint);
write_loudest(epv);
}
}
else if (fstatsearchtype == SCSEARCH) {
get_SC2F(&p);
PrintPoint(Nseg(),&p,"Targeted search:",BinarySearch());
p.stat = -p.stat / Nseg();
}
write_targeted_output(&p);
free_memory();
LALCheckMemoryLeaks();
// return FSTATNOMAD_ENORM;
} /* SEARCH_SP */
void FStatNomad::DoSignalSearch(INT4 fstatsearchtype, INT4 cnt, LALSegList *segList) {
init_SP(fstatsearchtype,segList);
REAL8 this_point_avg2F = 0.;
ExPoint p = SignalPoint;
if (fstatsearchtype == SCSEARCH) {
get_SC2F(&p);
this_point_avg2F = -p.stat/Nseg();
PrintPoint(INT4(Nseg()),&p,"DoSignalSearchSC",BinarySearch());
}
else if (fstatsearchtype == FCSEARCH) {
get_FC2F(&p);
this_point_avg2F = -p.stat;
PrintPoint(1.,&p,"DoSignalSearchFC",BinarySearch());
}
FILE *fpOut;
if ((fpOut = fopen(ControlFile(fstatsearchtype,cnt), "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", ControlFile(fstatsearchtype,cnt));
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
if (fstatsearchtype == SCSEARCH) {
fprintf ( fpOut, "%%%% number of segments: %d\n",(INT4)Nseg());
}
else if (fstatsearchtype == FCSEARCH) {
fprintf ( fpOut, "%%%% number of segments: 1\n");
}
if (BinarySearch()) {
fprintf ( fpOut, "%.16g %.16g %.16g %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16f\n",SignalPoint.freq, SignalPoint.alpha, SignalPoint.delta, SignalPoint.f1dot, SignalPoint.f2dot, SignalPoint.f3dot,XLALGPSGetREAL8(&SignalPoint.TpSSB),SignalPoint.argp,SignalPoint.asini,SignalPoint.ecc,SignalPoint.period, this_point_avg2F);
}
else {
fprintf ( fpOut, "%.16g %.16g %.16g %.16e %.16e %.16e %.16f\n",SignalPoint.freq, SignalPoint.alpha, SignalPoint.delta, SignalPoint.f1dot, SignalPoint.f2dot, SignalPoint.f3dot, this_point_avg2F);
}
fclose (fpOut);
free_memory();
}
void FStatNomad::DoStartPointSearch(INT4 fstatsearchtype, INT4 cnt, LALSegList *segList) {
init_SP(fstatsearchtype,segList);
ExPoint p;
SetPoint(&p);
REAL8 this_point_avg2F = 0.;
if (fstatsearchtype == SCSEARCH) {
get_SC2F(&p);
this_point_avg2F = -p.stat/Nseg();
PrintPoint(INT4(Nseg()),&p,"DoStartPointSearchSC",BinarySearch());
}
else if (fstatsearchtype == FCSEARCH) {
get_FC2F(&p);
this_point_avg2F = -p.stat;
PrintPoint(1.,&p,"DoStartPointSearchFC",BinarySearch());
}
FILE *fpOut;
if ((fpOut = fopen(StartPointFile(fstatsearchtype,cnt), "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", ControlFile(fstatsearchtype,cnt));
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
if (fstatsearchtype == SCSEARCH) {
fprintf ( fpOut, "%%%% number of segments: %d\n",(INT4)Nseg());
}
else if (fstatsearchtype == FCSEARCH) {
fprintf ( fpOut, "%%%% number of segments: 1\n");
}
if (BinarySearch()) {
fprintf ( fpOut, "%.16g %.16g %.16g %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16e %.16f\n",p.freq, p.alpha, p.delta, p.f1dot, p.f2dot, p.f3dot,p.TpSSB,p.argp,p.asini,p.ecc,p.period, this_point_avg2F);
}
else {
fprintf ( fpOut, "%.16g %.16g %.16g %.16e %.16e %.16e %.16f\n",p.freq, p.alpha, p.delta, p.f1dot, p.f2dot, p.f3dot, this_point_avg2F);
}
fclose (fpOut);
free_memory();
}
void FStatNomad::free_memory() {
/* free first stage memory */
int k;
// for ( k = 0; k < nStacks; k++) {
// XLALDestroyMultiSFTVector(stackMultiSFT.data[k]);
// XLALDestroyMultiNoiseWeights(stackMultiNoiseWeights.data[k]);
// XLALDestroyMultiDetectorStateSeries ( stackMultiDetStates.data[k] );
// }
// if (stackMultiSFT.data)
// LALFree(stackMultiSFT.data);
// if (stackMultiNoiseWeights.data)
// LALFree(stackMultiNoiseWeights.data);
// if (stackMultiDetStates.data)
// LALFree(stackMultiDetStates.data);
XLALDestroyTimestampVector(startTstack);
XLALDestroyTimestampVector(midTstack);
XLALDestroyTimestampVector(endTstack);
XLALDestroyFstatInputVector(Fstat_in_vec);
/* free dopplerscan stuff */
// for (k = 0; k < nStacks; k++)
// {
// if ( catalogSeq.data[k].length > 0 ) {
// LALFree(catalogSeq.data[k].data);
// } /* end if */
// } /* loop over stacks */
// LALFree( catalogSeq.data);
} /* FREE_MEMORY */
/** Set up stacks, read SFTs, calculate SFT noise weights and calculate
detector-state */
void FStatNomad::SetUpSFTs( LALStatus *status, /**< pointer to LALStatus structure */
FstatInputVector** p_Fstat_in_vec, /**< pointer to vector of Fstat input data structures for XLALComputeFstat(), one per stack */
UsefulStageVariables *in, /**< input params */
LALSegList *segList)
{
static SFTConstraints constraints;
REAL8 timebase, tObs, deltaFsft;
UINT4 k,numSFT;
REAL8 midTseg,startTseg,endTseg;
REAL8 doppWings;
INT4 extraBins;
INT4 sft_check_result = 0;
INITSTATUS(status);
ATTATCHSTATUSPTR (status);
/* get sft catalog */
constraints.minStartTime = &(in->minStartTimeGPS);
constraints.maxStartTime = &(in->maxEndTimeGPS);
if (catalog) {
LogPrintf(LOG_NORMAL,"Reuse SFT catalog of length: %d\n",catalog->length);
}
else {
TRY( LALSFTdataFind( status->statusPtr, &catalog, in->sftbasename, &constraints), status);
}
/* check CRC sums of SFTs */
TRY ( LALCheckSFTCatalog ( status->statusPtr, &sft_check_result, catalog ), status );
if (sft_check_result) {
LogPrintf(LOG_CRITICAL,"SFT validity check failed (%d)\n", sft_check_result);
ABORT ( status, FSTATFCNOMAD_ESFT, FSTATFCNOMAD_MSGESFT );
}
/* set some sft parameters */
deltaFsft = catalog->data[0].header.deltaF;
timebase = 1.0/deltaFsft;
/* get sft catalogs for each stack */
// if ( in->segmentList ) /* if segment list was given by user */
// {
SFTCatalogSequence *catalogSeq_p;
if ( (catalogSeq_p = XLALSetUpStacksFromSegmentList ( catalog, segList )) == NULL ) {
XLALPrintError ( "%s: XLALSetUpStacksFromSegmentList() failed to set up segments from given list.\n", __func__ );
ABORT ( status, FSTATFCNOMAD_ESUB, FSTATFCNOMAD_MSGESUB );
}
catalogSeq = (*catalogSeq_p);/* copy top-level struct */
XLALFree ( catalogSeq_p ); /* free alloc'ed top-level struct after copying (contents survive in catalogSeq!) */
/* we need to set tStack here:
* this will be used for setting Freq,f1dot resolution on segments, therefore we use the longest segment duration
*/
UINT4 iSeg;
REAL8 maxT = 0;
for ( iSeg=0; iSeg < segList->length; iSeg++)
{
REAL8 T = XLALGPSDiff ( &(segList->segs[iSeg].end), &(segList->segs[iSeg].start) );
maxT = HSMAX ( maxT, T );
}
in->tStack = maxT;
// }
// else /* set up nStacks segments of fixed span tStack */
// {
// TRY( SetUpStacks( status->statusPtr, &catalogSeq, in->tStack, catalog, in->nStacks), status);
// }
/* reset number of stacks */
UINT4 numSegments = catalogSeq.length;
in->nStacks = numSegments;
/* calculate start and end times and tobs from segmented catalog*/
tStartGPS = catalogSeq.data[0].data[0].header.epoch;
in->tStartGPS = tStartGPS;
SFTCatalog *LastSegmentCat = &(catalogSeq.data[numSegments - 1]);
UINT4 numSFTsInLastSeg = LastSegmentCat->length;
tEndGPS = LastSegmentCat->data[numSFTsInLastSeg-1].header.epoch;
XLALGPSAdd(&tEndGPS, timebase);
tObs = XLALGPSDiff(&tEndGPS, &tStartGPS);
in->tObs = tObs;
/* get timestamps of start, mid and end times of each stack */
/* set up vector containing mid times of stacks */
midTstack = XLALCreateTimestampVector ( in->nStacks );
/* set up vector containing start times of stacks */
startTstack = XLALCreateTimestampVector ( in->nStacks );
/* set up vector containing end times of stacks */
endTstack = XLALCreateTimestampVector ( in->nStacks );
/* now loop over stacks and get time stamps */
for (k = 0; k < in->nStacks; k++) {
if ( catalogSeq.data[k].length == 0 ) {
/* something is wrong */
ABORT ( status, FSTATFCNOMAD_EVAL, FSTATFCNOMAD_MSGEVAL );
}
/* start time of stack = time of first sft in stack */
startTstack->data[k] = catalogSeq.data[k].data[0].header.epoch;
/* end time of stack = time of last sft in stack */
numSFT = catalogSeq.data[k].length;
endTstack->data[k] = catalogSeq.data[k].data[numSFT - 1].header.epoch;
/* reference time for Fstat has to be midpoint of segment */
startTstackGPS = startTstack->data[k];
endTstackGPS = endTstack->data[k];
startTseg = XLALGPSGetREAL8( &startTstackGPS );
endTseg = XLALGPSGetREAL8( &endTstackGPS );
/*
TRY ( LALGPStoFloat( status->statusPtr, &startTseg, &startTstackGPS ), status);
TRY ( LALGPStoFloat( status->statusPtr, &endTseg, &endTstackGPS ), status);
*/
midTseg = startTseg + ((endTseg - startTseg + timebase)*0.5);
XLALGPSSetREAL8( &midTstackGPS, midTseg );
/*
TRY ( LALFloatToGPS( status->statusPtr, &midTstackGPS, &midTseg), status);
*/
midTstack->data[k] = midTstackGPS;
} /* loop over k */
/* set reference time for pulsar parameters */
/* first calculate the mid time of observation time span*/
{
REAL8 tStart8, tEnd8, tMid8;
tStart8 = XLALGPSGetREAL8( &tStartGPS );
tEnd8 = XLALGPSGetREAL8( &tEndGPS );
tMid8 = 0.5 * (tStart8 + tEnd8);
XLALGPSSetREAL8( &tMidGPS, tMid8 );
}
if ( in->refTime > 0 ) {
REAL8 refTime = in->refTime;
XLALGPSSetREAL8(&refTimeGPS, refTime);
}
else { /* set refTime to exact midtime of the total observation-time spanned */
refTimeGPS = tMidGPS;
}
/* get frequency and fdot bands at start time of sfts by extrapolating from reftime */
in->spinRange_refTime.refTime = refTimeGPS;
TRY( LALExtrapolatePulsarSpinRange( status->statusPtr, &in->spinRange_startTime, tStartGPS, &in->spinRange_refTime), status);
TRY( LALExtrapolatePulsarSpinRange( status->statusPtr, &in->spinRange_endTime, tEndGPS, &in->spinRange_refTime), status);
TRY( LALExtrapolatePulsarSpinRange( status->statusPtr, &in->spinRange_midTime, tMidGPS, &in->spinRange_refTime), status);
/* set Fstat spindown resolution (coarse grid) */
in->df1dot = HSMIN(in->df1dot, in->spinRange_midTime.fkdotBand[1]);
/* set Fstat 2nd spindown resolution (coarse grid) */
in->df2dot = HSMIN(in->df2dot, in->spinRange_midTime.fkdotBand[2]);
/* calculate number of bins for Fstat overhead due to residual spin-down */
in->extraBinsFstat = (UINT4)( 0.25*(in->tObs*in->df1dot + in->tObs*in->tObs*in->df2dot)/in->dFreqStack + 1e-6) + 1;
/* set wings of sfts to be read */
/* the wings must be enough for the Doppler shift and extra bins
for the running median block size and Dterms for Fstat calculation.
In addition, it must also include wings for the spindown correcting
for the reference time */
/* calculate Doppler wings at the highest frequency */
startTime_freqLo = in->spinRange_startTime.fkdot[0]; /* lowest search freq at start time */
startTime_freqHi = startTime_freqLo + in->spinRange_startTime.fkdotBand[0]; /* highest search freq. at start time*/
endTime_freqLo = in->spinRange_endTime.fkdot[0];
endTime_freqHi = endTime_freqLo + in->spinRange_endTime.fkdotBand[0];
freqLo = HSMIN ( startTime_freqLo, endTime_freqLo );
freqHi = HSMAX ( startTime_freqHi, endTime_freqHi );
LogPrintf(LOG_DEBUG,"freqLo: %f\n",freqLo);
LogPrintf(LOG_DEBUG,"freqHi: %f\n",freqHi);
if (BinarySearch()) {
doppWings = freqHi * GetMaxDopplerShift(GetUpperBoundFromDirection(PID_ASINI),GetLowerBoundFromDirection(PID_PERIOD),GetUpperBoundFromDirection(PID_ECC));
}
else {
doppWings = freqHi * in->dopplerMax; /* maximum Doppler wing -- probably larger than it has to be */
}
extraBins = HSMAX ( in->blocksRngMed/2 + 1, in->Dterms );
LogPrintf(LOG_DEBUG,"doppWings: %f\n",doppWings);
freqmin = freqLo - doppWings - extraBins * deltaFsft - in->extraBinsFstat * in->dFreqStack;
freqmax = freqHi + doppWings + extraBins * deltaFsft + in->extraBinsFstat * in->dFreqStack;
/* fill detector name vector with all detectors present in any data sements */
in->detectorIDs = NULL;
for (k = 0; k < in->nStacks; k++) {
if ( ( in->detectorIDs = XLALGetDetectorIDsFromSFTCatalog ( in->detectorIDs, catalogSeq.data + k ) ) == NULL ) {
ABORT ( status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
}
}
const UINT4 numDetectors = in->detectorIDs->length;
/* set up vector of Fstat input data structs */
(*p_Fstat_in_vec) = XLALCreateFstatInputVector( in->nStacks );
if ( (*p_Fstat_in_vec) == NULL ) {
ABORT ( status, FSTATFCNOMAD_EMEM, FSTATFCNOMAD_MSGEMEM );
}
/* loop over segments and read sfts */
in->nSFTs = 0;
for (UINT4 X = 0; X < numDetectors; X++) {
in->NSegmentsInvX[X] = 0;
}
FstatExtraParams XLAL_INIT_DECL(extraParams);
extraParams.SSBprec = in->SSBprec;
extraParams.Dterms = in->Dterms;
for (k = 0; k < in->nStacks; k++) {
MultiNoiseFloor assumeSqrtSX, *p_assumeSqrtSX;
p_assumeSqrtSX = NULL;
PulsarParamsVector *injectSources = NULL;
MultiNoiseFloor *injectSqrtSX = NULL;
(*p_Fstat_in_vec)->data[k] = XLALCreateFstatInput( &catalogSeq.data[k],freqmin, freqmax, injectSources, injectSqrtSX, p_assumeSqrtSX, in->blocksRngMed, in->edat, in->Fmethod, &extraParams );
if ( (*p_Fstat_in_vec)->data[k] == NULL ){
XLALPrintError("%s: XLALSetupFstatInput() failed with errno=%d", __func__, xlalErrno);
ABORT ( status, FSTATFCNOMAD_EXLAL, FSTATFCNOMAD_MSGEXLAL );
}
/* get SFT detectors and timestamps */
const MultiLALDetector *multiIFO = XLALGetFstatInputDetectors( (*p_Fstat_in_vec)->data[k] );
if ( multiIFO == NULL ) {
XLALPrintError("%s: XLALGetFstatInputDetectors() failed with errno=%d", __func__, xlalErrno);
ABORT ( status, FSTATFCNOMAD_EXLAL, FSTATFCNOMAD_MSGEXLAL );
}
const MultiLIGOTimeGPSVector *multiTS = XLALGetFstatInputTimestamps( (*p_Fstat_in_vec)->data[k] );
if ( multiTS == NULL ) {
XLALPrintError("%s: XLALGetFstatInputTimestamps() failed with errno=%d", __func__, xlalErrno);
ABORT ( status, FSTATFCNOMAD_EXLAL, FSTATFCNOMAD_MSGEXLAL );
}
/* ----- get effective inverse number of segments per detector (needed for correct averaging in single-IFO F calculation) ----- */
for (UINT4 X = 0; X < numDetectors; X++) {
/* for each detector, check if present in each segment, and save the number of segments where it is */
for (UINT4 Y = 0; Y < multiTS->length; Y++) {
if ( strcmp( multiIFO->sites[Y].frDetector.prefix, in->detectorIDs->data[X] ) == 0 )
in->NSegmentsInvX[X] += 1;
} /* for Y < numDetectors */
} /* for X < numDetectors */
/* ----- print debug info about SFTs in this stack ----- */
LogPrintf(LOG_DETAIL, "Segment %d ", k+1);
for ( UINT4 j = 0; j < multiIFO->length; j++) {
LogPrintfVerbatim(LOG_DETAIL, "%s: %d ", multiIFO->sites[j].frDetector.prefix, multiTS->data[j]->length);
}
LogPrintfVerbatim(LOG_DETAIL, "\n");
/* ----- count the total and per-segment number of SFTs used ----- */
UINT4 nSFTsInSeg = 0;
for ( UINT4 X = 0; X < multiTS->length; ++X ) {
nSFTsInSeg += multiTS->data[X]->length;
}
in->nSFTs += nSFTsInSeg;
/* ----- if we have a segment-list: double-check number of SFTs ----- */
if ( in->segmentList ) {
/* check the number of SFTs we found in this segment against the nominal value, stored in the segment list field 'id' */
UINT4 nSFTsExpected = in->segmentList->segs[k].id;
if ( nSFTsInSeg != nSFTsExpected ) {
XLALPrintError ("%s: Segment list seems inconsistent with data read: segment %d contains %d SFTs, should hold %d SFTs\n", __func__, k, nSFTsInSeg, nSFTsExpected );
ABORT ( status, FSTATFCNOMAD_EBAD, FSTATFCNOMAD_MSGEBAD );
}
} /* if have segmentList */
} /* loop over k */
for (UINT4 X = 0; X < numDetectors; X++) {
in->NSegmentsInvX[X] = 1.0 / in->NSegmentsInvX[X]; /* now it is the inverse number */
}
LogPrintf( LOG_NORMAL, "Number of segments: %d, total number of SFTs in segments: %d\n", in->nStacks, in->nSFTs );
DETATCHSTATUSPTR (status);
RETURN(status);
} /* SetUpSFTs */
/** \brief Breaks up input sft catalog into specified number of stacks
Loops over elements of the catalog, assigns a bin index and
allocates memory to the output catalog sequence appropriately. If
there are long gaps in the data, then some of the catalogs in the
output catalog sequence may be of zero length.
*/
void SetUpStacks(LALStatus *status, /**< pointer to LALStatus structure */
SFTCatalogSequence *out, /**< Output catalog of sfts -- one for each stack */
REAL8 tStack, /**< Output duration of each stack */
SFTCatalog *in, /**< Input sft catalog to be broken up into stacks (ordered in increasing time)*/
UINT4 nStacksMax ) /**< User specified number of stacks */
{
UINT4 j, stackCounter, length;
REAL8 tStart, thisTime;
REAL8 Tsft;
printf("SETUP STACK: %f\n",tStack);
INITSTATUS( status );
ATTATCHSTATUSPTR (status);
/* check input parameters */
ASSERT ( in != NULL, status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
ASSERT ( in->length > 0, status, FSTATFCNOMAD_EVAL, FSTATFCNOMAD_MSGEVAL );
ASSERT ( nStacksMax > 0, status, FSTATFCNOMAD_EVAL, FSTATFCNOMAD_MSGEVAL );
ASSERT ( in != NULL, status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
ASSERT ( tStack > 0, status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
ASSERT ( out != NULL, status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
/* set memory of output catalog sequence to maximum possible length */
out->length = nStacksMax;
out->data = (SFTCatalog *)LALCalloc( 1, nStacksMax * sizeof(SFTCatalog));
if ( out->data == NULL ) {
ABORT ( status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
}
Tsft = 1.0 / in->data[0].header.deltaF;
/* get first sft timestamp */
/* tStart will be start time of a given stack.
This initializes tStart to the first sft time stamp as this will
be the start time of the first stack */
tStart = XLALGPSGetREAL8(&(in->data[0].header.epoch));
/* loop over the sfts */
stackCounter = 0;
for ( j = 0; j < in->length; j++)
{
/* thisTime is current sft timestamp */
thisTime = XLALGPSGetREAL8(&(in->data[j].header.epoch));
/* if sft lies in stack duration then add
this sft to the stack. Otherwise move
on to the next stack */
if ( (thisTime - tStart + Tsft <= tStack) )
{
out->data[stackCounter].length += 1;
length = out->data[stackCounter].length;
/* realloc to increase length of catalog */
out->data[stackCounter].data = (SFTDescriptor *)LALRealloc( out->data[stackCounter].data, length * sizeof(SFTDescriptor));
if ( out->data[stackCounter].data == NULL ) {
ABORT ( status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
}
out->data[stackCounter].data[length - 1] = in->data[j];
}
else /* move onto the next stack */
{
if ( stackCounter + 1 == nStacksMax )
break;
stackCounter++;
/* reset start time of stack */
tStart = XLALGPSGetREAL8(&(in->data[j].header.epoch));
/* realloc to increase length of catalog and copy data */
out->data[stackCounter].length = 1; /* first entry in new stack */
out->data[stackCounter].data = (SFTDescriptor *)LALRealloc( out->data[stackCounter].data, sizeof(SFTDescriptor));
if ( out->data[stackCounter].data == NULL ) {
ABORT ( status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
}
out->data[stackCounter].data[0] = in->data[j];
} /* if new stack */
} /* loop over sfts */
/* realloc catalog sequence length to actual number of stacks */
out->length = stackCounter + 1;
out->data = (SFTCatalog *)LALRealloc( out->data, (stackCounter+1) * sizeof(SFTCatalog) );
if ( out->data == NULL ) {
ABORT ( status, FSTATFCNOMAD_ENULL, FSTATFCNOMAD_MSGENULL );
}
DETATCHSTATUSPTR (status);
RETURN(status);
} /* SetUpStacks() */
LALSegList *XLALSingleSegment(LALSegList* segments) {
const char *fn = __func__;
LALSegList *segList = NULL;
LALSeg thisSeg;
INT4 NSFT = 0;
LIGOTimeGPS start, end;
/* allocate and initialized segment list */
if ( (segList = (LALSegList*)XLALCalloc ( 1, sizeof(*segList) )) == NULL ) XLAL_ERROR_NULL ( XLAL_ENOMEM );
if ( XLALSegListInit ( segList ) != XLAL_SUCCESS ) XLAL_ERROR_NULL ( XLAL_EFUNC );
start = segments->segs[0].start;
end = segments->segs[segments->length - 1].end;
for (INT4 k = 0; k < segments->length; k++) {
NSFT += segments->segs[k].id;
}
if ( XLALSegSet ( &thisSeg, &start, &end, NSFT ) != XLAL_SUCCESS ) XLAL_ERROR_NULL ( XLAL_EFUNC );
if ( XLALSegListAppend ( segList, &thisSeg ) != XLAL_SUCCESS ) XLAL_ERROR_NULL ( XLAL_EFUNC );
if ( XLALSegListSort( segList ) != XLAL_SUCCESS ) XLAL_ERROR_NULL ( XLAL_EFUNC );
LogPrintf(LOG_DEBUG,"Number of single segments (must be 1!!!): %d\n",segList->length);
return segList;
}
LALSegList *XLALJoinSegments(LALSegList* segments, INT4 j) {
const char *fn = __func__;
LALSegList *segList = NULL;
LALSeg thisSeg;
INT4 jseg = 0;
INT4 c = 1;
INT4 NSFT = 0;
LIGOTimeGPS start, end;
BOOLEAN dumplastseg = TRUE;
/* allocate and initialized segment list */
if ( (segList = (LALSegList*)XLALCalloc ( 1, sizeof(*segList) )) == NULL )
XLAL_ERROR_NULL ( XLAL_ENOMEM );
if ( XLALSegListInit ( segList ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
for (INT4 k = 0; k < segments->length; k++) {
if ( c == 1) {
start = segments->segs[k].start;
end = segments->segs[k].end;
NSFT = segments->segs[k].id;
dumplastseg = TRUE;
}
else if ( c <= j) {
end = segments->segs[k].end;
NSFT += segments->segs[k].id;
}
c++;
if (c>j) {
/* we set number of SFTs as 'id' field, as we have no other use for it */
if ( XLALSegSet ( &thisSeg, &start, &end, NSFT ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
if ( XLALSegListAppend ( segList, &thisSeg ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
jseg++;
// printf("Append joined segment %d: %d %d %d\n",jseg,start.gpsSeconds,end.gpsSeconds,NSFT);
dumplastseg = FALSE;
NSFT = 0;
c = 1;
}
}
if (dumplastseg) {
if ( XLALSegSet ( &thisSeg, &start, &end, NSFT ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
if ( XLALSegListAppend ( segList, &thisSeg ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
jseg++;
// printf("Append joined segment %d: %d %d %d\n",jseg,start.gpsSeconds,end.gpsSeconds,NSFT);
}
/* sort final segment list in increasing GPS start-times */
if ( XLALSegListSort( segList ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
LogPrintf(LOG_DETAIL,"Joined number of segments: %d\n",jseg);
return segList;
}
/** Function to read a segment list from given filename, returns a *sorted* SegmentList
*
* The segment-list format parse here is consistent with Xavie's segment lists used previously
* and follows the format <repeated lines of form "startGPS endGPS duration[h] NumSFTs">,
* allowed comment-characters are '%' and '#'
*
* \note we (ab)use the integer 'id' field in LALSeg to carry the total number of SFTs
* contained in that segment. This will be used as a consistency check in
* XLALSetUpStacksFromSegmentList().
*
*/
LALSegList *
LocalXLALReadSegmentsFromFile ( const char *fname /**< name of file containing segment list */
)
{
const char *fn = __func__;
LALSegList *segList = NULL;
/** check input consistency */
if ( !fname ) {
XLALPrintError ( "%s: NULL input 'fname'", fn );
XLAL_ERROR_NULL ( XLAL_EINVAL );
}
/* read and parse segment-list file contents*/
LALParsedDataFile *flines = NULL;
if ( XLALParseDataFile ( &flines, fname ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
UINT4 numSegments = flines->lines->nTokens;
/* allocate and initialized segment list */
if ( (segList = (LALSegList*)XLALCalloc ( 1, sizeof(*segList) )) == NULL )
XLAL_ERROR_NULL ( XLAL_ENOMEM );
if ( XLALSegListInit ( segList ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
UINT4 iSeg;
for ( iSeg = 0; iSeg < numSegments; iSeg ++ )
{
REAL8 t0, t1, TspanHours;
INT4 NSFT;
LALSeg thisSeg;
int ret;
ret = sscanf ( flines->lines->tokens[iSeg], "%lf %lf %lf %d", &t0, &t1, &TspanHours, &NSFT );
if ( ret != 4 ) {
XLALPrintError ("%s: failed to parse data-line %d (%d) in segment-list %s: '%s'\n", fn, iSeg, ret, fname, flines->lines->tokens[iSeg] );
XLALSegListClear ( segList );
XLALFree ( segList );
XLALDestroyParsedDataFile ( flines );
XLAL_ERROR_NULL ( XLAL_ESYS );
}
/* check internal consistency of these numbers */
REAL8 hours = 3600.0;
if ( fabs ( t1 - t0 - TspanHours * hours ) >= 1.0 ) {
XLALPrintError ("%s: Inconsistent segment list, in line %d: t0 = %f, t1 = %f, Tspan = %f != t1 - t0 (to within 1s)\n", fn, iSeg, t0, t1, TspanHours );
XLAL_ERROR_NULL ( XLAL_EDOM );
}
LIGOTimeGPS start, end;
XLALGPSSetREAL8( &start, t0 );
XLALGPSSetREAL8( &end, t1 );
/* we set number of SFTs as 'id' field, as we have no other use for it */
if ( XLALSegSet ( &thisSeg, &start, &end, NSFT ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
if ( XLALSegListAppend ( segList, &thisSeg ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
} /* for iSeg < numSegments */
/* sort final segment list in increasing GPS start-times */
if ( XLALSegListSort( segList ) != XLAL_SUCCESS )
XLAL_ERROR_NULL ( XLAL_EFUNC );
/* free parsed segment file contents */
//if ( XLALDestroyParsedDataFile ( flines ) != XLAL_SUCCESS ) XLAL_ERROR_NULL ( XLAL_EFUNC );
return segList;
} /* XLALReadSegmentsFromFile() */
/** Set up 'segmented' SFT-catalogs for given list of segments and a total SFT-catalog.
*
* Note: this function does not allow 'empty' segments to be returned, i.e. if there is any
* segment that would contain no SFTs from the given SFT-catalog, an error is returned.
* These segment-lists are 'precomputed' and therefore one can assume that empty segments
* are not intended.
* However, the function will not complain if some SFTs from the catalog are 'unused', i.e.
* they didn't fit into any of the given segments.
*
* \note the input segment list must be sorted, otherwise an error is returned
*
*/
SFTCatalogSequence *
XLALSetUpStacksFromSegmentList ( const SFTCatalog *catalog, /**< complete list of SFTs read in */
const LALSegList *segList /**< pre-computed list of segments to split SFTs into */
)
{
const char *fn = __func__;
SFTCatalogSequence *stacks; /* output: segmented SFT-catalogs */
/* check input consistency */
if ( !catalog || !segList ) {
XLALPrintError ("%s: invalid NULL input\n", fn );
XLAL_ERROR_NULL ( XLAL_EINVAL );
}
/* check that segment list is sorted */
if ( ! segList->sorted ) {
XLALPrintError ("%s: input segment list must be sorted! -> Use XLALSegListSort()\n", fn );
XLAL_ERROR_NULL ( XLAL_EDOM );
}
UINT4 numSegments = segList->length;
UINT4 numSFTs = catalog->length;
/* set memory of output catalog sequence to maximum possible length */
if ( (stacks = (SFTCatalogSequence*)XLALCalloc ( 1, sizeof(*stacks) ) ) == NULL ) {
XLALPrintError ("%s: XLALCalloc(%d) failed.\n", fn, sizeof(*stacks) );
XLAL_ERROR_NULL ( XLAL_ENOMEM );
}
stacks->length = numSegments;
if ( (stacks->data = (SFTCatalog*)XLALCalloc( stacks->length, sizeof(*stacks->data) )) == NULL ) {
XLALPrintError ("%s: failed to allocate segmented SFT-catalog\n", fn );
XLAL_ERROR_NULL ( XLAL_ENOMEM );
}
REAL8 Tsft = 1.0 / catalog->data[0].header.deltaF;
/* Step through segment list:
* for every segment:
* - find earliest and last SFT fitting completely into this segment
* - copy this range of SFT-headers into the segmented SFT-catalog 'stacks'
*/
UINT4 iSeg;
INT4 iSFT0 = 0, iSFT1 = 0; /* indices of earliest and last SFT fitting into segment iSeg */
for ( iSeg = 0; iSeg < numSegments; iSeg ++ )
{
REAL8 tSeg0, tSeg1; /* boundaries of this segment */
LALSeg *thisSeg = &(segList->segs[iSeg]);
tSeg0 = XLALGPSGetREAL8 ( & thisSeg->start );
tSeg1 = XLALGPSGetREAL8 ( & thisSeg->end );
/* ----- find earliest SFT fitting into this segment */
iSFT0 = iSFT1; /* start from previous segment's last SFT */
while ( 1 )
{
LIGOTimeGPS gpsStart = catalog->data[iSFT0].header.epoch;
int cmp = XLALGPSInSeg ( &gpsStart, thisSeg );
if ( cmp < 0 ) /* iSFT0 lies *before* current segment => advance */
iSFT0 ++;
if ( cmp == 0 ) /* iSFT0 lies *inside* current segment ==> stop */
break;
/* if no more SFTs or iSFT0 lies *past* current segment => ERROR: empty segment! */
if ( cmp > 0 || iSFT0 == (INT4)numSFTs )
{
XLALPrintError ("%s: Empty segment! No SFTs fit into segment iSeg=%d\n", fn, iSeg );
XLAL_ERROR_NULL (XLAL_EDOM );
}
} /* while true */
/* ----- find last SFT completely fitting into this segment */
iSFT1 = iSFT0;
while ( 1 )
{
LIGOTimeGPS gpsEnd = catalog->data[iSFT1].header.epoch;
XLALGPSAdd( &gpsEnd, Tsft - 1e-3 ); /* subtract 1ms from end: segments are half-open intervals [t0, t1) */
int cmp = XLALGPSInSeg ( &gpsEnd, thisSeg );
if ( cmp < 0 ) { /* end of iSFT1 lies *before* current segment ==> something is screwed up! */
XLALPrintError ("%s: end of current SFT %d lies before current segment %d ==> code seems inconsistent!\n", fn, iSFT1, iSeg );
XLAL_ERROR_NULL (XLAL_EFAILED );
}
if ( cmp == 0 ) /* end of iSFT1 lies *inside* current segment ==> advance */
iSFT1 ++;
if ( cmp > 0 || iSFT1 == (INT4)numSFTs ) { /* last SFT reached or end of iSFT1 lies *past* current segment => step back once and stop */
iSFT1 --;
break;
}
} /* while true */
INT4 numSFTsInSeg = iSFT1 - iSFT0 + 1;
/* ----- allocate and copy this range of SFTs into the segmented catalog */
stacks->data[iSeg].length = (UINT4)numSFTsInSeg;
UINT4 size = sizeof(*stacks->data[iSeg].data);
if ( (stacks->data[iSeg].data = (SFTDescriptor*)XLALCalloc ( numSFTsInSeg, size)) == NULL ) {
XLALPrintError ("%s: failed to XLALCalloc(%d, %d)\n", fn, numSFTsInSeg, size );
XLAL_ERROR_NULL (XLAL_ENOMEM );
}
INT4 iSFT;
for ( iSFT = iSFT0; iSFT <= iSFT1; iSFT++ )
stacks->data[iSeg].data[iSFT - iSFT0] = catalog->data[iSFT];
} /* for iSeg < numSegments */
return stacks;
} /* XLALSetUpStacksFromSegmentList() */
REAL8 FStatNomad::compute_metric_distance(INT4 SType, gsl_matrix *g,ExPoint *sig,ExPoint *cand) {
REAL8 sumg = 0;
gsl_vector *c;
gsl_vector *s;
switch (SType) {
case STYPE_ISO_ALPHA_DELTA:
s = gsl_vector_calloc(3);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(sig->alpha,sig->delta));
gsl_vector_set(s,1,n3y_equ(sig->alpha,sig->delta));
gsl_vector_set(s,2,n3z_equ(sig->alpha,sig->delta));
c = gsl_vector_calloc(3);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(cand->alpha,cand->delta));
gsl_vector_set(c,1,n3y_equ(cand->alpha,cand->delta));
gsl_vector_set(c,2,n3z_equ(cand->alpha,cand->delta));
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT:
s = gsl_vector_calloc(5);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(sig->alpha,sig->delta));
gsl_vector_set(s,1,n3y_equ(sig->alpha,sig->delta));
gsl_vector_set(s,2,n3z_equ(sig->alpha,sig->delta));
gsl_vector_set(s,3,sig->freq);
gsl_vector_set(s,4,sig->f1dot);
c = gsl_vector_calloc(5);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(cand->alpha,cand->delta));
gsl_vector_set(c,1,n3y_equ(cand->alpha,cand->delta));
gsl_vector_set(c,2,n3z_equ(cand->alpha,cand->delta));
gsl_vector_set(c,3,cand->freq);
gsl_vector_set(c,4,cand->f1dot);
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT:
s = gsl_vector_calloc(6);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(sig->alpha,sig->delta));
gsl_vector_set(s,1,n3y_equ(sig->alpha,sig->delta));
gsl_vector_set(s,2,n3z_equ(sig->alpha,sig->delta));
gsl_vector_set(s,3,sig->freq);
gsl_vector_set(s,4,sig->f1dot);
gsl_vector_set(s,5,sig->f2dot);
c = gsl_vector_calloc(6);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(cand->alpha,cand->delta));
gsl_vector_set(c,1,n3y_equ(cand->alpha,cand->delta));
gsl_vector_set(c,2,n3z_equ(cand->alpha,cand->delta));
gsl_vector_set(c,3,cand->freq);
gsl_vector_set(c,4,cand->f1dot);
gsl_vector_set(c,5,cand->f2dot);
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT:
s = gsl_vector_calloc(7);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(sig->alpha,sig->delta));
gsl_vector_set(s,1,n3y_equ(sig->alpha,sig->delta));
gsl_vector_set(s,2,n3z_equ(sig->alpha,sig->delta));
gsl_vector_set(s,3,sig->freq);
gsl_vector_set(s,4,sig->f1dot);
gsl_vector_set(s,5,sig->f2dot);
gsl_vector_set(s,6,sig->f3dot);
c = gsl_vector_calloc(7);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(cand->alpha,cand->delta));
gsl_vector_set(c,1,n3y_equ(cand->alpha,cand->delta));
gsl_vector_set(c,2,n3z_equ(cand->alpha,cand->delta));
gsl_vector_set(c,3,cand->freq);
gsl_vector_set(c,4,cand->f1dot);
gsl_vector_set(c,5,cand->f2dot);
gsl_vector_set(c,6,cand->f3dot);
break;
case STYPE_ISO_FREQ_F1DOT:
s = gsl_vector_calloc(2);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,sig->freq);
gsl_vector_set(s,1,sig->f1dot);
c = gsl_vector_calloc(2);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,cand->freq);
gsl_vector_set(c,1,cand->f1dot);
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT:
s = gsl_vector_calloc(3);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,sig->freq);
gsl_vector_set(s,1,sig->f1dot);
gsl_vector_set(s,2,sig->f2dot);
c = gsl_vector_calloc(3);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,cand->freq);
gsl_vector_set(c,1,cand->f1dot);
gsl_vector_set(c,2,cand->f2dot);
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT_F3DOT:
s = gsl_vector_calloc(4);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,sig->freq);
gsl_vector_set(s,1,sig->f1dot);
gsl_vector_set(s,2,sig->f2dot);
gsl_vector_set(s,3,sig->f3dot);
c = gsl_vector_calloc(4);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,cand->freq);
gsl_vector_set(c,1,cand->f1dot);
gsl_vector_set(c,2,cand->f2dot);
gsl_vector_set(c,3,cand->f3dot);
break;
case STYPE_BNS_ALPHA_DELTA_NU_A_TASC_OMEGA_K_ETA:
s = gsl_vector_calloc(9);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(sig->alpha,sig->delta));
gsl_vector_set(s,1,n3y_equ(sig->alpha,sig->delta));
gsl_vector_set(s,2,n3z_equ(sig->alpha,sig->delta));
gsl_vector_set(s,3,sig->freq);
gsl_vector_set(s,4,sig->asini);
gsl_vector_set(s,5,bnsleltasc(XLALGPSGetREAL8(&sig->TpSSB),bnslelomega(sig->period),sig->argp));
gsl_vector_set(s,6,bnslelomega(sig->period));
gsl_vector_set(s,7,bnsk(sig->ecc,sig->argp));
gsl_vector_set(s,8,bnseta(sig->ecc,sig->argp));
c = gsl_vector_calloc(9);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(cand->alpha,cand->delta));
gsl_vector_set(c,1,n3y_equ(cand->alpha,cand->delta));
gsl_vector_set(c,2,n3z_equ(cand->alpha,cand->delta));
gsl_vector_set(c,3,cand->freq);
gsl_vector_set(c,4,cand->asini);
gsl_vector_set(c,5,bnsleltasc(XLALGPSGetREAL8(&cand->TpSSB),bnslelomega(cand->period),cand->argp));
gsl_vector_set(c,6,bnslelomega(cand->period));
gsl_vector_set(c,7,bnsk(cand->ecc,cand->argp));
gsl_vector_set(c,8,bnseta(cand->ecc,cand->argp));
break;
case STYPE_BNS_ALPHA_DELTA_NU_NU1DOT_A_TASC_OMEGA_K_ETA:
s = gsl_vector_calloc(10);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,n3x_equ(sig->alpha,sig->delta));
gsl_vector_set(s,1,n3y_equ(sig->alpha,sig->delta));
gsl_vector_set(s,2,n3z_equ(sig->alpha,sig->delta));
gsl_vector_set(s,3,sig->freq);
gsl_vector_set(s,4,sig->f1dot);
gsl_vector_set(s,5,sig->asini);
gsl_vector_set(s,6,bnsleltasc(XLALGPSGetREAL8(&sig->TpSSB),bnslelomega(sig->period),sig->argp));
gsl_vector_set(s,7,bnslelomega(sig->period));
gsl_vector_set(s,8,bnsk(sig->ecc,sig->argp));
gsl_vector_set(s,9,bnseta(sig->ecc,sig->argp));
c = gsl_vector_calloc(10);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,n3x_equ(cand->alpha,cand->delta));
gsl_vector_set(c,1,n3y_equ(cand->alpha,cand->delta));
gsl_vector_set(c,2,n3z_equ(cand->alpha,cand->delta));
gsl_vector_set(c,3,cand->freq);
gsl_vector_set(c,4,cand->f1dot);
gsl_vector_set(c,5,cand->asini);
gsl_vector_set(c,6,bnsleltasc(XLALGPSGetREAL8(&cand->TpSSB),bnslelomega(cand->period),cand->argp));
gsl_vector_set(c,7,bnslelomega(cand->period));
gsl_vector_set(c,8,bnsk(cand->ecc,cand->argp));
gsl_vector_set(c,9,bnseta(cand->ecc,cand->argp));
break;
case STYPE_BNS_NU_A_TASC_OMEGA_K_ETA:
s = gsl_vector_calloc(6);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,sig->freq);
gsl_vector_set(s,1,sig->asini);
gsl_vector_set(s,2,bnsleltasc(XLALGPSGetREAL8(&sig->TpSSB),bnslelomega(sig->period),sig->argp));
gsl_vector_set(s,3,bnslelomega(sig->period));
gsl_vector_set(s,4,bnsk(sig->ecc,sig->argp));
gsl_vector_set(s,5,bnseta(sig->ecc,sig->argp));
c = gsl_vector_calloc(6);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,cand->freq);
gsl_vector_set(c,1,cand->asini);
gsl_vector_set(c,2,bnsleltasc(XLALGPSGetREAL8(&cand->TpSSB),bnslelomega(cand->period),cand->argp));
gsl_vector_set(c,3,bnslelomega(cand->period));
gsl_vector_set(c,4,bnsk(cand->ecc,cand->argp));
gsl_vector_set(c,5,bnseta(cand->ecc,cand->argp));
break;
case STYPE_BNS_NU_NU1DOT_A_TASC_OMEGA_K_ETA:
s = gsl_vector_calloc(7);
gsl_vector_set_zero(s);
gsl_vector_set(s,0,sig->freq);
gsl_vector_set(s,1,sig->f1dot);
gsl_vector_set(s,2,sig->asini);
gsl_vector_set(s,3,bnsleltasc(XLALGPSGetREAL8(&sig->TpSSB),bnslelomega(sig->period),sig->argp));
gsl_vector_set(s,4,bnslelomega(sig->period));
gsl_vector_set(s,5,bnsk(sig->ecc,sig->argp));
gsl_vector_set(s,6,bnseta(sig->ecc,sig->argp));
c = gsl_vector_calloc(7);
gsl_vector_set_zero(c);
gsl_vector_set(c,0,cand->freq);
gsl_vector_set(c,1,cand->f1dot);
gsl_vector_set(c,2,cand->asini);
gsl_vector_set(c,3,bnsleltasc(XLALGPSGetREAL8(&cand->TpSSB),bnslelomega(cand->period),cand->argp));
gsl_vector_set(c,4,bnslelomega(cand->period));
gsl_vector_set(c,5,bnsk(cand->ecc,cand->argp));
gsl_vector_set(c,6,bnseta(cand->ecc,cand->argp));
break;
default:
LogPrintf(LOG_CRITICAL,"Unknown STYPE: %d! Abort!\n",SType);
exit(-1);
break;
}
for (INT4 i=0; i < g->size1; i++) {
for (INT4 j=0; j < g->size2; j++) {
sumg += gsl_matrix_get(g,i,j)*(gsl_vector_get(c,i)-gsl_vector_get(s,i))*(gsl_vector_get(c,j)-gsl_vector_get(s,j));
}
}
gsl_vector_free(s);
gsl_vector_free(c);
return sumg;
}
BOOLEAN FStatNomad::point_inside(INT4 t, gsl_matrix *g, gsl_vector *c, gsl_vector *s) {
REAL8 sumg = 0;
for (INT4 i=0; i < g->size1; i++) {
for (INT4 j=0; j < g->size2; j++) {
sumg += gsl_matrix_get(g,i,j)*(gsl_vector_get(c,i)-gsl_vector_get(s,i))*(gsl_vector_get(c,j)-gsl_vector_get(s,j));
}
}
#ifdef DEBUG
printf("SUMG: %.16e\n",sumg);
#endif
if (t == PITYPEMETRIC) {
if ( ( sumg >= 0 ) && ( sumg <= 1. ) ) {
return true;
}
else {
return false;
}
}
else if (t == PITYPEFISHER) {
if ( ( sumg >= 0 ) && ( sumg <= pow(PointInsideFactor(),2) ) ) {
return true;
}
else {
return false;
}
}
}
INT4 FStatNomad::compute_extends(gsl_vector *extends, INT4 ExType, INT4 SType, LALSegList *segments, ExPoint *point, REAL8 avg2F) {
gsl_matrix *m = gsl_matrix_calloc(DimMetric(SType,USETYPE_EXTENDS),DimMetric(SType,USETYPE_EXTENDS));
compute_metric(m,segments,SType,point,USETYPE_EXTENDS);
INT4 maxdim = m->size1;
INT4 signum;
gsl_matrix* im = gsl_matrix_calloc(maxdim,maxdim);
gsl_permutation * p = gsl_permutation_calloc (maxdim);
if (ExType == PITYPEFISHER) {
gsl_matrix_scale(m,(avg2F-4.)*segments->length);
}
#ifdef DEBUG
DoDisplayMetric(m);
#endif
gsl_linalg_LU_decomp(m, p, &signum);
gsl_linalg_LU_invert(m, p, im);
for (INT4 i = 0; i < maxdim; i++) {
gsl_vector_set(extends,i,sqrt(fabs(gsl_matrix_get(im,i,i))));
}
LogPrintf(LOG_LAST,"extends: ");
for (INT4 i = 0; i < maxdim; i++) {
LogPrintf(LOG_LAST," %e ",gsl_vector_get(extends,i));
}
LogPrintf(LOG_NORMAL,"\n");
return 0;
}
BOOLEAN FStatNomad::check_matrix(gsl_matrix *g) {
for (INT4 i = 0; i< g->size1; i++) {
for (INT4 j = 0; j < g->size2; j++) {
if (gsl_isnan(gsl_matrix_get(g,i,j)) || gsl_isinf(gsl_matrix_get(g,i,j))) {
LogPrintf(LOG_NORMAL,"check_matix finds nan or inf for element: %d,%d\n",i,j);
DoDisplayMetric(g);
return FALSE;
}
}
}
return TRUE;
}
INT4 FStatNomad::compute_fisher(gsl_matrix *m, REAL8 avg2F, REAL8 Nseg) {
if (fisher) gsl_matrix_free(fisher);
LogPrintf(LOG_DEBUG,"AVG2F: %f, NSEG: %f\n",avg2F,Nseg);
fisher = gsl_matrix_calloc(m->size1,m->size2);
for (INT4 i = 0; i< m->size1; i++) {
for (INT4 j = 0; j<m->size2; j++) {
gsl_matrix_set(fisher,i,j,gsl_matrix_get(m,i,j));
}
}
gsl_matrix_scale(fisher,Nseg*(avg2F-4));
return 0;
}
INT4 FStatNomad::compute_metric(gsl_matrix *g, LALSegList *seglist, INT4 SType, ExPoint *p, INT4 EType) {
DopplerMetric *metric;
if (seglist != NULL) {
XLALGPSSetREAL8( &refTimeGPS, usefulParams.refTime );
int Nseg = seglist->length;
#ifdef DEBUG
PrintPoint(Nseg,p,"compute metric in this point",BinarySearch());
#endif
LogPrintf(LOG_NORMAL,"for refTime: %d\n",p->f1dot,refTimeGPS.gpsSeconds);
int s;
ConfigVariables config = empty_ConfigVariables;
DopplerMetricParams XLAL_INIT_DECL(metricParams);
LALStringVector* coords;
switch (SType) {
case STYPE_ISO_ALPHA_DELTA:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", "freq", "f1dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", "f2dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", "freq", "f1dot", "f2dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_ALPHA_DELTA_FREQ_F1DOT_F2DOT_F3DOT:
if (EType == USETYPE_EXTENDS) {
if ( (coords = XLALCreateStringVector ( "alpha", "delta", "freq", "f1dot", "f2dot", "f3dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
else {
if ( (coords = XLALCreateStringVector ( "n3x_equ", "n3y_equ", "n3z_equ", "freq", "f1dot", "f2dot", "f3dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
}
break;
case STYPE_ISO_FREQ_F1DOT:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", "f2dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_ISO_FREQ_F1DOT_F2DOT_F3DOT:
if ( (coords = XLALCreateStringVector ( "freq", "f1dot", "f2dot", "f3dot", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_ALPHA_DELTA_NU_A_TASC_OMEGA_K_ETA:
LogPrintf (LOG_CRITICAL, "SType %d is not implemented yet\n", SType );
return -1;
break;
case STYPE_BNS_ALPHA_DELTA_NU_NU1DOT_A_TASC_OMEGA_K_ETA:
LogPrintf (LOG_CRITICAL, "SType %d is not implemented yet\n", SType );
return -1;
break;
case STYPE_BNS_NU_A_TASC_OMEGA_K_ETA:
if ( (coords = XLALCreateStringVector ( "bns_lel_nu", "bns_lel_a", "bns_lel_tasc", "bns_lel_omega", "bns_lel_k", "bns_lel_eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
case STYPE_BNS_NU_NU1DOT_A_TASC_OMEGA_K_ETA:
if ( (coords = XLALCreateStringVector ( "bns_lel_nu", "bns_lel_nu1dot", "bns_lel_a", "bns_lel_tasc", "bns_lel_omega", "bns_lel_k", "bns_lel_eta", NULL )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
return xlalErrno;
}
break;
default:
LogPrintf (LOG_CRITICAL, "Uknown SType: %d\n", SType );
return -1;
}
if ( XLALDopplerCoordinateNames2System ( &config.coordSys, coords ) ) {
LogPrintf (LOG_CRITICAL, "Call to XLALDopplerCoordinateNames2System() failed. errno = %d\n\n", xlalErrno );
return xlalErrno;
}
config.signalParams.Amp.h0 = 1.;
config.signalParams.Amp.cosi = 0.;
config.signalParams.Amp.psi = 0.;
config.signalParams.Amp.phi0 = 0.;
PulsarDopplerParams *dop = &(config.signalParams.Doppler);
XLAL_INIT_DECL(*dop);
dop->Alpha = p->alpha;
dop->Delta = p->delta;
dop->fkdot[0] = p->freq;
dop->fkdot[1] = p->f1dot;
dop->fkdot[2] = p->f2dot;
dop->fkdot[3] = p->f3dot;
if (BinarySearch()) {
dop->tp.gpsSeconds = p->TpSSB.gpsSeconds;
dop->tp.gpsNanoSeconds = p->TpSSB.gpsNanoSeconds;
dop->argp = p->argp;
dop->asini = p->asini;
dop->ecc = p->ecc;
dop->period = p->period;
}
dop->refTime = refTimeGPS;
if ( XLALParseMultiLALDetector ( &config.multiIFO, var_MetricIFO ) != XLAL_SUCCESS ) {
LogPrintf (LOG_CRITICAL, "XLALParseMultiDetectorInfo() failed to parse detector names and/or weights. errno = %d.\n\n", xlalErrno);
return xlalErrno;
}
metricParams.coordSys = config.coordSys;
metricParams.detMotionType = DETMOTION_ORBIT; /* FIXME spin + motion */
metricParams.metricType = (MetricType_t)0; /* phase metric */
metricParams.multiIFO = config.multiIFO;
metricParams.multiNoiseFloor = config.multiNoiseFloor;
metricParams.signalParams = config.signalParams;
metricParams.projectCoord = - 1; /* no projection */
metricParams.segmentList = *seglist;
/* ----- compute metric full metric + Fisher matrix ---------- */
if ( (metric = XLALDopplerFstatMetric ( &metricParams, edat )) == NULL ) {
LogPrintf (LOG_CRITICAL, "Something failed in XLALDopplerFstatMetric(). xlalErrno = %d\n\n", xlalErrno);
return xlalErrno;
}
else {
#ifdef DEBUG
DoDisplayMetric(metric->g_ij);
#endif
if (check_matrix(metric->g_ij)) {
gsl_matrix_add(g,metric->g_ij);
#ifdef DEBUG
DoDisplayMetric(g);
#endif
}
else {
LogPrintf(LOG_CRITICAL,"check_matrix failed for the metric of segment: %d. Abort!\n",s);
exit(-1);
}
}
if (check_matrix(g)) {
#ifdef DEBUG
DoDisplayMetric(g);
#endif
}
else {
LogPrintf(LOG_CRITICAL,"check_matrix failed for the rescaled metric. Abort!\n");
exit(-1);
}
}
else {
LogPrintf(LOG_NORMAL,"Cannot compute metric from NULL segment list!\n");
return -1;
}
}
/** write full 'PulsarCandidate' (i.e. Doppler params + Amplitude params + error-bars + Fa,Fb, F, + A,B,C,D
* RETURN 0 = OK, -1 = ERROR
*/
INT4 FStatNomad::write_PulsarCandidate_to_fp ( FILE *fp, const PulsarCandidate *pulsarParams, const FstatCandidate *Fcand, REAL8 *twoFi, REAL8 texp2F, REAL8 tsigma2F )
{
if ( !fp || !pulsarParams || !Fcand )
return -1;
fprintf (fp, "\n");
fprintf (fp, "refTime = %d.%d;\n", pulsarParams->Doppler.refTime.gpsSeconds,pulsarParams->Doppler.refTime.gpsNanoSeconds ); /* forget about ns... */
fprintf (fp, "\n");
/* Amplitude parameters with error-estimates */
fprintf (fp, "h0 = % .6g;\n", pulsarParams->Amp.h0 );
fprintf (fp, "dh0 = % .6g;\n", pulsarParams->dAmp.h0 );
fprintf (fp, "cosi = % .6g;\n", pulsarParams->Amp.cosi );
fprintf (fp, "dcosi = % .6g;\n", pulsarParams->dAmp.cosi );
fprintf (fp, "phi0 = % .6g;\n", pulsarParams->Amp.phi0 );
fprintf (fp, "dphi0 = % .6g;\n", pulsarParams->dAmp.phi0 );
fprintf (fp, "psi = % .6g;\n", pulsarParams->Amp.psi );
fprintf (fp, "dpsi = % .6g;\n", pulsarParams->dAmp.psi );
fprintf (fp, "\n");
/* Doppler parameters */
fprintf (fp, "Alpha = % .16g;\n", pulsarParams->Doppler.Alpha );
fprintf (fp, "Delta = % .16g;\n", pulsarParams->Doppler.Delta );
fprintf (fp, "Freq = % .16g;\n", pulsarParams->Doppler.fkdot[0] );
fprintf (fp, "f1dot = % .16g;\n", pulsarParams->Doppler.fkdot[1] );
fprintf (fp, "f2dot = % .16g;\n", pulsarParams->Doppler.fkdot[2] );
fprintf (fp, "f3dot = % .16g;\n", pulsarParams->Doppler.fkdot[3] );
fprintf (fp, "\n");
/* Binary parameters */
if (pulsarParams->Doppler.asini > 0)
{
fprintf (fp, "orbitPeriod = % .16g;\n", pulsarParams->Doppler.period );
fprintf (fp, "orbitasini = % .16g;\n", pulsarParams->Doppler.asini );
fprintf (fp, "orbitTpSSBsec = % .8d;\n", pulsarParams->Doppler.tp.gpsSeconds );
fprintf (fp, "orbitTpSSBnan = % .8d;\n", pulsarParams->Doppler.tp.gpsNanoSeconds );
fprintf (fp, "orbitArgp = % .16g;\n", pulsarParams->Doppler.argp );
fprintf (fp, "orbitEcc = % .16g;\n", pulsarParams->Doppler.ecc );
}
/* Amplitude Modulation Coefficients */
fprintf (fp, "Ad = % .6g;\n", Fcand->Mmunu.Ad );
fprintf (fp, "Bd = % .6g;\n", Fcand->Mmunu.Bd );
fprintf (fp, "Cd = % .6g;\n", Fcand->Mmunu.Cd );
fprintf (fp, "Ed = % .6g;\n", Fcand->Mmunu.Ed );
fprintf (fp, "Sinv_Tsft= % .6g;\n", Fcand->Mmunu.Sinv_Tsft );
fprintf (fp, "\n");
/* Fstat-values */
fprintf (fp, "Fa = % .6g %+.6gi;\n", creal(Fcand->Fa), cimag(Fcand->Fa) );
fprintf (fp, "Fb = % .6g %+.6gi;\n", creal(Fcand->Fb), cimag(Fcand->Fb) );
fprintf (fp, "twoF = % .6g;\n", Fcand->twoF );
fprintf (fp, "twoFi = %");
for (INT4 i = 0; i < Ndet(); i++) {
fprintf (fp, " %.6g",twoFi[i]);
}
fprintf (fp, ";\n");
fprintf (fp, "exp2F = % .6g;\n",texp2F );
fprintf (fp, "sigma2F = % .6g;\n",tsigma2F );
fprintf (fp, "\nAmpFisher = \\\n" );
XLALfprintfGSLmatrix ( fp, "%.9g",pulsarParams->AmpFisherMatrix );
return 0;
} /* write_PulsarCandidate_to_fp() */
INT4 FStatNomad::write_loudest(vector<ExPoint> epv) {
INT4 max_pos = 0;
REAL8 max_F = 0.0;
INT4 k;
for (k=0; k<epv.size(); k++) {
LogPrintf(LOG_DEBUG,"pos: %d, max_pos: %d, F: %.16f, maxF: %.16f a: %.16f, d: %.16f,f:%.16f, f1d: %.16e, f2d: %.16e f3d: %.16e\n",k,max_pos,epv.at(k).stat,max_F,epv.at(k).alpha,epv.at(k).delta,epv.at(k).freq,epv.at(k).f1dot,epv.at(k).f2dot,epv.at(k).f3dot);
if ((epv.at(k).stat - max_F) > 0) {
max_F = epv.at(k).stat;
max_pos = k;
}
}
FstatCandidate loudestFCand = empty_FstatCandidate;
PulsarDopplerParams XLAL_INIT_DECL(cPoint);
cPoint.Alpha = epv.at(max_pos).alpha;
cPoint.Delta = epv.at(max_pos).delta;
cPoint.fkdot[0] = epv.at(max_pos).freq;
cPoint.fkdot[1] = epv.at(max_pos).f1dot;
cPoint.fkdot[2] = epv.at(max_pos).f2dot;
cPoint.fkdot[3] = epv.at(max_pos).f3dot;
cPoint.refTime = refTimeGPS;
if (BinarySearch()) {
cPoint.argp = epv.at(max_pos).argp;
cPoint.asini = epv.at(max_pos).asini;
cPoint.ecc = epv.at(max_pos).ecc;
cPoint.period = epv.at(max_pos).period;
cPoint.tp = epv.at(max_pos).TpSSB;
LogPrintf(LOG_DEBUG,"orbital: argp: %f, asini: %f, ecc: %e, period: %f, tp: %f\n",cPoint.argp,cPoint.asini,cPoint.ecc,cPoint.period,XLALGPSGetREAL8(&cPoint.tp));
}
LogPrintf(LOG_DEBUG,"recompute for: %.16f %.16f %.16f %.16e %.16e %.16e\n",epv.at(max_pos).alpha,epv.at(max_pos).delta,epv.at(max_pos).freq,epv.at(max_pos).f1dot,epv.at(max_pos).f2dot,epv.at(max_pos).f3dot);
tFstat = NULL;
FstatQuantities Fstat_what = FSTATQ_2F;
Fstat_what |= FSTATQ_FAFB;
XLALComputeFstat(&tFstat,Fstat_in_vec->data[0],&cPoint,0,1,Fstat_what);
loudestFCand.twoF = tFstat->twoF[0];
loudestFCand.FaFb_refTime = refTimeGPS;
loudestFCand.doppler = cPoint;
loudestFCand.Mmunu = tFstat->Mmunu;
loudestFCand.Fa = tFstat->FaFb[0].Fa;
loudestFCand.Fb = tFstat->FaFb[0].Fb;
Fstat_what = FSTATQ_2F_PER_DET;
tFstat->numDetectors--;
XLALComputeFstat(&tFstat,Fstat_in_vec->data[0],&cPoint,0,1,Fstat_what);
for (INT4 i = 0; i < Ndet(); i++ ) {
twoFi[i] = tFstat->twoFPerDet[i][0];
}
FILE *fpLoudest;
PulsarCandidate XLAL_INIT_DECL(pulsarParams);
pulsarParams.Doppler = loudestFCand.doppler;
XLALEstimatePulsarAmplitudeParams (&pulsarParams, &loudestFCand.FaFb_refTime, loudestFCand.Fa, loudestFCand.Fb, &loudestFCand.Mmunu);
if ( (fpLoudest = fopen (LoudestFile(), "wb")) == NULL)
{
XLALPrintError ("\nError opening file '%s' for writing..\n\n", LoudestFile());
return FSTATFCNOMAD_EFILE;
}
/* write header with run-info */
// fprintf (fpLoudest, "%s", GV.logstring );
/* write this 'candidate' to disc */
if ( write_PulsarCandidate_to_fp ( fpLoudest, &pulsarParams, &loudestFCand, twoFi, exp2F(), sigma2F()) != XLAL_SUCCESS )
{
LogPrintf(LOG_CRITICAL, "call to write_PulsarCandidate_to_fp() failed!\n");
return FSTATFCNOMAD_EFILE;
}
fclose (fpLoudest);
gsl_matrix_free ( pulsarParams.AmpFisherMatrix );
}
INT4 FStatNomad::get_SC2F(ExPoint *p) {
// REAL8 timeStart2F = XLALGetTimeOfDay();
cPoint.Alpha = p->alpha;
cPoint.Delta = p->delta;
cPoint.fkdot[0] = p->freq;
cPoint.fkdot[1] = p->f1dot;// /scale_f1dot();
cPoint.fkdot[2] = p->f2dot;// /scale_f2dot();
cPoint.fkdot[3] = p->f3dot;// /scale_f3dot();
cPoint.refTime = refTimeGPS;
// LogPrintf(LOG_NORMAL,"refTime: %d, Alpha: %.16f, Delta: %.16f, Freq: %.16f f1dot: %.16e f2dot: %.16e f3dot: %.16e\n",cPoint.refTime.gpsSeconds,cPoint.Alpha,cPoint.Delta,cPoint.fkdot[0],cPoint.fkdot[1],cPoint.fkdot[2],cPoint.fkdot[3]);
if (BinarySearch())
{
cPoint.tp.gpsSeconds = p->TpSSB.gpsSeconds;
cPoint.tp.gpsNanoSeconds = p->TpSSB.gpsNanoSeconds;
cPoint.argp = p->argp;
cPoint.asini = p->asini;
cPoint.ecc = p->ecc;
cPoint.period = p->period;
// LogPrintf(LOG_NORMAL,"orbitalParams: tp: %d.%d, argp: %.16e, asini: %.16e\n, ecc: %.16e, period: %.16f\n",cPoint.orbit->tp.gpsSeconds,cPoint.orbit->tp.gpsNanoSeconds,cPoint.orbit->argp,cPoint.orbit->asini,cPoint.orbit->ecc,cPoint.orbit->period);
}
#ifdef DOTIMING
REAL8 timeStartFSTAT = XLALGetTimeOfDay();
#endif
REAL8 lF = 0;
FstatQuantities Fstat_what = FSTATQ_2F;
for (UINT4 k = 0; k < nStacks; k++) {
XLALComputeFstat(&tFstat,Fstat_in_vec->data[k],&cPoint,0,1,Fstat_what);
#ifdef DOTIMING
REAL8 timeEndFSTAT = XLALGetTimeOfDay();
printf("NSft: %d,%d\tTiming2F: %.16e,\tTimingFStat: %.16e\tc0: %e\n",stackMultiSFT.data[k]->data[0]->length,stackMultiSFT.data[k]->data[1]->length,timeEndFSTAT - timeStartFSTAT,timeEndFSTAT - timeStartFSTAT,(timeEndFSTAT - timeStartFSTAT)/(REAL8)(stackMultiSFT.data[k]->data[0]->length+stackMultiSFT.data[k]->data[1]->length));
#endif
lF += tFstat->twoF[0];
}
p->stat = -lF;
return 0;
}
INT4 FStatNomad::get_FC2F(ExPoint *p) {
cPoint.Alpha = p->alpha;
cPoint.Delta = p->delta;
cPoint.fkdot[0] = p->freq;
cPoint.fkdot[1] = p->f1dot;
cPoint.fkdot[2] = p->f2dot;
cPoint.fkdot[3] = p->f3dot;
cPoint.refTime = refTimeGPS;
if (BinarySearch())
{
cPoint.tp.gpsSeconds = p->TpSSB.gpsSeconds;
cPoint.tp.gpsNanoSeconds = p->TpSSB.gpsNanoSeconds;
cPoint.argp = p->argp;
cPoint.asini = p->asini;
cPoint.ecc = p->ecc;
cPoint.period = p->period;
}
FstatQuantities Fstat_what = FSTATQ_2F;
XLALComputeFstat(&tFstat,Fstat_in_vec->data[0],&cPoint,0,1,Fstat_what);
p->stat = -tFstat->twoF[0];
return 0;
}
REAL8 FStatNomad::Randomize(REAL8 s, REAL8 u, REAL8 l, INT4 p) {
REAL8 r;
INT4 c = 0;
REAL8 ar = 0;
while (c<10) {
ar = RandomFloat(-0.5*p/100.*(u-l),0.5*p/100.*(u-l));
r = s + ar;
if ((u - r) > 0 && (r - l) > 0) {
#ifdef DEBUG
printf("c: %d, l: %.16f, r: %.16f, u: %.16f, ar: %.16e, if: %d\n",c,l,r,u,ar,((u - r) > 0 && (r - l) > 0));
#endif
return r;
}
c++;
}
printf("return mid point: %.16g\n",l + 0.5*(u-l));
return l + 0.5*(u-l);
}
BOOLEAN FStatNomad::check_point_inside(INT4 t,gsl_matrix *g,NOMAD::Point p) {
gsl_vector *c = gsl_vector_calloc(g->size1);
gsl_vector *s = gsl_vector_calloc(g->size2);
ExPoint sp;
SetPoint(&sp);
if (BinarySearch()) {
INT4 offset = 0;
INT4 skyoffset = 0;
INT4 pskyoffset = 0;
if (BNSLEL()) {
if (!DirectedSearch()) {
gsl_vector_set(s,0,n3x_equ(sp.alpha,sp.delta));
gsl_vector_set(s,1,n3y_equ(sp.alpha,sp.delta));
gsl_vector_set(s,2,n3z_equ(sp.alpha,sp.delta));
gsl_vector_set(s,0,n3x_equ(sp.alpha,sp.delta));
gsl_vector_set(s,1,n3y_equ(sp.alpha,sp.delta));
gsl_vector_set(s,2,n3z_equ(sp.alpha,sp.delta));
skyoffset = 3;
pskyoffset = 2;
}
gsl_vector_set(s,0 + skyoffset,sp.freq);
gsl_vector_set(c,0 + skyoffset,p[0+pskyoffset].value());
if (g->size1 > 6 + skyoffset) {
gsl_vector_set(s,1 + skyoffset,sp.f1dot);
gsl_vector_set(c,1 + skyoffset,p[1+pskyoffset].value());
offset++;
}
gsl_vector_set(s,1 + skyoffset + offset,sp.asini);
gsl_vector_set(c,1 + skyoffset + offset,p[1 + pskyoffset + offset].value());
gsl_vector_set(s,2 + skyoffset + offset,bnsleltasc(XLALGPSGetREAL8(&sp.TpSSB),bnslelomega(sp.period),sp.argp));
gsl_vector_set(c,2 + skyoffset + offset,bnsleltasc(p[4 + pskyoffset + offset].value(),bnslelomega(p[5 + pskyoffset + offset].value()),p[2+pskyoffset+offset].value()));
gsl_vector_set(s,3 + skyoffset + offset,bnslelomega(sp.period));
gsl_vector_set(c,3 + skyoffset + offset,bnslelomega(p[5 + pskyoffset + offset].value()));
gsl_vector_set(s,4 + skyoffset + offset,bnsk(sp.ecc,sp.argp));
gsl_vector_set(c,4 + skyoffset + offset,bnsk(p[3 + pskyoffset + offset].value(),p[3 + pskyoffset + offset].value()));
gsl_vector_set(s,5 + skyoffset + offset,bnseta(sp.ecc,sp.argp));
gsl_vector_set(c,5 + skyoffset + offset,bnseta(p[3 + pskyoffset + offset].value(),p[3 + pskyoffset + offset].value()));
// [0:alpha,1:delta,2,0:freq,3,1:f1dot,4,2,asini,5,3,argp,6,4,ecc,7,5,tpssb,8,6,period]
}
else {
if (!DirectedSearch()) {
gsl_vector_set(c,0,n3x_equ(p[0].value(),p[1].value()));
gsl_vector_set(c,1,n3y_equ(p[0].value(),p[1].value()));
gsl_vector_set(c,2,n3z_equ(p[0].value(),p[1].value()));
gsl_vector_set(s,0,n3x_equ(sp.alpha,sp.delta));
gsl_vector_set(s,1,n3y_equ(sp.alpha,sp.delta));
gsl_vector_set(s,2,n3z_equ(sp.alpha,sp.delta));
skyoffset = 3;
}
for (INT4 i = skyoffset; i < g->size1; g++) {
gsl_vector_set(c,i,p[i].value());
}
gsl_vector_set(s,0 + skyoffset,sp.freq);
if (g->size1 > 6 + skyoffset) {
gsl_vector_set(s,1 + skyoffset,sp.f1dot);
offset++;
}
gsl_vector_set(s,1 + skyoffset + offset,sp.asini);
gsl_vector_set(s,2 + skyoffset + offset,XLALGPSGetREAL8(&sp.TpSSB));
gsl_vector_set(s,3 + skyoffset + offset,sp.period);
gsl_vector_set(s,4 + skyoffset + offset,sp.ecc);
gsl_vector_set(s,5 + skyoffset + offset,sp.argp);
}
}
else if (!DirectedSearch()) {
gsl_vector_set(c,0,n3x_equ(p[0].value(),p[1].value()));
gsl_vector_set(c,1,n3y_equ(p[0].value(),p[1].value()));
gsl_vector_set(c,2,n3z_equ(p[0].value(),p[1].value()));
gsl_vector_set(s,0,n3x_equ(sp.alpha,sp.delta));
gsl_vector_set(s,1,n3y_equ(sp.alpha,sp.delta));
gsl_vector_set(s,2,n3z_equ(sp.alpha,sp.delta));
for (INT4 i = 3; i < g->size1; i++) {
gsl_vector_set(c,i,p[i-1].value());
if (i == 3) {
gsl_vector_set(s,3,sp.freq);
}
else if (i == 4) {
gsl_vector_set(s,4,sp.f1dot);
}
else if (i == 5) {
gsl_vector_set(s,5,sp.f2dot);
}
else if (i == 6) {
gsl_vector_set(s,6,sp.f3dot);
}
}
}
else {
for (INT4 i = 0; i < g->size1; i++) {
gsl_vector_set(c,i,p[i].value());
if (i == 0) {
gsl_vector_set(s,0,sp.freq);
}
else if (i == 1) {
gsl_vector_set(s,1,sp.f1dot);
}
else if (i == 2) {
gsl_vector_set(s,2,sp.f2dot);
}
else if (i == 3) {
gsl_vector_set(s,3,sp.f3dot);
}
}
}
#ifdef DEBUG
printf("cand vec: \n");
gsl_vector_fprintf(stdout,c,"%.16e");
printf("sign vec: \n");
gsl_vector_fprintf(stdout,s,"%.16e");
#endif
return point_inside(t,g,c,s);
}
gsl_vector* FStatNomad::RandomPointInside(gsl_matrix *g, NOMAD::Point start_point,NOMAD::Point upper_bound, NOMAD::Point lower_bound, NOMAD::Point scale_point, INT4 p) {
/*
This is an example mapping for the directed case.
[0] nu <---> f <---> f [0]
[x=1] nu1dot <---> f1dot <---> f1dot [x=1]
[1+x] asini <---> asini <---> asini [1+x]
[2+x] tasc <---> (tp,omega(period),argp) <---> argp [2+x]
[3+x] omega <---> (period) <---> ecc [3+x]
[4+x] kappa <---> (ecc,argp) <---> tpssb [4+x]
[5+x] eta <---> (ecc,argp) <---> period [5+x]
*/
REAL8 alpha, delta, rnd_freq, rnd_f1dot, rnd_asini, rnd_argp, rnd_ecc, rnd_tpssb, rnd_period;
#ifdef DEBUG
DoDisplayMetric(g);
#endif
BOOLEAN search_point = TRUE;
INT4 cnt = 0;
INT4 spdoffset = 0;
gsl_vector *c = gsl_vector_calloc(g->size1);
gsl_vector *s = gsl_vector_calloc(g->size2);
while (search_point) {
if (BinarySearch()) {
spdoffset = 0;
if (DirectedSearch() && start_point.size() == 7) {
spdoffset = 1;
}
else if (!DirectedSearch() && start_point.size() == 9) {
spdoffset = 1;
}
if (BNSLEL()) {
INT4 offset = 0;
INT4 poffset = 0;
rnd_f1dot = 0;
if (!DirectedSearch()) {
alpha = Randomize(start_point[0].value(),upper_bound[0].value(),lower_bound[0].value(),p)/scale_point[0].value();
delta = Randomize(start_point[1].value(),upper_bound[1].value(),lower_bound[1].value(),p)/scale_point[1].value();
gsl_vector_set(c,0,n3x_equ(alpha/scale_point[0].value(),delta/scale_point[1].value()));
gsl_vector_set(c,1,n3y_equ(alpha/scale_point[0].value(),delta/scale_point[1].value()));
gsl_vector_set(c,2,n3z_equ(alpha/scale_point[0].value(),delta/scale_point[1].value()));
gsl_vector_set(s,0,n3x_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
gsl_vector_set(s,1,n3y_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
gsl_vector_set(s,2,n3z_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
offset = 3;
poffset = 2;
}
rnd_freq = Randomize(start_point[0 + poffset ].value(),upper_bound[0 + poffset].value(),lower_bound[0 + poffset].value(),p)/scale_point[0 + poffset].value();
if (spdoffset > 0) {
rnd_f1dot = Randomize(start_point[1 + poffset].value(),upper_bound[1 +
poffset].value(),lower_bound[1 + poffset].value(),p)/scale_point[0 + poffset].value();
}
rnd_asini = Randomize(start_point[1 + poffset + spdoffset].value(),upper_bound[1 + poffset + spdoffset].value(),lower_bound[1 + poffset + spdoffset].value(),p)/scale_point[1 + poffset + spdoffset].value();
rnd_argp = Randomize(start_point[2 + poffset + spdoffset + spdoffset].value(),upper_bound[2 + poffset + spdoffset].value(),lower_bound[2 + poffset + spdoffset].value(),p)/scale_point[2 + poffset + spdoffset].value();
rnd_ecc = Randomize(start_point[3 + poffset + spdoffset].value(),upper_bound[3 + poffset + spdoffset].value(),lower_bound[3 + poffset + spdoffset].value(),p)/scale_point[3 + poffset + spdoffset].value();
rnd_tpssb = Randomize(start_point[4 + poffset + spdoffset].value(),upper_bound[4 + poffset + spdoffset].value(),lower_bound[4 + poffset + spdoffset].value(),p)/scale_point[4 + poffset + spdoffset].value();
rnd_period = Randomize(start_point[5 + poffset + spdoffset].value(),upper_bound[5 + poffset + spdoffset].value(),lower_bound[5 + poffset + spdoffset].value(),p)/scale_point[5 + poffset + spdoffset].value();
gsl_vector_set(c,0 + offset,rnd_freq);
gsl_vector_set(s,0 + offset,start_point[0 + poffset].value()/scale_point[0 + poffset].value());
if (spdoffset > 0) {
gsl_vector_set(c,1 + offset,rnd_f1dot);
gsl_vector_set(s,1 + offset,start_point[1 + poffset].value());
}
gsl_vector_set(c,1 + spdoffset + offset,rnd_asini);
gsl_vector_set(s,1 + spdoffset + offset,start_point[1 + spdoffset + poffset].value()/scale_point[1 + spdoffset + poffset].value());
gsl_vector_set(c,2 + spdoffset + offset,bnsleltasc(rnd_tpssb,bnslelomega(rnd_period),rnd_argp));
gsl_vector_set(s,2 + spdoffset + offset,bnsleltasc(start_point[4 + spdoffset + poffset].value()/scale_point[4 + spdoffset + poffset].value(),bnslelomega(start_point[5 + spdoffset + poffset].value()/scale_point[5 + spdoffset + poffset].value()),start_point[2 + poffset + spdoffset].value()/scale_point[2 + poffset + spdoffset].value()));
gsl_vector_set(c,3 + spdoffset + offset,bnslelomega(rnd_period));
gsl_vector_set(s,3 + spdoffset + offset,bnslelomega(start_point[5 + spdoffset + poffset].value()/scale_point[5 + spdoffset + poffset].value()));
gsl_vector_set(c,4 + spdoffset + offset,bnsk(rnd_ecc,rnd_argp));
gsl_vector_set(s,4 + spdoffset + offset,bnsk(start_point[3 + spdoffset + poffset].value()/scale_point[3 + spdoffset + poffset].value(),start_point[2 + spdoffset + poffset].value()/scale_point[2 + spdoffset + poffset].value()));
gsl_vector_set(c,5 + spdoffset + offset,bnseta(rnd_ecc,rnd_argp));
gsl_vector_set(s,5 + spdoffset + offset,bnseta(start_point[3 + spdoffset + poffset].value()/scale_point[3 + spdoffset + poffset].value(),start_point[2 + spdoffset + poffset].value()/scale_point[2 + spdoffset + poffset].value()));
// [0:alpha,1:delta,2,0:freq,3,1:f1dot,4,2,asini,5,3,argp,6,4,ecc,7,5,tpssb,8,6,period]
} else {
if (!DirectedSearch()) {
alpha = Randomize(start_point[0].value(),upper_bound[0].value(),lower_bound[0].value(),p)/scale_point[0].value();
delta = Randomize(start_point[1].value(),upper_bound[1].value(),lower_bound[1].value(),p)/scale_point[1].value();
gsl_vector_set(c,0,n3x_equ(alpha,delta));
gsl_vector_set(c,1,n3y_equ(alpha,delta));
gsl_vector_set(c,2,n3z_equ(alpha,delta));
gsl_vector_set(s,0,n3x_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
gsl_vector_set(s,1,n3y_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
gsl_vector_set(s,2,n3z_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
for (INT4 i = 3; i < g->size1; i++ ) {
gsl_vector_set(c,i,Randomize(start_point[i].value(),upper_bound[i].value(),lower_bound[i].value(),p)/scale_point[i].value());
gsl_vector_set(s,i,start_point[i].value()/scale_point[i].value());
}
}
else {
for (INT4 i = 0; i < g->size1; i++) {
gsl_vector_set(c,i,Randomize(start_point[i].value(),upper_bound[i].value(),lower_bound[i].value(),p)/scale_point[i].value());
gsl_vector_set(s,i,start_point[i].value()/scale_point[i].value());
}
}
}
}
else if (!DirectedSearch()) {
alpha = Randomize(start_point[0].value(),upper_bound[0].value(),lower_bound[0].value(),p)/scale_point[0].value();
delta = Randomize(start_point[1].value(),upper_bound[1].value(),lower_bound[1].value(),p)/scale_point[1].value();
gsl_vector_set(c,0,n3x_equ(alpha,delta));
gsl_vector_set(c,1,n3y_equ(alpha,delta));
gsl_vector_set(c,2,n3z_equ(alpha,delta));
gsl_vector_set(s,0,n3x_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
gsl_vector_set(s,1,n3y_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
gsl_vector_set(s,2,n3z_equ(start_point[0].value()/scale_point[0].value(),start_point[1].value()/scale_point[1].value()));
for (INT4 i = 3; i < g->size1; i++ ) {
gsl_vector_set(c,i,Randomize(start_point[i-1].value(),upper_bound[i-1].value(),lower_bound[i-1].value(),p)/scale_point[i-1].value());
gsl_vector_set(s,i,start_point[i-1].value()/scale_point[i-1].value());
}
}
else {
for (INT4 i = 0; i < g->size1; i++) {
gsl_vector_set(c,i,Randomize(start_point[i].value(),upper_bound[i].value(),lower_bound[i].value(),p)/scale_point[i].value());
gsl_vector_set(s,i,start_point[i].value()/scale_point[i].value());
}
}
if (point_inside(RestrictionType(),g,c,s))
{
search_point = FALSE;
}
else {
cnt++;
if (cnt > 10) {
for (INT4 i = 0; i < g->size1; i++) {
gsl_vector_set(c,i,gsl_vector_get(s,i));
}
search_point = FALSE;
}
}
}
if (BinarySearch() && BNSLEL()) {
INT4 poffset = 0;
if (!DirectedSearch()) {
poffset = 2;
gsl_vector_set(c,0,alpha*scale_point[0].value());
gsl_vector_set(c,1,delta*scale_point[1].value());
}
gsl_vector_set(c,0 + poffset,rnd_freq*scale_point[0+poffset].value());
if (spdoffset > 0) {
gsl_vector_set(c,1 + poffset,rnd_f1dot*scale_point[1 + poffset].value());
}
gsl_vector_set(c,1 + poffset + spdoffset,rnd_asini*scale_point[1 + poffset + spdoffset].value());
gsl_vector_set(c,2 + poffset + spdoffset,rnd_argp*scale_point[2 + poffset + spdoffset].value());
gsl_vector_set(c,3 + poffset + spdoffset,rnd_ecc*scale_point[3 + poffset + spdoffset].value());
gsl_vector_set(c,4 + poffset + spdoffset,rnd_tpssb*scale_point[4 + poffset + spdoffset].value());
gsl_vector_set(c,5 + poffset + spdoffset,rnd_period*scale_point[5 + poffset + spdoffset].value());
}
else if (!DirectedSearch()) {
gsl_vector *tc = gsl_vector_calloc(c->size-1);
gsl_vector_set(tc,0,alpha*scale_point[0].value());
gsl_vector_set(tc,1,delta*scale_point[1].value());
for (INT4 i = 3; i < c->size; i++) {
gsl_vector_set(tc,i-1,gsl_vector_get(c,i)*scale_point[i-1].value());
}
c = tc;
}
else {
for (INT4 i = 0; i < c->size; i++) {
gsl_vector_set(c,i,gsl_vector_get(c,i)*scale_point[i].value());
}
}
#ifdef DEBUG
printf("signal...\n");
gsl_vector_fprintf(stdout,s,"%.16e");
printf("candidate...\n");
gsl_vector_fprintf(stdout,c,"%.16e");
printf("return...\n");
gsl_vector_fprintf(stdout,c,"%.16e");
#endif
return c;
}
INT4 FStatNomad::DoSearch(NOMAD::Display out,vector<ExPoint> &epv,INT4 fstatsearchtype,vector<INT4> searchtype, INT4 &stcounter, BOOLEAN randomstartpoint) {
NOMAD::Point lower_bound(directions.size());
NOMAD::Point upper_bound(directions.size());
NOMAD::Point start_point(directions.size());
NOMAD::Point scale_point(directions.size());
BOOLEAN dosearch = TRUE;
for (stcounter = 0; stcounter < searchtype.size() && dosearch; stcounter++) {
NOMAD::Parameters p ( out );
p.set_DIMENSION (directions.size()); // number of variables
vector<NOMAD::bb_output_type> bbot (1); // definition of
bbot[0] = NOMAD::OBJ; // output types
p.set_BB_OUTPUT_TYPE ( bbot );
p.set_DISPLAY_STATS ( "time - bbe ( %.16gsol ) obj" );
p.set_DISPLAY_DEGREE ( DisplayDegree() );
p.set_ADD_SEED_TO_FILE_NAMES(false);
p.set_MAX_BB_EVAL(MaxBBEval()); // the algorithm terminates after
p.set_MAX_MESH_INDEX(MaxMeshIndex());
p.set_DIRECTION_TYPE((NOMAD::direction_type)searchtype.at(stcounter));
p.set_MAX_ITERATIONS(MaxIterations());
p.set_MODEL_SEARCH(ModelSearch());
p.set_MODEL_EVAL_SORT(ModelEvalSort());
p.set_VNS_SEARCH(VnsSearch());
p.set_MESH_REFINING_EXPONENT(MeshRefiningExponent());
p.set_MESH_UPDATE_BASIS(MeshUpdateBasis());
p.set_MODEL_SEARCH_MAX_TRIAL_PTS(ModelSearchMaxTrialPts());
p.set_MODEL_QUAD_MIN_Y_SIZE(ModelQuadMinYSize());
p.set_MODEL_QUAD_MAX_Y_SIZE(ModelQuadMaxYSize());
p.set_MODEL_QUAD_RADIUS_FACTOR(ModelQuadRadiusFactor());
p.set_MODEL_SEARCH_OPTIMISTIC(ModelSearchOptimistic());
p.set_MODEL_QUAD_USE_WP(ModelQuadUseWP());
p.set_MODEL_SEARCH_PROJ_TO_MESH(ModelSearchProjToMesh());
p.set_OPPORTUNISTIC_EVAL(OpportunisticEval());
p.set_OPPORTUNISTIC_LUCKY_EVAL(OpportunisticLuckyEval());
NOMAD::Point init_mesh(directions.size());
for (INT4 i = 0; i < directions.size(); i++) {
init_mesh[i] = directions.at(i).startstep * directions.at(i).scale;
}
PrintNomadPoint(init_mesh,"init_mesh");
p.set_INITIAL_MESH_SIZE(init_mesh);
INT4 s;
for (s=0; s<numsteps; s++) {
if (SaveNomadStat()) {
p.set_HISTORY_FILE(string(HistoryFile(fstatsearchtype,0,countstep)));
p.set_SOLUTION_FILE(string(SolutionFile(fstatsearchtype,0,countstep)));
p.set_STATS_FILE(string(StatsFile(fstatsearchtype,0,countstep)),"time - bbe ( %.16gsol ) obj");
}
INT4 tMeshCoarseningExponent = MinMeshCoarseningExponent() + s*StepMeshCoarseningExponent();
LogPrintf(LOG_NORMAL,"STEP: %d of %d\tNOMAD STEP: %d\tNOMADMeshCoarseningExponent: %d\n",countstep,totalnumsteps,s,tMeshCoarseningExponent);
p.set_MESH_COARSENING_EXPONENT(tMeshCoarseningExponent);
// parameters validation:
p.set_EPSILON(Epsilon());
p.set_MIN_MESH_SIZE(MinMeshSize());
p.set_MIN_POLL_SIZE(MinPollSize());
for (INT4 i = 0; i < directions.size(); i++) {
start_point[i] = directions.at(i).startpoint * directions.at(i).scale;
lower_bound[i] = directions.at(i).min * directions.at(i).scale;
upper_bound[i] = directions.at(i).max * directions.at(i).scale;
scale_point[i] = directions.at(i).scale;
}
if ((stcounter == 1) && randomstartpoint) {
if (UseMetric()) {
gsl_vector* rpi = RandomPointInside(Semicohmetric(),start_point,upper_bound,lower_bound,scale_point,LTMADS());
for (INT4 i = 0; i < rpi->size; i++) {
#ifdef DEBUG
printf("UseMetric OLD start_point[%d]: %.16e\n",start_point[i].value());
printf("UseMetric NEW start point[%d]: %.16e\n",gsl_vector_get(rpi,i));
#endif
start_point[i] = gsl_vector_get(rpi,i);
}
}
else {
for (INT4 i = 0; i < start_point.size(); i++) {
#ifdef DEBUG
printf("OLD start_point[%d]: %.16e\n",start_point[i].value());
#endif
start_point[i] = Randomize(start_point[i].value(),upper_bound[i].value(),lower_bound[i].value(),LTMADS());
#ifdef DEBUG
printf("NEW start_point[%d]: %.16e\n",start_point[i].value());
#endif
}
}
}
PrintNomadPoint(lower_bound,"lower_bound");
PrintNomadPoint(start_point,"start_point");
PrintNomadPoint(upper_bound,"upper_bound");
Cell(lower_bound,upper_bound,start_point);
p.reset_X0();
p.set_X0 ( start_point ); // starting point
p.set_LOWER_BOUND ( lower_bound );
p.set_UPPER_BOUND ( upper_bound );
p.check();
// custom evaluator creation:
if (fstatsearchtype == SCSEARCH) {
SCFStatNomadEvaluator scev( p );
NOMAD::Mads mads ( p , &scev );
NOMAD::Evaluator_Control *ec;
ec = &mads.get_evaluator_control();
if (ClearCache()) {
NOMAD::Cache *cache = &ec->get_cache();
cache->clear();
}
NOMAD::Evaluator *pev;
pev = ec->get_evaluator();
FirstPoint(TRUE);
SkipSearch(false);
trials_init();
mads.run();
CacheIn(&ec->get_cache());
mads.reset(true,true);
// algorithm creation and execution:
countstep++;
const NOMAD::Eval_Point *bfp = mads.get_best_feasible();
ExPoint tExPoint;
if (bfp != NULL) {
SetPointFromSolution(&tExPoint,bfp);
}
else {
SetPoint(&tExPoint);
tExPoint.stat = 0;
}
tExPoint.stat = fabs(bfp->get_f().value());
epv.push_back(tExPoint);
trials_push();
}
else {
FCFStatNomadEvaluator fcev( p );
NOMAD::Mads mads ( p , &fcev );
NOMAD::Evaluator_Control *ec;
ec = &mads.get_evaluator_control();
if (ClearCache()) {
NOMAD::Cache *cache = &ec->get_cache();
cache->clear();
}
NOMAD::Evaluator *pev;
pev = ec->get_evaluator();
FirstPoint(true);
SkipSearch(false);
trials_init();
mads.run();
CacheIn(&ec->get_cache());
mads.reset(true,true);
// algorithm creation and execution:
countstep++;
const NOMAD::Eval_Point *bfp = mads.get_best_feasible();
ExPoint tExPoint;
if (bfp != NULL) {
SetPointFromSolution(&tExPoint,bfp);
}
else {
SetPoint(&tExPoint);
tExPoint.stat = 0;
}
epv.push_back(tExPoint);
trials_push();
if (BreakASAP()) {
if (((maxAccept2F() - tExPoint.stat) > 0) && ((tExPoint.stat - minAccept2F()) > 0)) {
LogPrintf(LOG_NORMAL,"Found point loud %f within the expectation [%f,%f], break the search!\n",tExPoint.stat,minAccept2F(),maxAccept2F());
var_ReSearch = FALSE;
var_FinalReSearch = FALSE;
dosearch = FALSE;
break;
}
else {
LogPrintf(LOG_NORMAL,"The found 2F %f is outside the expectation [%f,%f]!\n",tExPoint.stat,minAccept2F(),maxAccept2F());
}
}
}
}
} // end loop stcounter
}
INT4 FStatNomad::SaveFCSignalDistance(const char* fname, REAL8 number, const char* comment) {
LogPrintf(LOG_NORMAL,"Save the coherent distance %e from the signal to %s\n",number,fname);
FILE *fpOut;
if ((fpOut = fopen(fname, "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", fname);
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
fprintf(fpOut,"%%%% %s\n",comment);
fprintf(fpOut,"%.16e\n",number);
fclose(fpOut);
return 0;
}
INT4 FStatNomad::SaveSCSignalDistance(const char* fname, REAL8 numbera, REAL8 numberb, const char* comment) {
LogPrintf(LOG_NORMAL,"Save the coherent distance before %e and after %e from the signal to %s\n",numbera,numberb,fname);
FILE *fpOut;
if ((fpOut = fopen(fname, "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", fname);
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
fprintf(fpOut,"%%%% %s\n",comment);
fprintf(fpOut,"%.16e %.16e\n",numbera,numberb);
fclose(fpOut);
return 0;
}
INT4 FStatNomad::SaveMetric(const char* fname, gsl_matrix* g) {
LogPrintf(LOG_NORMAL,"Save the metric (%d,%d) to %s\n",g->size1,g->size2,fname);
FILE *fpOut;
if ((fpOut = fopen(fname, "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", fname);
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
for (INT4 i = 0; i < g->size1; i++) {
for (INT4 j = 0; j < g->size2; j++) {
fprintf(fpOut,"%8.16e\t",gsl_matrix_get(g,i,j));
}
fprintf(fpOut,"\n");
}
fclose(fpOut);
return 0;
}
void FStatNomad::DoDisplayMetric(gsl_matrix* m) {
for (INT4 i = 0; i < m->size1; i++) {
for (INT4 j = 0; j < m->size2; j++) {
printf("%8.16e\t",gsl_matrix_get(m,i,j));
}
printf("\n");
}
}
INT4 DumpSegmentsToFile(const char *fname, LALSegList *segments) {
FILE *fpOut;
if ((fpOut = fopen(fname, "wb")) == NULL) {
fprintf(stderr, "Unable to open file %s for writing\n", fname);
/*exit*/
return(FSTATFCNOMAD_EFILE);
}
for (INT4 k = 0; k < segments->length; k++) {
fprintf(fpOut,"%d %d %f %d\n",segments->segs[k].start.gpsSeconds,segments->segs[k].end.gpsSeconds,(segments->segs[k].end.gpsSeconds-segments->segs[k].start.gpsSeconds)/3600.,segments->segs[k].id);
}
fclose(fpOut);
return 0;
}
void FStatNomad::CacheIn(NOMAD::Cache *cache) {
LogPrintf(LOG_NORMAL,"CacheIn %d points.\n",cache->size());
LogPrintf(LOG_NORMAL,"Cache size before CacheIn: %d\n",eCache.size());
INT4 n;
for (INT4 i = 0; i < cache->size(); i++) {
const NOMAD::Eval_Point *x;
n = 0;
if (i == 0) {
x = cache->begin();
}
else {
x = cache->next();
}
if (x != NULL && !x->check_nan() && x->get_f().is_defined()) {
ExPoint p;
if (FixedDirection.at(PID_ALPHA)) {
p.alpha = FixedPoint.alpha;
}
else {
p.alpha = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_DELTA)) {
p.delta = FixedPoint.delta;
}
else {
p.delta = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_FREQ)) {
p.freq = FixedPoint.freq;
}
else {
p.freq = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_F1DOT)) {
p.f1dot = FixedPoint.f1dot;
}
else {
p.f1dot = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_F2DOT)) {
p.f2dot = FixedPoint.f2dot;
}
else {
p.f2dot = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_F3DOT)) {
p.f3dot = FixedPoint.f3dot;
}
else {
p.f3dot = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (BinarySearch()) {
if (FixedDirection.at(PID_ASINI)) {
p.asini = FixedPoint.asini;
}
else {
p.asini = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_ARGP)) {
p.argp = FixedPoint.argp;
}
else {
p.argp = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_ECC)) {
p.ecc = FixedPoint.ecc;
}
else {
p.ecc = x->operator[](n).value() / directions.at(n).scale;
n++;
}
if (FixedDirection.at(PID_TPSSB)) {
p.TpSSB = FixedPoint.TpSSB;
}
else {
XLALGPSSetREAL8(&p.TpSSB,x->operator[](n).value() / directions.at(n).scale);
// XLALGPSSetREAL8(&p.TpSSB,x[n].value());
n++;
}
if (FixedDirection.at(PID_PERIOD)) {
p.period = FixedPoint.period;
}
else {
p.period = x->operator[](n).value() / directions.at(n).scale;
n++;
}
}
p.stat = -x->get_f().value();
eCache.push_back(p);
}
}
LogPrintf(LOG_NORMAL,"Cache size after CacheIn: %d\n",eCache.size());
}
FStatNomad::FStatNomad() {
catalog = NULL;
freq_index = NULL;
for (INT4 i = 0; i < 6; i++) {
twoFi[i] = 0;
}
};
FStatNomad::~FStatNomad() {};
/**
* XLAL function to get a list of detector IDs from multi-segment multiSFT vectors
* returns all unique detector IDs for cases with some detectors switching on and off
*/
LALStringVector *
XLALGetDetectorIDs(
LALStringVector *IFOList, ///< IFO string vector for returning
const SFTCatalog *SFTcatalog ///< SFT catalog
)
{
/* check input parameters and report errors */
XLAL_CHECK_NULL( SFTcatalog != NULL, XLAL_EFAULT );
for (UINT4 X = 0; X < SFTcatalog->length; X++)
{
/* check if this segment k, IFO X contains a new detector */
const char *thisIFO = SFTcatalog->data[X].header.name;
if ( XLALFindStringInVector ( thisIFO, IFOList ) >= 0 )
{ // if already in list, do nothing
continue;
}
else
{ // otherwise, append to IFOList
if ( (IFOList = XLALAppendString2Vector ( IFOList, thisIFO )) == NULL )
XLAL_ERROR_NULL ( XLAL_EFUNC );
}
} /* for X < number of detectors */
/* sort final list by detector-name */
XLALSortStringVector ( IFOList );
if ( xlalErrno != 0 ) {
XLALPrintError ("\nError in function %s, line %d : Failed call to XLALSortStringVector().\n\n", __func__, __LINE__);
XLAL_ERROR_NULL ( XLAL_EFUNC );
}
return IFOList;
} /* XLALGetDetectorIDs() */