ComputeFstat_Demod.c 12.3 KB
Newer Older
1
//
2 3 4
// Copyright (C) 2012--2015 Karl Wette
// Copyright (C) 2005--2007, 2009, 2010, 2012, 2014 Reinhard Prix
// Copyright (C) 2007--2010, 2012 Bernd Machenschalk
5 6
// Copyright (C) 2007 Chris Messenger
// Copyright (C) 2006 John T. Whelan, Badri Krishnan
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
//
// 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
//

24 25 26
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
27

28 29 30
#include "ComputeFstat_internal.h"

#include <lal/Factorial.h>
31
#include <lal/LogPrintf.h>
32
#include <lal/SinCosLUT.h>
33

34 35 36
// ========== Demod internals ==========

// ----- local types ----------
37
typedef struct {
38 39 40 41 42 43
  int (*computefafb_func) (			// XLALComputeFaFb_...() function for the selected Demod hotloop
    COMPLEX8 *, COMPLEX8 *, FstatAtomVector **, const SFTVector *, const PulsarSpins, const SSBtimes *, const AMCoeffs *, const UINT4 Dterms
    );
  UINT4 Dterms;					// Number of terms to keep in Dirichlet kernel
  MultiSFTVector *multiSFTs;			// Input multi-detector SFTs
  REAL8 prevAlpha, prevDelta;			// buffering: previous skyposition computed
44
  LIGOTimeGPS prevRefTime;			// buffering: keep track of previous refTime for SSBtimes buffering
45 46
  MultiSSBtimes *prevMultiSSBtimes;		// buffering: previous multiSSB times, unique to skypos + SFTs
  MultiAMCoeffs *prevMultiAMcoef;		// buffering: previous AM-coeffs, unique to skypos + SFTs
47
} DemodMethodData;
48

49 50 51
// ----- local prototypes ----------

int XLALSetupFstatDemod  ( void **method_data, FstatCommon *common, FstatMethodFuncs* funcs, MultiSFTVector *multiSFTs, const FstatOptionalArgs *optArgs );
Karl Wette's avatar
Karl Wette committed
52

53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
int XLALComputeFaFb_Generic ( COMPLEX8 *Fa, COMPLEX8 *Fb, FstatAtomVector **FstatAtoms, const SFTVector *sfts,
                              const PulsarSpins fkdot, const SSBtimes *tSSB, const AMCoeffs *amcoe, const UINT4 Dterms );

int XLALComputeFaFb_OptC    ( COMPLEX8 *Fa, COMPLEX8 *Fb, FstatAtomVector **FstatAtoms, const SFTVector *sfts,
                              const PulsarSpins fkdot, const SSBtimes *tSSB, const AMCoeffs *amcoe, const UINT4 Dterms );

#ifdef HAVE_ALTIVEC
int XLALComputeFaFb_Altivec ( COMPLEX8 *Fa, COMPLEX8 *Fb, FstatAtomVector **FstatAtoms, const SFTVector *sfts,
                              const PulsarSpins fkdot, const SSBtimes *tSSB, const AMCoeffs *amcoe, const UINT4 Dterms );
#endif

#ifdef HAVE_SSE_COMPILER
int XLALComputeFaFb_SSE     ( COMPLEX8 *Fa, COMPLEX8 *Fb, FstatAtomVector **FstatAtoms, const SFTVector *sfts,
                              const PulsarSpins fkdot, const SSBtimes *tSSB, const AMCoeffs *amcoe, const UINT4 Dterms );
#endif

69
// ----- local function definitions ----------
Karl Wette's avatar
Karl Wette committed
70
static int
71 72 73 74
XLALComputeFstatDemod ( FstatResults* Fstats,
                        const FstatCommon *common,
                        void *method_data
                      )
75
{
76 77 78
  // Check input
  XLAL_CHECK(Fstats != NULL, XLAL_EFAULT);
  XLAL_CHECK(common != NULL, XLAL_EFAULT);
79 80
  XLAL_CHECK(method_data != NULL, XLAL_EFAULT);

81 82 83 84 85
#if COLLECT_TIMING
  // get internal timing info
  REAL8 tic, toc, tauBary, tauTotal;
#endif

86
  DemodMethodData *demod = (DemodMethodData*) method_data;
87 88 89 90

  // Get which F-statistic quantities to compute
  const FstatQuantities whatToCompute = Fstats->whatWasComputed;

91 92
  // handy shortcuts
  BOOLEAN returnAtoms = (whatToCompute & FSTATQ_ATOMS_PER_DET);
93 94
  PulsarDopplerParams thisPoint = Fstats->doppler;
  const REAL8 fStart = thisPoint.fkdot[0];
95
  const MultiSFTVector *multiSFTs = demod->multiSFTs;
96 97
  const MultiNoiseWeights *multiWeights = common->multiNoiseWeights;
  const MultiDetectorStateSeries *multiDetStates = common->multiDetectorStates;
98

99 100 101
  UINT4 numDetectors = multiSFTs->length;
  XLAL_CHECK ( multiDetStates->length == numDetectors, XLAL_EINVAL );
  XLAL_CHECK ( multiWeights==NULL || (multiWeights->length == numDetectors), XLAL_EINVAL );
102

103 104 105
#if COLLECT_TIMING
  tic = XLALGetCPUTime();
#endif
106 107 108
  MultiSSBtimes *multiSSB = NULL;
  MultiAMCoeffs *multiAMcoef = NULL;
  // ----- check if we have buffered SSB+AMcoef for current sky-position
109 110 111 112
  if ( (demod->prevAlpha == thisPoint.Alpha) && (demod->prevDelta == thisPoint.Delta ) &&
       (demod->prevMultiSSBtimes != NULL) && ( XLALGPSDiff(&demod->prevRefTime, &thisPoint.refTime) == 0 ) &&	// have SSB times for same reftime?
       (demod->prevMultiAMcoef != NULL)
       )
113 114 115 116 117 118 119 120 121 122 123 124
    { // if yes ==> reuse
      multiSSB    = demod->prevMultiSSBtimes;
      multiAMcoef = demod->prevMultiAMcoef;
    }
  else
    { // if not, compute SSB + AMcoef for this skyposition
      SkyPosition skypos;
      skypos.system = COORDINATESYSTEM_EQUATORIAL;
      skypos.longitude = thisPoint.Alpha;
      skypos.latitude  = thisPoint.Delta;
      XLAL_CHECK ( (multiSSB = XLALGetMultiSSBtimes ( multiDetStates, skypos, thisPoint.refTime, common->SSBprec )) != NULL, XLAL_EFUNC );
      XLAL_CHECK ( (multiAMcoef = XLALComputeMultiAMCoeffs ( multiDetStates, multiWeights, skypos )) != NULL, XLAL_EFUNC );
125

126 127 128
      // store these for possible later re-use in buffer
      XLALDestroyMultiSSBtimes ( demod->prevMultiSSBtimes );
      demod->prevMultiSSBtimes = multiSSB;
129
      demod->prevRefTime = thisPoint.refTime;
130 131 132 133 134
      XLALDestroyMultiAMCoeffs ( demod->prevMultiAMcoef );
      demod->prevMultiAMcoef = multiAMcoef;
      demod->prevAlpha = thisPoint.Alpha;
      demod->prevDelta = thisPoint.Delta;
    } // if could not reuse previously buffered quantites
135

136 137
  MultiSSBtimes *multiBinary = NULL;
  MultiSSBtimes *multiSSBTotal = NULL;
138
  // handle binary-orbital timing corrections, if applicable
139 140 141 142 143
  if ( thisPoint.asini > 0 )
    {
      // compute binary time corrections to the SSB time delays and SSB time derivitive
      XLAL_CHECK ( XLALAddMultiBinaryTimes ( &multiBinary, multiSSB, &thisPoint ) == XLAL_SUCCESS, XLAL_EFUNC );
      multiSSBTotal = multiBinary;
144
    }
145 146 147
  else
    {
      multiSSBTotal = multiSSB;
148
    }
149 150 151 152
#if COLLECT_TIMING
  toc = XLALGetCPUTime();
  tauBary = (toc - tic);
#endif
153

154
  // ----- compute final Fstatistic-value -----
155 156 157 158 159
  REAL4 Ad = multiAMcoef->Mmunu.Ad;
  REAL4 Bd = multiAMcoef->Mmunu.Bd;
  REAL4 Cd = multiAMcoef->Mmunu.Cd;
  REAL4 Ed = multiAMcoef->Mmunu.Ed;;
  REAL4 Dd_inv = 1.0 / multiAMcoef->Mmunu.Dd;
160

161 162
  // ---------- Compute F-stat for each frequency bin ----------
  for ( UINT4 k = 0; k < Fstats->numFreqBins; k++ )
163 164 165 166
    {
      // Set frequency to search at
      thisPoint.fkdot[0] = fStart + k * Fstats->dFreq;

167 168
      COMPLEX8 Fa = 0;       		// complex amplitude Fa
      COMPLEX8 Fb = 0;                 // complex amplitude Fb
169 170 171 172 173 174 175 176 177 178 179 180 181
      MultiFstatAtomVector *multiFstatAtoms = NULL;	// per-IFO, per-SFT arrays of F-stat 'atoms', ie quantities required to compute F-stat

      // prepare return of 'FstatAtoms' if requested
      if ( returnAtoms )
        {
          XLAL_CHECK ( (multiFstatAtoms = XLALMalloc ( sizeof(*multiFstatAtoms) )) != NULL, XLAL_ENOMEM );
          multiFstatAtoms->length = numDetectors;
          XLAL_CHECK ( (multiFstatAtoms->data = XLALMalloc ( numDetectors * sizeof(*multiFstatAtoms->data) )) != NULL, XLAL_ENOMEM );
        } // if returnAtoms

      // loop over detectors and compute all detector-specific quantities
      for ( UINT4 X=0; X < numDetectors; X ++)
        {
182
          COMPLEX8 FaX, FbX;
183 184 185
          FstatAtomVector *FstatAtoms = NULL;
          FstatAtomVector **FstatAtoms_p = returnAtoms ? (&FstatAtoms) : NULL;

186 187 188
          // call XLALComputeFaFb_...() function for the user-requested hotloop variant
          XLAL_CHECK ( (demod->computefafb_func) ( &FaX, &FbX, FstatAtoms_p, multiSFTs->data[X], thisPoint.fkdot,
                                                   multiSSBTotal->data[X], multiAMcoef->data[X], demod->Dterms ) == XLAL_SUCCESS, XLAL_EFUNC );
189 190 191 192 193 194 195

          if ( returnAtoms ) {
            multiFstatAtoms->data[X] = FstatAtoms;     // copy pointer to IFO-specific Fstat-atoms 'contents'
          }

          XLAL_CHECK ( isfinite(creal(FaX)) && isfinite(cimag(FaX)) && isfinite(creal(FbX)) && isfinite(cimag(FbX)), XLAL_EFPOVRFLW );

196 197
          if ( whatToCompute & FSTATQ_FAFB_PER_DET )
            {
198 199
              Fstats->FaPerDet[X][k] = FaX;
              Fstats->FbPerDet[X][k] = FbX;
200 201
            }

202
          // compute single-IFO F-stats, if requested
203
          if ( whatToCompute & FSTATQ_2F_PER_DET )
204
            {
205 206 207 208 209
              REAL4 AdX = multiAMcoef->data[X]->A;
              REAL4 BdX = multiAMcoef->data[X]->B;
              REAL4 CdX = multiAMcoef->data[X]->C;
              REAL4 EdX = 0;
              REAL4 DdX_inv = 1.0 / multiAMcoef->data[X]->D;
210

211
              // compute final single-IFO F-stat
212
              Fstats->twoFPerDet[X][k] = XLALComputeFstatFromFaFb ( FaX, FbX, AdX, BdX, CdX, EdX, DdX_inv );
213

214
            } // if FSTATQ_2F_PER_DET
215 216 217 218 219 220 221 222 223

          /* Fa = sum_X Fa_X */
          Fa += FaX;

          /* Fb = sum_X Fb_X */
          Fb += FbX;

        } // for  X < numDetectors

224 225
      if ( whatToCompute & FSTATQ_2F )
        {
226
          Fstats->twoF[k] = XLALComputeFstatFromFaFb ( Fa, Fb, Ad, Bd, Cd, Ed, Dd_inv );
227
        }
228

229
      // Return multi-detector Fa & Fb
230 231
      if ( whatToCompute & FSTATQ_FAFB )
        {
232 233
          Fstats->Fa[k] = Fa;
          Fstats->Fb[k] = Fb;
234 235 236
        }

      // Return F-atoms per detector
237 238 239 240 241
      if ( whatToCompute & FSTATQ_ATOMS_PER_DET )
        {
          XLALDestroyMultiFstatAtomVector ( Fstats->multiFatoms[k] );
          Fstats->multiFatoms[k] = multiFstatAtoms;
        }
242

243
    } // for k < Fstats->numFreqBins
244

245 246 247
  // this needs to be free'ed, as it's currently not buffered
  XLALDestroyMultiSSBtimes ( multiBinary );

248
  // Return amplitude modulation coefficients
249
  Fstats->Mmunu = demod->prevMultiAMcoef->Mmunu;
250

251 252 253 254 255 256 257
#if COLLECT_TIMING
  toc = XLALGetCPUTime();
  tauTotal = (toc - tic);
  Fstat_tauF1NoBuf = tauTotal / ( Fstats->numFreqBins * numDetectors );
  Fstat_tauF1Buf   = (tauTotal - tauBary) / ( Fstats->numFreqBins * numDetectors );
#endif

258 259
  return XLAL_SUCCESS;

260
} // XLALComputeFstatDemod()
261 262 263


static void
264
XLALDestroyDemodMethodData ( void* method_data )
265
{
266 267 268

  DemodMethodData *demod = (DemodMethodData*) method_data;

269 270 271 272 273
  XLALDestroyMultiSFTVector ( demod->multiSFTs);
  XLALDestroyMultiSSBtimes  ( demod->prevMultiSSBtimes );
  XLALDestroyMultiAMCoeffs  ( demod->prevMultiAMcoef );
  XLALFree ( demod );

274
} // XLALDestroyDemodMethodData()
275

276
int
277 278 279 280 281 282
XLALSetupFstatDemod ( void **method_data,
                      FstatCommon *common,
                      FstatMethodFuncs* funcs,
                      MultiSFTVector *multiSFTs,
                      const FstatOptionalArgs *optArgs
                    )
283 284
{
  // Check input
285
  XLAL_CHECK ( method_data != NULL, XLAL_EFAULT );
286 287
  XLAL_CHECK ( common != NULL, XLAL_EFAULT );
  XLAL_CHECK ( funcs != NULL, XLAL_EFAULT );
288
  XLAL_CHECK ( multiSFTs != NULL, XLAL_EFAULT );
289
  XLAL_CHECK ( optArgs != NULL, XLAL_EFAULT );
290

291 292 293
  // Allocate method data
  DemodMethodData *demod = *method_data = XLALCalloc( 1, sizeof(*demod) );
  XLAL_CHECK( demod != NULL, XLAL_ENOMEM );
294

295 296
  // Set method function pointers
  funcs->compute_func = XLALComputeFstatDemod;
297 298
  funcs->method_data_destroy_func = XLALDestroyDemodMethodData;
  funcs->workspace_destroy_func = NULL;
299

300 301
  // Save pointer to SFTs
  demod->multiSFTs = multiSFTs;
302

303
  // Save Dterms
304
  demod->Dterms = optArgs->Dterms;
305

306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328
  // Select XLALComputeFaFb_...() function for the user-requested hotloop variant
  switch ( optArgs->FstatMethod ) {
  case  FMETHOD_DEMOD_GENERIC:
    demod->computefafb_func = XLALComputeFaFb_Generic;
    break;
  case FMETHOD_DEMOD_OPTC:
    demod->computefafb_func = XLALComputeFaFb_OptC;
    break;
#ifdef HAVE_ALTIVEC
  case FMETHOD_DEMOD_ALTIVEC:
    demod->computefafb_func = XLALComputeFaFb_Altivec;
    break;
#endif
#ifdef HAVE_SSE_COMPILER
  case FMETHOD_DEMOD_SSE:
    demod->computefafb_func = XLALComputeFaFb_SSE;
    break;
#endif
  default:
    XLAL_ERROR ( XLAL_EINVAL, "Invalid Demod hotloop optArgs->FstatMethod='%d'", optArgs->FstatMethod );
    break;
  }

329
  return XLAL_SUCCESS;
330

331
} // XLALSetupFstatDemod()