Skip to content
GitLab
Commit 0c2cf112 authored by Miroslav Shaltev's avatar Miroslav Shaltev
Browse files

fixes related to hybrid search

parent 9b38e68d
Hide whitespace changes
Inline Side-by-side
src/FollowUp/Nomad/libFStatNomad.cpp
View file @ 0c2cf112
......@@ -1131,7 +1131,7 @@ void FStatNomad::SetUpSFTs( LALStatus *status, /**< pointer to LALStatus struc
XLALPrintError("%s: XLALSetupFstatInput() failed with errno=%d", __func__, xlalErrno);
ABORT ( status, FSTATFCNOMAD_EXLAL, FSTATFCNOMAD_MSGEXLAL );
}
LogPrintf (LOG_NORMAL, "FstatMethod used: '%s'\n", XLALGetFstatInputMethodName( (*p_Fstat_in_vec)->data[k] ) );
LogPrintf (LOG_DEBUG, "FstatMethod used: '%s'\n", XLALGetFstatInputMethodName( (*p_Fstat_in_vec)->data[k] ) );
/* get SFT detectors and timestamps */
const MultiLALDetector *multiIFO = XLALGetFstatInputDetectors( (*p_Fstat_in_vec)->data[k] );
if ( multiIFO == NULL ) {
......@@ -1891,7 +1891,7 @@ INT4 FStatNomad::compute_min_extent(gsl_vector *minext, const gsl_matrix *g) {
gsl_linalg_LU_invert(m, p, im);
for (INT4 i = 0; i < maxdim; i++) {
gsl_vector_set(minext,i,2*sqrt(fabs(gsl_matrix_get(im,i,i))));
gsl_vector_set(minext,i,MINEXTSCALE*2*sqrt(fabs(gsl_matrix_get(im,i,i))));
}
for (INT4 i = 0; i < maxdim; i++) {
......@@ -3406,7 +3406,6 @@ void SCHybridSearch::search( Mads & mads, int &nb_search_pts, bool &stop, stop_t
return;
}
// initial display:
const Display &out = _p.out();
dd_type display_degree = out.get_search_dd();
if ( display_degree == FULL_DISPLAY ){
......@@ -3415,20 +3414,14 @@ void SCHybridSearch::search( Mads & mads, int &nb_search_pts, bool &stop, stop_t
out << std::endl << open_block ( oss.str() ) << std::endl;
}
// the barriers:
Barrier & true_barrier = mads.get_true_barrier();
Barrier & sgte_barrier = mads.get_sgte_barrier();
const Barrier & active_barrier = mads.get_active_barrier();
// point x:
Double best_f;
bool x_feas = true;
const Eval_Point * x = active_barrier.get_best_feasible();
Double best_f;
bool x_feas = true;
const Eval_Point *x = mads.get_best_feasible();
if ( x ) {
best_f = x->get_f();
}
else {
x = active_barrier.get_best_infeasible();
x = mads.get_best_infeasible();
x_feas = false;
}
......@@ -3439,10 +3432,8 @@ void SCHybridSearch::search( Mads & mads, int &nb_search_pts, bool &stop, stop_t
return;
}
// get the signature:
Signature * signature = x->get_signature();
if ( !signature )
{
Signature *signature = x->get_signature();
if ( !signature ){
if ( display_degree == FULL_DISPLAY )
out.close_block ( "end of HybridSearch (no signature)" );
return;
......@@ -3454,26 +3445,27 @@ void SCHybridSearch::search( Mads & mads, int &nb_search_pts, bool &stop, stop_t
out.close_block ( "end of HybridSearch (incompatible signature)" );
return;
}
// template bank
LatticeTilingIterator *itr;
LatticeTilingLocator *tloc;
LatticeTiling *ptiling;
if (partiallattice){
LatticeTiling *ptiling = XLALCreateLatticeTiling(mdim);
ptiling = XLALCreateLatticeTiling(mdim);
XLAL_CHECK_MAIN(ptiling != NULL, XLAL_EFUNC);
// set bounds for each dimension using a rough estimate
std::vector<tDirection> pdirections;
for (INT4 i = 0; i < directions.size(); i++) {
pdirections.push_back(directions.at(i));
pdirections.at(i).min = (*x)[i].value()/directions.at(i).scale - gsl_vector_get(minext,i);
pdirections.at(i).max = (*x)[i].value()/directions.at(i).scale + gsl_vector_get(minext,i);
// printf("size in direction %d:%.16e,%.16e\n",i,gsl_matrix_get(g,i,i),pdirections.at(i).max-pdirections.at(i).min);
LogPrintf(LOG_DEBUG,"size in direction %d:%.16e,%.16e\n",i,gsl_matrix_get(g,i,i),pdirections.at(i).max-pdirections.at(i).min);
}
for (INT4 i = 0; i < directions.size(); i++) {
......@@ -3483,7 +3475,7 @@ void SCHybridSearch::search( Mads & mads, int &nb_search_pts, bool &stop, stop_t
const size_t ln = XLALTotalLatticeTilingDimensions(ptiling);
// Set lattice and metric
// printf("g dimensions [%d,%d]\n",g->size1,g->size2);
LogPrintf(LOG_DEBUG,"g dimensions [%d,%d]\n",g->size1,g->size2);
XLAL_CHECK_MAIN(XLALSetTilingLatticeAndMetric(ptiling, lattice, g, m) == XLAL_SUCCESS, XLAL_EFUNC);
......@@ -3496,19 +3488,17 @@ void SCHybridSearch::search( Mads & mads, int &nb_search_pts, bool &stop, stop_t
// Print number of templates
UINT8 ntemplates = XLALTotalLatticeTilingPoints(itr);
XLAL_CHECK_MAIN(ntemplates > 0, XLAL_EFUNC);
// printf("HybridSearch total number of partial templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
LogPrintf(LOG_DETAIL,"HybridSearch total number of partial templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
}
else {
itr = XLALCreateLatticeTilingIterator( tiling, n );
tloc = XLALCreateLatticeTilingLocator(tiling);
// Print number of templates
UINT8 ntemplates = XLALTotalLatticeTilingPoints(itr);
XLAL_CHECK_MAIN(ntemplates > 0, XLAL_EFUNC);
// printf("HybridSearch total number of templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
LogPrintf(LOG_DEBUG,"HybridSearch total number of templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
}
gsl_matrix *points = gsl_matrix_calloc(n,1);
for (INT4 i = 0; i < n; i++){
gsl_matrix_set(points,i,0,(*x)[i].value()/directions.at(i).scale);
......@@ -3518,59 +3508,54 @@ void SCHybridSearch::search( Mads & mads, int &nb_search_pts, bool &stop, stop_t
gsl_matrix *GAMAT( nearest_points, n, np );
/*
for (INT4 i = 0; i < total; i++){
gsl_vector *tpoint = gsl_vector_alloc (n);
INT4 r = XLALNextLatticeTilingPoint( itr, tpoint );
printf("r: %d\n",r);
for (INT4 j = 0; j < n; j++) {
printf("%.16f\t",gsl_vector_get(tpoint,j));
}
printf("\n");
gsl_vector_free(tpoint);
}
*/
XLALNearestLatticeTilingPoints (tloc, points, &nearest_points, NULL);
Evaluator_Control & ev_control = mads.get_evaluator_control();
Evaluator_Control & ev_control = mads.get_evaluator_control();
for (INT4 i = 0; i < nearest_points->size2; i++) {
gsl_vector *cv = gsl_vector_alloc(n);
gsl_matrix_get_col(cv,nearest_points,i);
Eval_Point * tk = new Eval_Point;
tk->set ( n , 1 );
tk->set_signature ( signature );
Eval_Point * tk = new Eval_Point;
tk->set ( n , 1 );
tk->set_signature ( signature );
for (INT4 j = 0; j < n; j++) {
(*tk)[j] = (directions.at(j).scale * gsl_vector_get(cv,j));
for (INT4 j = 0; j < n; j++) {
(*tk)[j] = (directions.at(j).scale * gsl_vector_get(cv,j));
}
if ( display_degree == NOMAD::FULL_DISPLAY ){
out << "candidate";
out << " (before projection)";
out << ": ( " << *tk << " )" << std::endl;
}
tk->project_to_mesh(*x,signature->get_mesh()->get_delta(),signature->get_lb(),signature->get_ub() );
if ( display_degree == NOMAD::FULL_DISPLAY ){
out << "candidate";
out << " (after projection)";
out << ": ( " << *tk << " )" << std::endl;
}
ev_control.add_eval_point ( tk, _p.out().get_search_dd(), false, Double(), Double(), Double(), Double());
gsl_matrix_free(points);
gsl_matrix_free(nearest_points);
gsl_vector_free(cv);
}
ev_control.add_eval_point ( tk, _p.out().get_search_dd(), false, Double(), Double(), Double(), Double());
nb_search_pts = 1;
gsl_matrix_free(points);
gsl_matrix_free(nearest_points);
gsl_vector_free(cv);
}
if (partiallattice){
XLALDestroyLatticeTiling(ptiling);
XLALDestroyLatticeTilingIterator(itr);
XLALDestroyLatticeTilingLocator(tloc);
}
new_feas_inc = new_infeas_inc = NULL;
ev_control.eval_list_of_points ( _type, mads.get_true_barrier(), mads.get_sgte_barrier(), mads.get_pareto_front(), stop, stop_reason, new_feas_inc ,new_infeas_inc, success);
nb_search_pts = 1;
// evaluation:
new_feas_inc = new_infeas_inc = NULL;
ev_control.eval_list_of_points ( _type ,
mads.get_true_barrier() ,
mads.get_sgte_barrier() ,
mads.get_pareto_front() ,
stop ,
stop_reason ,
new_feas_inc ,
new_infeas_inc ,
success );
}
void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_type &stop_reason, success_type &success, bool &count_search , const Eval_Point *& new_feas_inc, const Eval_Point *& new_infeas_inc ) {
......@@ -3584,7 +3569,6 @@ void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_t
return;
}
// initial display:
const Display &out = _p.out();
dd_type display_degree = out.get_search_dd();
if ( display_degree == FULL_DISPLAY ){
......@@ -3593,20 +3577,14 @@ void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_t
out << std::endl << open_block ( oss.str() ) << std::endl;
}
// the barriers:
Barrier & true_barrier = mads.get_true_barrier();
Barrier & sgte_barrier = mads.get_sgte_barrier();
const Barrier & active_barrier = mads.get_active_barrier();
// point x:
Double best_f;
bool x_feas = true;
const Eval_Point * x = active_barrier.get_best_feasible();
Double best_f;
bool x_feas = true;
const Eval_Point *x = mads.get_best_feasible();
if ( x ) {
best_f = x->get_f();
}
else {
x = active_barrier.get_best_infeasible();
x = mads.get_best_infeasible();
x_feas = false;
}
......@@ -3617,10 +3595,8 @@ void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_t
return;
}
// get the signature:
Signature * signature = x->get_signature();
if ( !signature )
{
if ( !signature ){
if ( display_degree == FULL_DISPLAY )
out.close_block ( "end of HybridSearch (no signature)" );
return;
......@@ -3632,26 +3608,24 @@ void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_t
out.close_block ( "end of HybridSearch (incompatible signature)" );
return;
}
// template bank
LatticeTilingIterator *itr;
LatticeTilingLocator *tloc;
LatticeTilingIterator *itr;
LatticeTilingLocator *tloc;
LatticeTiling *ptiling;
if (partiallattice){
LatticeTiling *ptiling = XLALCreateLatticeTiling(mdim);
ptiling = XLALCreateLatticeTiling(mdim);
XLAL_CHECK_MAIN(ptiling != NULL, XLAL_EFUNC);
// set bounds for each dimension using a rough estimate
std::vector<tDirection> pdirections;
for (INT4 i = 0; i < directions.size(); i++) {
pdirections.push_back(directions.at(i));
pdirections.at(i).min = (*x)[i].value()/directions.at(i).scale - gsl_vector_get(minext,i);
pdirections.at(i).max = (*x)[i].value()/directions.at(i).scale + gsl_vector_get(minext,i);
// printf("size in direction %d:%.16e,%.16e,%.16e\n",i,gsl_matrix_get(g,i,i),(*x)[i].value()/directions.at(i).scale,pdirections.at(i).max-pdirections.at(i).min);
LogPrintf(LOG_DEBUG,"size in direction %d:%.16e,%.16e,%.16e\n",i,gsl_matrix_get(g,i,i),(*x)[i].value()/directions.at(i).scale,pdirections.at(i).max-pdirections.at(i).min);
}
for (INT4 i = 0; i < directions.size(); i++) {
......@@ -3660,12 +3634,8 @@ void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_t
const size_t n = XLALTotalLatticeTilingDimensions(ptiling);
// Set lattice and metric
XLAL_CHECK_MAIN(XLALSetTilingLatticeAndMetric(ptiling, lattice, g, m) == XLAL_SUCCESS, XLAL_EFUNC);
// Create a lattice iterator
itr = XLALCreateLatticeTilingIterator(ptiling, n);
XLAL_CHECK_MAIN(itr != NULL, XLAL_EFUNC);
......@@ -3674,17 +3644,14 @@ void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_t
// Print number of templates
UINT8 ntemplates = XLALTotalLatticeTilingPoints(itr);
XLAL_CHECK_MAIN(ntemplates > 0, XLAL_EFUNC);
// printf("HybridSearch total number of partial templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
LogPrintf(LOG_DETAIL,"HybridSearch total number of partial templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
}
else {
itr = XLALCreateLatticeTilingIterator( tiling, n );
tloc = XLALCreateLatticeTilingLocator(tiling);
UINT8 ntemplates = XLALTotalLatticeTilingPoints( itr );
// printf("HybridSearch total number of templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
LogPrintf(LOG_DEBUG,"HybridSearch total number of templates: %" LAL_UINT8_FORMAT "\n", ntemplates);
}
gsl_matrix *points = gsl_matrix_calloc(n,1);
for (INT4 i = 0; i < n; i++){
......@@ -3695,47 +3662,52 @@ void FCHybridSearch::search ( Mads &mads, int &nb_search_pts, bool &stop, stop_t
gsl_matrix *GAMAT( nearest_points, n, np );
/*
for (INT4 i = 0; i < total; i++){
gsl_vector *tpoint = gsl_vector_alloc (n);
INT4 r = XLALNextLatticeTilingPoint( itr, tpoint );
printf("r: %d\n",r);
for (INT4 j = 0; j < n; j++) {
printf("%.16f\t",gsl_vector_get(tpoint,j));
}
printf("\n");
gsl_vector_free(tpoint);
}
*/
XLALNearestLatticeTilingPoints (tloc, points, &nearest_points, NULL);
Evaluator_Control & ev_control = mads.get_evaluator_control();
Evaluator_Control & ev_control = mads.get_evaluator_control();
for (INT4 i = 0; i < nearest_points->size2; i++) {
for (INT4 i = 0; i < nearest_points->size2; i++){
gsl_vector *cv = gsl_vector_alloc(n);
gsl_matrix_get_col(cv,nearest_points,i);
Eval_Point * tk = new Eval_Point;
tk->set ( n , 1 );
tk->set_signature ( signature );
Eval_Point * tk = new Eval_Point;
tk->set ( n , 1 );
tk->set_signature ( signature );
for (INT4 j = 0; j < n; j++) {
(*tk)[j] = (directions.at(j).scale * gsl_vector_get(cv,j));
}
for (INT4 j = 0; j < n; j++) {
(*tk)[j] = (directions.at(j).scale * gsl_vector_get(cv,j));
}
if ( display_degree == NOMAD::FULL_DISPLAY ){
out << "candidate";
out << " (before projection)";
out << ": ( " << *tk << " )" << std::endl;
}
tk->project_to_mesh(*x,signature->get_mesh()->get_delta(),signature->get_lb(),signature->get_ub() );
ev_control.add_eval_point ( tk, _p.out().get_search_dd(), false, Double(), Double(), Double(), Double());
if ( display_degree == NOMAD::FULL_DISPLAY ){
out << "candidate";
out << " (after projection)";
out << ": ( " << *tk << " )" << std::endl;
}
ev_control.add_eval_point ( tk, _p.out().get_search_dd(), false, Double(), Double(), Double(), Double());
gsl_matrix_free(points);
gsl_matrix_free(nearest_points);
gsl_vector_free(cv);
gsl_matrix_free(points);
gsl_matrix_free(nearest_points);
gsl_vector_free(cv);
}
nb_search_pts = 1;
nb_search_pts = 1;
new_feas_inc = new_infeas_inc = NULL;
ev_control.eval_list_of_points ( _type, mads.get_true_barrier(), mads.get_sgte_barrier(), mads.get_pareto_front(), stop, stop_reason, new_feas_inc, new_infeas_inc, success);
if (partiallattice){
XLALDestroyLatticeTiling(ptiling);
XLALDestroyLatticeTilingIterator(itr);
XLALDestroyLatticeTilingLocator(tloc);
}
new_feas_inc = new_infeas_inc = NULL;
ev_control.eval_list_of_points ( _type, mads.get_true_barrier(), mads.get_sgte_barrier(), mads.get_pareto_front(), stop, stop_reason, new_feas_inc, new_infeas_inc, success);
}
src/FollowUp/Nomad/libFStatNomad.h
View file @ 0c2cf112
......@@ -122,6 +122,9 @@ using namespace NOMAD;
#define MTPHASE 0
#define MTAVFSTAT 1
#define MINEXTSCALE 0.01
/******************************************************
* Protection against C++ name mangling
*/
......
Supports Markdown
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