test.cpp

/*--------------------------------------------------------------------------*/
/*-------------------------- File test.cpp ---------------------------------*/
/*--------------------------------------------------------------------------*/
/** @file
 * Main for testing LagBFunction
 *
 * Given three input parameters n, k and p:
 *
 * - k "random" PolyhedralFunction are constructed, each inside a
 *   PolyhedralFunctionBlock.
 *
 * - p random uncapacitated *max-cost* transportation problems
 *   (with *negative* costs) on complete bipartite n-graphs (n origins, n
 *   destinations) are constructed "by hand", each inside an AbstractBlock;
 *   the problems have balanced supply and demands, no directed cycles,
 *   and no capacities on a craftly chosen set of "crucial" arcs, so as to
 *   ensure that they surely attain finite optimal solutions.
 *
 * - The above p AbstractBlock are inserted as inner Block, each inside a
 *   LagBFunction. If the command line parameter dictates that the
 *   LagBFunction will actually be computed (the NDO Solver does not have
 *   the "easy components" feature), an appropriate TP Solver (typically
 *   an LP Solver) is registered to all the inner Block of the LagBFunction.
 *
 * - The above p LagBFunction are put each inside the FRealObjective of a
 *   new AbstractBlock, otherwise empty.
 *
 * - The k PolyhedralFunctionBlock (configured to use "natural"
 *   representation) and the p AbstractBlock are inserted inside a single
 *   AbstractBlock (NDOBlock), possibly with a linear function, to represent
 *   a problem of the form (for k = 1 and p = 2)
 *
 *   min { l x + f(x) + max { ( B x + c_1 ) z_1 : E z_1 = b_1 , 0 <= z_1 }
 *                    + max { ( B x + c_2 ) z_2 : E z_2 = b_2 , 0 <= z_2 } }
 *
 *   Assuming that arcs are ordered lexicographically (first all the ones
 *   outgoing node 0, ordered by tail node, then all the ones outgoing node
 *   1, ...), the matrix E has the form (ignoring bound constraints)
 *
 *             n^2
 *      | e^T  0  ...  0  |
 *   n  |  0  e^T ...  0  |
 *      |  :  :        :  |
 *      |  0  0   ... e^T |
 *      +-----------------+
 *   n  |  I  I   ...  I  |
 *
 *   corresponding to constraints (with I = J =  { 0 ... n - 1 })
 *
 *      \sum_{ j \in J } z[ i ][ j ] == s[ i ]   i \in I
 *
 *      \sum_{ i \in I } z[ i ][ j ] == d[ j ]   j \in J
 *
 *   with s[] and d[] being the vectors of supplies and demand. The matrix
 *   B has the form
 *
 *             n^2
 *   n  |  I  I   ...  I  |
 *
 *   corresponding to the fact that the cost of arc ( i , j ) is
 *
 *         c[ i ][ j ] + x[ j ]
 *
 *   i.e., x[ j ] is added to the cost of all arcs ingoing destination j.
 *
 * - Some arc ( i , j ) will have a finite upper bound u[ i ][ j ] >= 0,
 *   with the value 0 being possible (basically, fixing the variable).
 *   However, this immediately creates the risk that the transportation
 *   problem is unfeasible, which we don't want to handle. To avoid that,
 *   the following cunning plan has been devised:
 *
 *   = none of the arcs ( i , i ) will ever have finite upper bound;
 *
 *   = the i-th supply and demand will be equal: s[ i ] = d[ i ].
 *
 *   This guarantees that satisfying the i-th supply/demand pair via the
 *   direct arc ( i , j ) is always possible, albeit is may easily not be
 *   the best choice due to the costs being random. Note that one may have
 *   also put any finite upper bound >= s[ i ] = d[ i ] for this to work,
 *   but if bounds and capacities are randomly changed then one should be
 *   careful to guarantee that this always holds; by not having the bound
 *   at all we guarantee that this can never be a problem, since we will
 *   never create a finite upper bound when an infinite one was (nor
 *   vice-versa, for that matter).
 *
 * - An appropriate NDO Solver is attached to NDOBlock; this can in general
 *   be any Solver capable of solving it, but some specific provisions
 *   are done for BundleSolver, in particular when very verbose log is
 *   activated.
 *
 * Then, an LP equivalent of NDOBlock is constructed into a different
 * AbstractBlock (LPBlock) with the following steps:
 *
 * - The linear objective and the k PolyhedralFunctionBlock are just
 *   copied over, the latter using the R3Block.
 *
 * - For the p LagBFunction, an LP equivalent is constructed by the
 *   following derivation:
 *
 *   min { l x + f(x) + max { ( B x + c_1 ) z_1 : E z_1 = b_1  , 0 <= z_1 }
 *                    + max { ( B x + c_2 ) z_2 : E z_2 = b_2  , 0 <= z_2 }
 *         } =
 *
 *   min { l x + f(x) + min { y_1 b_1 : y_1 E >= B x + c_1 }
 *                    + min { y_2 b_2 : y_2 E >= B x + c_2 } } =
 *
 *   min { l x + f(x) + y_1 b_1 + y_2 b_2 :
 *         y_1 E >= B x + c_1 , y_2 E >= B x + c_2 }
 *
 *   Since the transpose of E has the form
 *
 *              n        n
 *         | e 0 ... 0 | I |
 *   n^2   | 0 e ... 0 | I |
 *         | : :     : | : |
 *         | 0 0 ... e | I |
 *
 *   this corresponds to variables yo[ i ] and yd[ j ] for each origin and
 *   destination, with costs s[ i ] and d[ j ] respectively, as well as
 *   constraints
 *
 *      ys[ i ] + yd[ j ] - x[ j ] >= c[ i ][ j ]
 *
 *   for all i \in I and j \in J. This works if the variable z[ i ][ j ] has
 *   *no* finite bound; if, instead, the constraint
 *
 *      z[ i ][ j ] <= u[ i ][ j ]
 *
 *   is present (with u[ i ][ j ] = 0 possible), then it has a dual variable
 *   w[ i ][ j ]; this means that a term u[ i ][ j ] * w[ i ][ j ] is added
 *   to the objective function, and that the constraint becomes
 *
 *      ys[ i ] + yd[ j ] + w[ i ][ j ] - x[ j ] >= c[ i ][ j ]
 *
 *      w[ i ][ j ] >= 0
 *
 *   All this, however, is only correct for the convex case; in the concave
 *   one, the problem is rather
 *
 *   max { l x + f(x) + min { ( B x + c_1 ) z_1 : E z_1 = b_1  , 0 <= z_1 }
 *                    + min { ( B x + c_2 ) z_2 : E z_2 = b_2  , 0 <= z_2 } }
 *   yelding
 *
 *   max { l x + f(x) + y_1 b_1 + y_2 b_2 :
 *         y_1 E <= B x + c_1 , y_2 E <= B x + c_2 }
 *
 *   and therefore the constraints are
 *
 *      ys[ i ] + yd[ j ] - w[ i ][ j ] - x[ j ] <= c[ i ][ j ]
 *
 *      w[ i ][ j ] >= 0
 *
 *   assuming u[ i ][ j ] is finite, else without the "- w[ i ][ j ]" term
 *   and the w[ i ][ j ] variable; note having put a "-" to keep w[ i ][ j ]
 *   non-negative, since the natural sign of the dual variable of a <=
 *   constraint in a minimization LP is <= 0. However, this means that the
 *   corresponding term in the objective, that would ordinarily be
 *
 *       + u[ i ][ j ] * w[ i ][ j ]
 *
 *   must then become
 *
 *       - u[ i ][ j ] * w[ i ][ j ]
 *
 *   because w[ i ][ j ] has changed sign (w[ i ][ j ] ==> - w[ i ][ j ]).
 *
 * - The variables ys[ i ] and yd[ j ], the objective function and the
 *   constraints (linking them with x[]) are constructed manually into an
 *   AbstractBlock for each p.
 *
 * - A LPSolver is attached to this AbstractBlock
 *
 * The PolyhedralFunction and/or the costs and demands (not supplies) of the
 * uncapacitated transportation problems are then repeatedly randomly
 * modified "in the same way", and re-solved several times; results of the
 * two solvers are compared.
 *
 * \author Antonio Frangioni \n
 *         Dipartimento di Informatica \n
 *         Universita' di Pisa \n
 *
 * \copyright &copy; by Antonio Frangioni
 */
/*--------------------------------------------------------------------------*/
/*-------------------------------- MACROS ----------------------------------*/
/*--------------------------------------------------------------------------*/

#define LOG_LEVEL 0
// 0 = only pass/fail
// 1 = result of each test
// 2 = + solver log
// 3 = + save LP file
// 4 = + save every LP for every iteration
// 5 = + print data

#if( LOG_LEVEL >= 1 )
 #define LOG1( x ) cout << x
 #define CLOG1( y , x ) if( y ) cout << x

 #if( LOG_LEVEL >= 2 )
  #define LOG_ON_COUT 1
  // if nonzero, the NDO Solver log is sent on cout rather than on a file
 #endif
#else
 #define LOG1( x )
 #define CLOG1( y , x )
#endif

/*--------------------------------------------------------------------------*/
// if HAVE_CONSTRAINTS == 1, then about 50% of the variables will have a
// non-negativity constraint implemented via ColVariable::is_positive()
// if HAVE_CONSTRAINTS == 2, then about 50% of the variables will have
// bound constraints; of these, 33% will only have 0 lower bound, 33% will
// only have random upper bound, and the rest will have both. of the
// remaining 50% of the variables, another 50% will have a non-negativity
// constraint implemented via ColVariable::is_positive()

#define HAVE_CONSTRAINTS 2

/*--------------------------------------------------------------------------*/
// if nonzero, the Solver attached to the LPBlock is detached and re-attached
// to it at all iterations

#define DETACH_LP 0

/*--------------------------------------------------------------------------*/
// if nonzero, the two Block are not solved at every round of changes, but
// only every SKIP_BEAT + 1 rounds. this allows changes to accumulate, and
// therefore puts more pressure on the Modification handling of the Solver
// (in case this tries to do "smart" things rather than dumbly processing
// each one in turn)
//
// note that the number of rounds of changes is them multiplied by
// SKIP_BEAT + 1, so that the input parameter still dictates the number of
// compute() calls per (Solver per) Block

#define SKIP_BEAT 3

/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
// if 1, the w variables are dynamic

#define DYNAMIC_w 0

/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
// if 1, the bc constraints are dynamic

#define DYNAMIC_bc 0

/*- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*/
// if 1, half of the variables are dynamic

#define DYNAMIC_VARS 0

/*--------------------------------------------------------------------------*/
/*------------------------------ INCLUDES ----------------------------------*/
/*--------------------------------------------------------------------------*/

#include <algorithm>
#include <chrono>
#include <cmath>
#include <fstream>
#include <sstream>
#include <iomanip>
#include <limits>
#include <map>
#include <memory>

#include <random>

#include "MILPSolver.h"

#include "common_utils.h"

#include "LagBFunction.h"

#include "PolyhedralFunctionBlock.h"

#include "UpdateSolver.h"

/*--------------------------------------------------------------------------*/
/*-------------------------------- USING -----------------------------------*/
/*--------------------------------------------------------------------------*/

using namespace std;

using namespace SMSpp_di_unipi_it;

/*--------------------------------------------------------------------------*/
/*-------------------------------- TYPES -----------------------------------*/
/*--------------------------------------------------------------------------*/

using Index = Block::Index;

using Range = Block::Range;

using Subset = Block::Subset;

using FunctionValue = Function::FunctionValue;
using c_FunctionValue = Function::c_FunctionValue;

using MultiVector = PolyhedralFunction::MultiVector;
using RealVector = PolyhedralFunction::RealVector;
using VarVector = PolyhedralFunction::VarVector;

using v_coeff_pair = LinearFunction::v_coeff_pair;

using p_AB = AbstractBlock *;
using p_PFB = PolyhedralFunctionBlock *;

using p_LF = LinearFunction *;
using p_PF = PolyhedralFunction *;
using p_LBF = LagBFunction *;

/*--------------------------------------------------------------------------*/
/*------------------------------- CONSTANTS --------------------------------*/
/*--------------------------------------------------------------------------*/

const double scale = 10;
const char * const logF = "log.bn";

c_FunctionValue INF = SMSpp_di_unipi_it::Inf< FunctionValue >();

/*--------------------------------------------------------------------------*/
/*------------------------------- GLOBALS ----------------------------------*/
/*--------------------------------------------------------------------------*/

p_AB LPBlock;              // the "partially dualised" LP representaion

p_AB NDOBlock;             // the "natural" NDO representation

bool convex = true;        // true if everything is convex

double bound = 1000;       // a tentative bound to detect unbounded instances

FunctionValue BND;         // the bound in the PolyhedralFunction (if any)

Index nvar = 10;           // number of variables
#if DYNAMIC_VARS > 0
 Index nsvar;              // number of static variables
 Index ndvar;              // number of dynamic variables
#else
 #define nsvar nvar        // all variables are static
#endif

std::mt19937 rg;           // base random generator
std::uniform_real_distribution<> dis( 0.0 , 1.0 );

MultiVector A;             // rows
RealVector b;              // constants

MultiVector C;             // arc costs
MultiVector U;             // arc capacities
RealVector s;              // supplies == demands

int SizeEasy = 0;          // how an easy transportation problem is sized

bool DetachNDO = false;    // if the Solver of NDOBlock are detached and
                           // re-attached to it before each compute()

/* The last command line parameter, size, says how the inner Block of an easy
 * LagBFunction carries the size Variable that the dual master problem of
 * BundleSolver scales it with [see MasterProblemBlock::configure()]: in the
 * dual of the level and doubly stabilized master problems, as in that of any
 * master problem with a global lower bound, the mass lambda of the lower
 * model may differ from one, and the feasible set of each easy component has
 * to be scaled with it. With a size Variable tau the balance and capacity
 * constraints read
 *
 *          sum_j f_ij - s_i tau = 0,
 *          sum_i f_ij - s_j tau = 0,
 *          f_ij - U_ij tau <= 0,
 *
 * and tau is
 *
 * - 0: absent, the balances keep their right-hand sides and the capacities
 *   are BoxConstraint; the inner Block is a plain AbstractBlock, which takes
 *   no size Variable, so that lambda is fixed to one, unless bit 5 of
 *   intDoEasy has the master put in its place a copy of it sized by lambda
 *   [see AbstractBlock::set_size_variable()], in which the master writes
 *   the changes the tester makes to it;
 *
 * - 1: a nonnegative ColVariable of the inner Block, declared with
 *   Block::set_owned_size_variable() and tied to lambda by the master;
 *
 * - 2: a Variable of the master given to the inner Block, a
 *   SizedTransportBlock, by Block::set_size_variable(), which writes it into
 *   the rows above as it gets it; the capacities are FRowConstraint, so that
 *   they can take it.
 *
 * Whether the master uses the size Variable, or the copy, is up to intDoEasy
 * of the BundleSolver [see BundleSolver::intDoEasy], and the stabilization
 * is up to its ComputeConfig. With a hard LagBFunction the parameter is
 * ignored.
 * The changes of supplies and capacities are written in the coefficients of
 * tau when it is there, in the right-hand sides otherwise. */

/*--------------------------------------------------------------------------*/
/*------------------------------- CLASSES ----------------------------------*/
/*--------------------------------------------------------------------------*/
/// the transportation problem of an easy LagBFunction that takes its size
/** An AbstractBlock with the transportation problem built by main(), whose
 * balances "sc" and "dc" and finite capacities "bc" are all FRowConstraint,
 * which takes the size Variable of its father [see
 * Block::set_size_variable()]: each of those rows, whose finite sides are
 * all equal to some c, gets the Variable with coefficient - c, and the
 * finite sides go to zero. The supplies and capacities are then the
 * coefficients of the Variable, which is how main() changes them [see
 * size_variable_of()]. It takes a single Variable, and the same one again
 * changes nothing; the rows that are there are rewritten, at any time, and
 * \p issueAMod says which Modification are issued. */

class SizedTransportBlock : public AbstractBlock
{
 public:

 bool set_size_variable( Variable * size_var ,
                         c_ModParam issueAMod = eNoBlck ) override {
  auto tau = dynamic_cast< ColVariable * >( size_var );
  if( ( ! tau ) || ( f_size && ( f_size != tau ) ) )
   return( false );
  if( f_size )  // there already
   return( true );

  f_size = tau;
  auto scale = [ & ]( FRowConstraint & row ) {
   const double c = std::isfinite( row.get_rhs() ) ? row.get_rhs()
                                                     : row.get_lhs();
   static_cast< p_LF >( row.get_function() )->add_variable( tau , - c ,
                                                            issueAMod );
   // the row registers itself with the Variable only when it receives the
   // Modification, which issueAMod may not ask for
   const auto & act = tau->active_stuff();
   ThinVarDepInterface * me = & row;
   if( ! std::binary_search( act.begin() , act.end() , me ) )
    tau->add_active( me );
   if( row.get_lhs() == row.get_rhs() )
    row.set_both( 0.0 , issueAMod );
   else
    if( std::isfinite( row.get_rhs() ) )
     row.set_rhs( 0.0 , issueAMod );
    else
     row.set_lhs( 0.0 , issueAMod );
   };

  for( auto & row : *get_static_constraint_v< FRowConstraint >( "sc" ) )
   scale( row );
  for( auto & row : *get_static_constraint_v< FRowConstraint >( "dc" ) )
   scale( row );
  #if DYNAMIC_bc
   if( auto bc = get_dynamic_constraint< FRowConstraint >( "bc" ) )
  #else
   if( auto bc = get_static_constraint_v< FRowConstraint >( "bc" ) )
  #endif
    for( auto & row : *bc )
     scale( row );

  return( true );
  }

 /// the size Variable this Block has taken, nullptr if none
 [[nodiscard]] ColVariable * get_given_size_variable( void ) const {
  return( f_size );
  }

 private:

 ColVariable * f_size = nullptr;  ///< the size Variable taken, if any
};

/*--------------------------------------------------------------------------*/
/*------------------------------ FUNCTIONS ---------------------------------*/
/*--------------------------------------------------------------------------*/
/// the size Variable the rows of the transportation problem \p tb carry
/** The one \p tb owns [see Block::get_size_variable()] or the one it has
 * taken [see SizedTransportBlock], and nullptr if its rows carry none. */

static ColVariable * size_variable_of( Block * tb )
{
 if( auto tau = tb->get_size_variable() )
  return( static_cast< ColVariable * >( tau ) );
 if( auto stb = dynamic_cast< SizedTransportBlock * >( tb ) )
  return( stb->get_given_size_variable() );
 return( nullptr );
 }

/*--------------------------------------------------------------------------*/
// convex ==> minimize ==> negative numbers

static double rs( double x ) { return( convex ? -x : x ); }

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

static double rndfctr( void )
{
 // a random number between 0.5 and 2, with 50% probability of being < 1
 double fctr = dis( rg ) - 0.5;
 return( fctr < 0 ? - fctr : fctr * 4 );
 }

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

static void GenerateAi( RealVector & Ai , Index nc )
{
 Ai.resize( nc );
 for( auto & aij : Ai )
  aij = scale * ( 2 * dis( rg ) - 1 );
 }

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

static void GenerateA( Index nr , Index nc )
{
 A.resize( nr );
 for( auto & Ai : A )
  GenerateAi( Ai , nc );
 }

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

static void Generateb( Index nr )
{
 b.resize( nr );
 for( auto & bj : b )
  bj = scale * nvar * ( 2 * dis( rg ) - 1 ) / 4;
 }

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

static void GenerateAb( Index nr , Index nc )
{
 // rationale: the solution x^* will be more or less the solution of some
 // square sub-system A_B x = b_B. We want x^* to be "well scaled", i.e.,
 // the entries to be ~= 1 (in absolute value). The average of each row A_i
 // is 0, the maximum (and minimum) expected value is something like
 // scale * nvar / 2. So we take each b_j in +- scale * nvar / 4

 GenerateA( nr , nc );
 Generateb( nr );
 }

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

static void GenerateBND( void )
{
 // rationale: we expect the solution x^* to have entries ~= 1 (in absolute
 // value, and the coefficients of A are <= scale (in absolute value), so
 // the LHS should be at most around - scale * nvar; the RHS can add it
 // a further - scale * nvar / 4, so we expect - (5/4) * scale * nvar to
 // be a "natural" LB. We therefore set the LB to a mean of 1/2 of that
 // (tight) 33% of the time, a mean of 2 times that (loose) 33% of the time,
 // and -INF the rest

 if( dis( rg ) <= 0.333 ) {  // "tight" bound
  BND = rs( dis( rg ) * 5 * scale * nvar / 4 );
  return;
  }

 if( dis( rg ) <= 0.333 ) {  // "loose" bound
  BND = rs( dis( rg ) * 5 * scale * nvar );
  return;
  }

 BND = INF;
 }

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

static Index GenerateCi( RealVector & Ci )
{
 // 10% of costs are zero, the rest random between:
 // -10 and 0 in the convex case (where the problem is max)
 // 0 and 10 in the concave case (where the problem is min)

 Index nzc = 0;
 if( convex ) {
  for( auto & cij : Ci )
   if( dis( rg ) < 0.1 )
    cij = 0;
   else {
    cij = - 10 * dis( rg );
    ++nzc;
    }
  }
 else
  for( auto & cij : Ci )
   if( dis( rg ) < 0.1 )
    cij = 0;
   else {
    cij = 10 * dis( rg );
    ++nzc;
    }

 return( nzc );
 }

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

static Index GenerateCosts( void )
{
 Index nzc = 0;
 for( auto & Ci : C )
  nzc += GenerateCi( Ci );

 return( nzc );
 }

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

static void GenerateSupplies( void )
{
 s.resize( nvar );
 for( auto & si : s )
  si = 10 * dis( rg );
 }

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

static void GenerateCapacities( RealVector & Ui )
{
 // in 85% of the cases the capacity is random between 0 and 5 (to make it
 // hopefully "byte" over a demand between 0 and 10), in the remaining 15%
 // of the cases it is 0 (which "surely bytes")
 for( auto & uij : Ui )
  uij = dis( rg ) < 0.15 ? 0 : 5 * dis( rg );
 }

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

static Index GenerateCapacities( void )
{
 // in 30% and in all arcs ( i , i ) the capacity is infinite, otherwise as
 // in GenerateCapacities( Ui )
 Index nic = 0;

 for( Index i = 0 ; i < nvar ; ++i )
  for( Index j = 0 ; j < nvar ; ++j )
   if( ( i == j ) || ( dis( rg ) < 0.3 ) )
    U[ i ][ j ] = INF;
   else {
    U[ i ][ j ] = dis( rg ) < 0.15 ? 0 : 5 * dis( rg );
    ++nic;
    }

 return( nic );
 }

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

static Subset GenerateSubset( Index m , Index k )
{
 // generate a sorted random k-vector of unique integers in 0 ... m - 1

 Subset rnd( m );
 std::iota( rnd.begin() , rnd.end() , 0 );
 std::shuffle( rnd.begin() , rnd.end() , rg );
 rnd.resize( k );
 sort( rnd.begin() , rnd.end() );

 return( std::move( rnd ) );
 }

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

static void printAb( const MultiVector & tA , const RealVector & tb ,
		     double bound )
{
 cout << "n = " << nvar << ", m = " << tA.size();
 if( convex )
  cout << " (convex)";
 else
  cout << " (concave)";
 if( std::abs( bound ) == INF )
  cout << " (no bound)" << endl;
 else
  cout << ", bound = " << bound << endl;

 for( Index i = 0 ; i < tA.size() ; ++i ) {
  cout << "A[ " << i << " ] = [ ";
  for( Index j = 0 ; j < nvar ; ++j )
   cout << tA[ i ][ j ] << " ";
  cout << "], b[ " << i << " ] = " << tb[ i ] << endl;
  }
 }

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

static void printC( void )
{
 for( Index i = 0 ; i < nvar ; ++i ) {
  cout << "C[ " << i << " ] = [ ";
  for( Index j = 0 ; j < nvar ; ++j )
   cout << C[ i ][ j ] << " ";
  cout << "]" << endl;
  }
 }

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

static void printU( void )
{
 for( Index i = 0 ; i < nvar ; ++i ) {
  cout << "U[ " << i << " ] = [ ";
  for( Index j = 0 ; j < nvar ; ++j )
   if( U[ i ][ j ] == INF )
    cout << "INF ";
   else
    cout << U[ i ][ j ] << " ";
  cout << "]" << endl;
  }
 }

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

static void printT( void )
{
 cout << "s = [ ";
 for( Index j = 0 ; j < nvar ; ++j )
  cout << s[ j ] << " ";
 cout << "]" << endl;

 printC();
 printU();
 }

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

static void ConstructObj( p_AB AB )
{
 // construct the Linear Objective (FRealObjective with a LinearFunction
 // inside) in the given AbstractBlock on the "x" static Variable; this is
 // only called if nf < 0

 auto x = AB->get_static_variable_v< ColVariable >( "x" );
 #if DYNAMIC_VARS > 0
  auto xd = AB->get_dynamic_variable< ColVariable >( "x" );
 #endif

 v_coeff_pair cp( nvar );
 Index i = 0;
 // static x
 for( ; i < nsvar ; ++i )
  cp[ i ] = std::make_pair( &((*x)[ i ] ) , A[ 0 ][ i ] );

 #if DYNAMIC_VARS > 0
  // dynamic x
  auto xdit = xd.begin();
  for( ; i < nvar ; ++i , ++xdit )
   cp[ j ] = std::make_pair( &(*xdit) , A[ 0 ][ i ] );
 #endif

 auto obj = new FRealObjective( AB , new LinearFunction( std::move( cp ) ) );
 obj->set_sense( convex ? Objective::eMin : Objective::eMax , eNoMod );
 AB->set_objective( obj , eNoMod );
 }

/*--------------------------------------------------------------------------*/
/// the value of a LagBFunction in the current value of its Variable, by copy
/** Returns the value of \p lbf in the current value of the multipliers y
 * without computing \p lbf, which is what a LagBFunction treated as easy
 * needs: its inner Block has no Solver, the NDO Solver handling it exactly
 * in its master, and computing it with a Solver lent for the purpose would
 * put into the Objective of its inner Block the Variable it keeps aside,
 * changing what the following rounds of changes pick. The inner Block is
 * copied [see AbstractBlock::mirror()], the copy is given the Lagrangian
 * objective c( x ) + < y , g( x ) > in the sense of the original, and the
 * copy is solved with the Solver of TPPar.txt; the value is NaN if anything
 * of this cannot be done. */

static double LBF_value_by_copy( LagBFunction * lbf )
{
 const double NaN = std::numeric_limits< double >::quiet_NaN();

 AbstractBlock copy;
 copy.mirror( lbf->get_inner_block() );
 if( ! copy.get_mirror_issues().empty() )
  return( NaN );

 // Evaluate the original transportation problem, irrespective of the
 // size chosen for this component in the last master solution: the size
 // Variable is fixed to 1 in the copy if it is a Variable of the inner
 // Block, and otherwise, being shared by the copy, taken out of its rows
 // with its coefficient moved into their sides
 if( auto tau = size_variable_of( lbf->get_inner_block() ) ) {
  if( auto tau_copy = copy.mirror_of( tau ) ) {
   tau_copy->is_fixed( false , eNoMod );
   tau_copy->set_value( 1.0 );
   tau_copy->is_fixed( true , eNoMod );
   }
  else {
   std::vector< FRowConstraint * > rows;
   for( auto stuff : tau->active_stuff() )
    if( auto row = dynamic_cast< FRowConstraint * >( stuff ) )
     if( row->get_Block() == & copy )
      rows.push_back( row );
   for( auto row : rows ) {
    auto lf = static_cast< p_LF >( row->get_function() );
    const Index ti = lf->is_active( tau );
    const double a = lf->get_coefficient( ti );
    row->remove_variable( ti , eNoMod );
    if( std::isfinite( row->get_lhs() ) )
     row->set_lhs( row->get_lhs() - a , eNoMod );
    if( std::isfinite( row->get_rhs() ) )
     row->set_rhs( row->get_rhs() - a , eNoMod );
    }
   }
  }

 auto obj = copy.get_objective< FRealObjective >();
 auto clf = obj ? dynamic_cast< LinearFunction * >( obj->get_function() )
                : nullptr;
 if( ! clf )
  return( NaN );

 // the Lagrangian costs: first the original costs of the Variable of the
 // Objective, in their order, which are those the LagBFunction keeps [see
 // LagBFunction::get_A_by_col()] and not those in the Objective, where a
 // LagBFunction that has been computed has written its Lagrangian ones
 LinearFunction::v_coeff_pair cost( clf->get_v_var() );
 std::map< ColVariable * , Index > pos;
 for( Index h = 0 ; h < cost.size() ; ++h ) {
  auto col = lbf->get_A_by_col( copy.mirrored_of( cost[ h ].first ) );
  if( ! col )
   return( NaN );
  cost[ h ].second = col->first;
  pos[ cost[ h ].first ] = h;
  }
 double ct = clf->get_constant_term();

 for( Index i = 0 ; i < lbf->get_num_active_var() ; ++i ) {
  const double yi =
   static_cast< ColVariable * >( lbf->get_active_var( i ) )->get_value();
  auto gi = static_cast< LinearFunction * >( lbf->get_Lagrangian_term( i ) );
  ct += yi * gi->get_constant_term();
  for( const auto & [ xj , aij ] : gi->get_v_var() ) {
   auto cj = copy.mirror_of( static_cast< ColVariable * >( xj ) );
   if( ! cj )
    return( NaN );
   auto it = pos.find( cj );
   if( it == pos.end() ) {
    pos[ cj ] = cost.size();
    cost.emplace_back( cj , yi * aij );
    }
   else
    cost[ it->second ].second += yi * aij;
   }
  }

 obj->set_function( new LinearFunction( std::move( cost ) , ct ) , eNoMod );

 auto cfg = Configuration::deserialize( "TPPar.txt" );
 auto tpc = dynamic_cast< BlockSolverConfig * >( cfg );
 if( ! tpc ) {
  delete( cfg );
  return( NaN );
  }

 tpc->apply( & copy );
 double value = NaN;
 if( ! copy.get_registered_solvers().empty() ) {
  auto slvr = copy.get_registered_solvers().front();
  const int rtrn = slvr->compute( false );
  if( ( rtrn >= Solver::kOK ) && ( rtrn < Solver::kError ) )
   value = slvr->get_var_value();
  }

 tpc->clear();  // take the Solver away again
 tpc->apply( & copy );
 delete( tpc );

 return( value );
 }

/*--------------------------------------------------------------------------*/
/// the value of the objective of NDOBlock in the current values of x
/** A LagBFunction whose inner Block has a size variable or no Solver is
 * valued by copy [see LBF_value_by_copy()], a PolyhedralFunctionBlock by
 * computing its PolyhedralFunction, since it may have no Objective, and any
 * other Function by computing it; a Block with no Objective, as NDOBlock
 * without the linear function, gives nothing. A sized inner Block must be
 * evaluated at unit size. */

static double NDO_value( void )
{
 double value = 0;
 if( auto obj = dynamic_cast< FRealObjective * >(
                                            NDOBlock->get_objective() ) ) {
  obj->compute();
  value += obj->value();
  }
 for( Index b = 0 ; b < NDOBlock->get_number_nested_Blocks() ; ++b ) {
  auto nb = NDOBlock->get_nested_Block( b );
  if( auto pfb = dynamic_cast< p_PFB >( nb ) ) {
   auto & pf = pfb->get_PolyhedralFunction();
   pf.compute();
   value += pf.get_value();
   continue;
   }
  if( auto obj = dynamic_cast< FRealObjective * >( nb->get_objective() ) ) {
   auto lbf = dynamic_cast< LagBFunction * >( obj->get_function() );
   if( lbf && ( size_variable_of( lbf->get_inner_block() ) ||
                lbf->get_inner_block()->get_registered_solvers().empty() ) )
    value += LBF_value_by_copy( lbf );
   else {
    obj->compute();
    value += obj->value();
    }
   }
  }
 return( value );
 }

/*--------------------------------------------------------------------------*/
/// true if the objective of NDOBlock improves along the ray of the LP
/** Called when the LP Solver has proved the LP unbounded, this reads the
 * unbounded direction it gives as the certificate and, from the point where
 * the NDO Solver has stopped, moves x along it by two steps s and 2 s, with s
 * large with respect to both the point and the direction: the objective of
 * NDOBlock has to improve at both, as it does along a direction of recession
 * of the LP, which is an exact reformulation of the NDO problem. A direction
 * that does not move x says nothing about NDOBlock, and is not taken as an
 * agreement. The values of the Variable of both Block are put back as they
 * were, and a LagBFunction treated as easy is valued by copy [see
 * NDO_value()]. */

static bool NDO_improves_along_LP_ray( Solver * slvrLP )
{
 if( ! slvrLP->has_var_direction() )
  return( false );

 auto LPx = LPBlock->get_static_variable_v< ColVariable >( "x" );
 auto NDOx = NDOBlock->get_static_variable_v< ColVariable >( "x" );
 const Index n = LPx->size();

 std::vector< double > x0( n );
 double xn = 0;
 for( Index i = 0 ; i < n ; ++i ) {
  x0[ i ] = (*NDOx)[ i ].get_value();
  xn = std::max( xn , std::abs( x0[ i ] ) );
  }

 // the ray is written in all the Variable of LPBlock, whose values are put
 // back as they were once the part of it along x has been read
 std::unique_ptr< Solution > saved( LPBlock->get_Solution( nullptr , false ) );
 slvrLP->get_var_direction();
 std::vector< double > d( n );
 double dn = 0;
 for( Index i = 0 ; i < n ; ++i ) {
  d[ i ] = (*LPx)[ i ].get_value();
  dn = std::max( dn , std::abs( d[ i ] ) );
  }
 if( saved )
  saved->write( LPBlock );

 if( dn == 0 )
  return( false );

 const double s = 1e+3 * ( 1 + xn ) / dn;
 auto value_at = [ & ]( double k ) {
  for( Index i = 0 ; i < n ; ++i )
   (*NDOx)[ i ].set_value( x0[ i ] + k * d[ i ] );
  return( convex ? NDO_value() : - NDO_value() );
  };

 const double f0 = value_at( 0 );
 const double f1 = value_at( s );
 const double f2 = value_at( 2 * s );

 for( Index i = 0 ; i < n ; ++i )
  (*NDOx)[ i ].set_value( x0[ i ] );

 const double eps = 1e-9 * std::max( { 1.0 , std::abs( f0 ) , std::abs( f1 ) } );
 return( ( f1 < f0 - eps ) && ( f2 < f1 - eps ) );
 }

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

static bool SolveEach( void )
{
 // LPBlock with one Solver and NDOBlock with as many as its BlockSolverConfig
 // attaches, each of the latter checked against the former
 try {
  // solve the LPBlock- - - - - - - - - - - - - - - - - - - - - - - - - - - -
  Solver * slvrLP = ( LPBlock->get_registered_solvers() ).front();
  #if DETACH_LP
   LPBlock->unregister_Solver( slvrLP );
   LPBlock->register_Solver( slvrLP , true );  // push it to the front
  #endif
  auto startLP = std::chrono::system_clock::now();
  int rtrnLP = slvrLP->compute( false );
  auto endLP = std::chrono::system_clock::now();
  double tLP = std::chrono::duration< double >( endLP - startLP ).count();
  bool hsLP = ( ( rtrnLP >= Solver::kOK ) && ( rtrnLP < Solver::kError ) )
              || ( rtrnLP == Solver::kLowPrecision );
  double foLP = hsLP ? ( convex ? slvrLP->get_ub() : slvrLP->get_lb() )
                     : ( convex ? INF : -INF );

  auto tok = []( bool hs , int rtrn , double fo ) -> std::string {
   if( hs )                              return( fmt_obj( fo ) );
   if( rtrn == Solver::kInfeasible )     return( "Unfeas" );
   if( rtrn == Solver::kUnbounded )      return( "Unbounded" );
   return( "Error!" );
   };
  std::vector< double > times = { tLP };
  std::vector< std::string > values = { tok( hsLP , rtrnLP , foLP ) };

  /* The two values are compared at the relative tolerance below, save that a
   * Solver answering kLowPrecision has declared that it did not reach the
   * accuracy it was asked for, and has said how far it is: the two bounds it
   * reports are around the value, and their distance is the accuracy it did
   * reach (they may well have crossed, which is what not reaching it looks
   * like). Holding such an answer to the tolerance of an exact one asks for
   * more than the Solver ever promised, so the declared gap is added to the
   * tolerance: the answer is off by the gap at worst, and the comparison
   * carries the ordinary tolerance on top of it. */

  auto declared_gap = []( Solver * slvr , int rtrn ) -> double {
   if( rtrn != Solver::kLowPrecision )
    return( 0 );
   const double lb = slvr->get_lb();
   const double ub = slvr->get_ub();
   if( ( lb <= -INF ) || ( ub >= INF ) )
    return( 0 );
   return( abs( ub - lb ) );
   };

  // solve the NDOBlock with each of its Solver- - - - - - - - - - - - - - - -
  // bespoke verdict of each against LPBlock, covering the conditional valid
  // bound, which is doubled once, after all of them, if any of them has
  // stopped at it; then the per-instance line of common_utils
  bool ok = true;
  std::string verdict;
  bool dbl_bound = false;
  const auto slvrs = NDOBlock->get_registered_solvers();  // a copy
  for( auto slvrNDO : slvrs ) {
   if( DetachNDO ) {
    NDOBlock->unregister_Solver( slvrNDO );
    NDOBlock->register_Solver( slvrNDO );
    }
   auto startNDO = std::chrono::system_clock::now();
   int rtrnNDO = slvrNDO->compute( false );
   auto endNDO = std::chrono::system_clock::now();
   times.push_back(
           std::chrono::duration< double >( endNDO - startNDO ).count() );
   bool hsNDO = ( ( rtrnNDO >= Solver::kOK ) && ( rtrnNDO < Solver::kError ) )
                || ( rtrnNDO == Solver::kLowPrecision );
   double foNDO = hsNDO ? ( convex ? slvrNDO->get_ub() : slvrNDO->get_lb() )
                        : ( convex ? INF : -INF );
   values.push_back( tok( hsNDO , rtrnNDO , foNDO ) );

   std::string v = "KO";

   const double gap = declared_gap( slvrLP , rtrnLP ) +
                      declared_gap( slvrNDO , rtrnNDO );
   const double tol = 2e-7 * max( double( 1 ) , abs( max( foLP , foNDO ) ) )
                      + gap;

   if( hsLP && hsNDO && ( abs( foLP - foNDO ) <= tol ) )
    v = ( gap > 0 ) ? "OK(f~)" : "OK(f)";

   if( ( v == "KO" ) && hsLP && ( rtrnNDO == Solver::kUnbounded ) ) {
    /* Weird case: the LP found an optimal solution but the NDO declared the
     * problem unbounded below, because the tentative lb is too high; if so,
     * accept the run and lower (double) the bound. */
    if( ( convex && ( foNDO <= bound * ( 1 + 1e-9 ) ) ) ||
        ( ( ! convex ) && ( foNDO >= bound * ( 1 + 1e-9 ) ) ) ) {
     dbl_bound = true;
     v = "OK(?bound?)";
     }
    }

   if( ( v == "KO" ) && ( rtrnLP == Solver::kInfeasible ) &&
       ( rtrnNDO == Solver::kInfeasible ) )
    v = "OK(?e?)";

   if( ( v == "KO" ) && ( rtrnLP == Solver::kUnbounded ) &&
       ( rtrnNDO == Solver::kUnbounded ) )
    v = "OK(u)";

   if( ( v == "KO" ) && ( rtrnLP == Solver::kUnbounded ) &&
       ( ( rtrnNDO == Solver::kStopIter ) ||
         ( rtrnNDO == Solver::kStopTime ) ) ) {
    /* The LP has been proved unbounded, and the NDO Solver has run out of
     * iterations or of time on its way down, which is what it does when the
     * objective decreases so slowly along the unbounded direction that the
     * tentative bound cannot be crossed in time. The two agree if the
     * objective of NDOBlock improves along the ray the LP gives as the
     * certificate [see NDO_improves_along_LP_ray()]. */
    if( NDO_improves_along_LP_ray( slvrLP ) )
     v = "OK(u~)";
    }

   // the verdict of the line is the first KO, or that of the last Solver
   if( v == "KO" )
    ok = false;
   if( verdict.empty() || ( verdict.compare( 0 , 2 , "KO" ) != 0 ) )
    verdict = v;
   }

  if( dbl_bound ) {
   bound *= 2;
   if( convex )
    NDOBlock->set_valid_lower_bound( -bound );
   else
    NDOBlock->set_valid_upper_bound( bound );
   }

  print_instance_line( times , values ,
                       std::numeric_limits< double >::quiet_NaN() , verdict );
  return( ok );
  }
 catch( exception &e ) {
  cerr << e.what() << endl;
  exit( 1 );
  }
 catch(...) {
  cerr << "Error: unknown exception thrown" << endl;
  exit( 1 );
  }
 }  // end( SolveEach )

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

int main( int argc , char **argv )
{
 // override the default terminate handler to print the exception message
 std::set_terminate( smspp_terminate );

 // reading command line parameters - - - - - - - - - - - - - - - - - - - - -
 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 assert( SKIP_BEAT >= 0 );

 round_opts.wchg = 511;
 static int nf = 1;
 static int nt = 1;
 static double dens = 3;

 docopt_desc = "SMS++ LagBFunction test; -S gives the BlockSolverConfig of\n"
               "the NDOBlock [NDOPar.txt, NDOPar-Hard.txt if nt > 0].\n";
 add_round_options(
  "                                    bit 0 add rows, bit 1 delete rows,\n"
  "                                    bit 2 modify rows, bit 3 modify\n"
  "                                    constants, bit 4 change global\n"
  "                                    lower/upper bound, bit 5 modify\n"
  "                                    costs, bit 6 modify demands,\n"
  "                                    bit 7 modify flow bounds, bit 8\n"
  "                                    change linear objective, bit 9\n"
  "                                    each LagBFunction dualises a\n"
  "                                    random subset, bit 10 remove one\n"
  "                                    dual_pair from a LagBFunction,\n"
  "                                    bit 11 add one dual_pair to a\n"
  "                                    LagBFunction, bit 12 the inner\n"
  "                                    Block of an easy LagBFunction\n"
  "                                    issues an NBModification\n"
 #if DYNAMIC_VARS > 0
  "                                    (with DYNAMIC_VARS, bit 7 add\n"
  "                                    variables, bit 8 delete variables)\n"
 #endif
  );
 short_opts += "N:f:t:d:z:X:";
 const std::vector< option > my_opts = {
   { "nvar" , required_argument , nullptr , 'N' } ,
   { "nf" , required_argument , nullptr , 'f' } ,
   { "nt" , required_argument , nullptr , 't' } ,
   { "dens" , required_argument , nullptr , 'd' } ,
   { "size" , required_argument , nullptr , 'z' } ,
   { "detach" , required_argument , nullptr , 'X' } };
 long_opts.insert( std::prev( long_opts.end() ) ,
                   my_opts.begin() , my_opts.end() );
 help +=
  "  -N, --nvar <n>                  number of variables [10]\n"
  "  -f, --nf <n>                    |nf| PolyhedralFunction (< 0:\n"
  "                                  linear function) [1]\n"
  "  -t, --nt <n>                    |nt| transportation problems\n"
  "                                  (< 0: easy components) [1]\n"
  "  -d, --dens <x>                  rows / variables [3]\n"
  "  -z, --size <0|1|2>              size Variable of the easy\n"
  "                                  components: 0 none, 1 owned,\n"
  "                                  2 given [0]\n"
  "  -X, --detach <0|1>              if the Solver of NDOBlock are\n"
  "                                  detached and re-attached [0]\n"
;

 // the instance is generated, so no input file is needed
 filename_optional = true;

 process_args( argc , argv , []( int opt ) {
   if( process_round_option( opt ) )
    return( true );
   switch( opt ) {
    case( 'N' ): Str2Sthg( optarg , nvar ); return( true );
    case( 'f' ): Str2Sthg( optarg , nf ); return( true );
    case( 't' ): Str2Sthg( optarg , nt ); return( true );
    case( 'd' ): Str2Sthg( optarg , dens ); return( true );
    case( 'z' ): Str2Sthg( optarg , SizeEasy ); return( true );
    case( 'X' ): Str2Sthg( optarg , DetachNDO ); return( true );
    default: return( false );
    }
   } );

 long int seed = round_opts.seed;
 Index wchg = round_opts.wchg;
 Index n_repeat = round_opts.n_repeat;
 Index n_change = round_opts.n_change;
 double p_change = round_opts.p_change;
 std::string ndo_par = sconf_file;

 if( nvar < 1 ) {
  cout << "error: nvar too small";
  exit( 1 );
  }

 if( ( SizeEasy < 0 ) || ( SizeEasy > 2 ) ) {
  cout << "error: size must be 0, 1 or 2";
  exit( 1 );
  }

 bool HasLin = ( nf < 0 );
 nf = std::abs( nf );
 bool HasEasy = ( nt < 0 );
 nt = std::abs( nt );

 if( ndo_par.empty() )
  ndo_par = ( nt && ( ! HasEasy ) ) ? "NDOPar-Hard.txt" : "NDOPar.txt";

 // wchg bit 9: each LagBFunction dualises only a random subset (about
 // half) of the NDOBlock x variables, instead of all of them. Different
 // LagBFunction get independent random subsets, so the union LamVcblr
 // covers all nvar x's but each v_c05f[ h ] exposes a strict subset (in
 // ascending order, as required by MPSolver::SetItemBse), which makes
 // the BundleSolver::set_Block auto-detect engage the sparse Lambda
 // path. To keep the LPBlock-vs-NDOBlock comparison mathematically
 // sound, the corresponding (i,j) constraint in the LPBlock
 // transportation block omits the xLP[ j ] term whenever j is not in
 // the subset: this mirrors the fact that x_NDO[ j ] does not shift
 // column j's cost in the inner LP of that LagBFunction.
 bool sparse_subset = ( wchg & ( 1 << 9 ) );

 if( ( ! nf ) && ( ! nt ) ) {
  cout << "error: no sub-Block";
  exit( 1 );
  }

 #if DYNAMIC_VARS > 0
  nsvar = nvar / 2;      // half of the variables are dynamic
  ndvar = nvar - nsvar;  // the other half are static
 #endif

 Index m = nvar * dens;  // number of rows
 if( m < 1 ) {
  cout << "error: dens too small";
  exit( 1 );
  }

 // adjust the bound depending on the number of components
 bound *= std::max( 1 , std::abs( nf ) );

 rg.seed( seed );  // seed the pseudo-random number generator

 cout.setf( ios::scientific, ios::floatfield );
 cout << setprecision( 10 );

 // construction and loading of the objects - - - - - - - - - - - - - - - - -
 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 // choosing whether convex or concave: toss a(n unbiased, two-sided) coin
 convex = ( dis( rg ) < 0.5 );

 LPBlock = new AbstractBlock();
 LPBlock->set_name( "LPBlock" );

 NDOBlock = new AbstractBlock();
 NDOBlock->set_name( "NDOBlock" );

 // immediately set the global conditional bound on NDOBlock
 if( convex )
  NDOBlock->set_valid_lower_bound( -bound , true );
 else
  NDOBlock->set_valid_upper_bound( bound , true );

 {
  // ensure all original pointers go out of scope immediately after that
  // the construction has finished

  // construct the Variable in LPBlock
  auto xLP = new std::vector< ColVariable >( nsvar );
  VarVector LPvars( nvar );
  auto vit = LPvars.begin();
  for( auto & xi : *xLP )
   *(vit++) = & xi;
  #if DYNAMIC_VARS > 0
   auto xLPd = new std::list< ColVariable >( ndvar );
   for( auto & xi : *xLPd )
    *(vit++) = & xi;
  #endif

  // now set the Variable in LPBlock
   LPBlock->add_static_variable( *xLP , "x" );
  #if DYNAMIC_VARS > 0
   LPBlock->add_dynamic_variable( *xLPd , "d" );
  #endif

   // construct the Variable in NDOBlock
  auto xNDO = new std::vector< ColVariable >( nsvar );
  VarVector NDOvars( nvar );
  vit = NDOvars.begin();
  for( auto & xi : *xNDO )
   *(vit++) = & xi;
  #if DYNAMIC_VARS > 0
   auto xNDOd = new std::list< ColVariable >( ndvar );
   for( auto & xi : *xNDOd )
    *(vit++) = & xi;
  #endif

  // now set the Variable in NDOBlock
   NDOBlock->add_static_variable( *xNDO , "x" );
  #if DYNAMIC_VARS > 0
   NDOBlock->add_dynamic_variable( *xNDOd , "x" );
  #endif

  // construct the linear objective (in both)
  if( HasLin ) {  // if any
   A.resize( 1 );
   GenerateAi( A[ 0 ] , nvar );

   ConstructObj( LPBlock );
   ConstructObj( NDOBlock );

   #if( LOG_LEVEL >= 5 )
    cout << "L = [ ";
    for( Index j = 0 ; j < nvar ; ++j )
     cout << A[ 0 ][ j ] << " ";
    cout << "]" << endl;
   #endif
   }

  // construct the PolyhedralFunctionBlocks (in both) - - - - - - - - - - - -

  for( Index k = 0 ; k < Index( nf ) ; ++k ) {
   // construct the PolyhedralFunctionBlock
   auto PFBLPk = new PolyhedralFunctionBlock( LPBlock );
   PFBLPk->set_name( "LP-PFB_" + std::to_string( k ) );

   // pass it to LPBlock
   LPBlock->add_nested_Block( PFBLPk );

   // construct the m x nvar matrix A, the m-vector b, and the bound
   GenerateAb( m , nvar );
    GenerateBND();

   #if( LOG_LEVEL >= 5 )
    cout << "PF[ " << k << " ] = " << endl;
    printAb( A , b , rs( BND ) );
   #endif

   // pass all the data of the PolyhedralFunction
   PFBLPk->get_PolyhedralFunction().set_PolyhedralFunction( std::move( A ) ,
							    std::move( b ) ,
							    rs( BND ) ,
							    convex );
  // copy the PolyhedralFunctionBlock
   auto PFBNDOk = static_cast< p_PFB >(
		      PFBLPk->get_R3_Block( nullptr , nullptr , NDOBlock ) );
   PFBNDOk->set_name( "NDO-PFB_" + std::to_string( k ) );

   // pass it to NDOBlock
   NDOBlock->add_nested_Block( PFBNDOk );

   // pass the Variable to the PolyhedralFunction in LP (copy the vector)
   PFBLPk->get_PolyhedralFunction().set_variables( VarVector( LPvars ) );

   // pass the Variable to the PolyhedralFunction in NDO (copy the vector)
   PFBNDOk->get_PolyhedralFunction().set_variables( VarVector( NDOvars ) );

   // configure the PolyhedralFunctionBlock in LPBlock to use the
   // "linearised" representation
   auto bc = new BlockConfig();
   bc->f_static_variables_Configuration = new SimpleConfiguration< int >( 1 );
   PFBLPk->set_BlockConfig( bc );

   }  // end( for( k ) )

  // now construct the transportation problems (and their duals)- - - - - - -

  // first allocate C and U memory once and for all
  C.resize( nvar );
  for( auto & Ci : C )
   Ci.resize( nvar );
  U.resize( nvar );
  for( auto & Ui : U )
   Ui.resize( nvar );

  // if LagBFunctions are treated as not-easy, load once and for all the
  // ComputeConfig containing the BlockSolverConfig that will be used to
  // have the appropriate Solver attached to them
  auto cfg = Configuration::deserialize( ( nt && ( ! HasEasy ) ) ?
					 "LBFTPPar-Hard.txt" :
					 "LBFTPPar-Easy.txt" );
  auto LBFC = dynamic_cast< ComputeConfig * >( cfg );
  if( ! LBFC ) {
   cout << "error loading Configuration file for "
	<< ( ( nt && ( ! HasEasy ) ) ? "hard" : "easy" )
	<< " LagBFunction" << endl;
   delete( cfg );
   exit( 1 );
   }

  for( Index p = 0 ; p < Index( nt ) ; ++p ) {
   // for all transportation sub-Block- - - - - - - - - - - - - - - - - - - -
   Index nzc = GenerateCosts();       // generate random costs
   GenerateSupplies();                // generate random supplies == demands
   Index nic = GenerateCapacities();  // generate random capacities

   // sparse-subset support: when wchg bit 9 is on, choose a random
   // ascending subset of [ 0 , nvar ) of size ~ nvar / 2 for this
   // LagBFunction; otherwise the subset is the full [ 0 , nvar ). The
   // subset is used on both the LPBlock side (to omit the corresponding
   // xLP[ j ] term from the (i, j) constraints) and the NDOBlock side
   // (to omit the dual_pair (xNDO[ j ], cost-shift_j)). LagBFunction p
   // ends up dualising only the x_j with j in subset[ p ]; the others
   // do not shift its inner-LP cost in either representation, which
   // keeps the LPBlock-vs-NDOBlock comparison mathematically sound.
   Subset subset;
   if( sparse_subset )
    subset = GenerateSubset( nvar , std::max( Index( 1 ) , Index( nvar / 2 ) ) );
   else {
    subset.resize( nvar );
    std::iota( subset.begin() , subset.end() , Index( 0 ) );
    }
   std::vector< bool > in_subset( nvar , false );
   for( auto j : subset )
    in_subset[ j ] = true;

   #if( LOG_LEVEL >= 5 )
    cout << "T[ " << p << " ] = " << endl;
    printT();
   #endif

   // construct the AbstractBlock for LP transportation - - - - - - - - - - -
   auto TLPp = new AbstractBlock( LPBlock );
   TLPp->set_name( "LP-TB_" + std::to_string( p ) );

   // construct the (dual) variables ys and yd
   auto ys = new std::vector< ColVariable >( nvar );
   auto yd = new std::vector< ColVariable >( nvar );

   // construct the (dual) objective
   v_coeff_pair docf( 2 * nvar + nic );

   for( Index i = 0 ; i < nvar ; ++i )
    docf[ i ] = std::make_pair( & (*ys)[ i ] , s[ i ] );

   for( Index j = 0 ; j < nvar ; ++j )
    docf[ nvar + j ] = std::make_pair( & (*yd)[ j ] , s[ j ] );

   // construct and initialise the (potential) constraints
   auto pc = new boost::multi_array< FRowConstraint , 2 >(
					    boost::extents[ nvar ][ nvar ] );

   // ... that link with the variables in the root LPBlock
   auto xLP = LPBlock->get_static_variable_v< ColVariable >( "x" );

   if( nic ) {
    // construct the (dual) variables w
    #if DYNAMIC_w
     auto w = new std::list< ColVariable >( nic );
    #else
     auto w = new std::vector< ColVariable >( nic );
    #endif
    auto wit = w->begin();
    auto docfit = docf.begin() + 2 * nvar;

    for( Index i = 0 ; i < nvar ; ++i )
     for( Index j = 0 ; j < nvar ; ++j ) {
      (*pc)[ i ][ j ].set_lhs( convex ? C[ i ][ j ] : -INF , eNoMod );
      (*pc)[ i ][ j ].set_rhs( convex ? INF : C[ i ][ j ] , eNoMod );

      // when j is in the subset for this LBF the (i, j) constraint
      // carries the -x[ j ] Lagrangian-shift term; when j is outside
      // the subset the term is omitted (this LBF does not dualise
      // x_j, so its inner LP cost for column j is the unshifted c)
      const bool with_x = in_subset[ j ];
      v_coeff_pair cf;
      if( U[ i ][ j ] == INF )
       // no bound ==> ys[ i ] + yd[ j ] - x[ j ] >= c[ i ][ j ]
       cf.resize( with_x ? 3 : 2 );
      else {
       cf.resize( with_x ? 4 : 3 );
       // bound ==> ys[ i ] + yd[ j ] + w[ i ][ j ] - x[ j ] >= c[ i ][ j ]
       wit->set_type( ColVariable::kNonNegative , eNoMod );
       cf[ with_x ? 3 : 2 ] = std::make_pair( & (*wit) , convex ? 1 : -1 );
       *(docfit++) = std::make_pair( & (*(wit++)) , convex ?   U[ i ][ j ]
				                           : - U[ i ][ j ] );
       }

     cf[ 0 ] = std::make_pair( & (*ys)[ i ] , 1 );
     cf[ 1 ] = std::make_pair( & (*yd)[ j ] , 1 );
     if( with_x )
      cf[ 2 ] = std::make_pair( & (*xLP)[ j ] , -1 );

     (*pc)[ i ][ j ].set_function( new LinearFunction( std::move( cf ) ) );
     }

    // pass the (dual) variables w to the AbstractBlock
    #if DYNAMIC_w
     TLPp->add_dynamic_variable( *w , "w" );
    #else
     TLPp->add_static_variable( *w , "w" );
    #endif
    }
   else
    for( Index i = 0 ; i < nvar ; ++i )
     for( Index j = 0 ; j < nvar ; ++j ) {
      // construct constraint ys[ i ] + yd[ j ] - x[ j ] >= c[ i ][ j ]
      // (-x[ j ] is dropped when j is not in this LBF's subset; see
      // the comment in the cap branch above)
      (*pc)[ i ][ j ].set_lhs( convex ? C[ i ][ j ] : -INF , eNoMod );
      (*pc)[ i ][ j ].set_rhs( convex ? INF : C[ i ][ j ] , eNoMod );

      const bool with_x = in_subset[ j ];
      v_coeff_pair cf( with_x ? 3 : 2 );
      cf[ 0 ] = std::make_pair( & (*ys)[ i ] , 1 );
      cf[ 1 ] = std::make_pair( & (*yd)[ j ] , 1 );
      if( with_x )
       cf[ 2 ] = std::make_pair( & (*xLP)[ j ] , -1 );

      (*pc)[ i ][ j ].set_function( new LinearFunction( std::move( cf ) ) );
      }

   // pass the (dual) variables to the AbstractBlock
   TLPp->add_static_variable( *ys , "ys" );
   TLPp->add_static_variable( *yd , "yd" );

   // pass the (potential) constraints to the AbstractBlock
   TLPp->add_static_constraint( *pc , "pc" );

   // construct the (dual) objective
   auto dobj = new FRealObjective( TLPp ,
				   new LinearFunction( std::move( docf ) ) );

   // set the right optimization sense (min if convex, max otherwise)
   dobj->set_sense( convex ? Objective::eMin : Objective::eMax , eNoMod );
   
   // pass the (dual) objective to the AbstractBlock
   TLPp->set_objective( dobj , eNoMod );

   // finally, pass the AbstractBlock to LPBlock
   LPBlock->add_nested_Block( TLPp );

   // construct the AbstractBlock for NDO transportation- - - - - - - - - - -
   // the AbstractBlock just has a FRealObjective with a LagBFunction inside;
   // the inner Block of the LagBFunction is the transportation probem
   
   auto TNDOp = new AbstractBlock( NDOBlock );
   TNDOp->set_name( "NDO-TB_" + std::to_string( p ) );

   // construct the inner Block for the LagBFunction, one that takes its
   // size Variable if it is to be given one [see SizeEasy]
   const bool own_size = HasEasy && ( SizeEasy == 1 );
   const bool row_caps = HasEasy && ( SizeEasy > 0 );
   p_AB IBNDOp = ( HasEasy && ( SizeEasy == 2 ) )
                 ? new SizedTransportBlock() : new AbstractBlock();
   IBNDOp->set_name( "IB-NDO-TB_" + std::to_string( p ) );

   // the size Variable owned by the inner Block, if any, which the master
   // ties to lambda
   ColVariable * tau_k = nullptr;
   if( own_size ) {
    tau_k = new ColVariable();

    tau_k->set_type( ColVariable::kNonNegative , eNoMod );
    tau_k->set_value( 1.0 );  // initial value only; tau_k remains free

    IBNDOp->add_static_variable( *tau_k , "tau" );
    IBNDOp->set_owned_size_variable( tau_k );
    }

   // construct the flow variables
   auto f = new boost::multi_array< ColVariable , 2 >(
				            boost::extents[ nvar ][ nvar ] );

   // set the sign of the flow variables
   for( Index i = 0 ; i < nvar ; ++i )
    for( Index j = 0 ; j < nvar ; ++j )
     (*f)[ i ][ j ].set_type( ColVariable::kNonNegative , eNoMod );

   // pass the flow variables to the inner Block
   IBNDOp->add_static_variable( *f , "f" );

   // construct the source constraints
   auto sc = new std::vector< FRowConstraint >( nvar );

   // initialize the source constraints
   for( Index i = 0 ; i < nvar ; ++i ) {
    // sum_j f[i][j] == s[i], or sum_j f[i][j] - s[i] * tau_k == 0
    auto lf = new LinearFunction();

    for( Index j = 0 ; j < nvar ; ++j )
     lf->add_variable( &(*f)[ i ][ j ] , 1.0 , eNoMod );

    if( own_size ) {
     lf->add_variable( tau_k , -s[ i ] , eNoMod );
     (*sc)[ i ].set_both( 0.0 , eNoMod );
     }
    else
     (*sc)[ i ].set_both( s[ i ] , eNoMod );

    (*sc)[ i ].set_function( lf , eNoMod );
    }

   // pass the source constraints to the inner Block
   IBNDOp->add_static_constraint( *sc , "sc" );

   // construct the destination constraints
   auto dc = new std::vector< FRowConstraint >( nvar );
 
   // initialize the destination constraints
   for( Index j = 0 ; j < nvar ; ++j ) {
    // sum_i f[i][j] == s[j], or sum_i f[i][j] - s[j] * tau_k == 0
    auto lf = new LinearFunction();

    for( Index i = 0 ; i < nvar ; ++i )
     lf->add_variable( &(*f)[ i ][ j ] , 1.0 , eNoMod );

    if( own_size ) {
     lf->add_variable( tau_k , -s[ j ] , eNoMod );
     (*dc)[ j ].set_both( 0.0 , eNoMod );
     }
    else
     (*dc)[ j ].set_both( s[ j ] , eNoMod );

    (*dc)[ j ].set_function( lf , eNoMod );
    }

   // pass the destination constraints to the inner Block
   IBNDOp->add_static_constraint( *dc , "dc" );

   if( nic ) {  // construct the finite capacity constraints (if any)
    if( row_caps ) {
     // f[i][j] <= U[i][j], or f[i][j] - U[i][j] * tau_k <= 0: a row, so
     // that it can take the size Variable
     #if DYNAMIC_bc
      auto bc = new std::list< FRowConstraint >( nic );
     #else
      auto bc = new std::vector< FRowConstraint >( nic );
     #endif
     auto bcit = bc->begin();

     for( Index i = 0 ; i < nvar ; ++i )
      for( Index j = 0 ; j < nvar ; ++j )
       if( U[ i ][ j ] < INF ) {
        auto lf = new LinearFunction();
        lf->add_variable( &(*f)[ i ][ j ] , 1.0 , eNoMod );
        if( own_size )
         lf->add_variable( tau_k , -U[ i ][ j ] , eNoMod );
        bcit->set_function( lf , eNoMod );
        bcit->set_lhs( -INF , eNoMod );
        (bcit++)->set_rhs( own_size ? 0.0 : U[ i ][ j ] , eNoMod );
        }

     #if DYNAMIC_bc
      IBNDOp->add_dynamic_constraint( *bc , "bc" );
     #else
      IBNDOp->add_static_constraint( *bc , "bc" );
     #endif
     }
    else {
     #if DYNAMIC_bc
      auto bc = new std::list< BoxConstraint >( nic );
     #else
      auto bc = new std::vector< BoxConstraint >( nic );
     #endif
     auto bcit = bc->begin();

     for( Index i = 0 ; i < nvar ; ++i )
      for( Index j = 0 ; j < nvar ; ++j )
       if( U[ i ][ j ] < INF ) {
        bcit->set_variable( & (*f)[ i ][ j ] , eNoMod );
        (bcit++)->set_rhs( U[ i ][ j ] , eNoMod );
        }

     #if DYNAMIC_bc
      IBNDOp->add_dynamic_constraint( *bc , "bc" );
     #else
      IBNDOp->add_static_constraint( *bc , "bc" );
     #endif
     }
    }
 
   // construct the objective
   v_coeff_pair ocf( nzc );

   auto it = ocf.begin();
   for( Index i = 0 ; i < nvar ; ++i )
    for( Index j = 0 ; j < nvar ; ++j )
     if( C[ i ][ j ] ) {
      it->first = & (*f)[ i ][ j ];
      (it++)->second = C[ i ][ j ];
      }

   assert( it == ocf.end() );

   // construct the objective of the inner Block
   auto ibo = new FRealObjective( IBNDOp ,
				  new LinearFunction( std::move( ocf ) ) );
   ibo->set_sense( Objective::eMax );  // ... to be *max*imized

   // set the right optimization sense (max if convex, min otherwise)
   ibo->set_sense( convex ? Objective::eMax : Objective::eMin , eNoMod );

   // pass the objective to the inner Block
   IBNDOp->set_objective( ibo , eNoMod );

   // construct the LagBFunction, passing it the inner Block
   auto LBF = new LagBFunction( IBNDOp );

   LBF->set_ComputeConfig( LBFC );  // ComputeConfig-ure it

   // construct the dual pairs only for j in this LBF's subset; if the
   // sparse-subset flag is off the subset is [ 0 , nvar ) and lp has
   // size nvar, matching the legacy dense behaviour. When the flag is
   // on, lp has size subset.size() < nvar and the local-to-global map
   // built by BundleSolver::set_Block will be a strict subset of
   // LamVcblr (in ascending order), which triggers the sparse path.
   LagBFunction::v_dual_pair lp;
   lp.reserve( subset.size() );

   // ... that use the variables in NDOBlock
   auto xNDO = NDOBlock->get_static_variable_v< ColVariable >( 0 );

   for( auto j : subset ) {
    v_coeff_pair cfj( nvar );

    for( Index i = 0 ; i < nvar ; ++i ) {
     cfj[ i ].first = & (*f)[ i ][ j ];
     cfj[ i ].second = 1;
     }

    lp.emplace_back( & (*xNDO)[ j ] , new LinearFunction( std::move( cfj ) ) );
    }

   // pass the dual pairs to the LagBFunction
   LBF->set_dual_pairs( std::move( lp ) );

   // construct the objective to the transportation Block, passing it the
   // LagBFunction
   auto tobj = new FRealObjective( TNDOp , LBF );

   // set the right optimization sense (min if convex, max otherwise)
   tobj->set_sense( convex ? Objective::eMin : Objective::eMax , eNoMod );

   // pass the objective to the transportation Block
   TNDOp->set_objective( tobj , eNoMod );

   // finally, pass the transportation Block to the NDOBlock
   NDOBlock->add_nested_Block( TNDOp );

   }  // end( for( p ) )

  delete( LBFC );
  }

 // define bound constraints- - - - - - - - - - - - - - - - - - - - - - - - -

 #if HAVE_CONSTRAINTS == 1
 {
  auto LPx = LPBlock->get_static_variable_v< ColVariable >( "x" );
  auto NDOx = NDOBlock->get_static_variable_v< ColVariable >( "x" );
  for( Index i = 0 ; i < nvar ; ++i )
   if( dis( rg ) < 0.5 ) {
    (*LPx)[ i ].is_positive( true , eNoMod );
    (*NDOx)[ i ].is_positive( true , eNoMod );
    }
  }
 #endif
 #if HAVE_CONSTRAINTS == 2
 {
  auto LPx = LPBlock->get_static_variable_v< ColVariable >( "x" );
  auto NDOx = NDOBlock->get_static_variable_v< ColVariable >( "x" );
  auto LPbnd = new std::list< BoxConstraint >;
  auto NDObnd = new std::list< BoxConstraint >;
  for( Index i = 0 ; i < nvar ; ++i )
   if( dis( rg ) < 0.5 ) {
    LPbnd->resize( LPbnd->size() + 1 );
    NDObnd->resize( NDObnd->size() + 1 );
    LPbnd->back().set_variable( & (*LPx)[ i ] );
    NDObnd->back().set_variable( & (*NDOx)[ i ] );
    auto p = dis( rg );
    auto lhs = p < 0.666 ? 0 : -INF;
    auto rhs = p > 0.333 ? dis( rg ) : INF;
    LPbnd->back().set_lhs( lhs , eNoMod );
    NDObnd->back().set_lhs( lhs , eNoMod );
    LPbnd->back().set_rhs( rhs , eNoMod );
    NDObnd->back().set_rhs( rhs , eNoMod );
    }
   else
    if( dis( rg ) < 0.5 ) {
     (*LPx)[ i ].is_positive( true , eNoMod );
     (*NDOx)[ i ].is_positive( true , eNoMod );
     }

  LPBlock->add_dynamic_constraint( *LPbnd );
  NDOBlock->add_dynamic_constraint( *NDObnd );
  }
 #endif

 // final setup- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 LPBlock->generate_abstract_variables();
 LPBlock->generate_abstract_constraints();
 LPBlock->generate_objective();

 NDOBlock->generate_abstract_variables();
 NDOBlock->generate_abstract_constraints();
 NDOBlock->generate_objective();

 // final checks- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 /**!!
 LPBlock->is_correct();
 NDOBlock->is_correct();
 !!*/

 // attach the Solver to the Block- - - - - - - - - - - - - - - - - - - - - -
 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 {
  // for LPBlock do this by reading an appropriate BlockSolverConfig from
  // file and apply() it to the LPBlock; 
  // BSC may be a plain BlockSolverConfig or a meta-config
  // SimpleConfiguration< std::map< std::string , Configuration * > >;
  // s_config_Block() dispatches on the runtime type and clears the config(s).
  auto lsc = Configuration::deserialize( "LPPar.txt" );
  if( ! lsc ) {
   cerr << "Error: cannot load BSC from LPPar.txt" << endl;
   return( 1 );
   }
  s_config_Block( LPBlock , lsc , "LPPar.txt" );
  delete( lsc );

  // furthermore, "manually" attach an UpdateSolver to (each
  // PolyhedralFunctionBlock in) LPBlock
  for( Index i = 0 ; i < Index( nf ) ; ++i )
   LPBlock->get_nested_Block( i )->register_Solver(
		       new UpdateSolver( NDOBlock->get_nested_Block( i ) ) );
  }

 {
  // for NDOBlock do this by reading the BlockSolverConfig, or the
  // meta-config, of ndo_par and apply() it to the NDOBlock; whether the
  // Solver treat the LagBFunction as easy components, and which of their
  // dual values they keep, is in that file (intDoEasy of BundleSolver),
  // which must agree with the sign of #nt
  auto cfg = Configuration::deserialize( ndo_par );
  if( ! cfg ) {
   cerr << "Error: cannot load BSC from " << ndo_par << endl;
   return( 1 );
   }
  s_config_Block( NDOBlock , cfg , ndo_par );
  delete( cfg );

  #if( LOG_LEVEL >= 4 )
   // in the extremely verbose mode, set an event that spits out the LPs
   // in the LagBFunctions every k iterations; k must be in the parameter
   // file via intEverykIt (0 by default == never)

   if( ! HasEasy ) {  // transportation problems are treated as "difficult"
    NDOBlock->get_registered_solvers().front()->set_event_handler(
     ThinComputeInterface::eEverykIteration ,
     [ & ] () {
      // have the inner TP Solver spit out the LP at this iteration;
      // this implies that a :MILPSolver is used, otherwise the Solver
      // will complain, but note that by the clever use of str_par_str2idx()
      // one does not need to include MILPSolver.h
      for( Index p = nf ; p < Index( nf + nt ) ; ++p ) {
       auto FRO =
	 NDOBlock->get_nested_Block( p )->get_objective< FRealObjective >();
       auto LBF = static_cast< p_LBF >( FRO->get_function() );
       auto inner = LBF->get_inner_block();
       if( ! inner )
	throw( std::logic_error( "LagBFunction without inner Block" ) );
       auto slv = inner->get_registered_solvers().front();
       slv->set_par( slv->str_par_str2idx( "strOutputFile" ) ,
		     "TB-" + std::to_string( p - nf ) + "-" +
		     std::to_string( slv->get_elapsed_calls() ) + "-" +
		     std::to_string( slv->get_elapsed_iterations() ) +
		     ".lp" );
       }

      return( ThinComputeInterface::eContinue );
      }  // end of lambda
								     );
    }
  #endif
  }

 // open log-file - - - - - - - - - - -  - - - - - - - - - - - - - - - - - -
 //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 #if( LOG_LEVEL >= 2 )
  #if( LOG_ON_COUT )
   for( auto slvr : NDOBlock->get_registered_solvers() )
    slvr->set_log( &cout );
  #else
   ofstream LOGFile( logF , ofstream::out );
   if( ! LOGFile.is_open() )
    cerr << "Warning: cannot open log file """ << logF << """" << endl;
   else {
    LOGFile.setf( ios::scientific, ios::floatfield );
    LOGFile << setprecision( 10 );
    for( auto slvr : NDOBlock->get_registered_solvers() )
     slvr->set_log( &LOGFile );
    }
  #endif

  #if( LOG_LEVEL >= 3 )
  {
   auto LPslv = ( LPBlock->get_registered_solvers() ).front();
   LPslv->set_par( LPslv->str_par_str2idx( "strOutputFile" ) ,
		   "LPBlock.lp" );
   }
  #endif
 #endif

 // first solver call - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 LOG1( "First call: " );

 bool AllPassed = SolveEach();
 
 // main loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 // now, for n_repeat times, one inner Block is selected and:
 // - if it is a PolyhedralFunctionBlock
 //   = up to n_change rows are added
 //   = up to n_change rows are deleted
 //   = up to n_change rows are modified
 //   = up to n_change constants are modified
 // - if it is an AbstractBlock (transportation problem)
 //   = costs are modified
 //   = demands are modified
 //   = flow bounds are modified
 // - in either case
 //   the linear objective, if any, is modified
 //
 // then the two problems are re-solved
 //
 // IMPORTANT NOTE: only LPBlock is changed, because UpdateSolver takes
 //                 care of intercepting all (physical) Modification and
 //                 map_forward them to NDOBlock
 //
 // by playing with SKIP_BEAT one can re-solve the two problems after having
 // changed an arbitrary number of inner Blocks

 for( Index rep = 0 ; rep < n_repeat * ( SKIP_BEAT + 1 ) ; ) {
  if( ! AllPassed )
   break;

  p_AB LPTr = nullptr;
  p_AB NDOTr = nullptr;
  p_PFB LPBr = nullptr;
  p_LBF LBF = nullptr;  // the LagBFunction of NDOTr, when bn is a TB

  Index bn = rep % ( nf + nt );  // which sub-Block to change

  if( bn < Index( nf ) ) {
   LPBr = static_cast< p_PFB >( LPBlock->get_nested_Block( bn ) );
   LOG1( rep << "[PFB " << bn << "]: ");
   }
  else {
   LPTr = static_cast< p_AB >( LPBlock->get_nested_Block( bn ) );
   NDOTr = static_cast< p_AB >( NDOBlock->get_nested_Block( bn ) );
   LBF = static_cast< p_LBF >(
		NDOTr->get_objective< FRealObjective >()->get_function() );
   NDOTr = static_cast< p_AB >( LBF->get_inner_block() );
   if( ! NDOTr )
    throw( std::logic_error( "LagBFunction without inner Block" ) );
   LOG1( rep << "[TB " << bn - nf << "]: ");
   }

  // add rows - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPBr && ( wchg & 1 ) && ( dis( rg ) <= p_change ) )
   if( Index tochange = Index( dis( rg ) * n_change ) ) {
    LOG1( "added " << tochange << " rows - " );

    GenerateAb( tochange , nvar );

    auto cnst = LPBr->get_dynamic_constraint< FRowConstraint >( 0 );

    if( tochange == 1 )
     LPBr->get_PolyhedralFunction().add_row( std::move( A[ 0 ] ) , b[ 0 ] );
    else
     LPBr->get_PolyhedralFunction().add_rows( std::move( A ) , b );
    }

  // delete rows- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPBr && ( wchg & 2 ) && ( dis( rg ) <= p_change ) ) {
   m = LPBr->get_PolyhedralFunction().get_nrows();
   if( Index tochange = min( m - 1 , Index( dis( rg ) * n_change ) ) ) {
    LOG1( "deleted " << tochange << " rows" );

    auto cnst = LPBr->get_dynamic_constraint< FRowConstraint >( 0 );

    if( dis( rg ) <= 0.5 ) {  // in 50% of the cases do a ranged change
    
     Index strt = dis( rg ) * ( m - tochange );
     Index stp = strt + tochange;

     LOG1( "(r) - " );
     if( tochange == 1 )
      LPBr->get_PolyhedralFunction().delete_row( strt );
     else
      LPBr->get_PolyhedralFunction().delete_rows( Range( strt , stp ) );
     }
    else {  // in the other 50% of the cases, do a sparse change
     Subset nms = GenerateSubset( m , tochange );

     LOG1( "(s) - " );
     if( tochange == 1 )
      LPBr->get_PolyhedralFunction().delete_row( nms[ 0 ] );
     else
      LPBr->get_PolyhedralFunction().delete_rows( std::move( nms ) );
     }
    }
   }

  // modify rows- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPBr && ( wchg & 4 ) && ( dis( rg ) <= p_change ) ) {
   m = LPBr->get_PolyhedralFunction().get_nrows();
   if( Index tochange = std::min( m , Index( dis( rg ) * n_change ) ) ) {
    LOG1( "modified " << tochange << " rows" );

    GenerateAb( tochange , nvar );

    if( dis( rg ) <= 0.5 ) {  // in 50% of the cases do a ranged change
     Index strt = dis( rg ) * ( m - tochange );
     Index stp = strt + tochange;

     LOG1( "(r) - " );
     if( tochange == 1 )
      LPBr->get_PolyhedralFunction().modify_row( strt , std::move( A[ 0 ] ) ,
						 b[ 0 ] );
     else
      LPBr->get_PolyhedralFunction().modify_rows( std::move( A ) , b ,
						  Range( strt , stp ) );
     }
    else {  // in the other 50% of the cases, do a sparse change
     Subset nms = GenerateSubset( m , tochange );

     LOG1( "(s) - " );
     if( tochange == 1 )
      LPBr->get_PolyhedralFunction().modify_row( nms[ 0 ] ,
						 std::move( A[ 0 ] ) ,
						 b[ 0 ] );
     else
      LPBr->get_PolyhedralFunction().modify_rows( std::move( A ) , b ,
						  std::move( nms ) , true );
     }
    }
   }

  // modify constants - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPBr && ( wchg & 8 ) && ( dis( rg ) <= p_change ) ) {
   m = LPBr->get_PolyhedralFunction().get_nrows();
   if( Index tochange = std::min( m , Index( dis( rg ) * n_change ) ) ) {
    LOG1( "modified " << tochange << " constants" );

    Generateb( tochange );

    if( dis( rg ) <= 0.5 ) {  // in 50% of the cases do a ranged change
     Index strt = dis( rg ) * ( m - tochange );
     Index stp = strt + tochange;

     LOG1( "(r) - " );
     if( tochange == 1 )
      LPBr->get_PolyhedralFunction().modify_constant( strt , b[ 0 ] );
     else
      LPBr->get_PolyhedralFunction().modify_constants( b ,
						       Range( strt , stp ) );
     }
    else {  // in the other 50% of the cases, do a sparse change
     Subset nms = GenerateSubset( m , tochange );

     LOG1( "(s) - " );
     if( tochange == 1 )
      LPBr->get_PolyhedralFunction().modify_constant( nms[ 0 ] , b[ 0 ] );
     else
      LPBr->get_PolyhedralFunction().modify_constants( b , std::move( nms ) ,
						       true );
     }
    }
   }

  // modify bound - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPBr && ( wchg & 16 ) && ( dis( rg ) <= p_change ) ) {
   LOG1( "modified bound - " );

   GenerateBND();

   LPBr->get_PolyhedralFunction().modify_bound( rs( BND ) );
   }

  // change costs - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPTr && ( wchg & 32 ) && ( dis( rg ) <= p_change ) ) {
   // about 30% of the coefficients are 0: the corresponding variables were
   // not originally in the objective, but have been (hopefully) silently
   // added back to it by LagBFunction, "at the end", because all costs do
   // have a Lagrangian term. when costs are changed, there is a fair chance
   // that some of these coefficients will be changed (and made nonzero).
   // in 50% of the cases we rather eliminate (==> set the coefficient to 0)
   // a bunch of variables from the objective to test LagBFunction's
   // capability to add them back to the end

   GenerateCosts();
   #if( LOG_LEVEL >= 5 )
    printC();
   #endif

   auto lf = static_cast< p_LF >(
		  NDOTr->get_objective< FRealObjective >()->get_function() );
   auto & cf = lf->get_v_var();

   // because the order of the variables in the objective is messed up with
   // reconstruct the information about which position corresponds to what
   auto f = NDOTr->get_static_variable< ColVariable , 2 >( "f" );
   if( ! f ) {
    LOG1( "something very bad happened!" );
    exit( 1 );
    }
   std::vector< std::pair< Index , Index > > dict( cf.size() );
   for( Index h = 0 ; h < cf.size() ; ++h ) {
    bool found = false;
    for( Index i = 0 ; ( i < nvar ) && ( ! found ) ; ++i )
     for( Index j = 0 ; ( j < nvar ) && ( ! found ) ; ++j )
      if( & (*f)[ i ][ j ] == cf[ h ].first ) {
       dict[ h ].first = i;
       dict[ h ].second = j;
       found = true;
       }
    }

   if( dis( rg ) < 0.5 ) {  // in 50% of the cases, modify the costs
    // change up to 50% of the variables and *no less* than n_change
    Index tochange = std::min( Index( cf.size() ) ,
			       std::max( n_change ,
					 Index( dis( rg ) * cf.size() / 2 )
					 ) );

    LOG1( "modified " << tochange << " costs" );

    if( dis( rg ) < 0.5 ) {  // in 25% of the cases, do a ranged change
     LOG1( "(r) - " );
     Index strt = dis( rg ) * ( cf.size() - tochange );
     Index stp = strt + tochange;

     // in the dual transportation problem these are the RHS of the
     // potential constraints: send all the corresponding Modification to
     // a new channel
     auto pc = LPTr->get_static_constraint< FRowConstraint , 2 >( "pc" );

     Observer::ChnlName chnl = LPTr->open_channel();
     const auto iAM = Observer::make_par( eModBlck , chnl );
     if( convex )
      for( Index h = strt ; h < stp ; ++h )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_lhs(
			    C[ dict[ h ].first ][ dict[ h ].second ] , iAM );
     else
      for( Index h = strt ; h < stp ; ++h )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_rhs(
			    C[ dict[ h ].first ][ dict[ h ].second ] , iAM );

     LPTr->close_channel( chnl );  // then close the chanel

     // in the transportation problem inside the LagBFunction, just do it
     RealVector newc( tochange );
     for( Index h = strt ; h < stp ; ++h )
      newc[ h - strt ] = C[ dict[ h ].first ][ dict[ h ].second ];

     lf->modify_coefficients( std::move( newc ) , Range( strt , stp ) );
     }
    else {                   // in 25% of the cases, do a subset change
     LOG1( "(s) - " );
     Subset nms = GenerateSubset( cf.size() , tochange );

     // in the dual transportation problem these are the RHS of the
     // potential constraints: send all the corresponding Modification to
     // a new channel
     auto pc = LPTr->get_static_constraint< FRowConstraint , 2 >( "pc" );

     Observer::ChnlName chnl = LPTr->open_channel();
     const auto iAM = Observer::make_par( eModBlck , chnl );
     if( convex )
      for( Index h : nms )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_lhs(
                            C[ dict[ h ].first ][ dict[ h ].second ] , iAM );
     else
      for( Index h : nms )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_rhs(
                            C[ dict[ h ].first ][ dict[ h ].second ] , iAM );

     LPTr->close_channel( chnl );  // then close the chanel

     // in the transportation problem inside the LagBFunction, just do it
     RealVector newc( tochange );
     for( Index h = 0 ; h < tochange ; ++h )
      newc[ h ] = C[ dict[ nms[ h ] ].first ][ dict[ nms[ h ] ].second ];

     lf->modify_coefficients( std::move( newc ) , std::move( nms ) );
     }
    }
   else {                    // in 50% of the cases, eliminate coefficients
    // change up to 25% of the variables and *no less* than 1
    Index tochange = std::min( Index( cf.size() ) ,
			       std::max( Index( 1 ) ,
					 Index( dis( rg ) * cf.size() / 4 )
					 ) );

    LOG1( "deleted " << tochange << " costs" );

    if( dis( rg ) < 0.5 ) {  // in 25% of the cases, do a ranged removal
     LOG1( "(r) - " );
     Index strt = dis( rg ) * ( cf.size() - tochange );
     Index stp = strt + tochange;

     // in the dual transportation problem these are the RHS of the
     // potential constraints: send all the corresponding Modification to
     // a new channel
     auto pc = LPTr->get_static_constraint< FRowConstraint , 2 >( "pc" );

     Observer::ChnlName chnl = LPTr->open_channel();
     const auto iAM = Observer::make_par( eModBlck , chnl );
     if( convex )
      for( Index h = strt ; h < stp ; ++h )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_lhs( 0 , iAM );
     else
      for( Index h = strt ; h < stp ; ++h )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_rhs( 0 , iAM );

     LPTr->close_channel( chnl );  // then close the chanel

     // in the transportation problem inside the LagBFunction, just do it
     lf->remove_variables( Range( strt , stp ) );
     }
    else {                   // in 25% of the cases, do a subset removal
     LOG1( "(s) - " );
     Subset nms = GenerateSubset( cf.size() , tochange );

     // in the dual transportation problem these are the RHS of the
     // potential constraints: send all the corresponding Modification to
     // a new channel
     auto pc = LPTr->get_static_constraint< FRowConstraint , 2 >( "pc" );

     Observer::ChnlName chnl = LPTr->open_channel();
     const auto iAM = Observer::make_par( eModBlck , chnl );
     if( convex )
      for( Index h : nms )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_lhs( 0 , iAM );
     else
      for( Index h : nms )
       (*pc)[ dict[ h ].first ][ dict[ h ].second ].set_rhs( 0 , iAM );

     LPTr->close_channel( chnl );  // then close the chanel

     // in the transportation problem inside the LagBFunction, just do it
     lf->remove_variables( std::move( nms ) );
     }
    }
   }

  // change demands - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPTr && ( wchg & 64 ) && ( dis( rg ) <= p_change ) ) {
   LOG1( "modified demands/supplies - " );

   GenerateSupplies();
   #if( LOG_LEVEL >= 5 )
    cout << endl << "s = [ ";
    for( auto el : s )
     cout << el << " ";
    cout << " ]" << endl;
   #endif

   // in the transportation problem inside the LagBFunction these are the
   // RHS of the demand and supply constraints, or the coefficients of
   // tau_k when scaled: send all the corresponding Modification to a
   // new channel
   Observer::ChnlName chnl = NDOTr->open_channel();
   const auto iAM = Observer::make_par( eModBlck , chnl );

   auto tau = size_variable_of( NDOTr );
   auto update_balances = [&]( std::vector< FRowConstraint > * rows ) {
    if( tau ) {
     for( Index j = 0 ; j < nvar ; ++j ) {
      auto lf = static_cast< p_LF >( (*rows)[ j ].get_function() );
      lf->modify_coefficient( lf->is_active( tau ) , -s[ j ] , iAM );
      }
     return;  // the scaled rows keep both sides equal to zero
     }
    for( Index j = 0 ; j < nvar ; ++j )
     (*rows)[ j ].set_both( s[ j ] , iAM );
    };

   update_balances( NDOTr->get_static_constraint_v< FRowConstraint >( "dc" ) );
   update_balances( NDOTr->get_static_constraint_v< FRowConstraint >( "sc" ) );

   NDOTr->close_channel( chnl );  // then close the chanel

   // in the dual transportation problem these are a slice of the coefficients
   // of the objective
   s.resize( 2 * nvar );
   std::copy( s.begin() , s.begin() + nvar , s.begin() + nvar );

   auto lf = static_cast< p_LF >(
		  LPTr->get_objective< FRealObjective >()->get_function() );
   lf->modify_coefficients( std::move( s ) , Range( 0 , 2 * nvar ) ); 
   }

  // change flow bounds - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( LPTr && ( wchg & 128 ) && ( dis( rg ) <= p_change ) ) {
   // The same modification sequence applies to either representation and
   // to both static vectors and dynamic lists of capacity constraints.
   auto change_bounds = [&]( auto * bc ) {
    if( ! bc )
     return;
    const Index tochange = std::min( Index( bc->size() ) ,
                                    Index( dis( rg ) * nvar * nvar / 7 ) );
    if( ! tochange )
     return;
     LOG1( "changed " << tochange << " flow bounds" );

     RealVector tmpU( tochange );
     GenerateCapacities( tmpU );

     auto lf = static_cast< p_LF >(
		  LPTr->get_objective< FRealObjective >()->get_function() );

     // In the inner transportation problem these are numeric box bounds,
     // or coefficients of tau_k in scaled rows. Send all the corresponding
     // Modification to a new channel
     Observer::ChnlName chnl = NDOTr->open_channel();
     const auto iAM = Observer::make_par( eModBlck , chnl );

     // important note: due to a limitation in OSIMPSolver, the "finite
     // nonzero" status of the bounds need be preserved: if the bound is
     // finite and nonzero it must remain so and vice-versa, although
     // a zero bound could in principle become +INF and vice-versa; yet,
     // since we can only change finite bounds (for otherwise we'd need to
     // create new dual variables in the dual), this basically means that
     // if a bound is zero it must remain so, and if it is nonzero it must
     // remain so

     auto tau = size_variable_of( NDOTr );
     auto update_bound = [&]( RowConstraint & bound , double & value ) {
      if( tau ) {
       auto blf = static_cast< p_LF >(
                             static_cast< FRowConstraint & >( bound ).
                             get_function() );
       const Index ti = blf->is_active( tau );
       if( blf->get_coefficient( ti ) != 0 ) {
        while( value == 0 )
         value = 5 * dis( rg );
        blf->modify_coefficient( ti , -value , iAM );
        }
       else
        value = 0;
       return;  // the scaled capacity row keeps its RHS at zero
       }
      if( bound.get_rhs() != 0 ) {
       while( value == 0 )
        value = 5 * dis( rg );
       bound.set_rhs( value , iAM );
       }
      else
       value = 0;
      };

     if( dis( rg ) < 0.5 ) {  // in 50% of the cases, do a ranged change
      LOG1( "(r) - " );
      Index strt = dis( rg ) * ( bc->size() - tochange );
      Index stp = strt + tochange;

      auto bcit = std::next( bc->begin() , strt );
      for( Index h = 0 ; h < tochange ; ++h , ++bcit )
       update_bound( *bcit , tmpU[ h ] );

      NDOTr->close_channel( chnl );  // then close the chanel

      // in the transportation problem, just do it
      if( ! convex )
       for( auto & tUi : tmpU )
	tUi = - tUi;

      lf->modify_coefficients( std::move( tmpU ) ,
			       Range( 2 * nvar + strt , 2 * nvar + stp ) ); 
      }
     else {                   // in 50% of the cases, do a subset change
      LOG1( "(s) - " );
      Subset nms = GenerateSubset( bc->size() , tochange );

      #if DYNAMIC_bc
       Index nh = 0;
       auto bcit = bc->begin();
       for( Index h = 0 ; h < tochange ; ++h ) {
        auto nnh = nms[ h ];
        std::advance( bcit , nnh - nh );
        update_bound( *bcit , tmpU[ h ] );
        nh = nnh;
        }
      #else
       for( Index h = 0 ; h < tochange ; ++h )
        update_bound( (*bc)[ nms[ h ] ] , tmpU[ h ] );
      #endif

      NDOTr->close_channel( chnl );  // then close the chanel

      // in the transportation problem, just do it
      for( auto & el : nms )
       el += 2 * nvar;
      if( ! convex )
       for( auto & tUi : tmpU )
	tUi = - tUi;

      lf->modify_coefficients( std::move( tmpU ) , std::move( nms ) , true ); 
      }
    };

   auto change_box_bounds = [&]() {
    #if DYNAMIC_bc
     change_bounds( NDOTr->get_dynamic_constraint< BoxConstraint >( "bc" ) );
    #else
     change_bounds( NDOTr->get_static_constraint_v< BoxConstraint >( "bc" ) );
    #endif
    };

   // the capacities are rows whenever there is a size Variable to take
   if( HasEasy && SizeEasy ) {
    #if DYNAMIC_bc
     change_bounds( NDOTr->get_dynamic_constraint< FRowConstraint >( "bc" ) );
    #else
     change_bounds( NDOTr->get_static_constraint_v< FRowConstraint >( "bc" ) );
    #endif
    }
   else
    change_box_bounds();
   }

  // the inner Block "changes everything" - - - - - - - - - - - - - - - - - -
  //
  // bit 12: the inner Block of an easy LagBFunction issues the "nuclear"
  // NBModification, which says that anything may have changed and which a
  // sized copy of it in the master cannot follow [see
  // AbstractBlock::mirror_forward_Modification()]: the master puts the inner
  // Block back in place of the copy, unscaled, and fixes its lambda to 1
  // [see MasterProblemBlock::EasyMirror]; nothing has changed in fact, so
  // that the two Block stay the same problem

  if( LPTr && HasEasy && ( wchg & 4096 ) && ( dis( rg ) <= p_change ) ) {
   LOG1( "everything changed - " );
   NDOTr->add_Modification( std::make_shared< NBModification >( NDOTr ) );
   }

  // per-LagBFunction dual_pair removal - - - - - - - - - - - - - - - - - - -
  //
  // bit 10: remove one dual_pair from THIS LBF and mirror by zeroing the
  // corresponding xLP[ j ] coefficient on the LPBlock side. The Mod is
  // issued on a single LBF — i.e. naked, not part of a GroupModification
  // covering all components — to exercise:
  //
  //  (a) BundleSolver's dense → sparse auto-promotion in the 4th
  //      Modification loop, when this is the first per-Function var Mod
  //      and the Solver was in dense mode;
  //
  //  (b) the per-Function sparse FunctionModVarsRngd handler that
  //      decrements v_ref_count on the affected LamVcblr slot and
  //      queues globally-dead slots for end-of-loop compaction.
  //
  // Mathematically, removing dual_pair k from the LBF means the inner LP
  // no longer dualises x_NDO[ j ] (where j is k's outer-variable index):
  // its inner cost for column j becomes the unshifted c. We mirror on
  // the LP side by zeroing the - xLP[ j ] coefficient in every (i, j)
  // constraint, which is the dual of "column j's inner cost is c".

  if( LPTr && ( wchg & 1024 ) && ( dis( rg ) <= p_change ) ) {
   const Index loc_NV = LBF->get_num_active_var();
   if( loc_NV > 1 ) {  // keep at least one dual_pair alive
    const Index k = Index( dis( rg ) * loc_NV );
    auto xNDO = NDOBlock->get_static_variable_v< ColVariable >( 0 );
    const auto * xp = static_cast< ColVariable * >(
                                       LBF->get_active_var( k ) );
    const Index j = Index( xp - xNDO->data() );

    LOG1( "removed dual_pair k=" << k << " (j=" << j << ") from LBF - " );

    // NDOBlock side: naked FunctionModVarsRngd[ k , k+1 )
    LBF->remove_variable( k , eModBlck );

    // LPBlock side: zero the - xLP[ j ] term in every (i, j) constraint
    // of LPTr's "pc" — equivalent to dropping xLP[ j ] from the
    // LinearFunction, but doesn't shift the other coefficients'
    // positions, which is friendlier to subsequent removals.
    auto pc = LPTr->get_static_constraint< FRowConstraint , 2 >( "pc" );
    auto xLP = LPBlock->get_static_variable_v< ColVariable >( "x" );
    for( Index i = 0 ; i < nvar ; ++i ) {
     auto lf = static_cast< p_LF >( (*pc)[ i ][ j ].get_function() );
     const auto idx = lf->is_active( & (*xLP)[ j ] );
     if( idx < lf->get_num_active_var() )
      lf->modify_coefficient( idx , 0.0 );
     }
    }
   }

  // per-LagBFunction dual_pair addition- - - - - - - - - - - - - - - - - - -
  //
  // bit 11: add to THIS LBF one dual_pair ( xNDO[ j ] , \sum_i f[ i ][ j ] )
  // on an x_j that it does not have, and mirror it by putting back the
  // - xLP[ j ] coefficient in every (i, j) constraint on the LPBlock side.
  // Such j exist after a removal by bit 10, or from setup with bit 9. If no
  // other component has x_j it comes in as a new coordinate of the master,
  // otherwise this LBF starts depending on a coordinate the master already
  // has. The C05FunctionModVarsAddd is issued by add_dual_pairs() with its
  // default ModParam, as a user adding a dual_pair to a LagBFunction with
  // the Solver attached does.

  if( LPTr && ( wchg & 2048 ) && ( dis( rg ) <= p_change ) ) {
   auto xNDO = NDOBlock->get_static_variable_v< ColVariable >( 0 );
   Subset cand;
   for( Index j = 0 ; j < nvar ; ++j )
    if( LBF->is_active( & (*xNDO)[ j ] ) >= LBF->get_num_active_var() )
     cand.push_back( j );

   if( ! cand.empty() ) {
    const Index j = cand[ Index( dis( rg ) * cand.size() ) ];

    LOG1( "added dual_pair (j=" << j << ") to LBF - " );

    // NDOBlock side: the dual_pair on the flow variables entering j
    auto f = NDOTr->get_static_variable< ColVariable , 2 >( "f" );
    v_coeff_pair cfj( nvar );
    for( Index i = 0 ; i < nvar ; ++i )
     cfj[ i ] = std::make_pair( & (*f)[ i ][ j ] , 1.0 );

    LagBFunction::v_dual_pair lp;
    lp.emplace_back( & (*xNDO)[ j ] , new LinearFunction( std::move( cfj ) ) );
    LBF->add_dual_pairs( std::move( lp ) );

    // LPBlock side: - xLP[ j ] back in every (i, j) constraint of "pc",
    // modifying the coefficient zeroed by bit 10 or adding the term that
    // bit 9 left out
    auto pc = LPTr->get_static_constraint< FRowConstraint , 2 >( "pc" );
    auto xLP = LPBlock->get_static_variable_v< ColVariable >( "x" );
    for( Index i = 0 ; i < nvar ; ++i ) {
     auto lf = static_cast< p_LF >( (*pc)[ i ][ j ].get_function() );
     const auto idx = lf->is_active( & (*xLP)[ j ] );
     if( idx < lf->get_num_active_var() )
      lf->modify_coefficient( idx , -1.0 );
     else
      lf->add_variable( & (*xLP)[ j ] , -1.0 );
     }
    }
   }

  // modify linear objective- - - - - - - - - - - - - - - - - - - - - - - - -
  // ... if there is any, of course

  if( ( nf < 0 ) && ( wchg & 256 ) && ( dis( rg ) <= p_change ) )
   if( Index tochange = Index( dis( rg ) * std::min( nvar , n_change ) ) ) {
    LOG1( "changed " << tochange << " objective coeff." );

    GenerateA( 1 , tochange );

    auto LPLF = static_cast< p_LF >(
	    ( LPBlock->get_objective< FRealObjective >() )->get_function() );
    auto NDOLF = static_cast< p_LF >(
	   ( NDOBlock->get_objective< FRealObjective >() )->get_function() );

    if( dis( rg ) <= 0.5 ) {  // in 50% of the cases do a ranged change
     Index strt = dis( rg ) * ( nvar - tochange );
     Index stp = strt + tochange;

     if( tochange == 1 ) {
      LPLF->modify_coefficient( strt , A[ 0 ][ 0 ] );
      NDOLF->modify_coefficient( strt , A[ 0 ][ 0 ] );
      }
     else {
      LPLF->modify_coefficients( RealVector( A[ 0 ] ) , Range( strt , stp ) );
      NDOLF->modify_coefficients( std::move( A[ 0 ] ) , Range( strt , stp ) );
      }
      
     LOG1( "(r) - " );
     }
    else {  // in the other 50% of the cases, do a sparse change
     Subset nms = GenerateSubset( nvar , tochange );

     if( tochange == 1 ) {
      LPLF->modify_coefficient( nms.front() , A[ 0 ][ 0 ] );
      NDOLF->modify_coefficient( nms.front() , A[ 0 ][ 0 ] );
      }
     else {
      LPLF->modify_coefficients( RealVector( A[ 0 ] ) , Subset( nms ) ,
				 true );
      NDOLF->modify_coefficients( std::move( A[ 0 ] ) , std::move( nms ) ,
				  true );
      }

     LOG1( "(s) - " );
     }
    }

  // add variables- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  #if DYNAMIC_VARS > 0
  if( ( wchg & 128 ) && ( dis( rg ) <= p_change ) ) {
   Index tochange = Index( dis( rg ) * n_change );
   if( tochange ) {
    LOG1( "added " << tochange << " variables - " );

    throw( std::logic_error( "adding variables not implemented yet" ) );

    GenerateA( m , tochange );

    // add them in the LP, *copying* the data
    std::list< ColVariable > nxLPd( tochange );
    std::vector< Variable * > nxpLP( tochange );
    auto nxit = nxLPd.begin();
    for( Index i = 0 ; i < tochange ; )
     nxpLP[ i++ ] = &(*(nxit++));

    LPBlock->add_dynamic_variables(
	      *(LPBlock->get_dynamic_variable< ColVariable >( 0 )) , nxLPd );

    if( tochange == 1 )
     LPBlock->get_PolyhedralFunction().add_variable( nxpLP[ 0 ] , A[ 0 ] );
    else
     LPBlock->get_PolyhedralFunction().add_variables( std::move( nxpLP ) ,
						      MultiVector( A ) );

    // add them in the NDO, letting the data go
    std::list< ColVariable > nxNDOd( tochange );
    std::vector< Variable * > nxpNDO( tochange );
     auto nxit = nxNDOd.begin();
    for( Index i = 0 ; i < tochange ; )
     nxpNDO[ i++ ] = &(*(nxit++));

    NDOBlock->add_dynamic_variables(
	    *(NDOBlock->get_dynamic_variable< ColVariable >( 0 )) , nxNDOd );

    if( tochange == 1 )
     NDOBlock->get_PolyhedralFunction().add_variable( nxpNDO[ 0 ] , A[ 0 ] );
    else
     NDOBlock->get_PolyhedralFunction().add_variables( std::move( nxpNDO ) ,
						       std::move( A ) );

    // update ndvar
    ndvar += tochange;
    }
   }

  // remove variables - - - - - - - - - - - - - - - - - - - - - - - - - - - -

  if( ( wchg & 256 ) && ( dis( rg ) <= p_change ) ) {
   Index tochange = min( ndvar , Index( dis( rg ) * n_change ) );
   if( tochange ) {
    LOG1( "removed " << tochange << " variables" );

    throw( std::logic_error( "removing variables not implemented yet" ) );

    // in 50% of the cases do a ranged change, in the others a sparse change
    if( dis( rg ) <= 0.5 ) {
     LOG1( "(r) - " );

     Index strt = dis( rg ) * ( ndvar - tochange );
     Index stp = strt + tochange;

     // remove them from the LP
     auto xLPd = NDOBlock->get_dynamic_variable< ColVariable >( 0 );
     if( tochange == 1 )
      LPBlock->get_PolyhedralFunction().remove_variable( strt );
     else
      LPBlock->get_PolyhedralFunction().remove_variables( Range( strt ,
								 stp ) );

     LPBlock->remove_dynamic_variables( *xLPd , Range( strt , stp ) );

     // remove them from the NDO
     auto xNDOd = NDOBlock->get_dynamic_variable< ColVariable >( 0 );
     if( tochange == 1 )
      NDOBlock->get_PolyhedralFunction().remove_variable( strt );
     else
      NDOBlock->get_PolyhedralFunction().remove_variables( Range( strt ,
								  stp ) );

     NDOBlock->remove_dynamic_variables( *xNDOd , Range( strt , stp ) );
     }
    else {
     LOG1( "(s) - " );
     Subset nms = GenerateSubset( ndvar , tochange );

     // remove them from the LP, *copying* names
     auto xLPd = LPBlock->get_dynamic_variable< ColVariable >( 0 );
     if( tochange == 1 ) {
      LPBlock->get_PolyhedralFunction().remove_variable( nms[ 0 ] );
      auto vp = &(*std::next( xLPd->begin() , nms[ 0 ] ));
      LPBlock->remove_dynamic_variable( *xLPd , vp );
      }
     else {
      LPBlock->get_PolyhedralFunction().remove_variables( Subset( nms ) );
      LPBlock->remove_dynamic_variables( *xLPd , Subset( nms ) );
      }

     // remove them from the NDO, finally letting names go
     auto xNDOd = NDOBlock->get_dynamic_variable< ColVariable >( 0 );
     if( tochange == 1 ) {
      NDOBlock->get_PolyhedralFunction().remove_variable( nms[ 0 ] );
      auto vp = &(*std::next( xNDOd->begin() , nms[ 0 ] ));
      NDOBlock->remove_dynamic_variable( *xNDOd , vp );
      }
     else {
      NDOBlock->get_PolyhedralFunction().remove_variables( Subset( nms ) );
      NDOBlock->remove_dynamic_variables( *xNDOd , std::move( nms ) );
      }
     }

    // update ndvar
    ndvar -= tochange;
    }
   }

  #endif  // DYNAMIC_VARS > 0

  // if verbose, print out stuff- - - - - - - - - - - - - - - - - - - - - - -

  #if( LOG_LEVEL >= 3 )
  {
   auto LPslv = ( LPBlock->get_registered_solvers() ).front();
   LPslv->set_par( LPslv->str_par_str2idx( "strOutputFile" ) ,
		   "LPBlock-" + std::to_string( rep ) + ".lp" );
   }

   #if( LOG_LEVEL >= 5 )
    if( bn < Index( nf ) ) {
     cout << endl << "LPBlock-PF: ";
     auto PF = & LPBr->get_PolyhedralFunction();
     printAb( PF->get_A() , PF->get_b() , convex
	      ? PF->get_global_lower_bound()
	      : PF->get_global_upper_bound() );
     }
   #endif
  #endif

  // finally, re-solve the problems- - - - - - - - - - - - - - - - - - - - -
  // ... every SKIP_BEAT + 1 rounds

  if( ! ( ++rep % ( SKIP_BEAT + 1 ) ) )
   AllPassed &= SolveEach();
  #if( LOG_LEVEL >= 1 )
  else
   cout << endl;
  #endif

  }  // end( main loop )- - - - - - - - - - - - - - - - - - - - - - - - - - -
     // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 if( AllPassed )
  cout << GREEN( All test passed!! ) << endl;
 else
  cout << RED( Shit happened!! ) << endl;
 
 // destroy objects and vectors - - - - - - - - - - - - - - - - - - - - - - - 
 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 // unregister (and delete) all Solvers attached to the Blocks
 if( HasEasy )
  for( Index p = nf ; p < Index( nf + nt ) ; ++p ) {
   auto FRO =
        NDOBlock->get_nested_Block( p )->get_objective< FRealObjective >();
   auto LBF = static_cast< LagBFunction * >( FRO->get_function() );
   if( auto inner = LBF->get_inner_block() )
    inner->unregister_Solvers();
   }

 NDOBlock->unregister_Solvers();

 for( Index i = 0 ; i < Index( nf ) ; )
  LPBlock->get_nested_Block( i++ )->unregister_Solvers();

 LPBlock->unregister_Solvers();

 // delete the Blocks
 delete( NDOBlock );
 delete( LPBlock );

 // terminate - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 return( AllPassed ? 0 : 1 );

 }  // end( main )

/*--------------------------------------------------------------------------*/
/*------------------------ End File test.cpp -------------------------------*/
/*--------------------------------------------------------------------------*/