ComputeFstatBenchmark.c 13.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
/*
*  Copyright (C) 2015 Reinhard Prix
*
*  This program is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as published by
*  the Free Software Foundation; either version 2 of the License, or
*  (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with with program; see the file COPYING. If not, write to the
*  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
*  MA  02111-1307  USA
*/

#include <lal/XLALError.h>
#include <lal/LALBarycenter.h>
#include <lal/LALInitBarycenter.h>
#include <lal/LogPrintf.h>
#include <lal/CWMakeFakeData.h>
#include <lal/LALConstants.h>
#include <lal/ExtrapolatePulsarSpins.h>
#include <lal/ComputeFstat.h>
#include <lal/DetectorStates.h>
#include <lal/LFTandTSutils.h>
#include <lal/LALString.h>
#include <lal/UserInput.h>

// benchmark ComputeFstat() functions for performance and memory usage
34
REAL8 XLALGetCurrentHeapUsageMB ( void );
35 36 37 38 39 40

typedef struct
{
  CHAR *FstatMethod;		//!< select which method/algorithm to use to compute the F-statistic
  REAL8 Freq;
  REAL8 f1dot;
41 42
  REAL8Vector *FreqResolution;
  INT4Vector *numFreqBins;
43 44
  REAL8 Tseg;
  INT4 numSegments;
45
  LALStringVector *IFOs;
46
  CHAR *outputInfo;
47 48 49 50
  INT4 numTrials;

  // ----- developer options
  REAL8 Tsft;
51
  BOOLEAN reuseInput;   // only useful for checking workspace management
52 53 54 55 56 57 58 59 60 61 62 63 64
} UserInput_t;

// ---------- main ----------
int
main ( int argc, char *argv[] )
{
  // ---------- handle user input ----------
  UserInput_t XLAL_INIT_DECL(uvar_s);
  UserInput_t *uvar = &uvar_s;

  uvar->FstatMethod = XLALStringDuplicate("ResampBest");
  uvar->Freq = 100;
  uvar->f1dot = -3e-9;
65
  uvar->FreqResolution = XLALCreateREAL8Vector ( 2 );
66 67
  uvar->FreqResolution->data[0] = 0.1;
  uvar->FreqResolution->data[1] = 1;
68 69 70
  uvar->numFreqBins = XLALCreateINT4Vector ( 2 );
  uvar->numFreqBins->data[0] = 1000;
  uvar->numFreqBins->data[1] = 100000;
71 72
  uvar->Tseg = 60 * 3600;
  uvar->numSegments = 90;
73 74
  uvar->numTrials = 1;

75
  uvar->Tsft = 1800;
76
  uvar->reuseInput = 1;
77 78 79

  XLAL_CHECK ( (uvar->IFOs = XLALCreateStringVector ( "H1", NULL )) != NULL, XLAL_EFUNC );
  uvar->outputInfo = NULL;
80

81
  XLAL_CHECK ( XLALRegisterUvarMember ( FstatMethod,    STRING,         0, OPTIONAL,  "F-statistic method to use. Available methods: %s", XLALFstatMethodHelpString() ) == XLAL_SUCCESS, XLAL_EFUNC );
82 83
  XLAL_CHECK ( XLALRegisterUvarMember ( Freq,           REAL8,          0, OPTIONAL,  "Search frequency in Hz" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( f1dot,          REAL8,          0, OPTIONAL,  "Search spindown f1dot in Hz/s" ) == XLAL_SUCCESS, XLAL_EFUNC );
84
  XLAL_CHECK ( XLALRegisterUvarMember ( FreqResolution, REAL8Vector,    0, OPTIONAL,  "Frequency resolution 'R' in natural units 1/Tseg such that: dFreq = R/Tseg) [range]" ) == XLAL_SUCCESS, XLAL_EFUNC );
85
  XLAL_CHECK ( XLALRegisterUvarMember ( Tseg,           REAL8,          0, OPTIONAL,  "Coherent segment length" ) == XLAL_SUCCESS, XLAL_EFUNC );
86 87 88
  XLAL_CHECK ( XLALRegisterUvarMember ( numSegments,    INT4,           0, OPTIONAL,  "Number of semi-coherent segments" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( numFreqBins,    INT4Vector,     0, OPTIONAL,  "Number of frequency bins to search [range]" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( IFOs,    	STRINGVector,   0, OPTIONAL,  "IFOs to use [list]" ) == XLAL_SUCCESS, XLAL_EFUNC );
89 90
  XLAL_CHECK ( XLALRegisterUvarMember ( numTrials,    	INT4,           0, OPTIONAL,  "Number of repeated trials to run (with potentially randomized parameters)" ) == XLAL_SUCCESS, XLAL_EFUNC );

91
  XLAL_CHECK ( XLALRegisterUvarMember ( outputInfo,     STRING,         0, OPTIONAL,  "Append Resampling internal info into this file") == XLAL_SUCCESS, XLAL_EFUNC );
92

93
  XLAL_CHECK ( XLALRegisterUvarMember ( Tsft,           REAL8,          0, DEVELOPER, "SFT length" ) == XLAL_SUCCESS, XLAL_EFUNC );
94
  XLAL_CHECK ( XLALRegisterUvarMember ( reuseInput,     BOOLEAN,        0, DEVELOPER, "Re-use FstatInput from previous setups (only useful for checking workspace management)" ) == XLAL_SUCCESS, XLAL_EFUNC );
95

96 97 98 99
  BOOLEAN should_exit = 0;
  XLAL_CHECK( XLALUserVarReadAllInput( &should_exit, argc, argv ) == XLAL_SUCCESS, XLAL_EFUNC );
  if ( should_exit ) {
    return EXIT_FAILURE;
100
  }
101
  // check user input
102
  XLAL_CHECK ( uvar->numSegments >= 1, XLAL_EINVAL );
103 104 105 106 107 108 109 110 111 112 113
  XLAL_CHECK ( (uvar->FreqResolution->length == 1) || (uvar->FreqResolution->length == 2), XLAL_EINVAL );
  XLAL_CHECK ( uvar->FreqResolution->data[0] > 0, XLAL_EINVAL );
  REAL8 FreqResolutionMin, FreqResolutionMax;
  FreqResolutionMin = FreqResolutionMax = uvar->FreqResolution->data[0];
  if ( uvar->FreqResolution->length == 2 )
    {
      XLAL_CHECK ( uvar->FreqResolution->data[1] > 0, XLAL_EINVAL );
      XLAL_CHECK ( uvar->FreqResolution->data[1] > uvar->FreqResolution->data[0], XLAL_EINVAL );
      FreqResolutionMax = uvar->FreqResolution->data[1];
    }

114 115 116
  XLAL_CHECK ( uvar->Freq > 0, XLAL_EINVAL );
  XLAL_CHECK ( uvar->Tseg > uvar->Tsft, XLAL_EINVAL );
  XLAL_CHECK ( uvar->Tsft > 1, XLAL_EINVAL );
117 118 119 120 121 122 123 124 125 126 127 128
  XLAL_CHECK ( (uvar->numFreqBins->length == 1) || (uvar->numFreqBins->length == 2), XLAL_EINVAL );
  XLAL_CHECK ( uvar->numFreqBins->data[0] > 0, XLAL_EINVAL );
  UINT4 numFreqBinsMax, numFreqBinsMin;
  numFreqBinsMin = numFreqBinsMax = uvar->numFreqBins->data[0];
  if ( uvar->numFreqBins->length == 2 )
    {
      XLAL_CHECK ( uvar->numFreqBins->data[1] > 0, XLAL_EINVAL );
      XLAL_CHECK ( uvar->numFreqBins->data[1] > uvar->numFreqBins->data[0], XLAL_EINVAL );
      numFreqBinsMax = uvar->numFreqBins->data[1];
    }

  XLAL_CHECK ( uvar->numTrials >= 1, XLAL_EINVAL );
129
  // ---------- end: handle user input ----------
130 131 132 133 134

  // common setup over repeated trials
  FstatMethodType FstatMethod;
  XLAL_CHECK ( XLALParseFstatMethodString ( &FstatMethod, uvar->FstatMethod ) == XLAL_SUCCESS, XLAL_EFUNC );

135 136
  EphemerisData *ephem;
  XLAL_CHECK ( (ephem = XLALInitBarycenter ( TEST_DATA_DIR "earth00-19-DE405.dat.gz", TEST_DATA_DIR "sun00-19-DE405.dat.gz" )) != NULL, XLAL_EFUNC );
137
  REAL8 memBase = XLALGetCurrentHeapUsageMB();
138

139
  UINT4 numDetectors = uvar->IFOs->length;
140 141 142 143 144 145 146 147 148 149 150 151 152
  // ----- setup injection and data parameters
  LIGOTimeGPS startTime = {711595934, 0};
  LIGOTimeGPS startTime_l = startTime;
  LIGOTimeGPS endTime_l;
  SFTCatalog **catalogs;
  XLAL_CHECK ( (catalogs = XLALCalloc ( uvar->numSegments, sizeof( catalogs[0] ))) != NULL, XLAL_ENOMEM );

  for ( INT4 l = 0; l < uvar->numSegments; l ++ )
    {
      endTime_l = startTime_l;
      XLALGPSAdd( &endTime_l, uvar->Tseg );
      MultiLIGOTimeGPSVector *multiTimestamps;
      XLAL_CHECK ( (multiTimestamps = XLALMakeMultiTimestamps ( startTime_l, uvar->Tseg, uvar->Tsft, 0, numDetectors )) != NULL, XLAL_EFUNC );
153
      XLAL_CHECK ( (catalogs[l] = XLALMultiAddToFakeSFTCatalog ( NULL, uvar->IFOs, multiTimestamps )) != NULL, XLAL_EFUNC );
154 155 156 157
      XLALDestroyMultiTimestamps ( multiTimestamps );
      startTime_l = endTime_l;
    } // for l < numSegments
  LIGOTimeGPS endTime = endTime_l;
158

159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
  PulsarSpinRange XLAL_INIT_DECL(spinRange);
  LIGOTimeGPS refTime = { startTime.gpsSeconds - 2.3 * uvar->Tseg, 0 };
  spinRange.refTime = refTime;
  spinRange.fkdot[0] = uvar->Freq;
  spinRange.fkdot[1] = uvar->f1dot;
  spinRange.fkdotBand[1] = 0;
  REAL8 asini = 0, Period = 0, ecc = 0;
  REAL8 minCoverFreq, maxCoverFreq;

  PulsarDopplerParams XLAL_INIT_DECL(Doppler);
  Doppler.refTime = refTime;
  Doppler.Alpha = 0.5;
  Doppler.Delta = 0.5;
  memcpy ( &Doppler.fkdot, &spinRange.fkdot, sizeof(Doppler.fkdot) );;
  Doppler.period = Period;
  Doppler.ecc = ecc;
  Doppler.asini = asini;

177 178
  // ----- setup optional Fstat arguments
  FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
179 180 181 182 183
  MultiNoiseFloor XLAL_INIT_DECL(injectSqrtSX);
  injectSqrtSX.length = numDetectors;
  for ( UINT4 X=0; X < numDetectors; X ++ ) {
    injectSqrtSX.sqrtSn[X] = 1;
  }
184 185
  optionalArgs.injectSqrtSX = &injectSqrtSX;
  optionalArgs.FstatMethod = FstatMethod;
186
  optionalArgs.collectTiming = 1;
187

188
  FILE *timingLogFILE = NULL;
189
  BOOLEAN printHeader = 0;
190
  if ( uvar->outputInfo != NULL )
191
    {
192 193 194 195 196 197
      FILE *tmp;
      if ( (tmp = fopen ( uvar->outputInfo, "r" )) == NULL ) {
        printHeader = 1;
      } else {
        fclose (tmp );
      }
198
      XLAL_CHECK ( (timingLogFILE = fopen (uvar->outputInfo, "ab")) != NULL, XLAL_ESYS, "Failed to open '%s' for appending\n", uvar->outputInfo );
199
    }
200
  FstatInputVector *inputs;
201 202
  FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_2F_PER_DET);
  FstatResults *results = NULL;
203 204
  REAL8 tauF1NoBuf = 0;
  REAL8 tauF1Buf = 0;
205 206
  // ---------- main loop over repeated trials ----------
  for ( INT4 i = 0; i < uvar->numTrials; i ++ )
207
    {
208 209 210
      // randomize numFreqBins
      UINT4 numFreqBins_i = numFreqBinsMin + (UINT4)round ( 1.0 * (numFreqBinsMax - numFreqBinsMin) * rand() / RAND_MAX );
      // randomize FreqResolution
211
      REAL8 FreqResolution_i = FreqResolutionMin + ( FreqResolutionMax - FreqResolutionMin ) * rand() / RAND_MAX;
212 213 214

      XLAL_CHECK ( (inputs = XLALCreateFstatInputVector ( uvar->numSegments )) != NULL, XLAL_EFUNC );

215
      REAL8 dFreq = FreqResolution_i / uvar->Tseg;
216 217

      REAL8 FreqBand = numFreqBins_i * dFreq;
218
      fprintf ( stderr, "trial %d/%d: Tseg = %.1f d, numSegments = %d, Freq = %.1f Hz, f1dot = %.1e Hz/s, FreqResolution R = %.2f, numFreqBins = %d [dFreq = %.2e Hz, FreqBand = %.2e Hz]\n",
219 220 221 222 223 224 225 226
                i+1, uvar->numTrials, uvar->Tseg / 86400.0, uvar->numSegments, uvar->Freq, uvar->f1dot, FreqResolution_i, numFreqBins_i, dFreq, FreqBand );

      spinRange.fkdotBand[0] = FreqBand;
      XLAL_CHECK ( XLALCWSignalCoveringBand ( &minCoverFreq, &maxCoverFreq, &startTime, &endTime, &spinRange, asini, Period, ecc ) == XLAL_SUCCESS, XLAL_EFUNC );

      // create per-segment input structs
      for ( INT4 l = 0; l < uvar->numSegments; l ++ )
        {
227
          if ( uvar->reuseInput && l > 0 ) {
228
            optionalArgs.prevInput = inputs->data[0];
229 230
          } else {
            optionalArgs.prevInput = NULL;
231
          }
232
          XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, dFreq, ephem, &optionalArgs )) != NULL, XLAL_EFUNC );
233 234 235
        }

      // ----- compute Fstatistics over segments
236 237
      REAL8 tauF1NoBuf_i = 0;
      REAL8 tauF1Buf_i = 0;
238 239 240 241
      for ( INT4 l = 0; l < uvar->numSegments; l ++ )
        {
          XLAL_CHECK ( XLALComputeFstat ( &results, inputs->data[l], &Doppler, numFreqBins_i, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );

242 243
          REAL8 Fstat_tauF1Buf, Fstat_tauF1NoBuf;
          XLAL_CHECK ( XLALGetFstatTiming ( inputs->data[l], &Fstat_tauF1Buf, &Fstat_tauF1NoBuf ) == XLAL_SUCCESS, XLAL_EFUNC );
244 245
          tauF1NoBuf_i += Fstat_tauF1NoBuf;
          tauF1Buf_i   += Fstat_tauF1Buf;
246 247
          // ----- output timing details to file if requested
          if ( timingLogFILE != NULL ) {
248
            XLAL_CHECK ( AppendFstatTimingInfo2File ( inputs->data[l], timingLogFILE, printHeader ) == XLAL_SUCCESS, XLAL_EFUNC );
249 250
          }

251
        } // for l < numSegments
252 253 254

      tauF1NoBuf_i /= uvar->numSegments;
      tauF1Buf_i   /= uvar->numSegments;
255 256
      tauF1Buf     += tauF1Buf_i;
      tauF1NoBuf   += tauF1NoBuf_i;
257

258 259
      REAL8 memEnd = XLALGetCurrentHeapUsageMB();
      const char *FmethodName = XLALGetFstatInputMethodName ( inputs->data[0] );
260

261 262
      fprintf (stderr, "%-15s: tauF1Buf = %8.2g s, tauF1NoBuf = %8.2g s, memResamp = %6.1f MB\n",
               FmethodName, tauF1Buf_i, tauF1NoBuf_i, memEnd - memBase );
263 264
      XLALDestroyFstatInputVector ( inputs );
    } // for i < numTrials
265

266 267 268 269 270
  tauF1Buf   /= uvar->numTrials;
  tauF1NoBuf /= uvar->numTrials;

  fprintf (stderr, "\nAveraged timings: <tauF1Buf> = %.2g s, <tauF1NoBuf> = %.2g s\n", tauF1Buf, tauF1NoBuf );

271
  // ----- free memory ----------
272 273
  if ( timingLogFILE != NULL ) {
    fclose ( timingLogFILE );
274
  }
275

276 277 278 279 280
  for ( INT4 l = 0; l < uvar->numSegments; l ++ )
    {
      XLALDestroySFTCatalog ( catalogs[l] );
    }
  XLALFree ( catalogs );
281
  XLALDestroyFstatResults ( results );
282 283 284 285 286 287 288 289
  XLALDestroyUserVars();
  XLALDestroyEphemerisData ( ephem );

  LALCheckMemoryLeaks();

  return XLAL_SUCCESS;

} // main()
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327


// --------------------------------------------------------------------------------
// code to read current process RSS memory usage from /proc, taken from
// https://stackoverflow.com/questions/63166/how-to-determine-cpu-and-memory-consumption-from-inside-a-process
static int parseLine(char* line);

REAL8
XLALGetCurrentHeapUsageMB ( void )
{
  FILE* file = fopen("/proc/self/status", "r");
  int result = -1;
  char line[128];
  if ( file == NULL ) {
    return result;
  }
  while ( fgets ( line, sizeof(line), file) != NULL )
    {
      if (strncmp(line, "VmRSS:", 6) == 0){
        result = parseLine(line);
        break;
      }
    }
  fclose(file);
  return (result / 1024.0);
} // XLALGetCurrentHeapUsageMB()

static int parseLine(char* line)
{
  // This assumes that a digit will be found and the line ends in " Kb".
  int i = strlen(line);
  const char* p = line;
  while (*p <'0' || *p > '9') p++;
  line[i-3] = '\0';
  i = atoi(p);
  return i;
}
// --------------------------------------------------------------------------------