LAL: BCVTemplatesFlatMesh.c Source File

00001 /*
00002 *  Copyright (C) 2007 Jolien Creighton, B.S. Sathyaprakash, Thomas Cokelaer
00003 *
00004 *  This program is free software; you can redistribute it and/or modify
00005 *  it under the terms of the GNU General Public License as published by
00006 *  the Free Software Foundation; either version 2 of the License, or
00007 *  (at your option) any later version.
00008 *
00009 *  This program is distributed in the hope that it will be useful,
00010 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012 *  GNU General Public License for more details.
00013 *
00014 *  You should have received a copy of the GNU General Public License
00015 *  along with with program; see the file COPYING. If not, write to the
00016 *  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
00017 *  MA  02111-1307  USA
00018 */
00019 
00020 /********************************* <lalVerbatim file="BCVTemplatesCV">
00021 Author: B.S. Sathyaprakash
00022 $Id: BCVTemplatesFlatMesh.c,v 1.6 2007/06/08 14:41:42 bema Exp $
00023 **************************************************** </lalVerbatim> */
00024 
00025 /********************************************************** <lalLaTeX>
00026 
00027 \subsection{Program \texttt{BCVTemplates.c}}
00028 \label{ss:BCVTemplates.c}
00029 
00030 Creates a template mesh for BCV (or, alternatively, for SPA but
00031 assuing a constant metric) using the mismatch metric.
00032 
00033 \subsubsection*{Usage}
00034 
00035 \subsubsection*{Description}
00036 
00037 \subsubsection*{Exit codes}
00038 ****************************************** </lalLaTeX><lalErrTable> */
00039 #define FLATMESHTESTC_ENORM 0
00040 #define FLATMESHTESTC_ESUB  1
00041 #define FLATMESHTESTC_EARG  2
00042 #define FLATMESHTESTC_EMEM  3
00043 #define FLATMESHTESTC_EDIM  4
00044 #define FLATMESHTESTC_ELEN  5
00045 #define FLATMESHTESTC_EFILE 6
00046 
00047 #define FLATMESHTESTC_MSGENORM "Normal exit"
00048 #define FLATMESHTESTC_MSGESUB  "Subroutine failed"
00049 #define FLATMESHTESTC_MSGEARG  "Error parsing arguments"
00050 #define FLATMESHTESTC_MSGEMEM  "Memory allocation error"
00051 #define FLATMESHTESTC_MSGEDIM  "Inconsistent parameter space dimension"
00052 #define FLATMESHTESTC_MSGELEN  "Too few points specified"
00053 #define FLATMESHTESTC_MSGEFILE "Could not open file"
00054 /******************************************** </lalErrTable><lalLaTeX>
00055 
00056 \subsubsection*{Algorithm}
00057 
00058 \subsubsection*{Uses}
00059 \begin{verbatim}
00060 lalDebugLevel
00061 LALPrintError()                 LALCheckMemoryLeaks()
00062 LALCalloc()                     LALFree()
00063 LALCreateFlatMesh()             LALSReadVectorSequence()
00064 LALSCreateVectorSequence()      LALSDestroyVectorSequence()
00065 LALSCreateVector()              LALSDestroyVector()
00066 \end{verbatim}
00067 
00068 \subsubsection*{Notes}
00069 
00070 \vfill{\footnotesize\input{BCVTemplatesCV}}
00071 
00072 ******************************************************* </lalLaTeX> */
00073 
00074 #include <math.h>
00075 #include <stdlib.h>
00076 #include <lal/LALInspiralBank.h>
00077 #include <lal/LALStdio.h>
00078 #include <lal/FileIO.h>
00079 #include <lal/LALStdlib.h>
00080 #include <lal/AVFactories.h>
00081 #include <lal/SeqFactories.h>
00082 #include <lal/FlatMesh.h>
00083 #include <lal/StreamInput.h>
00084 
00085 NRCSID(FLATMESHTESTC,"$Id: BCVTemplatesFlatMesh.c,v 1.6 2007/06/08 14:41:42 bema Exp $");
00086 
00087 /* Default parameter settings. */
00088 int lalDebugLevel = 0;
00089 #define MISMATCH 0.3
00090 #define DIM 2
00091 
00092 /* Usage format string. */
00093 #define USAGE "Usage: %s [-o outfile] [-d debuglevel] [-m mismatch]\n\
00094                     [eigenvectorfile inversefile rangefile]\n"
00095 
00096 /* Macros for printing errors and testing subroutines. */
00097 #define ERROR( code, msg, statement )                                \
00098 if ( lalDebugLevel & LALERROR )                                      \
00099 {                                                                    \
00100   LALPrintError( "Error[0] %d: program %s, file %s, line %d, %s\n"   \
00101                  "        %s %s\n", (code), *argv, __FILE__,         \
00102                  __LINE__, FLATMESHTESTC, statement ? statement : \
00103                  "", (msg) );                                        \
00104 }                                                                    \
00105 else (void)(0)
00106 
00107 #define INFO( statement )                                            \
00108 if ( lalDebugLevel & LALINFO )                                       \
00109 {                                                                    \
00110   LALPrintError( "Info[0]: program %s, file %s, line %d, %s\n"       \
00111                  "        %s\n", *argv, __FILE__, __LINE__,          \
00112                  FLATMESHTESTC, (statement) );                    \
00113 }                                                                    \
00114 else (void)(0)
00115 
00116 #define SUB( func, statusptr )                                       \
00117 if ( (func), (statusptr)->statusCode )                               \
00118 {                                                                    \
00119   ERROR( FLATMESHTESTC_ESUB, FLATMESHTESTC_MSGESUB,            \
00120          "Function call \"" #func "\" failed:" );                    \
00121   return FLATMESHTESTC_ESUB;                                      \
00122 }                                                                    \
00123 else (void)(0)
00124 
00125 /* A global pointer for debugging. */
00126 #ifndef NDEBUG
00127 char *lalWatch;
00128 #endif
00129 
00130 static void
00131 GetInspiralMoments (
00132                 LALStatus            *status,
00133                 InspiralMomentsEtc   *moments,
00134                 REAL8FrequencySeries *psd,
00135                 InspiralTemplate     *params );
00136 
00137 void 
00138 LALInspiralComputeBCVMetric(
00139    LALStatus            *status,
00140    InspiralMetric       *metric,
00141    REAL8FrequencySeries *shf,
00142    InspiralTemplate     *params
00143 );
00144 
00145 int
00146 main(int argc, char **argv)
00147 {
00148   INT4 arg;
00149   UINT4 dim;                 /* dimension of parameter space */
00150   static LALStatus status;     /* top-level status structure */
00151   REAL4 mismatch = MISMATCH; /* maximum mismatch level */
00152   REAL4VectorSequence *matrix = NULL;    /* tranformation matrix */
00153   REAL4VectorSequence *matrixInv = NULL; /* inverse tranformation */
00154   REAL4VectorSequence *corners = NULL;   /* corners of serach region */
00155   REAL4VectorSequence *mesh = NULL;      /* mesh of parameter values */
00156   FILE *fp;                  /* input/output file pointer */
00157 
00158   static InspiralMetric metric;
00159   static InspiralTemplate params;
00160   UINT4   nlist, numPSDpts=262144;
00161   REAL8FrequencySeries shf;
00162   REAL8 samplingRate;
00163   void (*noisemodel)(LALStatus*,REAL8*,REAL8) = LALLIGOIPsd;
00164   InspiralMomentsEtc moments;
00165 
00166 /* Number of templates is nlist */
00167 
00168   fp = fopen("BCVTemplatesFlatMesh.out", "w");
00169   dim = DIM;
00170   nlist = 0;
00171 
00172   params.OmegaS = 0.;
00173   params.Theta = 0.;
00174   params.ieta=1; 
00175   params.mass1=1.; 
00176   params.mass2=1.; 
00177   params.startTime=0.0; 
00178   params.startPhase=0.0;
00179   params.fLower=40.0; 
00180   params.fCutoff=2000.00;
00181   params.tSampling=4096.0;
00182   params.order=4;
00183   params.approximant=TaylorT3;
00184   params.signalAmplitude=1.0;
00185   params.nStartPad=0;
00186   params.nEndPad=1000;
00187   params.massChoice=m1Andm2;
00188   params.distance = 1.e8 * LAL_PC_SI/LAL_C_SI;
00189   LALInspiralParameterCalc(&status, &params);
00190 
00191   params.psi0 = 132250.;
00192   params.psi3 = -1314.2;
00193   /*
00194   params.alpha = 0.528;
00195   */
00196   params.alpha = 0.L;
00197   params.fFinal = 868.7;
00198   metric.space = Tau0Tau3;
00199 
00200   samplingRate = params.tSampling;
00201   memset( &(shf), 0, sizeof(REAL8FrequencySeries) );
00202   shf.f0 = 0;
00203   LALDCreateVector( &status, &(shf.data), numPSDpts );
00204   shf.deltaF = samplingRate / (2.*(REAL8) shf.data->length + 1.L);
00205   LALNoiseSpectralDensity (&status, shf.data, noisemodel, shf.deltaF );
00206 
00207   /* compute the metric at this point, update bankPars and add the params to the list */
00208           
00209   GetInspiralMoments (&status, &moments, &shf, &params);
00210   LALInspiralComputeMetric(&status, &metric, &params, &moments);
00211   /*
00212   LALInspiralComputeMetricBCV(&status, &metric, &shf, &params);
00213   */
00214 
00215   fprintf(fp, "%e %e %e\n", metric.G00, metric.G01, metric.G11);
00216   fprintf(fp, "%e %e %e\n", metric.g00, metric.g11, metric.theta);
00217   fprintf(fp, "dp0=%e dp1=%e\n", sqrt (mismatch/metric.G00), sqrt (mismatch/metric.G11));
00218   fprintf(fp, "dP0=%e dP1=%e\n", sqrt (mismatch/metric.g00), sqrt (mismatch/metric.g11));
00219   {
00220     CreateVectorSequenceIn in;
00221     in.length = dim;
00222     in.vectorLength = dim;
00223     SUB( LALSCreateVectorSequence( &status, &matrix, &in ), &status );
00224     SUB( LALSCreateVectorSequence( &status, &matrixInv, &in ), &status );
00225     in.length = 2;
00226     SUB( LALSCreateVectorSequence( &status, &corners, &in ), &status );
00227 
00228     corners->data[0] = 0.3;
00229     corners->data[1] = 0.15;
00230     corners->data[2] = 43.; 
00231     corners->data[3] = 1.25;
00232 
00233     /*
00234     corners->data[0] = 1.e5;
00235     corners->data[1] = -1.e3;
00236     corners->data[2] = 2.0e5;
00237     corners->data[3] = 1.e3;
00238     */
00239 
00240   }
00241   {
00242           REAL4 det;
00243           UINT4 i;
00244 
00245     
00246                   
00247           matrix->data[0] = -metric.G01/sqrt(pow(metric.G01,2.) + pow(metric.G00-metric.g00,2.))
00248           /sqrt(metric.g11);
00249           matrix->data[1] = -(metric.G00-metric.g00)/sqrt(pow(metric.G01,2.) + pow(metric.G00-metric.g00,2.))
00250           /sqrt(metric.g00);
00251           matrix->data[2] = -(metric.G11-metric.g11)/sqrt(pow(metric.G01,2.) + pow(metric.G11-metric.g11,2.))
00252           /sqrt(metric.g11);
00253           matrix->data[3] = -metric.G01/sqrt(pow(metric.G01,2.) + pow(metric.G00-metric.g00,2.))
00254           /sqrt(metric.g00);
00255 
00256           det = matrix->data[0]*matrix->data[3] - matrix->data[1]*matrix->data[2];
00257 
00258           matrixInv->data[0] = matrix->data[3]/det;
00259           matrixInv->data[1] = -matrix->data[1]/det;
00260           matrixInv->data[2] = -matrix->data[2]/det;
00261           matrixInv->data[3] = matrix->data[0]/det;
00262 
00263           for (i=0; i<2*dim; i+=2) fprintf(fp, "%e\t%e\n", matrix->data[i], matrix->data[i+1]);
00264           fprintf(fp, "Det=%e\n", det);
00265           for (i=0; i<2*dim; i+=2) fprintf(fp, "%e\t%e\n", matrixInv->data[i], matrixInv->data[i+1]);
00266   }
00267 
00268           
00269   /* Apply mismatch threshold to the transformation matrices. */
00270   {
00271     
00272           UINT4 i;
00273           REAL4 adjust = 2.0*mismatch/sqrt( (REAL4)(dim) );
00274           REAL4 *data;  /* pointer to matrix data */
00275     
00276           i = matrix->length*matrix->vectorLength;
00277           data = matrix->data;
00278           while ( i-- )
00279                   *(data++) *= adjust;
00280     
00281           adjust = 1.0/adjust;
00282           i = matrixInv->length*matrixInv->vectorLength;
00283           data = matrixInv->data;
00284           while ( i-- )
00285                   *(data++) *= adjust;
00286   }
00287 
00288   /* Extend the range boundary to ensure edge coverage. */
00289   {
00290           UINT4 i, j;    /* indecies */
00291           INT2 direct;  /* sign of direction from first corner to second */
00292           REAL4 *data;  /* pointer to matrix data */
00293           REAL4 *width; /* maximum width of a patch in each dimension */
00294     
00295           /* Allocate local memory. */
00296           width = (REAL4 *)LALCalloc( dim, sizeof(REAL4) );
00297           if ( !width ) {
00298                   ERROR( FLATMESHTESTC_EMEM, FLATMESHTESTC_MSGEMEM, 0 );
00299                   return FLATMESHTESTC_EMEM;
00300           }
00301     
00302           /* Determine patch width. */
00303           for ( data = matrix->data, i = 0; i < dim; i++ )
00304                   for ( j = 0; j < dim; j++, data++ )
00305                           width[j] += fabs( *data );
00306     
00307           /* Extend each corner by 0.5*width in the appropriate
00308              direction. */
00309           for ( data = corners->data, i = 0; i < dim; i++, data++ ) {
00310                   direct = ( data[0] < data[dim] ) ? -1 : 1;
00311                   data[0] += 0.5*direct*width[i];
00312                   data[dim] -= 0.5*direct*width[i];
00313           }
00314     
00315           /* Free local memory. */
00316           LALFree( width );
00317   }
00318 
00319   /* Generate mesh using LALFlatMesh() and LALRectIntersect(). */
00320   {
00321     /* Set up parameter structure for LALFlatMesh. */
00322     static FlatMeshParamStruc flatmesh;
00323     flatmesh.matrix = matrix;
00324     flatmesh.matrixInv = matrixInv;
00325     flatmesh.controlPoints = corners;
00326     /*
00327     flatmesh.intersection = LALRectIntersect;
00328     */
00329     flatmesh.intersection = NULL;
00330     SUB( LALSCreateVector( &status, &(flatmesh.xMin), dim ), &status );
00331     SUB( LALSCreateVector( &status, &(flatmesh.xMax), dim ), &status );
00332   
00333     flatmesh.xMin->data[0] = 0.3; 
00334     flatmesh.xMin->data[1] = 0.15;
00335     flatmesh.xMax->data[0] = 43.; 
00336     flatmesh.xMax->data[1] = 1.25;
00337     /*
00338     flatmesh.xMin->data[0] = 1.e5;
00339     flatmesh.xMin->data[1] = -1.e3;
00340     flatmesh.xMax->data[0] = 2.0e5;
00341     flatmesh.xMax->data[1] = 1.e3;
00342     */
00343 
00344 
00345     /* Compute the mesh, and clean up local memory. */
00346     SUB( LALCreateFlatMesh( &status, &mesh, &flatmesh ), &status );
00347     SUB( LALSDestroyVector( &status, &(flatmesh.xMin) ), &status );
00348     SUB( LALSDestroyVector( &status, &(flatmesh.xMax) ), &status );
00349     SUB( LALSDestroyVectorSequence( &status, &matrix ), &status );
00350     SUB( LALSDestroyVectorSequence( &status, &matrixInv ), &status );
00351     SUB( LALSDestroyVectorSequence( &status, &corners ), &status );
00352   }
00353 
00354   /* Prepare to print result. */
00355   {
00356     UINT4 i;
00357     InspiralBankParams   bankParams; 
00358     InspiralCoarseBankIn coarseIn;
00359     INT4 valid;
00360     UINT4 k=0;
00361     REAL4 *data;
00362 
00363     coarseIn.mmCoarse = 0.70;
00364     coarseIn.mmFine = 0.97;
00365     coarseIn.fLower = 40.L;
00366     coarseIn.fUpper = 2000.L;
00367     coarseIn.iflso = 0.0L;
00368     coarseIn.tSampling = 4096.L;
00369     coarseIn.order = twoPN;
00370     coarseIn.space = Tau0Tau3;
00371     coarseIn.approximant = TaylorT1;
00372 
00373     coarseIn.mMin = 1.0;
00374     coarseIn.mMax = 20.0;
00375     coarseIn.MMax = coarseIn.mMax * 2.;
00376 
00377     coarseIn.massRange = MinMaxComponentMass; 
00378     /* coarseIn.massRange = MinComponentMassMaxTotalMass;*/
00379 
00380     /* minimum value of eta */
00381     coarseIn.etamin = coarseIn.mMin * ( coarseIn.MMax - coarseIn.mMin) / pow(coarseIn.MMax,2.);
00382 
00383     /* Print out the template parameters */
00384     i = mesh->length;
00385     data = mesh->data;
00386     while ( i--) {
00387         /*
00388         Retain only those templates that have meaningful masses:
00389         */
00390         bankParams.x0 = (REAL8) data[k++];
00391         bankParams.x1 = (REAL8) data[k++];
00392         LALInspiralValidParams(&status, &valid, bankParams, coarseIn);
00393         if (valid) fprintf(fp, "%10.3e %10.3e\n", bankParams.x0, bankParams.x1);
00394     }
00395   }
00396 
00397 
00398   fclose(fp);
00399   /* Free the mesh, and exit. */
00400   SUB( LALSDestroyVectorSequence( &status, &mesh ), &status );
00401   LALDDestroyVector(&status, &shf.data);
00402   LALCheckMemoryLeaks();
00403   INFO( FLATMESHTESTC_MSGENORM );
00404   return FLATMESHTESTC_ENORM;
00405 }
00406 
00407 
00408 static void
00409 GetInspiralMoments (
00410                 LALStatus            *status,
00411                 InspiralMomentsEtc   *moments,
00412                 REAL8FrequencySeries *psd,
00413                 InspiralTemplate     *params )
00414 {
00415 
00416    UINT4 k;
00417    InspiralMomentsIn in;
00418 
00419    INITSTATUS (status, "GetInspiralMoments", FLATMESHTESTC);
00420    ATTATCHSTATUSPTR(status);
00421   
00422    ASSERT (params, status, LALINSPIRALBANKH_ENULL, LALINSPIRALBANKH_MSGENULL);
00423    ASSERT (params->fLower>0, status, LALINSPIRALBANKH_ENULL, LALINSPIRALBANKH_MSGENULL);
00424    ASSERT (moments, status, LALINSPIRALBANKH_ENULL, LALINSPIRALBANKH_MSGENULL);
00425    ASSERT (psd, status, LALINSPIRALBANKH_ENULL, LALINSPIRALBANKH_MSGENULL);
00426 
00427    moments->a01 = 3.L/5.L;
00428    moments->a21 = 11.L * LAL_PI/12.L;
00429    moments->a22 = 743.L/2016.L * pow(25.L/(2.L*LAL_PI*LAL_PI), 1.L/3.L);
00430    moments->a31 = -3.L/2.L;
00431    moments->a41 = 617.L * LAL_PI * LAL_PI / 384.L;
00432    moments->a42 = 5429.L/5376.L * pow ( 25.L * LAL_PI/2.L, 1.L/3.L);
00433    moments->a43 = 1.5293365L/1.0838016L * pow(5.L/(4.L*pow(LAL_PI,4.L)), 1.L/3.L);
00434    
00435    /* setup the input structure needed in the computation of the moments */
00436 
00437    in.shf = psd;
00438    in.shf->f0 /= params->fLower;
00439    in.shf->deltaF /= params->fLower;
00440    in.xmin = params->fLower/params->fLower;
00441    in.xmax = params->fCutoff/params->fLower;
00442            
00443    /* First compute the norm */
00444 
00445    in.norm = 1.L;
00446    in.ndx = 7.L/3.L; 
00447    LALInspiralMoments(status->statusPtr, &moments->j[7], in); 
00448    CHECKSTATUSPTR(status);
00449    in.norm = moments->j[7];
00450 
00451    if (lalDebugLevel & LALINFO)
00452    {
00453            fprintf (stderr, "a01=%e a21=%e a22=%e a31=%e a41=%e a42=%e a43=%e \n", 
00454                            moments->a01, moments->a21, moments->a22, moments->a31, 
00455                            moments->a41, moments->a42, moments->a43);
00456    
00457            fprintf(stderr, "j7=%e\n", moments->j[7]);
00458    }
00459 
00460    /* Normalised moments of the noise PSD from 1/3 to 17/3. */
00461 
00462    for (k=1; k<=17; k++)
00463    {
00464            in.ndx = (REAL8) k /3.L; 
00465            LALInspiralMoments(status->statusPtr,&moments->j[k],in);  
00466            CHECKSTATUSPTR(status);
00467            if (lalDebugLevel==1) fprintf(stderr, "j%1i=%e\n", k,moments->j[k]);
00468    }
00469    in.shf->deltaF *= params->fLower;
00470    in.shf->f0 *= params->fLower;
00471   
00472    DETATCHSTATUSPTR(status);
00473    RETURN (status);
00474 }