Commit 04fd4c0e authored by Matthew David Pitkin's avatar Matthew David Pitkin

LALInferenceGenerateROQ.c: check normalisation of new basis is not zero

 - if a new basis is entirely encompases within the current basis within
   numerical precision it can consist entirely of zeros. When normalising
   a basis vector of zeros it leads to a basis of NaNs, which gets added
   to the current basis set. This will give NaN values in likelihoods
   calculated using the ROQ.
 - this patch checks if the normalisation is zero and if so will not add
   the new basis to the current basis as it is unnecessary.
parent 4a6e0a04
......@@ -792,7 +792,7 @@ REAL8 LALInferenceGenerateREAL8OrthonormalBasis(REAL8Array **RBin,
REAL8Array *RB = NULL;
UINT4Vector *dims = NULL;
gsl_matrix_view RBview;
/* this memory should be freed here */
ts_el = gsl_vector_alloc(cols);
last_rb = gsl_vector_alloc(cols);
......@@ -889,13 +889,16 @@ REAL8 LALInferenceGenerateREAL8OrthonormalBasis(REAL8Array **RBin,
}
//////////// end check ////////////
gpts->data[dim_RB] = worst_app;
//greedy_err[dim_RB] = worst_err; // These aren't used or output at all
gpts->data[dim_RB] = worst_app;
/* add worst approximated solution to basis set */
gsl_matrix_get_row(ortho_basis, &TSview.matrix, worst_app);
iterated_modified_gm(ru, ortho_basis, &RBview.matrix, &deltaview.vector, dim_RB); /* use IMGS */ // TODO: need REAL8 version of this
iterated_modified_gm(ru, ortho_basis, &RBview.matrix, &deltaview.vector, dim_RB); /* use IMGS */
/* check normalisation of generated orthogonal basis is not NaN (cause by a new orthogonal basis
having zero residual with the current basis) - if this is the case do not add the new basis. */
if ( gsl_isnan(gsl_vector_get(ru, dim_RB)) ){ break; }
/* add to reduced basis */
dims->data[0] = dim_RB+1; /* add row */
RB = XLALResizeREAL8Array( RB, dims );
......@@ -977,7 +980,6 @@ REAL8 LALInferenceGenerateCOMPLEX16OrthonormalBasis(COMPLEX16Array **RBin,
size_t max_RB = rows;
//REAL8 greedy_err[max_RB]; /* approximate error */
UINT4Vector *gpts = NULL;
gpts = XLALCreateUINT4Vector(max_RB); /* selected greedy points (row selection) */
*greedypoints = gpts;
......@@ -1036,7 +1038,6 @@ REAL8 LALInferenceGenerateCOMPLEX16OrthonormalBasis(COMPLEX16Array **RBin,
gpts->data[0] = 0;
UINT4 dim_RB = 1;
//greedy_err[0] = 1.0;
/* loop to find reduced basis */
while( 1 ){
......@@ -1093,12 +1094,15 @@ REAL8 LALInferenceGenerateCOMPLEX16OrthonormalBasis(COMPLEX16Array **RBin,
//////////// end check ////////////
gpts->data[dim_RB] = worst_app;
//greedy_err[dim_RB] = worst_err; // These aren't used or output at all
/* add worst approximated solution to basis set */
gsl_matrix_complex_get_row(ortho_basis, &TSview.matrix, worst_app);
iterated_modified_gm_complex(ru, ortho_basis, &RBview.matrix, &deltaview.vector, dim_RB); /* use IMGS */
/* check normalisation of generated orthogonal basis is not NaN (cause by a new orthogonal basis
having zero residual with the current basis) - if this is the case do not add the new basis. */
if ( gsl_isnan(GSL_REAL(gsl_vector_complex_get(ru, dim_RB))) ){ break; }
/* add to reduced basis */
dims->data[0] = dim_RB+1; /* add row */
RB = XLALResizeCOMPLEX16Array( RB, dims );
......@@ -1345,7 +1349,7 @@ REAL8 LALInferenceEnrichREAL8Basis(const REAL8Vector *delta,
if ( keepcounter == 0 ){ return maxprojerr; } // nothing was above the tolerance
// add vectors used to produce the orgial reduced basis
// add vectors used to produce the original reduced basis
UINT4Vector *dims = XLALCreateUINT4Vector( 2 );
dims->data[0] = RBin->dimLength->data[0] + keepcounter;
dims->data[1] = dlength;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment