Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org on Tuesday 7th July 2020 starting at approximately 10am PDT and lasting for around 15 minutes. There will be a short period of downtime towards the end of the maintenance window. Please direct any comments, questions, or concerns to uwm-help@cgca.uwm.edu.

ComputeFstatBenchmark.c 13.6 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
  REAL8 f2dot;
42 43
  REAL8Vector *FreqResolution;
  INT4Vector *numFreqBins;
44 45
  REAL8 Tseg;
  INT4 numSegments;
46
  LALStringVector *IFOs;
47
  CHAR *outputInfo;
48 49 50 51
  INT4 numTrials;

  // ----- developer options
  REAL8 Tsft;
52
  BOOLEAN reuseInput;   // only useful for checking workspace management
53 54 55 56 57 58 59 60 61 62 63 64 65
} 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;
66
  uvar->f2dot = 0;
67
  uvar->FreqResolution = XLALCreateREAL8Vector ( 2 );
68 69
  uvar->FreqResolution->data[0] = 0.1;
  uvar->FreqResolution->data[1] = 1;
70 71 72
  uvar->numFreqBins = XLALCreateINT4Vector ( 2 );
  uvar->numFreqBins->data[0] = 1000;
  uvar->numFreqBins->data[1] = 100000;
73 74
  uvar->Tseg = 60 * 3600;
  uvar->numSegments = 90;
75 76
  uvar->numTrials = 1;

77
  uvar->Tsft = 1800;
78
  uvar->reuseInput = 1;
79 80 81

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

83
  XLAL_CHECK ( XLALRegisterUvarMember ( FstatMethod,    STRING,         0, OPTIONAL,  "F-statistic method to use. Available methods: %s", XLALFstatMethodHelpString() ) == XLAL_SUCCESS, XLAL_EFUNC );
84
  XLAL_CHECK ( XLALRegisterUvarMember ( Freq,           REAL8,          0, OPTIONAL,  "Search frequency in Hz" ) == XLAL_SUCCESS, XLAL_EFUNC );
85 86
  XLAL_CHECK ( XLALRegisterUvarMember ( f1dot,          REAL8,          0, OPTIONAL,  "Search 1st spindown f1dot in Hz/s" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( f2dot,          REAL8,          0, OPTIONAL,  "Search 2nd spindown f2dot in Hz/s^2" ) == XLAL_SUCCESS, XLAL_EFUNC );
87
  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 );
88
  XLAL_CHECK ( XLALRegisterUvarMember ( Tseg,           REAL8,          0, OPTIONAL,  "Coherent segment length" ) == XLAL_SUCCESS, XLAL_EFUNC );
89 90 91
  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 );
92 93
  XLAL_CHECK ( XLALRegisterUvarMember ( numTrials,    	INT4,           0, OPTIONAL,  "Number of repeated trials to run (with potentially randomized parameters)" ) == XLAL_SUCCESS, XLAL_EFUNC );

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

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

99 100 101 102
  BOOLEAN should_exit = 0;
  XLAL_CHECK( XLALUserVarReadAllInput( &should_exit, argc, argv ) == XLAL_SUCCESS, XLAL_EFUNC );
  if ( should_exit ) {
    return EXIT_FAILURE;
103
  }
104
  // check user input
105
  XLAL_CHECK ( uvar->numSegments >= 1, XLAL_EINVAL );
106 107 108 109 110 111 112 113 114 115 116
  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];
    }

117 118 119
  XLAL_CHECK ( uvar->Freq > 0, XLAL_EINVAL );
  XLAL_CHECK ( uvar->Tseg > uvar->Tsft, XLAL_EINVAL );
  XLAL_CHECK ( uvar->Tsft > 1, XLAL_EINVAL );
120 121 122 123 124 125 126 127 128 129 130 131
  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 );
132
  // ---------- end: handle user input ----------
133 134 135 136 137

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

138 139
  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 );
140
  REAL8 memBase = XLALGetCurrentHeapUsageMB();
141

142
  UINT4 numDetectors = uvar->IFOs->length;
143 144 145 146 147 148 149 150 151 152 153 154 155
  // ----- 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 );
156
      XLAL_CHECK ( (catalogs[l] = XLALMultiAddToFakeSFTCatalog ( NULL, uvar->IFOs, multiTimestamps )) != NULL, XLAL_EFUNC );
157 158 159 160
      XLALDestroyMultiTimestamps ( multiTimestamps );
      startTime_l = endTime_l;
    } // for l < numSegments
  LIGOTimeGPS endTime = endTime_l;
161

162 163
  // ----- setup optional Fstat arguments
  FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
164 165 166 167 168
  MultiNoiseFloor XLAL_INIT_DECL(injectSqrtSX);
  injectSqrtSX.length = numDetectors;
  for ( UINT4 X=0; X < numDetectors; X ++ ) {
    injectSqrtSX.sqrtSn[X] = 1;
  }
169 170
  optionalArgs.injectSqrtSX = &injectSqrtSX;
  optionalArgs.FstatMethod = FstatMethod;
171
  optionalArgs.collectTiming = 1;
172

173 174
  FILE *timingLogFILE = NULL;
  if ( uvar->outputInfo != NULL )
175
    {
176 177 178 179 180
      FILE *tmp;
      if ( (tmp = fopen ( uvar->outputInfo, "r" )) == NULL ) {
      } else {
        fclose (tmp );
      }
181
      XLAL_CHECK ( (timingLogFILE = fopen (uvar->outputInfo, "ab")) != NULL, XLAL_ESYS, "Failed to open '%s' for appending\n", uvar->outputInfo );
182
    }
183
  FstatInputVector *inputs;
184 185
  FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_2F_PER_DET);
  FstatResults *results = NULL;
186 187
  REAL8 tauF1NoBuf = 0;
  REAL8 tauF1Buf = 0;
188 189
  // ---------- main loop over repeated trials ----------
  for ( INT4 i = 0; i < uvar->numTrials; i ++ )
190
    {
191
      PulsarSpinRange XLAL_INIT_DECL(spinRange);
192
      LIGOTimeGPS refTime = { startTime.gpsSeconds + 0.5 * uvar->Tseg, 0 };
193 194 195
      spinRange.refTime = refTime;
      spinRange.fkdot[0] = uvar->Freq;
      spinRange.fkdot[1] = uvar->f1dot;
196
      spinRange.fkdot[2] = uvar->f2dot;
197 198 199 200 201 202 203 204 205 206 207 208 209
      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;

210 211 212
      // randomize numFreqBins
      UINT4 numFreqBins_i = numFreqBinsMin + (UINT4)round ( 1.0 * (numFreqBinsMax - numFreqBinsMin) * rand() / RAND_MAX );
      // randomize FreqResolution
213
      REAL8 FreqResolution_i = FreqResolutionMin + ( FreqResolutionMax - FreqResolutionMin ) * rand() / RAND_MAX;
214 215 216

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

217
      REAL8 dFreq = FreqResolution_i / uvar->Tseg;
218 219

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

      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 ++ )
        {
229
          if ( uvar->reuseInput && l > 0 ) {
230
            optionalArgs.prevInput = inputs->data[0];
231 232
          } else {
            optionalArgs.prevInput = NULL;
233
          }
234
          XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, dFreq, ephem, &optionalArgs )) != NULL, XLAL_EFUNC );
235 236 237
        }

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

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

253
        } // for l < numSegments
254 255 256

      tauF1NoBuf_i /= uvar->numSegments;
      tauF1Buf_i   /= uvar->numSegments;
257 258
      tauF1Buf     += tauF1Buf_i;
      tauF1NoBuf   += tauF1NoBuf_i;
259

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

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

268 269 270 271 272
  tauF1Buf   /= uvar->numTrials;
  tauF1NoBuf /= uvar->numTrials;

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

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

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

  LALCheckMemoryLeaks();

  return XLAL_SUCCESS;

} // main()
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 328 329


// --------------------------------------------------------------------------------
// 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;
}
// --------------------------------------------------------------------------------