1#include <lal/LALInferenceGenerateROQ.h>
2#include <lal/LALConstants.h>
3#include <lal/XLALGSL.h>
4#include <gsl/gsl_randist.h>
20#define TOLERANCE 10e-12
26double calc_phase(
double frequency,
double Mchirp);
29double real_model(
double frequency,
double Mchirp,
double modperiod);
38double real_model(
double frequency,
double Mchirp,
double modperiod){
47 REAL8Array *TS = NULL, *TSquad = NULL, *cTSquad = NULL;
53 size_t k = 0,
j = 0,
i = 0;
55 REAL8Array *RBlinear = NULL, *RBquad = NULL, *cRBquad = NULL;
68 TSdims->
data[0] = TSsize;
75 gsl_matrix_view TSview, TSviewquad, cTSviewquad;
76 TSview = gsl_matrix_view_array( TS->
data, TSdims->
data[0], TSdims->
data[1] );
77 TSviewquad = gsl_matrix_view_array( TSquad->data, TSdims->
data[0], TSdims->
data[1] );
78 cTSviewquad = gsl_matrix_view_array( cTSquad->data, TSdims->
data[0], TSdims->
data[1] );
79 gsl_matrix_complex_view cTSview;
80 cTSview = gsl_matrix_complex_view_array( (
double *)cTS->
data, TSdims->
data[0], TSdims->
data[1] );
85 double fmin0 = 48, fmax0 = 256, f0 = 0., m0 = 0.;
86 double df = (fmax0-fmin0)/(wl-1.);
87 double Mcmax = 2., Mcmin = 1.5, Mc = 0.;
88 double periodmax = 1./99.995, periodmin = 1./100., modperiod = 0.;
96 const gsl_rng_type *
T;
100 r = gsl_rng_alloc(
T);
103 for (
k=0;
k < TSsize;
k++ ){
104 Mc = pow(pow(Mcmin, 5./3.) + (
double)
k*(pow(Mcmax, 5./3.)-pow(Mcmin, 5./3.))/((
double)TSsize-1), 3./5.);
105 modperiod = gsl_ran_flat(
r, periodmin, periodmax);
107 for (
j=0;
j < wl;
j++ ){
108 f0 = fmin0 + (double)
j*(fmax0-fmin0)/((double)wl-1.);
114 GSL_SET_COMPLEX(&gctmp, creal(ctmp), cimag(ctmp));
116 gsl_matrix_set(&TSview.matrix,
k,
j, m0);
117 gsl_matrix_set(&TSviewquad.matrix,
k,
j, m0*m0);
118 gsl_matrix_complex_set(&cTSview.matrix,
k,
j, gctmp);
119 gsl_matrix_set(&cTSviewquad.matrix,
k,
j, creal(ctmp*conj(ctmp)));
124 REAL8 maxprojerr = 0.;
133 fprintf(stderr,
"No. quadratic nodes (real) = %d, %d x %d; Maximum projection err. = %le\n", RBquad->dimLength->data[0], RBquad->dimLength->data[0], RBquad->dimLength->data[1], maxprojerr);
136 fprintf(stderr,
"No. quadratic nodes (complex) = %d, %d x %d; Maximum projection err. = %le\n", cRBquad->dimLength->data[0], cRBquad->dimLength->data[0], cRBquad->dimLength->data[1], maxprojerr);
172 for (
i=0;
i<wl;
i++ ){
173 data->data[
i] = gsl_ran_gaussian(
r, 1.0);
174 cdata->
data[
i] = gsl_ran_gaussian(
r, 1.0) + I*gsl_ran_gaussian(
r, 1.0);
184 double randMc = 1.873;
185 double randmp = 1./99.9989;
187 gsl_vector *modelfull = gsl_vector_alloc(wl);
190 gsl_vector_complex *cmodelfull = gsl_vector_complex_alloc(wl);
195 for (
i=0;
i<wl;
i++ ){
201 GSL_SET_COMPLEX(&gcval, creal(cval), cimag(cval));
202 gsl_vector_complex_set(cmodelfull,
i, gcval);
206 for (
i=0;
i<modelreduced->
length;
i++ ){
208 modelreduced->
data[
i] = rm;
210 for (
i=0;
i<modelreducedquad->
length;
i++ ){
212 modelreducedquad->
data[
i] = rm*rm;
214 for (
i=0;
i<cmodelreduced->
length;
i++ ){
216 cmodelreduced->
data[
i] = crm;
218 for (
i=0;
i<cmodelreducedquad->
length;
i++ ){
220 cmodelreducedquad->
data[
i] = creal(crm*conj(crm));
226 struct timeval t1, t2, t3, t4;
232 gettimeofday(&t1, NULL);
233 XLAL_CALLGSL( gsl_blas_ddot(modelfull, modelfull, &mmfull) );
234 gettimeofday(&t2, NULL);
237 gettimeofday(&t3, NULL);
239 gettimeofday(&t4, NULL);
241 dt1 = (double)((t2.tv_sec + t2.tv_usec*1.e-6) - (t1.tv_sec + t1.tv_usec*1.e-6));
242 dt2 = (double)((t4.tv_sec + t4.tv_usec*1.e-6) - (t3.tv_sec + t3.tv_usec*1.e-6));
243 fprintf(stderr,
"Real Signal:\n - M dot M (full) = %le [%.9lf s], M dot M (reduced) = %le [%.9lf s], time ratio = %lf\n", mmfull, dt1, mmred, dt2, dt1/dt2);
247 gsl_vector_view dataview = gsl_vector_view_array(
data->data, wl);
248 gettimeofday(&t1, NULL);
249 XLAL_CALLGSL( gsl_blas_ddot(&dataview.vector, modelfull, &dmfull) );
250 gettimeofday(&t2, NULL);
253 gettimeofday(&t3, NULL);
255 gettimeofday(&t4, NULL);
257 dt1 = (double)((t2.tv_sec + t2.tv_usec*1.e-6) - (t1.tv_sec + t1.tv_usec*1.e-6));
258 dt2 = (double)((t4.tv_sec + t4.tv_usec*1.e-6) - (t3.tv_sec + t3.tv_usec*1.e-6));
259 fprintf(stderr,
" - D dot M (full) = %le [%.9lf s], D dot M (reduced) = %le [%.9lf s], time ratio = %lf\n", dmfull, dt1, dmred, dt2, dt1/dt2);
262 double Lfull, Lred, Lfrac;
264 Lfull = mmfull - 2.*dmfull;
265 Lred = mmred - 2.*dmred;
266 Lfrac = 100.*fabs(Lfull-Lred)/fabs(Lfull);
268 fprintf(stderr,
" - Fractional difference in log likelihoods = %lf%%\n", Lfrac);
271 gsl_vector_free(modelfull);
278 if ( Lfrac >
LTOL ) {
return 1; }
282 REAL8 cmmred, cmmfull;
283 gsl_complex cmmfulltmp;
284 gettimeofday(&t1, NULL);
285 XLAL_CALLGSL( gsl_blas_zdotc(cmodelfull, cmodelfull, &cmmfulltmp) );
286 cmmfull = GSL_REAL(cmmfulltmp);
287 gettimeofday(&t2, NULL);
289 gettimeofday(&t3, NULL);
291 gettimeofday(&t4, NULL);
292 dt1 = (double)((t2.tv_sec + t2.tv_usec*1.e-6) - (t1.tv_sec + t1.tv_usec*1.e-6));
293 dt2 = (double)((t4.tv_sec + t4.tv_usec*1.e-6) - (t3.tv_sec + t3.tv_usec*1.e-6));
294 fprintf(stderr,
"Complex Signal:\n - M dot M (full) = %le [%.9lf s], M dot M (reduced) = %le [%.9lf s], time ratio = %lf\n", cmmfull, dt1, cmmred, dt2, dt1/dt2);
297 gsl_complex cdmfulltmp;
298 gsl_vector_complex_view cdataview = gsl_vector_complex_view_array((
double*)cdata->
data, wl);
299 gettimeofday(&t1, NULL);
300 XLAL_CALLGSL( gsl_blas_zdotc(&cdataview.vector, cmodelfull, &cdmfulltmp) );
301 cdmfull = GSL_REAL(cdmfulltmp) + I*GSL_IMAG(cdmfulltmp);
302 gettimeofday(&t2, NULL);
304 gettimeofday(&t3, NULL);
306 gettimeofday(&t4, NULL);
308 dt1 = (double)((t2.tv_sec + t2.tv_usec*1.e-6) - (t1.tv_sec + t1.tv_usec*1.e-6));
309 dt2 = (double)((t4.tv_sec + t4.tv_usec*1.e-6) - (t3.tv_sec + t3.tv_usec*1.e-6));
310 fprintf(stderr,
" - D dot M (full) = %le [%.9lf s], D dot M (reduced) = %le [%.9lf s], time ratio = %lf\n", creal(cdmfull), dt1, creal(cdmred), dt2, dt1/dt2);
313 Lfull = cmmfull - 2.*creal(cdmfull);
314 Lred = cmmred - 2.*creal(cdmred);
315 Lfrac = 100.*fabs(Lfull-Lred)/fabs(Lfull);
317 fprintf(stderr,
" - Fractional difference in log likelihoods = %lf%%\n", Lfrac);
320 gsl_vector_complex_free(cmodelfull);
332 if ( Lfrac >
LTOL ) {
return 1; }
COMPLEX16 LALInferenceROQCOMPLEX16DotProduct(COMPLEX16Vector *weights, COMPLEX16Vector *model)
Calculate the dot product of two complex vectors using the ROQ iterpolant.
REAL8 LALInferenceGenerateCOMPLEX16OrthonormalBasis(COMPLEX16Array **RBin, const REAL8Vector *delta, REAL8 tolerance, COMPLEX16Array **TS, UINT4Vector **greedypoints)
Create a orthonormal basis set from a training set of complex waveforms.
void LALInferenceRemoveCOMPLEXROQInterpolant(LALInferenceCOMPLEXROQInterpolant *a)
Free memory for a LALInferenceCOMPLEXROQInterpolant.
COMPLEX16Vector * LALInferenceGenerateCOMPLEX16LinearWeights(COMPLEX16Array *B, COMPLEX16Vector *data, REAL8Vector *vars)
Create the weights for the ROQ interpolant for the complex data and model dot product <d|h>
REAL8Vector * LALInferenceGenerateREAL8LinearWeights(REAL8Array *B, REAL8Vector *data, REAL8Vector *vars)
Create the weights for the ROQ interpolant for the linear data and model dot product <d|h>
REAL8 LALInferenceROQREAL8DotProduct(REAL8Vector *weights, REAL8Vector *model)
Calculate the dot product of two vectors using the ROQ iterpolant.
REAL8Vector * LALInferenceGenerateQuadraticWeights(REAL8Array *B, REAL8Vector *vars)
Create the weights for the ROQ interpolant for the model quadratic model term real(<h|h>).
LALInferenceCOMPLEXROQInterpolant * LALInferenceGenerateCOMPLEXROQInterpolant(COMPLEX16Array *RB)
Create a complex empirical interpolant from a set of orthonormal basis functions.
REAL8 LALInferenceGenerateREAL8OrthonormalBasis(REAL8Array **RBin, const REAL8Vector *delta, REAL8 tolerance, REAL8Array **TS, UINT4Vector **greedypoints)
Create a orthonormal basis set from a training set of real waveforms.
void LALInferenceRemoveREALROQInterpolant(LALInferenceREALROQInterpolant *a)
Free memory for a LALInferenceREALROQInterpolant.
LALInferenceREALROQInterpolant * LALInferenceGenerateREALROQInterpolant(REAL8Array *RB)
Create a real empirical interpolant from a set of orthonormal basis functions.
double real_model(double frequency, double Mchirp, double modperiod)
COMPLEX16 imag_model(double frequency, double Mchirp, double modperiod)
double calc_phase(double frequency, double Mchirp)
void LALCheckMemoryLeaks(void)
#define XLAL_CALLGSL(statement)
void XLALDestroyREAL8Array(REAL8Array *)
COMPLEX16Array * XLALCreateCOMPLEX16Array(UINT4Vector *)
void XLALDestroyCOMPLEX16Array(COMPLEX16Array *)
REAL8Array * XLALCreateREAL8Array(UINT4Vector *)
void XLALDestroyUINT4Vector(UINT4Vector *vector)
REAL8Vector * XLALCreateREAL8Vector(UINT4 length)
COMPLEX16Vector * XLALCreateCOMPLEX16Vector(UINT4 length)
void XLALDestroyREAL8Vector(REAL8Vector *vector)
void XLALDestroyCOMPLEX16Vector(COMPLEX16Vector *vector)
UINT4Vector * XLALCreateUINT4Vector(UINT4 length)
A structure to hold a complex (double precision) interpolant matrix and interpolation node indices.
COMPLEX16Array * B
The interpolant matrix.
UINT4 * nodes
The nodes (indices) for the interpolation.
A structure to hold a real (double precision) interpolant matrix and interpolation node indices.
UINT4 * nodes
The nodes (indices) for the interpolation.
REAL8Array * B
The interpolant matrix.