ComputeFstatBenchmark.c 12.7 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 34 35 36 37 38 39 40
/*
*  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

typedef struct
{
  BOOLEAN help;			//!< output help-string */
  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 51
  INT4 numTrials;

  // ----- developer options
  REAL8 Tsft;
  BOOLEAN runBuffered;	// only useful for double-checking and Demod timing
52 53
} UserInput_t;

54 55
int XLALAppendResampInfo2File ( FILE *fp, const FstatInput *input );

56 57 58 59 60 61 62 63 64 65 66
// ---------- 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;
67 68 69 70 71 72
  uvar->FreqResolution = XLALCreateREAL8Vector ( 2 );
  uvar->FreqResolution->data[0] = 1;
  uvar->FreqResolution->data[1] = 10;
  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 79 80 81
  uvar->runBuffered = 0;

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

  XLAL_CHECK ( XLALRegisterUvarMember ( help,           BOOLEAN,        'h', HELP,    "Print help message" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( FstatMethod,    STRING,         0, OPTIONAL,  XLALFstatMethodHelpString() ) == XLAL_SUCCESS, XLAL_EFUNC );
  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 );
87
  XLAL_CHECK ( XLALRegisterUvarMember ( FreqResolution, REAL8Vector,    0, OPTIONAL,  "Range of frequency resolution factor 'r' (st dFreq = 1/(r*T)) [2-number range input]" ) == XLAL_SUCCESS, XLAL_EFUNC );
88 89 90 91 92 93 94
  XLAL_CHECK ( XLALRegisterUvarMember ( Tseg,           REAL8,          0, OPTIONAL,  "Coherent segment length" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( numSegments,    INT4,           0, OPTIONAL,  "Number of semi-coherent segment" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( numFreqBins,    INT4Vector,     0, OPTIONAL,  "Range of number of frequency bins to search [2-number range input]" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( IFOs,    	STRINGVector,   0, OPTIONAL,  "IFOs to use" ) == XLAL_SUCCESS, XLAL_EFUNC );
  XLAL_CHECK ( XLALRegisterUvarMember ( numTrials,    	INT4,           0, OPTIONAL,  "Number of repeated trials to run (with potentially randomized parameters)" ) == XLAL_SUCCESS, XLAL_EFUNC );

  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 ( runBuffered,    BOOLEAN,        0, DEVELOPER, "Explicitly time buffered Fstat call (only useful for double-checking and Demod timing)" ) == XLAL_SUCCESS, XLAL_EFUNC );
98 99 100 101 102

  XLAL_CHECK ( XLALUserVarReadAllInput(argc, argv) == XLAL_SUCCESS, XLAL_EFUNC );
  if (uvar->help) {	// if help was requested, we're done here
    return XLAL_SUCCESS;
  }
103
  // check user input
104
  XLAL_CHECK ( uvar->numSegments >= 1, XLAL_EINVAL );
105 106 107 108 109 110 111 112 113 114 115
  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];
    }

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

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

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

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

  FILE *fpInfo = NULL;
  if ( (uvar->outputInfo != NULL) && (FstatMethod == FMETHOD_RESAMP_GENERIC) )
    {
      XLAL_CHECK ( (fpInfo = fopen (uvar->outputInfo, "ab")) != NULL, XLAL_ESYS, "Failed to open '%s' for appending\n", uvar->outputInfo );
      XLALAppendResampInfo2File ( fpInfo, NULL ); // create header comment line
    }
168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186

  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;

187 188
  // ----- setup optional Fstat arguments
  FstatOptionalArgs optionalArgs = FstatOptionalArgsDefaults;
189 190 191 192 193
  MultiNoiseFloor XLAL_INIT_DECL(injectSqrtSX);
  injectSqrtSX.length = numDetectors;
  for ( UINT4 X=0; X < numDetectors; X ++ ) {
    injectSqrtSX.sqrtSn[X] = 1;
  }
194 195
  optionalArgs.injectSqrtSX = &injectSqrtSX;
  optionalArgs.FstatMethod = FstatMethod;
196

197
  FstatWorkspace *sharedWorkspace = NULL;
198
  FstatInputVector *inputs;
199 200 201 202
  FstatQuantities whatToCompute = (FSTATQ_2F | FSTATQ_2F_PER_DET);
  FstatResults *results = NULL;
  // ---------- main loop over repeated trials ----------
  for ( INT4 i = 0; i < uvar->numTrials; i ++ )
203
    {
204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
      // randomize numFreqBins
      UINT4 numFreqBins_i = numFreqBinsMin + (UINT4)round ( 1.0 * (numFreqBinsMax - numFreqBinsMin) * rand() / RAND_MAX );
      // randomize FreqResolution
      REAL8 FreqResolution_i = FreqResolutionMin + 1.0 * ( FreqResolutionMax - FreqResolutionMin ) * rand() / RAND_MAX;

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

      REAL8 dFreq = 1.0 / ( FreqResolution_i * uvar->Tseg );

      REAL8 FreqBand = numFreqBins_i * dFreq;
      fprintf ( stderr, "trial %d/%d: Tseg = %.1f d, numSegments = %d, Freq = %.1f Hz, f1dot = %.1e Hz/s, FreqResolution r = %f, numFreqBins = %d [dFreq = %.2e Hz, FreqBand = %.2e Hz]\n",
                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 );

      UINT4 numBinsSFT = ceil ( (maxCoverFreq - minCoverFreq) * uvar->Tsft + 2 * 8 );
      REAL8 memSFTs = uvar->numSegments * numSFTsPerSeg * ( sizeof(SFTtype) + numBinsSFT * sizeof(COMPLEX8)) / 1e6;

      // create per-segment input structs
      for ( INT4 l = 0; l < uvar->numSegments; l ++ )
        {
          XLAL_CHECK ( (inputs->data[l] = XLALCreateFstatInput ( catalogs[l], minCoverFreq, maxCoverFreq, dFreq, ephem, &optionalArgs )) != NULL, XLAL_EFUNC );
          if ( l == 0 ) {
            sharedWorkspace = XLALGetSharedFstatWorkspace ( inputs->data[0] );
          }
          optionalArgs.sharedWorkspace = sharedWorkspace;
        }

      // ----- compute Fstatistics over segments
      REAL8 tauFSumUnbuffered = 0;
      REAL8 tauFSumBuffered = 0;
      for ( INT4 l = 0; l < uvar->numSegments; l ++ )
        {
238
          REAL8 tic = XLALGetTimeOfDay();
239
          XLAL_CHECK ( XLALComputeFstat ( &results, inputs->data[l], &Doppler, numFreqBins_i, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
240
          REAL8 toc = XLALGetTimeOfDay();
241 242 243 244 245 246 247 248 249
          tauFSumUnbuffered += ( toc - tic );

          // ----- output more details if requested [only from first segment]
          if ( (l == 0) && (fpInfo != NULL) ) {
            XLALAppendResampInfo2File ( fpInfo, inputs->data[0] );
          }

          if ( uvar->runBuffered )
            {
250
              tic = XLALGetTimeOfDay();
251
              XLAL_CHECK ( XLALComputeFstat ( &results, inputs->data[l], &Doppler, numFreqBins_i, whatToCompute ) == XLAL_SUCCESS, XLAL_EFUNC );
252
              toc = XLALGetTimeOfDay();
253 254
              tauFSumBuffered += ( toc - tic );
            }
255

256 257 258 259
        } // for l < numSegments
      REAL8 tauF1Unbuffered = tauFSumUnbuffered / ( uvar->numSegments * numFreqBins_i * numDetectors );
      REAL8 tauF1Buffered   = tauFSumBuffered   / ( uvar->numSegments * numFreqBins_i * numDetectors );
      REAL8 memMaxCompute = XLALGetPeakHeapUsageMB() - memBase;
260

261 262
      fprintf (stderr, "%-15s: tauF1Buffered = %.2g s, tauF1Unbuffered = %.2g s, memSFTs = %.1f MB, memMaxCompute = %.1f MB\n",
               XLALGetFstatMethodName ( FstatMethod ), (tauF1Buffered > 0) ? tauF1Buffered : -1, tauF1Unbuffered, memSFTs, memMaxCompute  );
263

264 265 266 267
      XLALDestroyFstatInputVector ( inputs );
      XLALDestroyFstatWorkspace ( sharedWorkspace );
      optionalArgs.sharedWorkspace = NULL;
    } // for i < numTrials
268

269
  // ----- free memory ----------
270 271 272
  if ( fpInfo != NULL ) {
    fclose ( fpInfo );
  }
273

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

  LALCheckMemoryLeaks();

  return XLAL_SUCCESS;

} // main()