common_utils.cpp

/*--------------------------------------------------------------------------*/
/*-------------------------- common_utils.cpp ------------------------------*/
/*--------------------------------------------------------------------------*/
/** @file
 * Implementations of the utility functions declared in common_utils.h.
 *
 * \author Antonio Frangioni \n
 *         Dipartimento di Informatica \n
 *         Universita' di Pisa \n
 *
 * \author Donato Meoli \n
 *         Dipartimento di Informatica \n
 *         Universita' di Pisa \n
 *
 * \copyright © by Antonio Frangioni, Donato Meoli
 */
/*--------------------------------------------------------------------------*/

#include "common_utils.h"

#include <algorithm>
#include <cctype>
#include <chrono>
#include <cmath>
#include <cstdlib>
#include <exception>
#include <fstream>
#include <iomanip>
#include <limits>
#include <sstream>
#include <list>
#include <map>
#include <numeric>
#include <typeinfo>

#include <SMSTypedefs.h>

#include <Configuration.h>

#include <DQuadFunction.h>

#include <FRowConstraint.h>

#include <FRealObjective.h>

#include <Objective.h>

#include <QuadFunction.h>

/*--------------------------------------------------------------------------*/
/*--------------------- MPI / UCX SAFE-DEFAULTS ----------------------------*/
/*--------------------------------------------------------------------------*/
/* Some SMS++ targets (InvestmentBlock, and the tools using SDDPBlock) pull
 * in libboost_mpi / libmpi transitively even when they never call
 * MPI_Init(). Where Open MPI / UCX are installed but no usable transport is
 * available (no IB, missing UCX vfs.sock, ...), their runtime can hang on
 * startup spinning on futex / X11 sockets, so TCP-only defaults are
 * pre-seeded here, before main(), with setenv( ... , 0 ): the third
 * argument is "overwrite = false", hence whoever has exported UCX_TLS /
 * OMPI_MCA_* already, as on a cluster with a real fabric, keeps control. */

namespace {

void set_default_env( const char * name , const char * value ) {
#ifdef _WIN32
 if( std::getenv( name ) == nullptr )
  _putenv_s( name , value );
#else
 setenv( name , value , 0 );
#endif
}

struct SmsppMpiSafeEnvInit {
 SmsppMpiSafeEnvInit() {
  set_default_env( "UCX_TLS"     , "tcp,self" );
  set_default_env( "OMPI_MCA_btl", "tcp,self" );
  set_default_env( "OMPI_MCA_pml", "ob1"      );
  // the hwloc GL component probes the GPU topology via XOpenDisplay(),
  // which may hang inside MPI_Init(); no SMS++ target has a use for it
  set_default_env( "HWLOC_COMPONENTS", "-gl"  );
  }
 };

static SmsppMpiSafeEnvInit smspp_mpi_safe_env_init_;

}  // anonymous namespace

/*--------------------------------------------------------------------------*/
/*----------------------- OUTPUT AND ERROR HANDLING ------------------------*/
/*--------------------------------------------------------------------------*/

// print an objective value, or the reason why there is none

void PrintResults( bool hs , int rtrn , double fo )
{
 if( hs ) {
  std::cout.setf( std::ios::scientific , std::ios::floatfield );
  std::cout << def << fo;
  }
 else
  if( rtrn == Solver::kInfeasible )
   std::cout << "    Unfeas";
  else
   if( rtrn == Solver::kUnbounded )
    std::cout << "      Unbounded";
   else
    std::cout << "      Error!";
 }

/*--------------------------------------------------------------------------*/
// what a Solver that failed says of itself
/* A Solver that returns anything from kError up delivered nothing, and the
 * line says so; which of the errors it is, however, is the only thing left
 * to go by when the failure does not repeat, as one of a parallel Solver may
 * well not, hence the code, and its name where the core gives it one. */

std::string error_token( int status )
{
 std::string what = "Error!";
 if( status == Solver::kError )            what += " (kError)";
 else if( status == Solver::kBlockLocked ) what += " (Block locked)";
 else                                      what += " (status " +
                                            std::to_string( status ) + ")";
 return( what );
 }

/*--------------------------------------------------------------------------*/
// print the exception that reached std::terminate(), then abort

void smspp_terminate( void )
{
 std::cerr << "Uncaught exception in executing SMS++:\n";
 try {
  std::rethrow_exception( std::current_exception() );
  }
 catch( const std::exception & e ) {
  std::cerr << "\tException type: " << typeid( e ).name() << "\n";
  std::cerr << "\tException message: " << e.what() << "\n";
  }
 catch( ... ) {
  std::cerr << "\tUnknown exception" << std::endl;
  }
 std::abort();  // or exit( 1 )
 }

/*--------------------------------------------------------------------------*/
// load a Block out of a text file, or stop with a clear message

void load_Block_or_exit( Block * block , const std::string & fn , char frmt )
{
 std::ifstream file( Block::get_filename_prefix() + fn , std::fstream::in );
 if( ! file.is_open() ) {
  std::cerr << "Error: cannot open " << Block::get_filename_prefix() + fn
            << std::endl;
  std::exit( 1 );
  }

 block->load( file , frmt );
 }

/*--------------------------------------------------------------------------*/

std::string first_Solver_of( const std::vector< std::string > & names ,
			     const std::string & what )
{
 for( const auto & name : names )
  if( Solver::has_Solver( name ) )
   return( name );

 std::string all;
 for( const auto & name : names )
  all += ( all.empty() ? "" : ", " ) + name;

 std::cerr << ANSI_YELLOW << "[WARNING] none of " << all
	   << " is in this build, so " << what << " is skipped"
	   << ANSI_RESET << std::endl;

 return( std::string() );

 }  // end( first_Solver_of )

/*--------------------------------------------------------------------------*/
// a random subset of k distinct indices out of 0 ... m - 1

Block::Subset GenerateRand( Block::Index m , Block::Index k ,
                            std::mt19937 & rg , bool ord )
{
 if( k > m ) {
  std::cerr << "error: GenerateRand( " << m << " , " << k << " )"
            << std::endl;
  exit( 1 );
  }

 Block::Subset sbst( m );
 std::iota( sbst.begin() , sbst.end() , 0 );
 std::shuffle( sbst.begin() , sbst.end() , rg );
 sbst.resize( k );
 if( ord )
  std::sort( sbst.begin() , sbst.end() );

 return( sbst );
 }

/*--------------------------------------------------------------------------*/
/*----------------------------- CONFIGURATION ------------------------------*/
/*--------------------------------------------------------------------------*/

// apply a BlockConfig, or dispatch a meta one to the sub-Block it names

void b_config_Block( Block * block , Configuration * b_config ,
                     const std::string & fn )
{
 // handle the special case of a "meta" BlockConfig: BlockConfig-ure all
 // Block whose classname() matches, or the "*" entry (if any) the others
 if( auto * mb =
     dynamic_cast< SimpleConfiguration< std::map< std::string ,
                                                  Configuration * > >
                                        * >( b_config ) ) {
  for_each_by_classname( block , mb->f_value ,
			 [ & ]( Block * b , Configuration * c ) {
			  auto bc = dynamic_cast< BlockConfig * >( c );
			  if( ! bc ) {
			   std::cerr << "Error: meta-Configuration for :Block "
				     << b->classname() << " in file " << fn
				     << " is not a BlockConfig" << std::endl;
			   exit( 1 );
			   }
			  auto cbc = bc->clone();
			  cbc->apply( b );
			  delete cbc;
			  } );

  return;  // all done
  }

 // apply() moves the Configuration of the BlockConfig into the Block, hence a
 // clone() is applied, as for the meta one: @p b_config stays whole, and
 // can be applied again, e.g., to the next copy of the same Block
 if( auto * bc = dynamic_cast< BlockConfig * >( b_config ) ) {
  auto cbc = bc->clone();
  cbc->apply( block );
  delete cbc;
  return;  // all done
  }

 std::cerr << "Error: " << fn
           << " does not contain a valid [meta]BlockConfig" << std::endl;
 exit( 1 );

 }  // end( b_config_Block )

/*--------------------------------------------------------------------------*/
// the list of the (meta-)BlockConfig of the formulations of a battery

SimpleConfiguration< std::vector< Configuration * > > *
 read_formulations( const std::string & fn )
{
 auto c = Configuration::deserialize( fn );
 auto vc = dynamic_cast< SimpleConfiguration<
                          std::vector< Configuration * > > * >( c );
 if( ! vc ) {
  std::cerr << "Error: " << fn << " does not contain a "
            << "SimpleConfiguration< std::vector< Configuration * > >"
            << std::endl;
  delete c;
  exit( 1 );
  }

 if( vc->f_value.empty() ) {
  std::cerr << "Error: " << fn << " lists no formulation" << std::endl;
  exit( 1 );
  }

 for( std::size_t k = 0 ; k < vc->f_value.size() ; ++k ) {
  auto e = vc->f_value[ k ];
  if( auto p = dynamic_cast< SimpleConfiguration< std::pair< Configuration * ,
                                       Configuration * > > * >( e ) )
   e = p->f_value.first;
  else
   if( auto q = dynamic_cast< SimpleConfiguration< std::vector<
                                         Configuration * > > * >( e ) )
    e = q->f_value.empty() ? nullptr : q->f_value.front();
  if( ! e ) {
   std::cerr << "Error: formulation " << k << " of " << fn << " is empty"
             << std::endl;
   exit( 1 );
   }
  }

 return( vc );
 }

/*--------------------------------------------------------------------------*/
// the Solver that an option of a tester names by its position, as given and
// where it is now [see name_Solver()]

namespace {

struct NamedSolver {
 std::string option;
 std::size_t given;
 std::size_t at;
 };

std::vector< NamedSolver > named_solvers;

}  // end( anonymous namespace )

/*--------------------------------------------------------------------------*/
/* Takes out of a BlockSolverConfig the Solver that this build does not have.
 * A configuration names every Solver that makes sense on the instances of a
 * battery, so that a run solves with all of them and compares; but which
 * Solver are there depends on the modules the build has and on the external
 * libraries each of them has found, and asking the factory for one that is
 * not there throws [see Solver::new_Solver()], which would kill the whole
 * run rather than that one comparison. What is left out is said, name by
 * name, so that the log of a run tells what it has actually solved with; it
 * is returned, last first, with its index in the configuration. */

using Dropped = std::vector< std::pair< Block::Index , std::string > >;

static Dropped drop_missing_Solvers( BlockSolverConfig * bsc ,
				     const std::string & fn )
{
 Dropped dropped;
 if( ! bsc )
  return( dropped );

 const auto asked = bsc->num_ComputeConfig();

 /* A BlockSolverConfig that names no Solver at all is how a Solver is
  * detached rather than a configuration that cannot be applied: there is
  * nothing to look for in the factory and nothing to skip. */
 if( ! asked )
  return( dropped );

 for( Block::Index i = asked ; i-- ; ) {
  auto name = bsc->get_SolverName( i );
  bool there = Solver::has_Solver( name );

  // a LagrangianDualSolver is only there if its inner Solver is
  if( there )
   if( auto cc = bsc->get_SolverConfig( i ) )
    for( const auto & pair : cc->str_pars )
     if( ( pair.first == "str_LDSlv_ISName" ) &&
	 ( ! Solver::has_Solver( pair.second ) ) ) {
      name += " with " + pair.second;
      there = false;
      }

  if( there )
   continue;

  std::cerr << ANSI_YELLOW << "[WARNING] " << name << " is not in this build"
	    << ", so " << fn << " solves without it" << ANSI_RESET << std::endl;

  dropped.emplace_back( i , name );
  bsc->remove_ComputeConfig( i );
  }

 const auto left = bsc->num_ComputeConfig();

 if( ! left ) {
  /* A configuration none of whose Solver is there cannot be applied to
   * anything: a run that went through the motions would report that nothing
   * failed, which is worse than saying that it cannot be done. The status is
   * the one of a test that has nothing to do rather than of one that failed,
   * which is what the batteries and ctest read it as
   * [see SKIP_RETURN_CODE and run_test of batch_common.sh]. */
  std::cerr << ANSI_YELLOW << "[SKIPPED] none of the " << asked
	    << " Solver that " << fn << " names is in this build"
	    << ANSI_RESET << std::endl;
  exit( 77 );
  }

 if( ( asked > 1 ) && ( left == 1 ) )
  std::cerr << ANSI_YELLOW << "[WARNING] " << fn << " names " << asked
	    << " Solver and this build has one of them, "
	    << bsc->get_SolverName( 0 )
	    << ": there is no Solver to cross-check it against, and what this "
	       "run says is what that one Solver says" << ANSI_RESET
	    << std::endl;

 return( dropped );
 }

/*--------------------------------------------------------------------------*/
/* What -E, -R and the options that name a Solver declare is positional on the
 * order in which the Solver are registered to the Block, hence the entry of
 * a Solver that is not there has to go with it: what is left would otherwise
 * be read off the Solver that takes its place, i.e. a tolerance, a "this one
 * solves a relaxation" or a role landing on another Solver. A single value
 * of -E or -R applies to every Solver and stays as it is [see solver_eps and
 * solver_relaxation], and an option that names a Solver that is not there
 * names none, which is said. The Solver of @p bsc sit from @p first on,
 * @p dropped is what drop_missing_Solvers() returned, last first. */

static void drop_positions( const Dropped & dropped , std::size_t first ,
			    const std::string & fn )
{
 for( const auto & [ i , name ] : dropped ) {
  const auto p = first + i;
  if( solver_eps.size() > 1 )
   if( p < solver_eps.size() )
    solver_eps.erase( solver_eps.begin() + p );
  if( solver_relaxation.size() > 1 )
   if( p < solver_relaxation.size() )
    solver_relaxation.erase( solver_relaxation.begin() + p );

  for( auto & n : named_solvers ) {
   if( ( n.at == std::string::npos ) || ( n.at < p ) )
    continue;
   if( n.at > p ) {
    --n.at;
    continue;
    }
   std::cerr << ANSI_YELLOW << "[WARNING] " << n.option << " names " << name
	     << ", which is not in this build, so " << fn << " solves as "
	        "without " << n.option << ANSI_RESET << std::endl;
   n.at = std::string::npos;
   }
  }
 }

/*--------------------------------------------------------------------------*/
/* Where the Solver of a BlockSolverConfig sit among those of the Block once
 * it is apply()-ed: after the @p before ones already there if it adds to them,
 * from the first otherwise, since it replaces them in their order. Its mode
 * is not public, but adding is the only one that leaves @p before + @p kept
 * Solver when neither is zero: setting leaves @p kept, and the differential
 * one the largest of the two. */

static std::size_t first_position( std::size_t before , std::size_t kept ,
				   std::size_t after )
{
 return( ( before && kept && ( after == before + kept ) ) ? before : 0 );
 }

/*--------------------------------------------------------------------------*/
// apply() @p bsc to @p block, taking out first the Solver that are not there

static void apply_present( BlockSolverConfig * bsc , Block * block ,
			   const std::string & fn , bool positional )
{
 const auto dropped = drop_missing_Solvers( bsc , fn );
 const auto before = block->get_registered_solvers().size();
 bsc->apply( block );
 if( positional && ( ! dropped.empty() ) )
  drop_positions( dropped , first_position( before ,
                                            bsc->num_ComputeConfig() ,
                                            block->get_registered_solvers()
                                            .size() ) , fn );
 }

/*--------------------------------------------------------------------------*/
// the same for a BlockSolverConfig, which is also clear()-ed for the final
// cleanup, see common_utils.h

void s_config_Block( Block * block , Configuration * s_config ,
                     const std::string & fn ,
                     bool clear_after )
{
 // handle the special case of a "meta" BlockSolverConfig: the same, father
 // first, as the Solver of a Block finds those of its sub-Block when it is
 // compute()-d, not when it is registered
 if( auto * mb =
     dynamic_cast< SimpleConfiguration< std::map< std::string ,
                                                  Configuration * > >
                                        * >( s_config ) ) {
  for_each_by_classname( block , mb->f_value ,
			 [ & ]( Block * b , Configuration * c ) {
			  auto bsc = dynamic_cast< BlockSolverConfig * >( c );
			  if( ! bsc ) {
			   std::cerr << "Error: meta-Configuration for :Block "
				     << b->classname() << " in file " << fn
				     << " is not a BlockSolverConfig"
				     << std::endl;
			   exit( 1 );
			   }
			  // only the Solver of the Block itself are
			  // positional, not those of its sub-Block
			  apply_present( bsc , b , fn , b == block );
			  } );

  // finally, clear() all the BlockSolverConfig for final cleanup
  if( clear_after )
   mb->clear();

  return;  // all done
  }

 if( auto * bsc = dynamic_cast< BlockSolverConfig * >( s_config ) ) {
  apply_present( bsc , block , fn , true );  // just apply() it
  if( clear_after )
   bsc->clear();                // clear() it for final cleanup
  return;                       // all done
  }

 std::cerr << "Error: " << fn
           << " does not contain a valid [meta]BlockSolverConfig"
           << std::endl;
 exit( 1 );

 }  // end( s_config_Block )

/*--------------------------------------------------------------------------*/
/*------------------------------ CROSS-CHECK -------------------------------*/
/*--------------------------------------------------------------------------*/

// where a Solver publishes its answer

double get_obj_value( Solver * slvr , ObjGetter g )
{
 switch( g ) {
  case ObjGetter::VarValue:   return( slvr->get_var_value() );
  case ObjGetter::LowerBound: return( slvr->get_lb() );
  case ObjGetter::UpperBound: return( slvr->get_ub() );
  }
 return( std::numeric_limits< double >::quiet_NaN() );  // unreachable
 }

/*--------------------------------------------------------------------------*/
// the canonical format of every value the tests print

std::string fmt_obj( double v )
{
 std::ostringstream os;
 os.setf( std::ios::scientific , std::ios::floatfield );
 os << std::setprecision( 7 ) << v;
 return( os.str() );
 }

/*--------------------------------------------------------------------------*/
// the tolerance Solver k is held to

double eps_of( std::size_t k , Solver * s , double dflt )
{
 // a single value of -E is the tolerance of every Solver, a list is
 // positional and its empty fields declare nothing
 double e = std::numeric_limits< double >::quiet_NaN();
 if( solver_eps.size() == 1 )
  e = solver_eps.front();
 else if( k < solver_eps.size() )
  e = solver_eps[ k ];

 if( ! std::isnan( e ) )
  return( e );

 // else the accuracy the Solver was asked for in its ComputeConfig
 return( s ? s->get_dbl_par( Solver::dblRelAcc ) : dflt );
 }

/*--------------------------------------------------------------------------*/
// the status the last SolveAll() got from each of its Solver

static std::vector< int > last_statuses;

int last_status( std::size_t k )
{
 return( k < last_statuses.size() ? last_statuses[ k ] : Solver::kUnEval );
 }

/*--------------------------------------------------------------------------*/
// whether Solver k was declared to be solving a relaxation

bool is_relaxation( std::size_t k )
{
 // a single value of -R applies to every Solver, a list is positional
 if( solver_relaxation.size() == 1 )
  return( solver_relaxation.front() );

 return( ( k < solver_relaxation.size() ) && solver_relaxation[ k ] );
 }

/*--------------------------------------------------------------------------*/
// the default reading of a Solver: the interval of its base contract

SolverReading read_bounds( Solver * s , std::size_t k )
{
 return( SolverReading{ get_obj_value( s , ObjGetter::LowerBound ) ,
                        get_obj_value( s , ObjGetter::UpperBound ) ,
                        eps_of( k , s ) } );
 }

/*--------------------------------------------------------------------------*/
// the reading of a Solver that publishes an optimum and no interval

SolverClassifier exact_getter( ObjGetter g )
{
 return( [ g ]( Solver * s , std::size_t k ) -> SolverReading {
  return( SolverReading::exact( get_obj_value( s , g ) , eps_of( k , s ) ) );
  } );
 }

/*--------------------------------------------------------------------------*/

// whether the extended output is on, either via -v or via the `verbose`
// environment variable (`verbose=1 ./batch ...`, `verbose=1 ctest ...`),
// which the tests that parse their arguments positionally, and so do not
// understand -v, inherit as well
static bool tests_verbose()
{
 static const bool env_on = []() {
   const char * e = std::getenv( "verbose" );
   return e && ( e[ 0 ] != '\0' ) && ( std::string( e ) != "0" );
   }();
 return ( verbosity_level >= 1 ) || env_on;
}

/*--------------------------------------------------------------------------*/
// the ColVariable of a Block the father Objective is built over

void collect_vars( Block * b , const std::vector< std::string > & groups ,
                   std::vector< ColVariable * > & vars )
{
 if( ! groups.empty() ) {
  for( const auto & name : groups ) {
   auto grp = b->get_static_variable_v< ColVariable >( name );
   if( ! grp )
    throw( std::invalid_argument( "collect_vars: no variable group named '" +
                                  name + "' in the sub-Block" ) );
   for( auto & v : *grp )
    vars.push_back( & v );
   }
  return;
  }

 auto obj = dynamic_cast< FRealObjective * >( b->get_objective() );
 if( ! obj )
  throw( std::invalid_argument(
                            "collect_vars: child has no FRealObjective" ) );
 auto f = obj->get_function();
 const Block::Index n = f->get_num_active_var();
 for( Block::Index i = 0 ; i < n ; ++i )
  vars.push_back( static_cast< ColVariable * >( f->get_active_var( i ) ) );
 }

/*--------------------------------------------------------------------------*/
// a father AbstractBlock with k copies of the Block a file gives

AbstractBlock * build_father( const std::string & filename , int k ,
                              const std::string & bconf ,
                              const std::vector< std::string > & groups ,
                              std::vector< ColVariable * > & vars )
{
 auto father = new AbstractBlock();
 vars.clear();
 for( int j = 0 ; j < k ; ++j ) {
  Block * child = Block::deserialize( filename , father );
  if( ! child )
   throw( std::invalid_argument( "build_father: cannot read Block from " +
                                 filename ) );
  if( ! bconf.empty() ) {
   Configuration * bc = Configuration::deserialize( bconf );
   b_config_Block( child , bc , bconf );
   delete bc;
   }
  child->generate_abstract_variables();
  child->generate_abstract_constraints();
  child->generate_objective();
  father->add_nested_Block( child );
  collect_vars( child , groups , vars );
  }
 return( father );
 }

/*--------------------------------------------------------------------------*/
// random data of a convex PolyhedralFunction

void generate_poly( Block::Index nv , int poly_rows , double scale ,
                    std::mt19937 & rg ,
                    PolyhedralFunction::MultiVector & A ,
                    PolyhedralFunction::RealVector & b )
{
 using Index = Block::Index;
 Index nr = poly_rows > 0 ? Index( poly_rows ) : nv + 1;
 A.assign( nr , PolyhedralFunction::RealVector( nv ) );
 b.assign( nr , 0 );
 for( Index r = 0 ; r < nr ; ++r ) {
  for( Index i = 0 ; i < nv ; ++i )
   A[ r ][ i ] = scale * rnd( rg );
  b[ r ] = scale * nv * rnd( rg ) / 4;
  }
 }

/*--------------------------------------------------------------------------*/
// a random Function over the given Variable, for the father Objective

Function * make_father_objective( std::vector< ColVariable * > & vars ,
                                  int obj_type , double scale ,
                                  int poly_rows , std::mt19937 & rg )
{
 using Index = Block::Index;
 using Coefficient = DQuadFunction::Coefficient;
 const Index nv = Index( vars.size() );

 // DQuadFunction: sum_i ( a_i x_i^2 + b_i x_i ), a_i > 0
 if( obj_type == 0 ) {
  DQuadFunction::v_coeff_triple tr( nv );
  for( Index i = 0 ; i < nv ; ++i )
   tr[ i ] = std::make_tuple( vars[ i ] , Coefficient( scale * rnd( rg ) ) ,
                              Coefficient( scale * ( 0.5 + pos( rg ) ) ) );
  return( new DQuadFunction( std::move( tr ) ) );
  }

 // QuadFunction: off-diagonal terms (i+1,i), kept PSD by Gershgorin diagonal
 // dominance ( 2 a_i >= sum |q| )
 if( obj_type == 1 ) {
  QuadFunction::v_off_diag_term nd;
  std::vector< double > rowabs( nv , 0.0 );
  for( Index i = 0 ; i + 1 < nv ; ++i ) {
   double q = scale * 0.3 * rnd( rg );
   nd.push_back( std::make_tuple( i + 1 , i , Coefficient( q ) ) );
   rowabs[ i ]     += std::abs( q );
   rowabs[ i + 1 ] += std::abs( q );
   }
  DQuadFunction::v_coeff_triple tr( nv );
  for( Index i = 0 ; i < nv ; ++i )
   tr[ i ] = std::make_tuple( vars[ i ] , Coefficient( scale * rnd( rg ) ) ,
                              Coefficient( 0.5 * rowabs[ i ] +
                                           scale * ( 0.5 + pos( rg ) ) ) );
  return( new QuadFunction( std::move( tr ) , std::move( nd ) ) );
  }

 // obj_type == 2: convex PolyhedralFunction = max_r ( A_r . x + b_r )
 PolyhedralFunction::MultiVector A;
 PolyhedralFunction::RealVector b;
 generate_poly( nv , poly_rows , scale , rg , A , b );
 PolyhedralFunction::VarVector vv( vars.begin() , vars.end() );
 auto pf = new PolyhedralFunction( std::move( vv ) , std::move( A ) ,
                                   std::move( b ) ,
                                   - Inf< Function::FunctionValue >() );
 pf->set_is_convex( true , eNoMod );
 return( pf );
 }

/*--------------------------------------------------------------------------*/
// make every Solver of the Block log at the verbosity the -v option asks for

void apply_solver_verbosity( Block * b )
{
 if( verbosity_level <= 0 )
  return;

 for( auto s : b->get_registered_solvers() ) {
  s->set_log( & std::cout );
  s->set_par( Solver::intLogVerb , verbosity_level );
  }
 }

/*--------------------------------------------------------------------------*/
// print the one line that reports an instance: timings, Solver values,
// reference and verdict

void print_instance_line( const std::vector< double > & times ,
                          const std::vector< std::string > & value_tokens ,
                          double ref ,
                          const std::string & verdict ,
                          double diff ,
                          bool always ,
                          const std::vector< std::string > & names )
{
 // the detailed per-round line (times, solver values, verdict) of a test
 // that re-solves in a loop of modification rounds is "extended" output:
 // print it only when verbose, keeping the default output terse. The
 // one-per-instance cross-check line of SolveAll() (@p always) and any
 // failing comparison (KO) are always shown instead, so that what was
 // compared, and a failure, are visible without re-running in verbose mode
 if( ( ! always ) && ( ! tests_verbose() ) &&
     ( verdict.compare( 0 , 2 , "KO" ) != 0 ) )
  return;

 // with the names of the Solver, one line each, the values aligned one
 // under the other so that two that do not agree are seen at a glance
 if( ! names.empty() ) {
  std::size_t w = 3;   // "Ref"
  for( const auto & n : names )
   w = std::max( w , n.size() );

  // and the times as well, one under the other, whatever the values: the
  // values are padded to the widest of them (the reference included)
  const std::string reftok = std::isnan( ref ) ? "" : fmt_obj( ref );
  std::size_t v = reftok.size();
  for( const auto & t : value_tokens )
   v = std::max( v , t.size() );

  // the times right-aligned on the decimal point, all having the digits of
  // fixd()
  std::vector< std::string > tt( times.size() );
  std::size_t u = 0;
  for( std::size_t k = 0 ; k < times.size() ; ++k ) {
   std::ostringstream os;
   os << fixd << times[ k ];
   tt[ k ] = os.str();
   u = std::max( u , tt[ k ].size() );
   }

  // the block starts on a line of its own, since whatever the test printed
  // before it (e.g., "First call: ") would shift its first line
  std::cout << std::endl;

  for( std::size_t k = 0 ; k < value_tokens.size() ; ++k ) {
   std::cout << "  " << std::left << std::setw( int( w ) )
             << ( k < names.size() ? names[ k ] : std::string() )
             << " = " << std::setw( int( v ) ) << value_tokens[ k ]
             << std::right;
   if( k < times.size() )
    std::cout << "   " << std::setw( int( u ) ) << tt[ k ] << " s";
   std::cout << std::endl;
   }

  if( ! std::isnan( ref ) ) {
   std::cout << "  " << std::left << std::setw( int( w ) ) << "Ref"
             << " = " << std::setw( int( v ) ) << reftok << std::right;
   if( ! std::isnan( diff ) )
    std::cout << "   (|diff| = " << fmt_obj( diff ) << ")";
   std::cout << std::endl;
   }

  std::cout << "  -> " << verdict << std::endl;
  return;
  }

 // without them, the Solver are numbered and the report is one line
 for( std::size_t k = 0 ; k < times.size() ; ++k )
  std::cout << ( k ? " - " : "" ) << fixd << times[ k ];

 std::cout << " | ";
 for( std::size_t k = 0 ; k < value_tokens.size() ; ++k )
  std::cout << ( k ? "  " : "" ) << "S" << k << " = " << value_tokens[ k ];

 if( ! std::isnan( ref ) ) {
  std::cout << "  ~ Ref = " << fmt_obj( ref );
  if( ! std::isnan( diff ) )
   std::cout << " (|diff| = " << fmt_obj( diff ) << ")";
  }

 std::cout << "  -> " << verdict << std::endl;
 }

/*--------------------------------------------------------------------------*/
// how a reading appears in that line

std::string reading_token( const SolverReading & r )
{
 const std::string lb = fmt_obj( r.lb ) , ub = fmt_obj( r.ub );

 // the interval of a relaxation is marked, since it is an interval around
 // the optimum of another problem [see SolverReading::Valid]
 const std::string tag =
  ( r.valid == SolverReading::kBoth ) ? "" : " (rel)";

 // a Solver that closed the gap prints the optimum rather than an interval
 // whose two ends read the same
 if( lb == ub )
  return( lb + tag );

 return( "[ " + lb + " , " + ub + " ]" + tag );
 }

/*--------------------------------------------------------------------------*/
// A cross-check where every Solver reports infeasible, or every Solver
// reports unbounded, is a pass: they do agree. A run made mostly of such
// comparisons, however, passes without ever comparing a solution, so they
// are counted here and a warning is issued at exit if they dominate; no
// verdict and no exit code is touched, since a test can be legitimately
// infeasible.

namespace {
 struct MutualInfWatchdog {
  std::size_t n_total = 0;  // multi-Solver cross-checks performed
  std::size_t n_inf   = 0;  // ... of which unanimously infeasible (OK(e))
  std::size_t n_unb   = 0;  // ... of which unanimously unbounded  (OK(u))
  ~MutualInfWatchdog() {
   if( ( n_total >= 4 ) && ( 2 * ( n_inf + n_unb ) > n_total ) )
    std::cerr << ANSI_YELLOW << "[WARNING] " << ( n_inf + n_unb )
              << " of " << n_total << " multi-Solver comparisons were "
                 "unanimously "
              << ( n_unb > n_inf ? "unbounded" : "infeasible" )
              << ": this run passed mostly via solver agreement on the "
                 "absence of a solution, not by comparing solutions."
              << ANSI_RESET << "\n";
   }
  };
 MutualInfWatchdog mutual_inf_watchdog;
}

/*--------------------------------------------------------------------------*/
// the verdict on one instance, see common_utils.h for the criterion

bool cross_check( const std::vector< SolverReading > & rd ,
                  const std::vector< bool > & has_solution ,
                  const std::vector< int > & status ,
                  double ref , double tol ,
                  std::string & verdict_out , double & diff_out )
{
 diff_out = std::numeric_limits< double >::quiet_NaN();

 // a <= b up to the relative tolerance t
 auto le = []( double a , double b , double t ) {
  return( a - b <=
          t * std::max( double( 1 ) ,
                        std::max( std::abs( a ) , std::abs( b ) ) ) );
  };

 const std::size_t M = has_solution.size();
 std::size_t nFeas = 0 , nInf = 0 , nUnb = 0;
 for( std::size_t k = 0 ; k < M ; ++k ) {
  if( has_solution[ k ] )                        ++nFeas;
  else if( status[ k ] == Solver::kInfeasible )  ++nInf;
  else if( status[ k ] == Solver::kUnbounded )   ++nUnb;
  }

 // count genuine multi-Solver comparisons for the mutual-infeasibility
 // watchdog (see above)
 if( M >= 2 )
  ++mutual_inf_watchdog.n_total;

 /* A Solver that says it has a solution has to say what it is worth, and a
  * value that is not a number is not a bound. The check has to be explicit:
  * every comparison against a NaN is false, so it does not fail any of the
  * agreement tests below, and with a single Solver and no reference there is
  * nothing it could be compared with at all. */

 for( std::size_t k = 0 ; k < M ; ++k )
  if( has_solution[ k ] &&
      ( std::isnan( rd[ k ].lb ) || std::isnan( rd[ k ].ub ) ) ) {
   verdict_out = "KO";
   return( false );
   }

 // single Solver, no reference: just "did it find a solution?" (with the
 // lb <= ub sanity check when the reading is a bracket)
 if( ( M == 1 ) && std::isnan( ref ) ) {
  bool ok = has_solution[ 0 ];
  if( ok && std::isfinite( rd[ 0 ].lb ) && std::isfinite( rd[ 0 ].ub ) )
   ok = le( rd[ 0 ].lb , rd[ 0 ].ub , tol );
  verdict_out = ok ? "OK" : "KO";
  return( ok );
  }

 /* A Solver of a relaxation says nothing about whether the problem is
  * feasible: a lower bound of a minimization problem is a lower bound of
  * +INF as well, and the Lagrangian dual of an infeasible problem whose
  * subproblems are feasible grows without ever proving it. Nor does one
  * whose interval is the whole line, such as a heuristic that has found
  * nothing. Hence, when all the other Solver say infeasible, and these have
  * a one-sided bound or no bound at all, the infeasibility is unanimous. */
 std::size_t nOneSided = 0;
 for( std::size_t k = 0 ; k < M ; ++k )
  if( has_solution[ k ] &&
      ( ( rd[ k ].valid != SolverReading::kBoth ) ||
        ( ( rd[ k ].lb == - std::numeric_limits< double >::infinity() ) &&
          ( rd[ k ].ub == std::numeric_limits< double >::infinity() ) ) ) )
   ++nOneSided;
 if( ( nInf > 0 ) && ( nFeas == nOneSided ) && ( nInf + nFeas == M ) &&
     std::isnan( ref ) ) {
  if( M >= 2 ) ++mutual_inf_watchdog.n_inf;
  verdict_out = "OK(e)"; return( true );
  }

 // unanimous infeasible / unbounded is a pass only when NO reference was given
 if( nFeas == 0 ) {
  if( std::isnan( ref ) && nInf == M ) {
   if( M >= 2 ) ++mutual_inf_watchdog.n_inf;
   verdict_out = "OK(e)"; return( true );
   }
  if( std::isnan( ref ) && nUnb == M ) {
   if( M >= 2 ) ++mutual_inf_watchdog.n_unb;
   verdict_out = "OK(u)"; return( true );
   }
  verdict_out = "KO";
  return( false );
  }

 // some feasible, some not: disagreement (an errored Solver counts as
 // disagreeing with the feasible ones, too)
 if( nFeas < M ) { verdict_out = "KO"; return( false ); }

 // all feasible: the agreement check of common_utils.h. The optimum is not
 // known, so what every Solver is measured against is the best knowledge
 // the whole set of them provides, i.e. the largest of the lower bounds and
 // the smallest of the upper bounds (the reference value, when given, being
 // one more Solver that has both)
 constexpr double INF = std::numeric_limits< double >::infinity();
 double best_lb = - INF , best_ub = INF;
 for( std::size_t k = 0 ; k < M ; ++k ) {
  // the end of a relaxation's interval that bounds its own optimum and not
  // this one is no knowledge about this problem [see SolverReading::Valid]
  if( rd[ k ].valid != SolverReading::kUpper )
   best_lb = std::max( best_lb , rd[ k ].lb );
  if( rd[ k ].valid != SolverReading::kLower )
   best_ub = std::min( best_ub , rd[ k ].ub );
  }
 if( ! std::isnan( ref ) ) {
  best_lb = std::max( best_lb , ref );
  best_ub = std::min( best_ub , ref );
  }

 bool ok = true;
 double zstar = std::numeric_limits< double >::quiet_NaN();  // 1st claimed z*
 for( std::size_t k = 0 ; k < M ; ++k ) {
  const double lb = rd[ k ].lb , ub = rd[ k ].ub;
  const bool lb_is_ours = ( rd[ k ].valid != SolverReading::kUpper );
  const bool ub_is_ours = ( rd[ k ].valid != SolverReading::kLower );

  // correctness, which every Solver owes whatever it promised: the bounds it
  // gives on THIS problem cannot contradict the bounds of the others
  if( ub_is_ours && std::isfinite( ub ) && ( ! le( best_lb , ub , tol ) ) )
   ok = false;
  if( lb_is_ours && std::isfinite( lb ) && ( ! le( lb , best_ub , tol ) ) )
   ok = false;

  // quality, which only the Solver that delivered what it was asked owes:
  // each of its bounds has to be within its declared tolerance of the best
  // opposite one. A Solver that returns kLowPrecision, or that stopped on
  // any other condition, promised nothing and owes nothing here. The
  // tolerance is never taken below tol: under that the comparison would be
  // measuring the noise of the cross-check itself, and a Solver asked for
  // an accuracy that tight stops just outside it anyway, its own stopping
  // criterion not being normalized exactly like this one
  const double e = std::max( std::isnan( rd[ k ].eps ) ? tol : rd[ k ].eps ,
                             tol );
  if( ( status[ k ] == Solver::kOK ) && ( ! std::isinf( e ) ) ) {
   if( rd[ k ].valid == SolverReading::kBoth ) {
    if( std::isfinite( ub ) && std::isfinite( best_lb )
        && ( ! le( ub , best_lb , e ) ) )
     ok = false;
    if( std::isfinite( lb ) && std::isfinite( best_ub )
        && ( ! le( best_ub , lb , e ) ) )
     ok = false;
    }
   else
    // a relaxation is asked to solve exactly the problem it does solve, so
    // what its tolerance is about is the width of its own interval; against
    // this problem it owes correctness alone, since the distance between the
    // two optima is the gap of the relaxation and not an error of its
    if( std::isfinite( lb ) && std::isfinite( ub ) && ( ! le( ub , lb , e ) ) )
     ok = false;
   }

  // the optimum is pinned by the first Solver that bounds it on both sides,
  // which a relaxation never does whatever its own interval looks like
  if( std::isnan( zstar ) && ( rd[ k ].valid == SolverReading::kBoth )
      && std::isfinite( lb ) && std::isfinite( ub ) )
   zstar = rd[ k ].claimed();
  }

 if( ( ! std::isnan( zstar ) ) && ( ! std::isnan( ref ) ) )
  diff_out = std::abs( zstar - ref );

 verdict_out = ok ? ( std::isnan( zstar ) ? "OK" : "OK(f)" ) : "KO";
 return( ok );
 }

/*--------------------------------------------------------------------------*/
// solve an instance with every Solver, cross-check them, report

void print_solver_parameters( Block * block )
{
 if( ( verbosity_level < 2 ) || ( ! block ) )
  return;

 // one entry per distinct ( Solver class , Block class , parameters ), in
 // the order they are first met, with the number of Solvers sharing it
 struct Entry {
  std::string solver , owner , pars;
  std::size_t count;
  };
 std::vector< Entry > entries;

 std::function< void( Block * ) > visit = [ & ]( Block * b ) {
  for( auto solver : b->get_registered_solvers() ) {
   std::ostringstream pars;
   solver->print_parameters( pars );
   Entry e{ solver->classname() , b->classname() , pars.str() , 1 };
   auto it = std::find_if( entries.begin() , entries.end() ,
                           [ & ]( const Entry & o ) {
                            return( ( o.solver == e.solver ) &&
                                    ( o.owner == e.owner ) &&
                                    ( o.pars == e.pars ) ); } );
   if( it == entries.end() )
    entries.push_back( std::move( e ) );
   else
    ++it->count;
   }
  if( verbosity_level >= 3 )
   for( auto sb : b->get_nested_Blocks() )
    visit( sb );
  };
 visit( block );

 for( const auto & e : entries ) {
  std::cout << std::endl << "--- parameters of " << e.solver << " on "
            << e.owner;
  if( e.count > 1 )
   std::cout << " (" << e.count << " Solvers)";
  std::cout << std::endl << e.pars;
  }
 std::cout << std::endl;
 }

/*--------------------------------------------------------------------------*/

bool SolveAll( Block * block ,
               const SolverClassifier & classify ,
               double ref ,
               double tol ,
               double * out_fo1 ,
               bool   * out_hs1 ,
               double * out_time1 ,
               long   * out_it1 ,
               BlockSolverConfig * bsc ,
               Solution * start ,
               SolveRecord * record )
{
 constexpr double INF = std::numeric_limits< double >::has_infinity
                        ? std::numeric_limits< double >::infinity()
                        : std::numeric_limits< double >::max();

 try {
  // the Solver the BlockSolverConfig attached to this Block, in its order,
  // which is the order the -E tolerances are positional on; the Solver that
  // these in turn attach to the inner Block are theirs, not part of the
  // comparison [see BlockSolverConfig::get_Solvers()]. Without a
  // BlockSolverConfig the comparison is over all the Solver registered to
  // the Block
  std::vector< Solver * > S;
  if( bsc ) {
   const auto & mine = bsc->get_Solvers( block );
   S.assign( mine.begin() , mine.end() );
   }
  else {
   const auto & reg = block->get_registered_solvers();
   S.assign( reg.begin() , reg.end() );
   }
  const std::size_t M = S.size();
  if( M == 0 ) {
   std::cout << "no Solver registered to the Block!" << std::endl;
   return( false );
   }

  // with -v 2, before solving, print what each Solver was actually given
  print_solver_parameters( block );

  // solve every Solver, timing each, then read the feasible ones - - - - - - -
  std::vector< int >    status( M );
  std::vector< double > times( M );
  std::vector< long >   iters( M );
  std::vector< bool >   hs( M );
  std::vector< SolverReading > rd( M );
  std::vector< std::string > tok( M );
  for( std::size_t k = 0 ; k < M ; ++k ) {
   // the same starting point for every Solver, if one is given
   if( start )
    start->write( block );
   // a Solver that throws has failed on this instance, and is reported as
   // such in the line of the others rather than ending them all
   bool thrown = false;
   auto t_start = std::chrono::system_clock::now();
   try {
    status[ k ] = S[ k ]->compute( false );
    }
   catch( std::exception & e ) {
    std::cerr << S[ k ]->classname() << ": " << e.what() << std::endl;
    status[ k ] = Solver::kError;
    thrown = true;
    }
   auto t_end = std::chrono::system_clock::now();
   times[ k ] = std::chrono::duration< double >( t_end - t_start ).count();
   iters[ k ] = S[ k ]->get_elapsed_iterations();
   hs[ k ] = ( ( ( status[ k ] >= Solver::kOK )
                 && ( status[ k ] < Solver::kError )
                 && ( status[ k ] != Solver::kUnbounded )
                 && ( status[ k ] != Solver::kInfeasible ) )
               || ( status[ k ] == Solver::kLowPrecision ) );
   if( hs[ k ] ) {
    rd[ k ] = classify( S[ k ] , k );
    // a Solver declared a relaxation by -R bounds the optimum of the
    // relaxation: of this problem it bounds only the side the relaxation
    // is on, which the sense of the Objective decides (that of the Block,
    // which one without an Objective of its own, its objective being the
    // sum of these of its sub-Block, gives all the same)
    if( is_relaxation( k ) )
     rd[ k ].valid = ( block->get_objective_sense() == Objective::eMax )
                     ? SolverReading::kUpper : SolverReading::kLower;
    tok[ k ] = reading_token( rd[ k ] );
    // a Solver that did not return kOK did not deliver what it was asked
    // and is therefore only held to correctness: say so in the line, since
    // it is what its numbers are worth
    if( status[ k ] == Solver::kLowPrecision )  tok[ k ] += " (lowP)";
    else if( status[ k ] != Solver::kOK )       tok[ k ] += " (stop)";
    }
   else if( status[ k ] == Solver::kInfeasible )  tok[ k ] = "Unfeas";
   else if( status[ k ] == Solver::kUnbounded )   tok[ k ] = "Unbounded";
   else if( thrown )                              tok[ k ] = "Error! (throw)";
   else                                           tok[ k ] = error_token(
							      status[ k ] );
   }

  last_statuses = status;

  // out-params from the first Solver - - - - - - - - - - - - - - - - - - - -
  if( out_fo1 )   *out_fo1   = hs[ 0 ] ? rd[ 0 ].claimed() : -INF;
  if( out_hs1 )   *out_hs1   = hs[ 0 ];
  if( out_time1 ) *out_time1 = times[ 0 ];
  if( out_it1 )   *out_it1   = iters[ 0 ];

  // cross-check + uniform per-instance line - - - - - - - - - - - - - - - - -
  std::string verdict;
  double diff;
  bool ok = cross_check( rd , hs , status , ref , tol , verdict , diff );

  // what each Solver is called, so that the report says who returned what
  std::vector< std::string > names( M );
  for( std::size_t k = 0 ; k < M ; ++k )
   names[ k ] = S[ k ]->classname();

  print_instance_line( times , tok , ref , verdict , diff , true , names );

  if( record )
   *record = { rd , hs , status , times , tok , names };

  return( ok );
  }
 catch( std::exception & e ) {
  std::cerr << e.what() << std::endl;
  exit( 1 );
  }
 catch( ... ) {
  std::cerr << "Error: unknown exception thrown" << std::endl;
  exit( 1 );
  }
 }

/*--------------------------------------------------------------------------*/

bool SolveAll( Block * block ,
               double ref ,
               double tol ,
               double * out_fo1 ,
               bool   * out_hs1 ,
               double * out_time1 ,
               long   * out_it1 ,
               BlockSolverConfig * bsc )
{
 return( SolveAll( block , read_bounds , ref , tol ,
                   out_fo1 , out_hs1 , out_time1 , out_it1 , bsc ) );
 }

/*--------------------------------------------------------------------------*/
// SolveAll() of each copy of a Block, then the copies against each other

bool SolveFormulations( const std::vector< Block * > & blocks ,
                        const std::vector< std::string > & names ,
                        const SolverClassifier & classify ,
                        double ref , double tol ,
                        const std::vector< BlockSolverConfig * > & bscs ,
                        const std::vector< std::vector< bool > > & relax )
{
 auto bsc_of = [ & bscs ]( std::size_t k ) {
  return( k < bscs.size() ? bscs[ k ] : nullptr );
  };

 // SolveAll() of the copy k, with its own -R if it has one
 auto solve = [ & ]( std::size_t k , SolveRecord * rec ) {
  const auto global = solver_relaxation;
  if( k < relax.size() )
   solver_relaxation = relax[ k ];
  const bool ok = SolveAll( blocks[ k ] , classify , ref , tol , nullptr ,
                            nullptr , nullptr , nullptr , bsc_of( k ) ,
                            nullptr , rec );
  solver_relaxation = global;
  return( ok );
  };

 if( blocks.size() == 1 )
  return( solve( 0 , nullptr ) );

 bool ok = true;
 SolveRecord all;
 for( std::size_t k = 0 ; k < blocks.size() ; ++k ) {
  std::cout << "formulation " << names[ k ] << ":";
  SolveRecord rec;
  ok &= solve( k , & rec );
  all.rd.insert( all.rd.end() , rec.rd.begin() , rec.rd.end() );
  all.hs.insert( all.hs.end() , rec.hs.begin() , rec.hs.end() );
  all.status.insert( all.status.end() , rec.status.begin() ,
                     rec.status.end() );
  all.times.insert( all.times.end() , rec.times.begin() , rec.times.end() );
  all.tok.insert( all.tok.end() , rec.tok.begin() , rec.tok.end() );
  for( const auto & n : rec.names )
   all.names.push_back( names[ k ] + " " + n );
  }

 // the Solver of all the copies against each other: they bound the same
 // optimum, whatever the formulation of the copy they are attached to
 std::string verdict;
 double diff;
 const bool cross = cross_check( all.rd , all.hs , all.status , ref , tol ,
                                 verdict , diff );
 std::cout << "formulations against each other:";
 print_instance_line( all.times , all.tok , ref , verdict , diff , true ,
                      all.names );

 return( ok && cross );
 }

/*--------------------------------------------------------------------------*/
// the options of the rounds of random changes

RoundOptions round_opts;

void add_round_options( const std::string & wchg_help )
{
 const auto & r = round_opts;
 short_opts += "e:k:n:m:q:";
 const std::vector< option > my_opts = {
   { "seed"   , required_argument , nullptr , 'e' } ,
   { "wchg"   , required_argument , nullptr , 'k' } ,
   { "rounds" , required_argument , nullptr , 'n' } ,
   { "nchng"  , required_argument , nullptr , 'm' } ,
   { "pchng"  , required_argument , nullptr , 'q' } };
 long_opts.insert( std::prev( long_opts.end() ) ,
                   my_opts.begin() , my_opts.end() );
 std::ostringstream p;
 p << r.p_change;
 help += "  -e, --seed <n>                  seed of the random generator [" +
         std::to_string( r.seed ) + "]\n"
         "  -k, --wchg <bits>               what to change, bit-wise [" +
         std::to_string( r.wchg ) + "]:\n" + wchg_help +
         "  -n, --rounds <n>                number of rounds of changes [" +
         std::to_string( r.n_repeat ) + "]\n"
         "  -m, --nchng <n>                 number of changes [" +
         std::to_string( r.n_change ) + "]\n"
         "  -q, --pchng <p>                 probability of a change [" +
         p.str() + "]\n";
 }

/*--------------------------------------------------------------------------*/

bool process_round_option( int opt )
{
 auto & r = round_opts;
 switch( opt ) {
  case( 'e' ): Str2Sthg( optarg , r.seed );     return( true );
  case( 'k' ): Str2Sthg( optarg , r.wchg );     return( true );
  case( 'n' ): Str2Sthg( optarg , r.n_repeat ); return( true );
  case( 'm' ): Str2Sthg( optarg , r.n_change ); return( true );
  case( 'q' ): Str2Sthg( optarg , r.p_change ); return( true );
  default:                                      return( false );
  }
 }

/*--------------------------------------------------------------------------*/
// the main of a tester that loads a Block out of a file and cross-checks
// its Solver

int load_solve_test( int argc , char ** argv , const std::string & desc ,
                     const std::string & classname , bool any_class ,
                     double tol )
{
 std::set_terminate( smspp_terminate );

 static std::string block_class;
 block_class = classname;
 static bool b_option;
 b_option = any_class;

 docopt_desc = desc;
 short_opts += any_class ? "b:r:" : "r:";
 std::vector< option > my_opts;
 if( any_class )
  my_opts.push_back( { "block" , required_argument , nullptr , 'b' } );
 my_opts.push_back( { "ref" , required_argument , nullptr , 'r' } );
 long_opts.insert( std::prev( long_opts.end() ) ,
                   my_opts.begin() , my_opts.end() );
 if( any_class )
  help += "  -b, --block <classname>         the Block the instance is "
          "loaded into [" + classname + "]\n";
 help += "  -r, --ref <value>               the optimum of the instance "
         "[none]\n";

 process_args( argc , argv , []( int opt ) {
   switch( opt ) {
    case( 'b' ):
     if( ! b_option )
      return( false );
     block_class = optarg;
     return( true );
    case( 'r' ): Str2Sthg( optarg , RefObjective ); return( true );
    default:                                        return( false );
    }
   } );

 auto forms = get_formulations( true );

 std::vector< Block * > blocks;
 std::vector< Configuration * > bscs;
 std::vector< BlockSolverConfig * > bsolvers;
 std::vector< std::vector< bool > > relax;
 auto clean = [ & ]( void ) {
  for( std::size_t k = 0 ; k < blocks.size() ; ++k ) {
   // detach the Solver by applying the cleared BlockSolverConfig, then free
   bscs[ k ]->clear();
   s_config_Block( blocks[ k ] , bscs[ k ] , sconf_file );
   delete bscs[ k ];
   delete blocks[ k ];
   }
  delete forms.owner;
  };

 for( std::size_t k = 0 ; k < forms.configs.size() ; ++k ) {
  auto block = Block::new_Block( block_class );
  if( ! block ) {
   std::cerr << "Error: " << block_class << " not present in Block factory"
             << std::endl;
   exit( 1 );
   }

  block->load( filename );

  if( forms.configs[ k ] )
   b_config_Block( block , forms.configs[ k ] , forms.file );

  auto bsc = solver_config_of( forms , k );
  relax.push_back( attach_formulation_Solvers( block , bsc , forms , k ,
                                               false ) );
  blocks.push_back( block );
  bscs.push_back( bsc );
  bsolvers.push_back( dynamic_cast< BlockSolverConfig * >( bsc ) );
  if( block->get_registered_solvers().empty() ) {
   // none of the Solver it names is in this build: nothing to do
   clean();
   return( 77 );
   }
  }

 const bool ok = SolveFormulations( blocks , forms.names , read_bounds ,
                                    RefObjective , tol , bsolvers , relax );
 clean();
 return( ok ? 0 : 1 );
 }

/*--------------------------------------------------------------------------*/
// the (meta-)BlockConfig of a run: the list of -F, or the one of -B

Formulations get_formulations( bool optional )
{
 Formulations v;

 if( optional && fconf_file.empty() && bconf_file.empty() ) {
  v.configs = { nullptr };
  v.solvers = { nullptr };
  v.names = { "" };
  return( v );
  }

 if( fconf_file.empty() == bconf_file.empty() ) {
  std::cerr << "Error: give either a BlockConfig (-B) or a list of them (-F)"
            << std::endl;
  exit( 1 );
  }

 if( ! fconf_file.empty() ) {
  auto vc = read_formulations( fconf_file );
  v.owner = vc;
  for( std::size_t k = 0 ; k < vc->f_value.size() ; ++k ) {
   auto e = vc->f_value[ k ];
   if( auto p = dynamic_cast< SimpleConfiguration< std::pair<
                 Configuration * , Configuration * > > * >( e ) ) {
    v.configs.push_back( p->f_value.first );
    v.solvers.push_back( p->f_value.second );
    }
   else
    if( auto q = dynamic_cast< SimpleConfiguration< std::vector<
                                      Configuration * > > * >( e ) ) {
     // the BlockConfig, the BlockSolverConfig and which of its Solver solve
     // a relaxation
     const auto & w = q->f_value;
     auto r = w.size() > 2 ? dynamic_cast< SimpleConfiguration<
                                  std::vector< int > > * >( w[ 2 ] ) : nullptr;
     if( ( w.size() < 2 ) || ( w.size() > 3 ) ||
         ( ( w.size() == 3 ) && ( ! r ) ) ) {
      std::cerr << "Error: formulation " << k << " of " << fconf_file
                << " is not a BlockConfig, a BlockSolverConfig and the"
                << " SimpleConfiguration< std::vector< int > > of its -R"
                << std::endl;
      exit( 1 );
      }
     v.configs.push_back( w[ 0 ] );
     v.solvers.push_back( w[ 1 ] );
     if( r ) {
      v.relax.resize( k + 1 );
      v.relax[ k ] = r->f_value;
      }
     }
    else {
     v.configs.push_back( e );
     v.solvers.push_back( nullptr );
     }
   v.names.push_back( "F" + std::to_string( k ) );
   }
  v.file = fconf_file;
  return( v );
  }

 v.owner = Configuration::deserialize( bconf_file );
 if( ! v.owner ) {
  std::cerr << "Error: cannot load BlockConfig from " << bconf_file
            << std::endl;
  exit( 1 );
  }
 v.configs = { v.owner };
 v.solvers = { nullptr };
 v.names = { bconf_file };
 v.file = bconf_file;
 return( v );
 }

/*--------------------------------------------------------------------------*/

std::size_t name_Solver( const std::string & option , std::size_t k )
{
 named_solvers.push_back( { option , k , k } );
 return( named_solvers.size() - 1 );
 }

/*--------------------------------------------------------------------------*/

std::size_t named_Solver( std::size_t h )
{
 return( h < named_solvers.size() ? named_solvers[ h ].at
                                  : std::string::npos );
 }

/*--------------------------------------------------------------------------*/
// -E as given, which each copy of a Block filters anew

static std::vector< double > solver_eps_given;

/*--------------------------------------------------------------------------*/
// the BlockSolverConfig of a formulation applied, with its -R

std::vector< bool > attach_formulation_Solvers( Block * block ,
                                                Configuration * bsc ,
                                                const Formulations & v ,
                                                std::size_t k ,
                                                bool clear_after )
{
 // the -R of the formulation goes through s_config_Block(), which takes out of
 // it the Solver that are not there, as it does of -R; -E and the Solver
 // named by an option are filtered anew from what was given, the previous
 // copy having filtered them already, and left so for the solve, every copy
 // having the same Solver
 const auto global = solver_relaxation;
 const bool own = ( k < v.relax.size() ) && ( ! v.relax[ k ].empty() );
 if( own )
  solver_relaxation.assign( v.relax[ k ].begin() , v.relax[ k ].end() );
 solver_eps = solver_eps_given;
 for( auto & n : named_solvers )
  n.at = n.given;
 const bool mine_bsc = ( k < v.solvers.size() ) && v.solvers[ k ];
 s_config_Block( block , bsc , mine_bsc ? v.file : sconf_file ,
                 clear_after );
 auto mine = solver_relaxation;
 solver_relaxation = global;
 return( mine );
 }

/*--------------------------------------------------------------------------*/
// the (meta-)BlockSolverConfig of a formulation

Configuration * solver_config_of( const Formulations & v ,
                                  std::size_t k )
{
 if( ( k < v.solvers.size() ) && v.solvers[ k ] )
  return( v.solvers[ k ]->clone() );

 if( sconf_file.empty() ) {
  std::cerr << "Error: formulation " << v.names[ k ] << " of " << v.file
            << " has no BlockSolverConfig and none is given with -S"
            << std::endl;
  exit( 1 );
  }

 auto c = Configuration::deserialize( sconf_file );
 if( ! c ) {
  std::cerr << "Error: cannot load the BlockSolverConfig " << sconf_file
            << std::endl;
  exit( 1 );
  }
 return( c );
 }

/*--------------------------------------------------------------------------*/
// SolveAll() of one or two Solver read via the given getters

bool SolveBoth( Block * block ,
                ObjGetter g1 ,
                ObjGetter g2 ,
                bool one_sided_le ,
                double tol ,
                double * out_fo1 ,
                bool   * out_hs1 ,
                double * out_time1 ,
                long   * out_it1 )
{
 // SolveBoth is the M <= 2 special case of SolveAll: the first Solver is read
 // via g1, the second via g2. In the default (two-sided) mode both claim an
 // optimum that must agree; in the ProxHeur one-sided mode the first is a
 // lower bound and the second an upper bound, so the verdict becomes LB <= UB.
 SolverClassifier classify =
  [ g1 , g2 , one_sided_le ]( Solver * s , std::size_t k ) -> SolverReading {
   const double v = get_obj_value( s , k ? g2 : g1 );
   if( ! one_sided_le )
    return( SolverReading::exact( v , eps_of( k , s ) ) );

   // the one-sided bounds claim nothing unless -E says how tight they are
   const double e = eps_of( k , nullptr ,
                            std::numeric_limits< double >::infinity() );
   return( k ? SolverReading::upper_bound( v , e )
             : SolverReading::lower_bound( v , e ) );
   };

 return( SolveAll( block , classify ,
                   std::numeric_limits< double >::quiet_NaN() , tol ,
                   out_fo1 , out_hs1 , out_time1 , out_it1 ) );
 }

/*--------------------------------------------------------------------------*/

double tree_objective_value( Block * block )
{
 double value = 0;

 if( auto obj = dynamic_cast< RealObjective * >( block->get_objective() ) ) {
  obj->compute();
  value += obj->value();
  }

 for( auto sb : block->get_nested_Blocks() )
  value += tree_objective_value( sb );

 return( value );

 }  // end( tree_objective_value )

/*--------------------------------------------------------------------------*/

bool check_var_solutions( Block * block , double tol )
{
 bool allok = true;
 std::size_t k = 0;

 for( auto slvr : block->get_registered_solvers() ) {
  const std::size_t h = k++;

  if( is_relaxation( h ) )
   /* A Solver that solves a relaxation reports the value of the
    * relaxation, while what it writes in the Variable is a point of the
    * relaxed problem: the two are the same number only if the objective
    * is linear on it. A Lagrangian dual is the obvious case, its point
    * being a convex combination of the solutions of the sub-problems: as
    * soon as the objective is convex and not linear, the combination
    * costs strictly less than the combination of the costs, and the
    * difference is not an error of anybody's. Hence the check does not
    * apply here. */
   continue;

  if( ! slvr->has_var_solution() )
   continue;                    // nothing to read, hence nothing to check

  const double reported = slvr->get_var_value();
  if( ! std::isfinite( reported ) )
   continue;                    // nothing to compare the solution with

  double value;
  try {
   slvr->get_var_solution();
   value = tree_objective_value( block );
   }
  catch( std::exception & e ) {
   std::cerr << "Error: Solver " << h << " has a solution it cannot give: "
             << e.what() << std::endl;
   allok = false;
   continue;
   }

  if( verbosity_level > 0 )
   std::cout << "solution of Solver " << h << " is worth "
             << fmt_obj( value ) << ", reported "
             << fmt_obj( reported ) << std::endl;

  if( std::abs( value - reported ) >
      tol * std::max( { 1.0 , std::abs( value ) , std::abs( reported ) } ) ) {
   std::cerr << "Error: the solution of Solver " << h << " is worth "
             << fmt_obj( value ) << ", but it reports "
             << fmt_obj( reported ) << std::endl;
   allok = false;
   }
  }

 return( allok );

 }  // end( check_var_solutions )

/*--------------------------------------------------------------------------*/

double own_rows_violation( Block * block )
{
 double viol = 0;

 auto see = [ & viol ]( FRowConstraint & cnst ) {
  if( cnst.is_relaxed() )
   return;
  if( auto ret = cnst.compute() ;
      ( ret <= FRowConstraint::kUnEval ) || ( ret > FRowConstraint::kOK ) ) {
   viol = Inf< double >();
   return;
   }
  viol = std::max( viol , double( cnst.rel_viol() ) );
  };

 block->for_each_constraint_group( [ & see ]( const BaseGroup & group ) {
   group.for_each_as< FRowConstraint >( see ); } );

 return( viol );

 }  // end( own_rows_violation )

/*--------------------------------------------------------------------------*/

bool check_relaxation_solutions( Block * block , double tol , double ref ,
                                 double ref_tol )
{
 bool allok = true;
 std::size_t k = 0;

 for( auto slvr : block->get_registered_solvers() ) {
  const std::size_t h = k++;

  if( ! is_relaxation( h ) )
   continue;                    // check_var_solutions() covers these

  if( ! slvr->has_var_solution() )
   continue;                    // nothing to read, hence nothing to check

  /* A Solver that stopped before converging (a time limit, a low precision,
   * i.e. anything but kOK) has not driven the residual to zero, and the
   * combination it has built is only as good as where it stopped: it is held
   * to correctness alone, as in the cross-check [see SolveAll()]. */
  if( last_status( h ) != Solver::kOK )
   continue;

  /* What the reconstruction satisfies the dualised rows is the convex
   * combination of the answers of the components, and that is a point of the
   * original problem only where the relaxation is exact. Where it is not,
   * i.e. where the bound the relaxation reports is strictly better than the
   * optimum (a unit with a commitment is enough for that), the combination
   * violates the rows by the very gap, and holding it to them would call a
   * duality gap an error. */
  if( ( ! std::isnan( ref ) ) && ( std::abs( slvr->get_var_value() - ref ) >
                                   ref_tol * std::max( double( 1 ) ,
                                                       std::abs( ref ) ) ) ) {
   if( verbosity_level > 0 )
    std::cout << "relaxation of Solver " << h << " is not exact here, its "
              << "reconstruction is not held to the dualised rows"
              << std::endl;
   continue;
   }

  double viol;
  try {
   slvr->get_var_solution();
   viol = own_rows_violation( block );
   }
  catch( std::exception & e ) {
   std::cerr << "Error: Solver " << h << " has a solution it cannot give: "
             << e.what() << std::endl;
   allok = false;
   continue;
   }

  if( verbosity_level > 0 )
   std::cout << "relaxation of Solver " << h << " violates the dualised rows"
             << " by " << viol << std::endl;

  if( viol > tol ) {
   std::cerr << "Error: the solution Solver " << h << " reconstructs violates"
             << " the dualised rows by " << viol << " > " << tol << std::endl;
   allok = false;
   }
  }

 return( allok );

 }  // end( check_relaxation_solutions )

/*--------------------------------------------------------------------------*/
// compare a value against the reference one and report

bool CheckRefValue( double fo , double ref ,
                    double rel_tol ,
                    double time1 , long iters )
{
 double maxv = std::max( double( 1 ) ,
                         std::max( std::abs( fo ) , std::abs( ref ) ) );
 double diff = std::abs( fo - ref );
 double tol = rel_tol * maxv;

 bool OK = ( diff <= tol );

 std::cout << fixd << time1 << "\t" << iters << "\t"
           << def << fo
           << " ~ Ref = " << def << ref
           << " (|diff| = " << def << diff
           << ( OK ? ", OK" : ", KO" ) << ")" << std::endl;

 return( OK );
 }

/*--------------------------------------------------------------------------*/
// solve with the only Solver and compare against the reference value

bool SolveAndCheckRef( Block * block , double ref ,
                       ObjGetter g ,
                       double rel_tol )
{
 // single-Solver solve + reference check, expressed as SolveAll() of the only
 // registered Solver (read via g) against ref
 return( SolveAll( block , exact_getter( g ) , ref , rel_tol ) );
 }

/*--------------------------------------------------------------------------*/
/*------------------------------ COMMAND LINE ------------------------------*/
/*--------------------------------------------------------------------------*/
// the globals declared extern in common_utils.h

std::string docopt_desc;

std::string exe;
std::string filename;
std::string bconf_file;
std::string sconf_file;
std::string fconf_file;
std::string block_prefix;
std::string conf_prefix;

bool sol_verbose = false;
bool dryrun      = false;

int verbosity_level = 0;

double RefObjective = std::numeric_limits< double >::quiet_NaN();

std::vector< double > solver_eps;

std::vector< bool > solver_relaxation;

// the getopt baseline shared by the tests that opt in; those that need
// extra switches override short_opts / long_opts / help in their main()
// before calling process_args()
std::string short_opts = "B:F:S:p:c:E:R:Dv::h";

std::vector< option > long_opts = {
 { "help"            , no_argument       , nullptr , 'h' } ,
 { "blockcfg"        , required_argument , nullptr , 'B' } ,
 { "formulations"    , required_argument , nullptr , 'F' } ,
 { "solvercfg"       , required_argument , nullptr , 'S' } ,
 { "prefix"          , required_argument , nullptr , 'p' } ,
 { "configdir"       , required_argument , nullptr , 'c' } ,
 { "eps"             , required_argument , nullptr , 'E' } ,
 { "relaxation"      , required_argument , nullptr , 'R' } ,
 { "dryrun"          , no_argument       , nullptr , 'D' } ,
 { "verbose"         , optional_argument , nullptr , 'v' } ,
 { nullptr           , no_argument       , nullptr , 0   }
 };

std::string help =
 "  -h, --help                      print this help\n"
 "  -B, --blockcfg <file>           Block Configuration\n"
 "  -F, --formulations <file>       list of Block Configurations, in place\n"
 "                                  of -B: the Block is read once per\n"
 "                                  element, the copies are solved by the\n"
 "                                  same Solver and compared to each other\n"
 "  -S, --solvercfg <file>          Solver Configuration\n"
 "  -p, --prefix <path>             the prefix for all Block filenames\n"
 "  -c, --configdir <path>          the prefix for all Config filenames\n"
 "  -E, --eps <e[,e,...]>           optimality tolerance of each Solver, in\n"
 "                                  the order of the BlockSolverConfig, the\n"
 "                                  empty field being its own dblRelAcc; one\n"
 "                                  value applies to all, inf claims nothing\n"
 "                                  beyond a valid [ get_lb() , get_ub() ]\n"
 "  -R, --relaxation <r[,r,...]>    which Solver solve a relaxation of the\n"
 "                                  problem, r in the same positional order,\n"
 "                                  the empty field being an ordinary one:\n"
 "                                  a relaxation is held to solve its own\n"
 "                                  problem exactly, and only the bound on\n"
 "                                  its side is one on this problem\n"
 "  -D, --dryrun                    skip the compute() call\n"
 "  -v, --verbose[=N]               verbose output (0 = silent, 1 = basic,\n"
 "                                  2 = debug, with the parameters of the\n"
 "                                  Solvers of the Block, 3 = those of the\n"
 "                                  sub-Blocks' Solvers too)\n";

/*--------------------------------------------------------------------------*/
// open an SMS++ nc4 file and check that it is one

int read_open_netCDF( netCDF::NcFile & f , std::string fn )
{
 if( ! block_prefix.empty() && ! fn.empty()
     && fn.front() != '/' && fn.front() != '\\' )
  fn = block_prefix + fn;

 try {
  f.open( fn , netCDF::NcFile::read );
  }
 catch( netCDF::exceptions::NcException & ) {
  std::cerr << exe << ": cannot open nc4 file " << fn << std::endl;
  exit( 1 );
  }

 netCDF::NcGroupAtt gtype = f.getAtt( "SMS++_file_type" );
 if( gtype.isNull() ) {
  std::cerr << exe << ": " << fn << " is not an SMS++ nc4 file" << std::endl;
  exit( 1 );
  }

 int type;
 gtype.getValues( &type );

 if( ( type != eProbFile ) && ( type != eBlockFile ) ) {
  std::cerr << exe << ": " << fn << " is not a valid SMS++ file" << std::endl;
  exit( 1 );
  }

 return( type );
 }

/*--------------------------------------------------------------------------*/
// print the usage

void docopt( void )
{
 // http://docopt.org
 std::cout << docopt_desc << std::endl;
 std::cout << "Usage:" << std::endl
           << "  " << exe << " [options] <file>" << std::endl
           << "  " << exe << " -h | --help" << std::endl << std::endl
           << "Options:" << std::endl << help << std::endl;
 }

/*--------------------------------------------------------------------------*/
// read the tolerances of -E, a comma-separated list whose fields are either
// a number or, case-insensitively, "inf"
static void parse_eps_list( const std::string & arg )
{
 solver_eps.clear();
 std::size_t pos = 0;
 while( pos <= arg.size() ) {
  const std::size_t next = arg.find( ',' , pos );
  std::string field = arg.substr( pos , next == std::string::npos
                                        ? std::string::npos : next - pos );
  pos = ( next == std::string::npos ) ? arg.size() + 1 : next + 1;

  // trim the blanks that a quoted argument may carry
  const auto b = field.find_first_not_of( " \t" );
  const auto e = field.find_last_not_of( " \t" );
  field = ( b == std::string::npos ) ? "" : field.substr( b , e - b + 1 );

  // an empty field declares nothing for that Solver, but holds its place
  if( field.empty() ) {
   solver_eps.push_back( std::numeric_limits< double >::quiet_NaN() );
   continue;
   }

  std::string lc;
  for( auto c : field )
   lc += char( std::tolower( (unsigned char) c ) );

  if( ( lc == "inf" ) || ( lc == "+inf" ) || ( lc == "infinity" ) ) {
   solver_eps.push_back( std::numeric_limits< double >::infinity() );
   continue;
   }

  try {
   solver_eps.push_back( std::stod( field ) );
   }
  catch( ... ) {
   std::cerr << exe << ": invalid tolerance '" << field << "' in -E"
             << std::endl;
   exit( 1 );
   }
  }

 solver_eps_given = solver_eps;
 }

/*--------------------------------------------------------------------------*/
// read the roles of -R, a comma-separated list whose fields are either empty
// or, case-insensitively, "r"

static void parse_relaxation_list( const std::string & arg )
{
 solver_relaxation.clear();
 std::size_t pos = 0;
 while( pos <= arg.size() ) {
  const std::size_t next = arg.find( ',' , pos );
  std::string field = arg.substr( pos , next == std::string::npos
                                        ? std::string::npos : next - pos );
  pos = ( next == std::string::npos ) ? arg.size() + 1 : next + 1;

  // trim the blanks that a quoted argument may carry
  const auto b = field.find_first_not_of( " \t" );
  const auto e = field.find_last_not_of( " \t" );
  field = ( b == std::string::npos ) ? "" : field.substr( b , e - b + 1 );

  std::string lc;
  for( auto c : field )
   lc += char( std::tolower( (unsigned char) c ) );

  if( lc.empty() ) {          // an ordinary Solver, holding its place
   solver_relaxation.push_back( false );
   continue;
   }

  if( ( lc == "r" ) || ( lc == "rel" ) || ( lc == "relaxation" ) ) {
   solver_relaxation.push_back( true );
   continue;
   }

  std::cerr << exe << ": invalid role '" << field << "' in -R" << std::endl;
  exit( 1 );
  }
 }

/*--------------------------------------------------------------------------*/
// the command line, so that an option can look at the argument that follows
// it even when getopt does not hand it over [see the 'v' case below]

static int f_argc = 0;
static char ** f_argv = nullptr;

// whether the whole string is made of digits, so that it is the level of -v
// and not the argument that follows it

static bool is_number( const char * s )
{
 if( ( ! s ) || ( ! *s ) )
  return( false );
 for( ; *s ; ++s )
  if( ( *s < '0' ) || ( *s > '9' ) )
   return( false );
 return( true );
 }

/*--------------------------------------------------------------------------*/
// erase a standard short option, and its argument marker, from short_opts

void override_short_opt( char opt )
{
 const auto pos = short_opts.find( opt );
 if( pos == std::string::npos )
  return;

 auto end = pos + 1;
 while( ( end < short_opts.size() ) && ( short_opts[ end ] == ':' ) )
  ++end;

 short_opts.erase( pos , end - pos );
 }

/*--------------------------------------------------------------------------*/
// one standard option, or false if it is not one

bool process_standard_arg( int opt )
{
 switch( opt ) {
  case 'B': bconf_file = std::string( optarg ); break;
  case 'F': fconf_file = std::string( optarg ); break;
  case 'S': sconf_file = std::string( optarg ); break;
  case 'p': {
   block_prefix = normalize_prefix( std::string( optarg ) );
   Block::set_filename_prefix( std::string( block_prefix ) );
   break;
   }
  case 'c': conf_prefix = normalize_prefix( std::string( optarg ) );
            // also hand the -c prefix to Configuration, so filenames
            // referenced from inside a Configuration file (the "*filename"
            // includes and the strInnerBSC / str_LagBF_BSCfg meta-config
            // chains) are resolved against it, exactly as
            // Block::set_filename_prefix() does for -p
            Configuration::set_filename_prefix( std::string( conf_prefix ) );
            break;
  case 'E': parse_eps_list( std::string( optarg ) ); break;
  case 'R': parse_relaxation_list( std::string( optarg ) ); break;
  case 'D': dryrun = true; break;
  case 'v': {
   sol_verbose = true;
   // the level can be attached (-v2) or, since that is what one naturally
   // writes, separate (-v 2): getopt only gives the attached form, so the
   // next argument is taken when it is a number
   if( optarg )
    verbosity_level = std::atoi( optarg );
   else if( f_argv && ( optind < f_argc ) && is_number( f_argv[ optind ] ) )
    verbosity_level = std::atoi( f_argv[ optind++ ] );
   else
    verbosity_level = 1;
   break;
   }
  case 'h': docopt(); exit( 0 );
  case '?':
  default:  return( false );
  }
 return( true );
 }

/*--------------------------------------------------------------------------*/

bool filename_optional = false;

/*--------------------------------------------------------------------------*/
// the -S and -B that a test cannot do without

void require_solver_config( void )
{
 if( sconf_file.empty() )
  throw( std::invalid_argument(
   "a BlockSolverConfig must be provided (did you forget the -S option?)" ) );
 }

/*--------------------------------------------------------------------------*/

void require_block_config( void )
{
 if( bconf_file.empty() )
  throw( std::invalid_argument(
   "a BlockConfig must be provided (did you forget the -B option?)" ) );
 }

/*--------------------------------------------------------------------------*/
// parse the command line, with or without test-specific options

void process_args( int argc , char ** argv )
{
 process_args( argc , argv , nullptr );
 }

/*--------------------------------------------------------------------------*/

void process_args( int argc , char ** argv , bool ( *custom_arg )( int opt ) )
{
 exe = get_filename( argv[ 0 ] );
 f_argc = argc;
 f_argv = argv;

 // the `verbose` environment variable carries the level as well, so that a
 // battery, which has no way of passing -v to the test it runs, can ask for
 // the log of the Solver with `verbose=2 ./batch ...`; -v, where the test
 // understands it, is read afterwards and wins
 if( const char * e = std::getenv( "verbose" ) ; e && is_number( e ) )
  verbosity_level = std::atoi( e );

 while( true ) {  // options
  const auto opt = getopt_long( argc , argv , short_opts.data() ,
                                long_opts.data() , nullptr );
  if( opt == -1 ) break;

  // test-specific options are processed first: a test that re-defines one
  // of the standard letters means its own
  if( custom_arg && custom_arg( opt ) )  // test-specific option
   continue;                             // next

  if( process_standard_arg( opt ) )      // if it is a standard one
   continue;                             // next

  std::cout << "Try '" << exe << " --help' for more information"
            << std::endl;
  exit( 1 );
  }

 if( optind < argc )  // last positional argument == [Block] instance file
  filename = std::string( argv[ optind ] );
 else if( ! filename_optional ) {
  std::cout << exe << ": no input file" << std::endl
            << "Try '" << exe << " --help' for more information" << std::endl;
  exit( 1 );
  }
 }

/*--------------------------------------------------------------------------*/
/*------------------------ End common_utils.cpp ----------------------------*/
/*--------------------------------------------------------------------------*/