tests_NetCDF.cpp

/** @file
 * Unit tests for the netCDF format of the components of the core.
 *
 * Each component is written to a netCDF group with serialize() and read back
 * with deserialize() (or with the factory, as Block::new_Block() does), and
 * what is read is compared with the original: the sizes, the coefficients,
 * the bounds, the sides, the values. The components are the rows
 * (FRowConstraint), the bounds (the OneVarConstraint family) and the
 * FRealObjective of an AbstractBlock, which travel as the LP file of the
 * model, the PolyhedralFunctionBlock, the BendersBFunction with its
 * sub-Block and the LagBFunction; the empty cases, zero rows or zero
 * variables, are there too. For the :Function that have a State, the State
 * goes the same way, both the one of get_State() and the one written by
 * serialize_State(), and put_State() of what is read back gives a Function
 * with the same global pool. A Block with nested Block also goes through
 * the image of a netCDF file in memory (NcMemFile), before and after its
 * data are changed, and is compared with the one read from a file on disk;
 * copy_group() and write_via_memory() are checked against writing in place,
 * and write_stacked() and read_stacked() against the parts written alone.
 *
 * All of this is done with the core alone: no Solver is needed, since the
 * global pools are filled by compute() for a PolyhedralFunction and by
 * reading a State written here for the LagBFunction and the
 * BendersBFunction, which would need a Solver of their sub-Block to compute
 * anything.
 *
 * \author Donato Meoli \n
 *         Dipartimento di Informatica \n
 *         Universita' di Pisa \n
 *
 * \copyright © by Donato Meoli
 */

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

#include "AbstractBlock.h"
#include "BendersBFunction.h"
#include "ColVariableSolution.h"
#include "FRealObjective.h"
#include "FRowConstraint.h"
#include "LagBFunction.h"
#include "LinearFunction.h"
#include "OneVarConstraint.h"
#include "PolyhedralFunctionBlock.h"
#include "SMSnetCDF.h"

#include <algorithm>
#include <chrono>
#include <cmath>
#include <filesystem>
#include <fstream>
#include <iterator>
#include <iostream>
#include <limits>
#include <map>
#include <memory>
#include <sstream>
#include <stdexcept>
#include <tuple>
#include <vector>

// last, so that the headers above are read as the library was compiled
#include "TestAssert.h"

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

using namespace SMSpp_di_unipi_it;

using Index = Block::Index;
using MultiVector = PolyhedralFunction::MultiVector;
using RealVector = PolyhedralFunction::RealVector;

/*--------------------------------------------------------------------------*/
/*------------------------------ FUNCTIONS ---------------------------------*/
/*--------------------------------------------------------------------------*/

/// what the tests give to a Block that the Block does not own, deleted at the
/// end of main(), after every Block, the last given the first deleted

static std::vector< std::shared_ptr< void > > given;

template< class T >
static T * kept( T * t )
{
 given.emplace_back( t );
 return( t );
 }

/// an FRealObjective of b with Function f, to be given to b

static FRealObjective * objective( Block * b , Function * f )
{
 return( kept( new FRealObjective( b , f ) ) );
 }

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

/// the name of the netCDF file the tests write, removed at the end

static const std::string & nc_file( void )
{
 static const std::string f = ( std::filesystem::temp_directory_path() /
  ( "smspp_NetCDF_test_" + std::to_string(
       std::chrono::steady_clock::now().time_since_epoch().count() ) +
    ".nc4" ) ).string();
 return( f );
 }

/*--------------------------------------------------------------------------*/
/// the netCDF round trip of a Block through the factory

static Block * round_trip( const Block & b )
{
 {
  netCDF::NcFile f( nc_file() , netCDF::NcFile::replace );
  auto g = f.addGroup( "B" );
  b.serialize( g );
  }
 netCDF::NcFile f( nc_file() , netCDF::NcFile::read );
 return( Block::new_Block( f.getGroup( "B" ) ) );
 }

/*--------------------------------------------------------------------------*/
/// writes a State with the given writer and reads it back with the factory
/** \p write gets the group, and has to write in it a State. */

template< class W >
static State * state_round_trip( W write )
{
 {
  netCDF::NcFile f( nc_file() , netCDF::NcFile::replace );
  auto g = f.addGroup( "S" );
  write( g );
  }
 netCDF::NcFile f( nc_file() , netCDF::NcFile::read );
 return( State::new_State( f.getGroup( "S" ) ) );
 }

/*--------------------------------------------------------------------------*/
/// true if a and b are the same number, the infinite ones included

static bool near( double a , double b )
{
 if( a == b )
  return( true );
 return( std::abs( a - b ) <= 1e-12 * std::max( 1.0 , std::abs( a ) ) );
 }

/*--------------------------------------------------------------------------*/
/// the model an AbstractBlock describes, independent of how it is grouped
/** The columns in the order of the groups of Variable, each with its being
 * integer and its bounds, the ones the ColVariable has of its own and those
 * of the OneVarConstraint on it together; the rows as the sorted list of
 * the one-sided or equality rows, a row with two different finite sides
 * being two of them, as the LP format writes it; the sense, the
 * coefficients and the constant of the Objective. */

struct Model {
 using Row = std::tuple< std::vector< std::pair< Index , double > > ,
			 double , double >;

 std::vector< bool > integer;
 std::vector< double > lb;
 std::vector< double > ub;
 std::vector< Row > rows;
 int sense = Objective::eMin;
 std::vector< double > obj;
 double obj_const = 0;
 };

static Model model_of( const Block & b )
{
 Model m;
 std::map< const Variable * , Index > idx;

 b.for_each_variable_group( [ & ]( const BaseGroup & g ) {
   for( Index i = 0 ; i < g.get_num_elements() ; ++i ) {
    auto v = dynamic_cast< const ColVariable * >( g.get_Variable( i ) );
    assert( v );
    idx[ v ] = m.lb.size();
    m.integer.push_back( v->is_integer() );
    m.lb.push_back( v->is_positive() ? 0 : - Inf< double >() );
    m.ub.push_back( v->is_negative() ? 0 : Inf< double >() );
    }
   } );

 auto coefficients = [ & idx ]( const Function * f ) {
  std::vector< std::pair< Index , double > > c;
  auto lf = dynamic_cast< const LinearFunction * >( f );
  assert( lf );
  for( const auto & [ var , coeff ] : lf->get_v_var() )
   if( coeff != 0 )
    c.emplace_back( idx.at( var ) , coeff );
  std::sort( c.begin() , c.end() );
  return( c );
  };

 b.for_each_constraint_group( [ & ]( const BaseGroup & g ) {
   for( Index i = 0 ; i < g.get_num_elements() ; ++i ) {
    auto c = g.get_Constraint( i );
    if( auto o = dynamic_cast< const OneVarConstraint * >( c ) ) {
     auto k = idx.at( o->get_active_var( 0 ) );
     m.lb[ k ] = std::max( m.lb[ k ] , o->get_lhs() );
     m.ub[ k ] = std::min( m.ub[ k ] , o->get_rhs() );
     continue;
     }
    auto r = dynamic_cast< const FRowConstraint * >( c );
    assert( r );
    auto coeffs = coefficients( r->get_function() );
    const double lhs = r->get_lhs();
    const double rhs = r->get_rhs();
    if( lhs == rhs )
     m.rows.emplace_back( coeffs , lhs , rhs );
    else {
     if( lhs > - Inf< double >() )
      m.rows.emplace_back( coeffs , lhs , Inf< double >() );
     if( rhs < Inf< double >() )
      m.rows.emplace_back( coeffs , - Inf< double >() , rhs );
     }
    }
   } );
 std::sort( m.rows.begin() , m.rows.end() );

 m.obj.assign( m.lb.size() , 0 );
 if( auto o = dynamic_cast< const FRealObjective * >( b.get_objective() ) ) {
  m.sense = o->get_sense();
  for( auto & [ k , c ] : coefficients( o->get_function() ) )
   m.obj[ k ] = c;
  m.obj_const = static_cast< const LinearFunction * >( o->get_function() )
   ->get_constant_term();
  }

 return( m );
 }

/*--------------------------------------------------------------------------*/
/// true if the two models are the same

static bool same( const Model & a , const Model & b )
{
 if( ( a.integer != b.integer ) || ( a.lb.size() != b.lb.size() ) ||
     ( a.rows.size() != b.rows.size() ) || ( a.sense != b.sense ) ||
     ( ! near( a.obj_const , b.obj_const ) ) )
  return( false );
 for( std::size_t j = 0 ; j < a.lb.size() ; ++j )
  if( ( ! near( a.lb[ j ] , b.lb[ j ] ) ) ||
      ( ! near( a.ub[ j ] , b.ub[ j ] ) ) ||
      ( ! near( a.obj[ j ] , b.obj[ j ] ) ) )
   return( false );
 for( std::size_t i = 0 ; i < a.rows.size() ; ++i ) {
  const auto & [ ca , la , ra ] = a.rows[ i ];
  const auto & [ cb , lb , rb ] = b.rows[ i ];
  if( ( ca.size() != cb.size() ) || ( ! near( la , lb ) ) ||
      ( ! near( ra , rb ) ) )
   return( false );
  for( std::size_t k = 0 ; k < ca.size() ; ++k )
   if( ( ca[ k ].first != cb[ k ].first ) ||
       ( ! near( ca[ k ].second , cb[ k ].second ) ) )
    return( false );
  }
 return( true );
 }

/*--------------------------------------------------------------------------*/
/// a LinearFunction with the given coefficients on the given ColVariable

static LinearFunction * linear( std::vector< ColVariable > & x ,
				const std::vector< double > & c ,
				double constant = 0 )
{
 LinearFunction::v_coeff_pair p;
 for( std::size_t j = 0 ; j < c.size() ; ++j )
  if( c[ j ] != 0 )
   p.push_back( { & x[ j ] , c[ j ] } );
 return( new LinearFunction( std::move( p ) , constant ) );
 }

/*--------------------------------------------------------------------------*/
/*------------------------------- THE TESTS --------------------------------*/
/*--------------------------------------------------------------------------*/
/* The rows, the bounds and the Objective of an AbstractBlock: an equality
 * row, a one-sided one and a two-sided one, a BoxConstraint, an
 * LBConstraint and an UBConstraint, an integer column, a free one, and an
 * Objective to be maximised. */

static void test_abstract_block( void )
{
 AbstractBlock block;

 auto x = kept( new std::vector< ColVariable >( 4 ) );
 ( *x )[ 1 ].set_type( ColVariable::kNatural );
 block.add_static_variable( *x , "x" );

 auto rows = kept( new std::vector< FRowConstraint >( 3 ) );
 ( *rows )[ 0 ].set_function( linear( *x , { 1 , 2 , 0 , 0 } ) );
 ( *rows )[ 0 ].set_both( 3 );
 ( *rows )[ 1 ].set_function( linear( *x , { 0 , 0 , -3 , 0.5 } ) );
 ( *rows )[ 1 ].set_lhs( - Inf< double >() );
 ( *rows )[ 1 ].set_rhs( 7 );
 ( *rows )[ 2 ].set_function( linear( *x , { 1 , 0 , -1 , 0 } ) );
 ( *rows )[ 2 ].set_lhs( 1 );
 ( *rows )[ 2 ].set_rhs( 4 );
 block.add_static_constraint( *rows , "r" );

 auto box = kept( new BoxConstraint( & block , & ( *x )[ 0 ] , -1 , 5 ) );
 block.add_static_constraint( *box , "box" );
 auto lb = kept( new LBConstraint( & block , & ( *x )[ 2 ] , -2 ) );
 block.add_static_constraint( *lb , "lb" );
 auto ub = kept( new UBConstraint( & block , & ( *x )[ 3 ] , 8 ) );
 block.add_static_constraint( *ub , "ub" );
 auto lb3 = kept( new LBConstraint( & block , & ( *x )[ 3 ] , -4 ) );
 block.add_static_constraint( *lb3 , "lb3" );

 auto obj = objective( & block , linear( *x , { 5 , -1 , 2 , 1 } ) );
 obj->set_sense( Objective::eMax , eNoMod );
 block.set_objective( obj , eNoMod );

 const auto original = model_of( block );
 assert( original.lb.size() == 4 );
 assert( original.rows.size() == 4 );

 auto read = round_trip( block );
 assert( dynamic_cast< AbstractBlock * >( read ) );
 const auto copy = model_of( *read );
 assert( same( copy , original ) );
 assert( copy.integer[ 1 ] && ( copy.lb[ 1 ] == 0 ) );
 assert( ( copy.lb[ 0 ] == -1 ) && ( copy.ub[ 0 ] == 5 ) );
 assert( ( copy.lb[ 2 ] == -2 ) && ( copy.ub[ 2 ] == Inf< double >() ) );
 assert( ( copy.lb[ 3 ] == -4 ) && ( copy.ub[ 3 ] == 8 ) );
 assert( copy.sense == Objective::eMax );

 // and a second trip gives the same again
 auto again = round_trip( *read );
 assert( again && same( model_of( *again ) , original ) );
 delete again;
 delete read;

 // a column with no lower bound and a finite upper bound u, which the LP
 // format has to be told is -inf <= x <= u, a lower bound that is not
 // written being 0 to it
 {
  AbstractBlock only_ub;
  auto y = kept( new std::vector< ColVariable >( 1 ) );
  only_ub.add_static_variable( *y , "y" );
  auto u = kept( new UBConstraint( & only_ub , & ( *y )[ 0 ] , 8 ) );
  only_ub.add_static_constraint( *u , "u" );
  only_ub.set_objective( objective( & only_ub , linear( *y , { 1 } ) ) ,
			 eNoMod );
  auto r = round_trip( only_ub );
  assert( r && same( model_of( *r ) , model_of( only_ub ) ) );
  delete r;
  }

 std::cout << "rows, bounds and Objective of an AbstractBlock: OK"
	   << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* The empty cases of an AbstractBlock: nothing at all, columns and no row,
 * an Objective and no row. */

static void test_empty_abstract_block( void )
{
 {
  AbstractBlock block;
  auto read = round_trip( block );
  assert( dynamic_cast< AbstractBlock * >( read ) );
  const auto m = model_of( *read );
  assert( m.lb.empty() && m.rows.empty() && ( ! read->get_objective() ) );
  assert( ! read->get_number_nested_Blocks() );
  delete read;
  }
 {
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );
  auto obj = objective( & block , linear( *x , { 1 , -1 } ) );
  block.set_objective( obj , eNoMod );

  auto read = round_trip( block );
  assert( read && same( model_of( *read ) , model_of( block ) ) );
  delete read;
  }
 {
  // an empty group of rows
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );
  auto rows = kept( new std::vector< FRowConstraint >() );
  block.add_static_constraint( *rows , "none" );
  block.set_objective( objective( & block , linear( *x , { 2 , 3 } ) ) ,
		       eNoMod );

  auto read = round_trip( block );
  assert( read && same( model_of( *read ) , model_of( block ) ) );
  delete read;
  }

 // a column that is in no row and has a zero coefficient in the Objective
 {
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );
  block.set_objective( objective( & block , linear( *x , { 0 , 3 } ) ) ,
		       eNoMod );

  auto read = round_trip( block );
  assert( read && same( model_of( *read ) , model_of( block ) ) );
  delete read;
  }

 // columns and no Objective
 {
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );

  auto read = round_trip( block );
  assert( read && same( model_of( *read ) , model_of( block ) ) );
  delete read;
  }

 // an Objective whose coefficients are all zero, the columns being in a row
 // whose coefficients are all zero too: both are written as the constant 0
 {
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );
  auto rows = kept( new std::vector< FRowConstraint >( 2 ) );
  ( *rows )[ 0 ].set_function( linear( *x , { 1 , 1 } ) );
  ( *rows )[ 0 ].set_lhs( - Inf< double >() );
  ( *rows )[ 0 ].set_rhs( 4 );
  ( *rows )[ 1 ].set_function( linear( *x , { 0 , 0 } ) );
  ( *rows )[ 1 ].set_lhs( -1 );
  ( *rows )[ 1 ].set_rhs( Inf< double >() );
  block.add_static_constraint( *rows , "r" );
  block.set_objective( objective( & block , linear( *x , { 0 , 0 } ) ) ,
		       eNoMod );

  auto read = round_trip( block );
  assert( read && same( model_of( *read ) , model_of( block ) ) );
  delete read;
  }

 // an LP file that ends before its End section, or has a word out of its
 // place, is refused rather than read forever
 for( const char * lp : { "Minimize\n obj: 0\n" ,
			  "Minimize\n obj: x\nSubject To\n c: x >= 1\n" ,
			  "Minimize\n obj: x\nSubject To\n c: x >= 1\n"
			  "Bounds\n x <= 3\n" ,
			  "Minimize\n obj: x\nSubject To\n c: x >= 1\n"
			  "Bounds\n x <= <= 3\nEnd\n" ,
			  "Minimize\n obj: x\nSubject To\n x >= 1\nEnd\n" ,
			  "Maximise\n obj: x\nSubject To\nEnd\n" } ) {
  AbstractBlock block;
  std::istringstream in( lp );
  bool refused = false;
  try { block.load( in , 'L' ); }
  catch( const std::invalid_argument & ) { refused = true; }
  assert( refused );
  }

 // what read_lp() takes that write_lp() does not write: a constant in the
 // Objective and in a row, a column named in the Bounds section alone, the
 // bounds with the sense turned, the default bounds of a column
 {
  AbstractBlock block;
  std::istringstream in( "Maximize\n obj: 2 x + 3 + y\nSubject To\n"
			 " c: x + 1 <= 5\n d: 0 >= -2\nBounds\n"
			 " 4 >= x >= -1\n z <= 6\nEnd\n" );
  block.load( in , 'L' );
  const auto m = model_of( block );
  assert( m.lb.size() == 3 );
  assert( ( m.obj[ 0 ] == 2 ) && ( m.obj[ 1 ] == 1 ) && ( m.obj[ 2 ] == 0 ) );
  assert( ( m.obj_const == 3 ) && ( m.sense == Objective::eMax ) );
  assert( ( m.lb[ 0 ] == -1 ) && ( m.ub[ 0 ] == 4 ) );
  assert( ( m.lb[ 1 ] == 0 ) && ( m.ub[ 1 ] == Inf< double >() ) );
  assert( ( m.lb[ 2 ] == 0 ) && ( m.ub[ 2 ] == 6 ) );
  assert( m.rows.size() == 2 );
  // sorted: the row with no coefficient comes first
  assert( std::get< 0 >( m.rows[ 0 ] ).empty() &&
	  ( std::get< 1 >( m.rows[ 0 ] ) == -2 ) );
  assert( ( std::get< 0 >( m.rows[ 1 ] ).size() == 1 ) &&
	  ( std::get< 2 >( m.rows[ 1 ] ) == 4 ) );
  }

 std::cout << "empty AbstractBlock: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* PolyhedralFunctionBlock: the matrix, the constants, the verse and the
 * bound of the PolyhedralFunction, and the empty cases. */

static void check_pf( PolyhedralFunction & a , PolyhedralFunction & b )
{
 assert( a.get_A() == b.get_A() );
 assert( a.get_b() == b.get_b() );
 assert( a.get_nrows() == b.get_nrows() );
 assert( a.is_convex() == b.is_convex() );
 assert( a.is_bound_set() == b.is_bound_set() );
 if( a.is_bound_set() )
  assert( a.get_global_bound() == b.get_global_bound() );
 }

static void test_polyhedral_function_block( void )
{
 std::vector< ColVariable > x( 3 );

 // a convex max with a bound, a concave min with none
 for( bool convex : { true , false } ) {
  auto pfb = new PolyhedralFunctionBlock();
  auto & pf = pfb->get_PolyhedralFunction();
  pf.set_variables( { & x[ 0 ] , & x[ 1 ] , & x[ 2 ] } );
  pf.set_PolyhedralFunction( { { 1 , 2 , 0 } , { -3 , 0.5 , 4 } } ,
			     { 0 , 1.5 } ,
			     convex ? -10 : Inf< double >() , convex ,
			     eNoMod );

  auto read = dynamic_cast< PolyhedralFunctionBlock * >( round_trip( *pfb ) );
  assert( read );
  auto & rpf = read->get_PolyhedralFunction();
  check_pf( rpf , pf );
  assert( rpf.get_A()[ 0 ].size() == 3 );

  // what has been read, with no Variable given to it yet, is written as it
  // was read: the columns of A are those of the matrix
  auto again = dynamic_cast< PolyhedralFunctionBlock * >(
						      round_trip( *read ) );
  assert( again );
  check_pf( again->get_PolyhedralFunction() , pf );
  assert( again->get_PolyhedralFunction().get_A()[ 0 ].size() == 3 );
  delete again;

  // the Variable are not part of the format: those given afterwards must
  // be as many as the columns of A
  rpf.set_variables( { & x[ 0 ] , & x[ 1 ] , & x[ 2 ] } );
  assert( rpf.get_num_active_var() == 3 );
  delete read;
  delete pfb;
  }

 // zero rows
 {
  auto pfb = new PolyhedralFunctionBlock();
  auto & pf = pfb->get_PolyhedralFunction();
  pf.set_variables( { & x[ 0 ] , & x[ 1 ] } );
  pf.set_PolyhedralFunction( {} , {} , -1 , true , eNoMod );

  auto read = dynamic_cast< PolyhedralFunctionBlock * >( round_trip( *pfb ) );
  assert( read );
  check_pf( read->get_PolyhedralFunction() , pf );
  assert( ! read->get_PolyhedralFunction().get_nrows() );
  delete read;
  delete pfb;
  }

 // zero variables: a constant function, the max of the constants
 {
  auto pfb = new PolyhedralFunctionBlock();
  auto & pf = pfb->get_PolyhedralFunction();
  pf.set_PolyhedralFunction( { {} , {} } , { 2 , 3 } , - Inf< double >() ,
			     true , eNoMod );

  auto read = dynamic_cast< PolyhedralFunctionBlock * >( round_trip( *pfb ) );
  assert( read );
  check_pf( read->get_PolyhedralFunction() , pf );
  assert( read->get_PolyhedralFunction().get_nrows() == 2 );
  delete read;
  delete pfb;
  }

 // which rows are vertical linearizations goes in the file too
 {
  auto pfb = new PolyhedralFunctionBlock();
  auto & pf = pfb->get_PolyhedralFunction();
  pf.set_variables( { & x[ 0 ] , & x[ 1 ] } );
  pf.set_PolyhedralFunction( { { 1 , 0 } , { 0 , 1 } } , { 0 , -1 } ,
			     - Inf< double >() , true , eNoMod ,
			     { false , true } );
  auto read = dynamic_cast< PolyhedralFunctionBlock * >( round_trip( *pfb ) );
  assert( read );
  assert( read->get_PolyhedralFunction().get_is_vert() == pf.get_is_vert() );
  delete read;
  delete pfb;
  }

 std::cout << "PolyhedralFunctionBlock: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* The State of a PolyhedralFunction: a global pool with two linearizations
 * and a combination of them, and the important linearization; the State of
 * get_State() and the one of serialize_State() both give it back. */

/// a PolyhedralFunction on x with the data used for the State

static void fill_pf( PolyhedralFunction & pf , std::vector< ColVariable > & x )
{
 pf.set_variables( { & x[ 0 ] , & x[ 1 ] } );
 pf.set_PolyhedralFunction( { { 1 , 2 } , { -1 , 1 } , { 0.5 , -3 } } ,
			    { 0 , 1 , 2 } , - Inf< double >() , true ,
			    eNoMod );
 pf.set_par( PolyhedralFunction::intGPMaxSz , 4 );
 }

/// true if the two PolyhedralFunction have the same global pool

static bool same_pool( PolyhedralFunction & a , PolyhedralFunction & b )
{
 if( a.get_int_par( PolyhedralFunction::intGPMaxSz ) !=
     b.get_int_par( PolyhedralFunction::intGPMaxSz ) )
  return( false );
 const Index n = a.get_int_par( PolyhedralFunction::intGPMaxSz );
 for( Index i = 0 ; i < n ; ++i ) {
  if( a.is_linearization_there( i ) != b.is_linearization_there( i ) )
   return( false );
  if( ! a.is_linearization_there( i ) )
   continue;
  if( a.is_linearization_vertical( i ) != b.is_linearization_vertical( i ) )
   return( false );
  if( ! near( a.get_linearization_constant( i ) ,
	      b.get_linearization_constant( i ) ) )
   return( false );
  std::vector< double > ga( 2 ) , gb( 2 );
  a.get_linearization_coefficients( ga.data() ,
				    PolyhedralFunction::Range( 0 , 2 ) , i );
  b.get_linearization_coefficients( gb.data() ,
				    PolyhedralFunction::Range( 0 , 2 ) , i );
  if( ga != gb )
   return( false );
  }
 return( a.get_important_linearization_coefficients() ==
	 b.get_important_linearization_coefficients() );
 }

static void test_polyhedral_function_state( void )
{
 // a new PolyhedralFunction has the default tolerance on the multipliers
 {
  PolyhedralFunction fresh;
  assert( fresh.get_dbl_par( PolyhedralFunction::dblAAccMlt ) ==
	  fresh.get_dflt_dbl_par( PolyhedralFunction::dblAAccMlt ) );
  }

 std::vector< ColVariable > x( 2 );
 PolyhedralFunction pf;
 fill_pf( pf , x );

 // an empty pool first
 {
  std::unique_ptr< State > s( pf.get_State() );
  std::unique_ptr< State > r( state_round_trip( [ & s ]( auto & g ) {
			      s->serialize( g ); } ) );
  assert( dynamic_cast< PolyhedralFunctionState * >( r.get() ) );
  PolyhedralFunction other;
  fill_pf( other , x );
  other.put_State( *r );
  assert( same_pool( other , pf ) );
  }

 // two linearizations, in x = ( 1 , 1 ) and in x = ( -1 , 0 ), and their
 // combination with weights 1/2
 x[ 0 ].set_value( 1 );
 x[ 1 ].set_value( 1 );
 pf.compute();
 pf.store_linearization( 0 );
 x[ 0 ].set_value( -1 );
 x[ 1 ].set_value( 0 );
 pf.compute();
 pf.store_linearization( 1 );
 pf.store_combination_of_linearizations( { { 0 , 0.5 } , { 1 , 0.5 } } , 2 );
 pf.set_important_linearization( { { 0 , 0.25 } , { 2 , 0.75 } } );
 assert( pf.is_linearization_there( 0 ) && pf.is_linearization_there( 2 ) );
 assert( ! pf.is_linearization_there( 3 ) );

 // the State of get_State()
 {
  std::unique_ptr< State > s( pf.get_State() );
  std::unique_ptr< State > r( state_round_trip( [ & s ]( auto & g ) {
			      s->serialize( g ); } ) );
  assert( r );
  PolyhedralFunction other;
  fill_pf( other , x );
  other.put_State( *r );
  assert( same_pool( other , pf ) );

  // and the one moved in
  PolyhedralFunction moved;
  fill_pf( moved , x );
  moved.put_State( std::move( *r ) );
  assert( same_pool( moved , pf ) );
  }

 // the State written by serialize_State()
 {
  std::unique_ptr< State > r( state_round_trip( [ & pf ]( auto & g ) {
			      pf.serialize_State( g ); } ) );
  assert( dynamic_cast< PolyhedralFunctionState * >( r.get() ) );
  PolyhedralFunction other;
  fill_pf( other , x );
  other.put_State( *r );
  assert( same_pool( other , pf ) );
  }

 std::cout << "State of a PolyhedralFunction: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* BendersBFunction: the mapping, the sides, the paths to the RowConstraint
 * of the sub-Block and the sub-Block itself; zero rows and zero variables
 * too. */

/// the sub-Block: m rows on two ColVariable, and an Objective

static AbstractBlock * sub_Block( Index m ,
				  std::vector< RowConstraint * > & rows )
{
 auto b = new AbstractBlock();
 auto y = kept( new std::vector< ColVariable >( 2 ) );
 ( *y )[ 0 ].is_positive( true , eNoMod );
 b->add_static_variable( *y , "y" );

 auto r = kept( new std::vector< FRowConstraint >( m ) );
 for( Index i = 0 ; i < m ; ++i ) {
  ( *r )[ i ].set_function( linear( *y , { 1.0 + i , -1 } ) );
  ( *r )[ i ].set_lhs( - double( i ) );
  ( *r )[ i ].set_rhs( double( i ) + 1 );
  }
 b->add_static_constraint( *r , "r" );
 rows.clear();
 for( auto & row : *r )
  rows.push_back( & row );

 b->set_objective( objective( b , linear( *y , { 1 , 2 } ) ) , eNoMod );
 return( b );
 }

/// the BendersBFunction and the data it was built with

struct Benders {
 std::vector< ColVariable > x;
 BendersBFunction * f;

 Benders( MultiVector A , RealVector b ,
	  BendersBFunction::ConstraintSideVector sides , Index nx = 2 )
  : x( nx ) {
  std::vector< RowConstraint * > rows;
  auto inner = sub_Block( A.size() , rows );
  BendersBFunction::VarVector vars;
  for( auto & v : x )
   vars.push_back( & v );
  f = new BendersBFunction( inner , std::move( vars ) , std::move( A ) ,
			    std::move( b ) , std::move( rows ) ,
			    std::move( sides ) );
  }

 ~Benders() { delete f; }

 /// a BendersBFunction on the same x, read out of the netCDF format of f
 BendersBFunction * read( void ) {
  {
   netCDF::NcFile file( nc_file() , netCDF::NcFile::replace );
   auto g = file.addGroup( "B" );
   f->serialize( g );
   }
  auto r = new BendersBFunction();
  BendersBFunction::VarVector vars;
  for( auto & v : x )
   vars.push_back( & v );
  r->set_variables( std::move( vars ) );
  netCDF::NcFile file( nc_file() , netCDF::NcFile::read );
  r->deserialize( file.getGroup( "B" ) );
  return( r );
  }
 };

static void check_benders( BendersBFunction & r , BendersBFunction & f )
{
 assert( r.get_A() == f.get_A() );
 assert( r.get_b() == f.get_b() );
 assert( r.get_sides() == f.get_sides() );
 assert( r.get_constraints().size() == f.get_constraints().size() );
 assert( r.get_num_active_var() == f.get_num_active_var() );
 for( Index j = 0 ; j < f.get_num_active_var() ; ++j )
  assert( r.get_active_var( j ) == f.get_active_var( j ) );
 assert( r.get_inner_block() && f.get_inner_block() );
 assert( same( model_of( *r.get_inner_block() ) ,
	       model_of( *f.get_inner_block() ) ) );
 }

static void test_benders_function( void )
{
 // A dense, both kinds of sides
 {
  Benders b( { { 1 , 2 } , { 3 , -4 } } , { 5 , 6 } ,
	     { BendersBFunction::eLHS , BendersBFunction::eBoth } );
  auto r = b.read();
  check_benders( *r , *b.f );
  delete r;
  }

 // zero rows
 {
  Benders b( {} , {} , {} );
  auto r = b.read();
  check_benders( *r , *b.f );
  assert( r->get_A().empty() );
  delete r;
  }

 // zero variables and zero rows
 {
  Benders b( {} , {} , {} , 0 );
  auto r = b.read();
  check_benders( *r , *b.f );
  delete r;
  }

 // A sparse enough to be written in the sparse format, with an empty row
 {
  Benders b( { { 1 , 0 , 0 , 0 } , { 0 , 0 , 0 , 0 } , { 0 , 0 , 2 , 0 } ,
	       { 0 , 0 , 0 , 3 } } , { 1 , 2 , 3 , 4 } ,
	     { BendersBFunction::eRHS , BendersBFunction::eRHS ,
	       BendersBFunction::eRHS , BendersBFunction::eRHS } , 4 );
  auto r = b.read();
  check_benders( *r , *b.f );
  delete r;
  }

 // the factory, which builds a BendersBFunction with no active Variable:
 // their number is then the one of the netCDF group, and set_variables()
 // afterwards has to give as many
 {
  Benders b( { { 1 , 2 } } , { 5 } , { BendersBFunction::eBoth } );
  auto r = dynamic_cast< BendersBFunction * >( round_trip( *b.f ) );
  assert( r && ( r->get_A() == b.f->get_A() ) );
  assert( r->get_num_active_var() == 0 );
  bool refused = false;
  try { r->set_variables( { & b.x[ 0 ] } ); }
  catch( const std::logic_error & ) { refused = true; }
  assert( refused );
  r->set_variables( { & b.x[ 0 ] , & b.x[ 1 ] } );
  check_benders( *r , *b.f );
  delete r;
  }

 std::cout << "BendersBFunction: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* The State of a BendersBFunction: a global pool of two places, the first
 * holding a linearization and the second not, and the important
 * linearization. Filling the pool asks a Solver of the sub-Block, hence the
 * State is written here in the format BendersBFunctionState reads. */

static void write_benders_state( netCDF::NcGroup & g )
{
 g.putAtt( "type" , "BendersBFunctionState" );
 auto d = g.addDim( "BendersBFunction_MaxGlob" , 2 );
 std::vector< signed char > type = { 1 , 0 };
 g.addVar( "BendersBFunction_Type" , netCDF::NcByte() , d ).putVar(
							      type.data() );
 std::vector< double > constants = {
  1.5 , std::numeric_limits< double >::quiet_NaN() };
 g.addVar( "BendersBFunction_Constants" , netCDF::NcDouble() , d ).putVar(
							 constants.data() );
 auto c = g.addDim( "BendersBFunction_ImpCoeffNum" , 1 );
 g.addVar( "BendersBFunction_ImpCoeffInd" , netCDF::NcInt() , c ).putVar(
								  { 0 } , 0 );
 g.addVar( "BendersBFunction_ImpCoeffVal" , netCDF::NcDouble() , c ).putVar(
								{ 0 } , 1.0 );
 }

static void check_benders_pool( BendersBFunction & f )
{
 assert( f.is_linearization_there( 0 ) );
 assert( ! f.is_linearization_there( 1 ) );
 assert( ! f.is_linearization_vertical( 0 ) );
 assert( f.get_linearization_constant( 0 ) == 1.5 );
 assert( ( f.get_important_linearization_coefficients() ==
	   C05Function::LinearCombination( { { 0 , 1.0 } } ) ) );
 }

static void test_benders_state( void )
{
 Benders a( { { 1 , 2 } } , { 5 } , { BendersBFunction::eBoth } );
 std::unique_ptr< State > s( state_round_trip( write_benders_state ) );
 assert( dynamic_cast< BendersBFunctionState * >( s.get() ) );
 a.f->put_State( *s );
 check_benders_pool( *a.f );

 // the State of get_State()
 {
  std::unique_ptr< State > g( a.f->get_State() );
  std::unique_ptr< State > r( state_round_trip( [ & g ]( auto & grp ) {
			      g->serialize( grp ); } ) );
  Benders b( { { 1 , 2 } } , { 5 } , { BendersBFunction::eBoth } );
  b.f->put_State( *r );
  check_benders_pool( *b.f );
  }

 // the State written by serialize_State()
 {
  std::unique_ptr< State > r( state_round_trip( [ & a ]( auto & grp ) {
			      a.f->serialize_State( grp ); } ) );
  Benders b( { { 1 , 2 } } , { 5 } , { BendersBFunction::eBoth } );
  b.f->put_State( *r );
  check_benders_pool( *b.f );
  }

 // an empty pool
 {
  Benders e( { { 1 , 2 } } , { 5 } , { BendersBFunction::eBoth } );
  std::unique_ptr< State > g( e.f->get_State() );
  std::unique_ptr< State > r( state_round_trip( [ & g ]( auto & grp ) {
			      g->serialize( grp ); } ) );
  a.f->put_State( *r );
  assert( a.f->get_important_linearization_coefficients().empty() );
  // the places of the pool past those of the State, which it keeps, hold
  // no linearization any more
  assert( a.f->get_int_par( C05Function::intGPMaxSz ) == 2 );
  assert( ! a.f->is_linearization_there( 0 ) );
  assert( ! a.f->is_linearization_there( 1 ) );
  }

 // a State moved in a BendersBFunction whose pool is smaller than the one
 // of the State, and one moved in a BendersBFunction whose pool is larger
 {
  std::unique_ptr< State > r( state_round_trip( write_benders_state ) );
  Benders b( { { 1 , 2 } } , { 5 } , { BendersBFunction::eBoth } );
  b.f->put_State( std::move( *r ) );
  check_benders_pool( *b.f );

  Benders e( { { 1 , 2 } } , { 5 } , { BendersBFunction::eBoth } );
  std::unique_ptr< State > g( e.f->get_State() );
  b.f->put_State( std::move( *g ) );
  assert( b.f->get_int_par( C05Function::intGPMaxSz ) == 2 );
  assert( ! b.f->is_linearization_there( 0 ) );
  }

 std::cout << "State of a BendersBFunction: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* LagBFunction: the inner Block and the Lagrangian term, and its State.
 * Filling the global pool asks a Solver of the inner Block, hence the State
 * is written here in the format LagBFunctionState reads, with a Solution of
 * the inner Block in the first place of the pool. */

/// a LagBFunction and its inner Block
/** With terms == 1, one Lagrangian multiplier y and g( x ) = 2 x_0 + x_1;
 * with terms == 2, y with g_0( x ) = -4 x_1 + 1.5 and y2 with g_1( x ) = -2,
 * a function with no term; with terms == 0, none. */

struct Lagrangian {
 ColVariable y;
 ColVariable y2;
 AbstractBlock * inner;
 std::vector< ColVariable > * x;
 LagBFunction * f;

 Lagrangian( int terms = 1 ) {
  inner = new AbstractBlock();
  x = kept( new std::vector< ColVariable >( 2 ) );
  inner->add_static_variable( *x , "x" );
  inner->set_objective( objective( inner , linear( *x , { 1 , -1 } ) ) ,
			eNoMod );
  f = new LagBFunction( inner );
  LagBFunction::v_dual_pair dp;
  if( terms == 1 )
   dp.emplace_back( & y , linear( *x , { 2 , 1 } ) );
  if( terms == 2 ) {
   dp.emplace_back( & y , linear( *x , { 0 , -4 } , 1.5 ) );
   dp.emplace_back( & y2 , linear( *x , { 0 , 0 } , -2 ) );
   }
  f->set_dual_pairs( std::move( dp ) );
  }

 ~Lagrangian() { delete f; }
 };

static void write_lagrangian_state( netCDF::NcGroup & g , Block * inner ,
				    bool with_value = true )
{
 g.putAtt( "type" , "LagBFunctionState" );
 auto d = g.addDim( "LagBFunction_MaxGlob" , 2 );
 std::vector< signed char > type = { 1 , 0 };
 g.addVar( "LagBFunction_Type" , netCDF::NcByte() , d ).putVar( type.data() );
 if( with_value ) {
  std::vector< double > value = { 2.5 , 0 };
  g.addVar( "LagBFunction_Value" , netCDF::NcDouble() , d ).putVar(
							       value.data() );
  std::vector< signed char > convexified = { 1 , 0 };
  g.addVar( "LagBFunction_Convexified" , netCDF::NcByte() , d ).putVar(
							 convexified.data() );
  }
 ColVariableSolution sol;
 sol.read( inner );
 auto sg = g.addGroup( "LagBFunction_Sol_0" );
 sol.serialize( sg );
 auto c = g.addDim( "LagBFunction_ImpCoeffNum" , 1 );
 g.addVar( "LagBFunction_ImpCoeffInd" , netCDF::NcInt() , c ).putVar(
								  { 0 } , 0 );
 g.addVar( "LagBFunction_ImpCoeffVal" , netCDF::NcDouble() , c ).putVar(
								{ 0 } , 1.0 );
 }

static void check_lagrangian_pool( LagBFunction & f )
{
 assert( f.is_linearization_there( 0 ) );
 assert( ! f.is_linearization_there( 1 ) );
 assert( f.get_linearization_constant( 0 ) == 2.5 );
 assert( ( f.get_important_linearization_coefficients() ==
	   C05Function::LinearCombination( { { 0 , 1.0 } } ) ) );
 }

static void test_lagrangian_function( void )
{
 Lagrangian a;
 ( *a.x )[ 0 ].set_value( 1 );
 ( *a.x )[ 1 ].set_value( 2 );

 std::unique_ptr< State > s( state_round_trip( [ & a ]( auto & g ) {
			     write_lagrangian_state( g , a.inner ); } ) );
 assert( dynamic_cast< LagBFunctionState * >( s.get() ) );
 a.f->put_State( *s );
 check_lagrangian_pool( *a.f );

 // the State of get_State()
 {
  std::unique_ptr< State > g( a.f->get_State() );
  std::unique_ptr< State > r( state_round_trip( [ & g ]( auto & grp ) {
			      g->serialize( grp ); } ) );
  Lagrangian b;
  b.f->put_State( *r );
  check_lagrangian_pool( *b.f );
  Lagrangian c;
  c.f->put_State( std::move( *r ) );
  check_lagrangian_pool( *c.f );
  }

 // the State written by serialize_State(): the pool is the same
 {
  std::unique_ptr< State > r( state_round_trip( [ & a ]( auto & grp ) {
			      a.f->serialize_State( grp ); } ) );
  Lagrangian b;
  b.f->put_State( *r );
  check_lagrangian_pool( *b.f );
  }

 // a State with no LagBFunction_Value and LagBFunction_Convexified, which
 // are optional: the constants are 0
 {
  std::unique_ptr< State > r( state_round_trip( [ & a ]( auto & g ) {
			      write_lagrangian_state( g , a.inner , false );
			      } ) );
  Lagrangian b;
  b.f->put_State( *r );
  assert( b.f->is_linearization_there( 0 ) );
  assert( ! b.f->is_linearization_there( 1 ) );
  assert( b.f->get_linearization_constant( 0 ) == 0 );
  }

 // an empty pool
 {
  Lagrangian e;
  std::unique_ptr< State > g( e.f->get_State() );
  std::unique_ptr< State > r( state_round_trip( [ & g ]( auto & grp ) {
			      g->serialize( grp ); } ) );
  Lagrangian b;
  b.f->put_State( *r );
  assert( b.f->get_important_linearization_coefficients().empty() );
  }

 // a LagBFunction whose pool is larger than the one of the State, and has
 // been made so first: the constants are those of the State as well
 {
  Lagrangian b;
  b.f->set_par( LagBFunction::intGPMaxSz , 3 );
  std::unique_ptr< State > r( state_round_trip( [ & b ]( auto & g ) {
			      write_lagrangian_state( g , b.inner ); } ) );
  b.f->put_State( *r );
  check_lagrangian_pool( *b.f );
  assert( ! b.f->is_linearization_there( 2 ) );
  }

 // the netCDF format: the inner Block and the functions g_i( x ) of the
 // Lagrangian term, whose y are not in the format and are given afterwards
 {
  auto r = dynamic_cast< LagBFunction * >( round_trip( *a.f ) );
  assert( r && ( r->get_num_active_var() == 0 ) );
  assert( same( model_of( *r->get_inner_block() ) , model_of( *a.inner ) ) );

  ColVariable z;
  bool refused = false;
  try { r->set_variables( { & z , & a.y } ); }
  catch( const std::invalid_argument & ) { refused = true; }
  assert( refused );

  r->set_variables( { & z } );
  assert( ( r->get_num_active_var() == 1 ) &&
	  ( r->get_active_var( 0 ) == & z ) );
  auto g = static_cast< const LinearFunction * >(
					       r->get_Lagrangian_term( 0 ) );
  assert( g->get_num_active_var() == 2 );
  for( const auto & [ x , c ] : g->get_v_var() )
   assert( x->get_Block() == r->get_inner_block() );
  assert( ( g->get_v_var()[ 0 ].second == 2 ) &&
	  ( g->get_v_var()[ 1 ].second == 1 ) );
  assert( g->get_constant_term() == 0 );

  // and a second trip gives the same again
  auto again = dynamic_cast< LagBFunction * >( round_trip( *r ) );
  assert( again && same( model_of( *again->get_inner_block() ) ,
			 model_of( *a.inner ) ) );
  again->set_variables( { & z } );
  assert( again->get_num_active_var() == 1 );
  delete again;
  delete r;
  }

 // two Lagrangian terms, one with a constant and one with no term at all,
 // and none at all
 {
  Lagrangian t( 2 );
  auto r = dynamic_cast< LagBFunction * >( round_trip( *t.f ) );
  assert( r && ( r->get_num_active_var() == 0 ) );
  r->set_variables( { & t.y , & t.y2 } );
  auto g0 = static_cast< const LinearFunction * >(
					       r->get_Lagrangian_term( 0 ) );
  auto g1 = static_cast< const LinearFunction * >(
					       r->get_Lagrangian_term( 1 ) );
  assert( ( g0->get_num_active_var() == 1 ) &&
	  ( g0->get_v_var()[ 0 ].second == -4 ) &&
	  ( g0->get_constant_term() == 1.5 ) );
  assert( ( g1->get_num_active_var() == 0 ) &&
	  ( g1->get_constant_term() == -2 ) );
  delete r;

  Lagrangian e( 0 );
  auto n = dynamic_cast< LagBFunction * >( round_trip( *e.f ) );
  assert( n && ( n->get_num_active_var() == 0 ) );
  n->set_variables( {} );
  assert( n->get_num_active_var() == 0 );
  delete n;
  }

 std::cout << "LagBFunction: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* The constant terms of an AbstractBlock: that of the Objective, positive,
 * negative and alone, comes back as it is; that of a row comes back moved to
 * its sides, the row being the same. */

static void test_constants( void )
{
 for( double c0 : { 0.25 , -3.5 } ) {
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );
  block.set_objective( objective( & block , linear( *x , { 1 , -2 } , c0 ) ) ,
		       eNoMod );
  auto read = round_trip( block );
  assert( read && same( model_of( *read ) , model_of( block ) ) );
  assert( model_of( *read ).obj_const == c0 );
  delete read;
  }

 // a constant and no term: every cost is zero
 {
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );
  block.set_objective( objective( & block , linear( *x , { 0 , 0 } , 7 ) ) ,
		       eNoMod );
  auto read = round_trip( block );
  assert( read && same( model_of( *read ) , model_of( block ) ) );
  assert( model_of( *read ).obj_const == 7 );
  delete read;
  }

 // 1 <= x0 + x1 + 2 <= 5 is -1 <= x0 + x1 <= 3, and x0 - 0.5 = 1 is x0 = 1.5
 {
  AbstractBlock block;
  auto x = kept( new std::vector< ColVariable >( 2 ) );
  block.add_static_variable( *x , "x" );
  auto rows = kept( new std::vector< FRowConstraint >( 2 ) );
  ( *rows )[ 0 ].set_function( linear( *x , { 1 , 1 } , 2 ) );
  ( *rows )[ 0 ].set_lhs( 1 );
  ( *rows )[ 0 ].set_rhs( 5 );
  ( *rows )[ 1 ].set_function( linear( *x , { 1 , 0 } , -0.5 ) );
  ( *rows )[ 1 ].set_both( 1 );
  block.add_static_constraint( *rows , "r" );
  block.set_objective( objective( & block , linear( *x , { 1 , 1 } ) ) ,
		       eNoMod );

  auto read = round_trip( block );
  assert( read );
  const auto m = model_of( *read );
  assert( m.rows.size() == 3 );
  using Row = Model::Row;
  const std::vector< std::pair< Index , double > > both{ { 0 , 1 } ,
							 { 1 , 1 } };
  const std::vector< std::pair< Index , double > > first{ { 0 , 1 } };
  std::vector< Row > expected{ { both , -1 , Inf< double >() } ,
			       { both , - Inf< double >() , 3 } ,
			       { first , 1.5 , 1.5 } };
  std::sort( expected.begin() , expected.end() );
  assert( m.rows == expected );
  read->for_each_constraint_group( [ & ]( const BaseGroup & g ) {
    for( Index i = 0 ; i < g.get_num_elements() ; ++i )
     if( auto r = dynamic_cast< const FRowConstraint * >(
						    g.get_Constraint( i ) ) )
      assert( static_cast< const LinearFunction * >( r->get_function() )
	      ->get_constant_term() == 0 );
    } );
  delete read;
  }

 std::cout << "constant terms of an AbstractBlock: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* serialize( std::string , int , false ) appends to the file, creating it if
 * it is not there, and throws if it is there and cannot be appended to,
 * leaving it as it was rather than writing over it. */

static void test_append_to_file( void )
{
 const auto fn = ( std::filesystem::temp_directory_path() /
		   ( "smspp_NetCDF_test_append_" + std::to_string(
       std::chrono::steady_clock::now().time_since_epoch().count() ) +
		     ".nc4" ) ).string();

 AbstractBlock a;
 auto x = kept( new std::vector< ColVariable >( 1 ) );
 a.add_static_variable( *x , "x" );
 a.set_objective( objective( & a , linear( *x , { 2 } ) ) , eNoMod );
 AbstractBlock b;
 auto y = kept( new std::vector< ColVariable >( 2 ) );
 b.add_static_variable( *y , "y" );
 b.set_objective( objective( & b , linear( *y , { 1 , 3 } ) ) , eNoMod );

 // not there: created
 std::filesystem::remove( fn );
 a.Block::serialize( fn , eBlockFile , false );
 b.Block::serialize( fn , eBlockFile , false );
 {
  netCDF::NcFile f( fn , netCDF::NcFile::read );
  auto r0 = Block::deserialize( f , 0 );
  auto r1 = Block::deserialize( f , 1 );
  assert( r0 && same( model_of( *r0 ) , model_of( a ) ) );
  assert( r1 && same( model_of( *r1 ) , model_of( b ) ) );
  assert( ! Block::deserialize( f , 2 ) );
  delete r0;
  delete r1;
  }

 // there and not a netCDF file: an exception, and the file as it was
 {
  std::ofstream out( fn , std::ios::trunc );
  out << "not a netCDF file" << std::endl;
  }
 bool thrown = false;
 try { a.Block::serialize( fn , eBlockFile , false ); }
 catch( std::runtime_error & e ) {
  thrown = ( std::string( e.what() ).rfind( "Block::serialize: " , 0 ) == 0 );
  }
 assert( thrown );
 {
  std::ifstream in( fn );
  std::string line;
  std::getline( in , line );
  assert( line == "not a netCDF file" );
  }

 std::filesystem::remove( fn );
 std::cout << "appending to a file: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* A Block written to the image of a netCDF file in memory and read back
 * from it: an AbstractBlock with its rows, bounds and Objective and two
 * nested Block, an AbstractBlock and a PolyhedralFunctionBlock. The image
 * gives the same Block as the file serialize( std::string ) writes, it is a
 * netCDF file that can be written to disk and read from there, and the
 * bytes of a file read from disk are an image; writing the same Block twice
 * gives the same bytes; a Block
 * changed after it was built is read back with the changes, which is how a
 * copy of a Block as it is now is made; the Block of an image of type
 * eBlockFile are appended and read by position; and what is not an image,
 * or not the image of an SMS++ file, gives no Block. */

/// the AbstractBlock of the in-memory tests, with its two nested Block

static AbstractBlock * memory_block( std::vector< ColVariable > & y )
{
 auto block = new AbstractBlock;

 auto x = kept( new std::vector< ColVariable >( 4 ) );
 ( *x )[ 1 ].set_type( ColVariable::kNatural );
 block->add_static_variable( *x , "x" );

 auto rows = kept( new std::vector< FRowConstraint >( 3 ) );
 ( *rows )[ 0 ].set_function( linear( *x , { 1 , 2 , 0 , 0 } ) );
 ( *rows )[ 0 ].set_both( 3 );
 ( *rows )[ 1 ].set_function( linear( *x , { 0 , 0 , -3 , 0.5 } ) );
 ( *rows )[ 1 ].set_lhs( - Inf< double >() );
 ( *rows )[ 1 ].set_rhs( 7 );
 ( *rows )[ 2 ].set_function( linear( *x , { 1 , 0 , -1 , 0 } ) );
 ( *rows )[ 2 ].set_lhs( 1 );
 ( *rows )[ 2 ].set_rhs( 4 );
 block->add_static_constraint( *rows , "r" );

 auto box = kept( new BoxConstraint( block , & ( *x )[ 0 ] , -1 , 5 ) );
 block->add_static_constraint( *box , "box" );

 auto obj = objective( block , linear( *x , { 5 , -1 , 2 , 1 } ) );
 obj->set_sense( Objective::eMax , eNoMod );
 block->set_objective( obj , eNoMod );

 auto inner = new AbstractBlock( block );
 auto z = kept( new std::vector< ColVariable >( 2 ) );
 inner->add_static_variable( *z , "z" );
 auto irows = kept( new std::vector< FRowConstraint >( 1 ) );
 ( *irows )[ 0 ].set_function( linear( *z , { 1 , 1 } ) );
 ( *irows )[ 0 ].set_lhs( - Inf< double >() );
 ( *irows )[ 0 ].set_rhs( 4 );
 inner->add_static_constraint( *irows , "r" );
 inner->set_objective( objective( inner , linear( *z , { 0 , 1 } ) ) ,
		       eNoMod );
 block->add_nested_Block( inner );

 auto pfb = new PolyhedralFunctionBlock( block );
 auto & pf = pfb->get_PolyhedralFunction();
 pf.set_variables( { & y[ 0 ] , & y[ 1 ] , & y[ 2 ] } );
 pf.set_PolyhedralFunction( { { 1 , 2 , 0 } , { -3 , 0.5 , 4 } } ,
			    { 0 , 1.5 } , -10 , true , eNoMod );
 block->add_nested_Block( pfb );

 return( block );
 }

/// the models of an AbstractBlock and of its nested AbstractBlock

static std::vector< Model > models_of( const Block & b )
{
 std::vector< Model > m{ model_of( b ) };
 for( auto nb : b.get_nested_Blocks() )
  if( dynamic_cast< const AbstractBlock * >( nb ) &&
      ( ! dynamic_cast< const PolyhedralFunctionBlock * >( nb ) ) )
   m.push_back( model_of( *nb ) );
 return( m );
 }

/// true if a and b have the same models and the same PolyhedralFunction

static bool same_block( Block & a , Block & b )
{
 const auto ma = models_of( a );
 const auto mb = models_of( b );
 if( ma.size() != mb.size() )
  return( false );
 for( std::size_t i = 0 ; i < ma.size() ; ++i )
  if( ! same( ma[ i ] , mb[ i ] ) )
   return( false );
 if( a.get_number_nested_Blocks() != b.get_number_nested_Blocks() )
  return( false );
 auto pa = dynamic_cast< PolyhedralFunctionBlock * >(
					       a.get_nested_Blocks().back() );
 auto pb = dynamic_cast< PolyhedralFunctionBlock * >(
					       b.get_nested_Blocks().back() );
 if( ( ! pa ) || ( ! pb ) )
  return( false );
 auto & fa = pa->get_PolyhedralFunction();
 auto & fb = pb->get_PolyhedralFunction();
 return( ( fa.get_A() == fb.get_A() ) && ( fa.get_b() == fb.get_b() ) &&
	 ( fa.is_convex() == fb.is_convex() ) &&
	 ( fa.get_global_bound() == fb.get_global_bound() ) );
 }

/// the bytes of a file

static NcMemFile::Image bytes_of( const std::string & name )
{
 std::ifstream in( name , std::ios::binary );
 return( NcMemFile::Image( std::istreambuf_iterator< char >( in ) ,
			   std::istreambuf_iterator< char >() ) );
 }

static void test_in_memory( void )
{
 std::vector< ColVariable > y( 3 );
 auto block = memory_block( y );

 // the round trip through memory gives the same Block
 NcMemFile::Image image;
 block->Block::serialize( image );
 assert( ! image.empty() );
 const auto image0 = image;

 auto read = Block::deserialize( image );
 assert( dynamic_cast< AbstractBlock * >( read ) );
 assert( same_block( *read , *block ) );
 assert( read->get_nested_Blocks()[ 0 ]->get_f_Block() == read );

 // the bytes are read where they are, and left as they were
 assert( image == image0 );
 delete read;

 // writing the same Block again gives the same bytes
 NcMemFile::Image twice;
 block->Block::serialize( twice );
 assert( twice == image );

 // the file serialize( std::string ) writes gives the same Block
 block->Block::serialize( nc_file() , eProbFile );
 auto from_file = Block::deserialize( nc_file() );
 assert( from_file && same_block( *from_file , *block ) );
 delete from_file;

 // the image is a netCDF file: written to disk, it is read from there
 {
  std::ofstream out( nc_file() , std::ios::binary | std::ios::trunc );
  out.write( reinterpret_cast< const char * >( image.data() ) ,
	     image.size() );
  }
 auto from_disk = Block::deserialize( nc_file() );
 assert( from_disk && same_block( *from_disk , *block ) );
 delete from_disk;
 block->Block::serialize( nc_file() , eProbFile );

 // the bytes of a file written on disk are read from memory as well
 auto from_bytes = Block::deserialize( bytes_of( nc_file() ) );
 assert( from_bytes && same_block( *from_bytes , *block ) );
 delete from_bytes;

 // a Block changed after it was built: its image is a copy of it as it is
 // now, with the changes and with no Modification to be copied
 const auto before = models_of( *block );
 auto & x = *block->get_static_variable_v< ColVariable >( "x" );
 auto & r = *block->get_static_constraint_v< FRowConstraint >( "r" );
 r[ 1 ].set_rhs( 2 );
 static_cast< LinearFunction * >( r[ 2 ].get_function() )
					       ->modify_coefficient( 1 , -5 );
 auto box = block->get_static_constraint< BoxConstraint >( "box" );
 assert( box );
 box->set_rhs( 3 );
 block->get_objective()->set_sense( Objective::eMin );
 auto & ir = *block->get_nested_Blocks()[ 0 ]->
		      get_static_constraint_v< FRowConstraint >( "r" );
 ir[ 0 ].set_rhs( 9 );
 auto & pf = static_cast< PolyhedralFunctionBlock * >(
		  block->get_nested_Blocks()[ 1 ] )->get_PolyhedralFunction();
 pf.modify_constant( 0 , -7 );
 pf.add_row( { 0 , 0 , 1 } , 2 );
 x[ 3 ].set_type( ColVariable::kInteger );

 const auto after = models_of( *block );
 assert( ! same( after[ 0 ] , before[ 0 ] ) );
 assert( ! same( after[ 1 ] , before[ 1 ] ) );

 NcMemFile::Image changed;
 block->Block::serialize( changed );
 assert( changed != image0 );
 auto copy = Block::deserialize( changed );
 assert( copy && same_block( *copy , *block ) );
 const auto cm = models_of( *copy );
 assert( cm[ 0 ].integer[ 3 ] && ( cm[ 0 ].ub[ 0 ] == 3 ) );
 assert( cm[ 0 ].sense == Objective::eMin );
 auto & cpf = static_cast< PolyhedralFunctionBlock * >(
		   copy->get_nested_Blocks()[ 1 ] )->get_PolyhedralFunction();
 assert( ( cpf.get_nrows() == 3 ) && ( cpf.get_b()[ 0 ] == -7 ) );

 // the copy is a Block of its own: changing it leaves the original alone
 cpf.modify_constant( 1 , 5 );
 assert( pf.get_b()[ 1 ] == 1.5 );
 assert( same( models_of( *block )[ 0 ] , after[ 0 ] ) );
 delete copy;

 // and writing the changed Block again gives the same bytes again
 NcMemFile::Image changed2;
 block->Block::serialize( changed2 );
 assert( changed2 == changed );

 // two images opened at the same time
 {
  NcMemFile f1( image0 );
  NcMemFile f2( changed );
  auto b1 = Block::deserialize( f1 );
  auto b2 = Block::deserialize( f2 );
  assert( b1 && b2 && same_block( *b2 , *block ) );
  assert( same( models_of( *b1 )[ 0 ] , before[ 0 ] ) );
  delete b1;
  delete b2;
  }

 // eBlockFile: several Block written in the same NcMemFile, and read from
 // the image by their position
 {
  AbstractBlock small;
  auto v = kept( new std::vector< ColVariable >( 2 ) );
  small.add_static_variable( *v , "v" );
  small.set_objective( objective( & small , linear( *v , { 1 , -1 } ) ) ,
		       eNoMod );

  NcMemFile f;
  f.putAtt( "SMS++_file_type" , netCDF::NcInt() , eBlockFile );
  block->Block::serialize( f , eBlockFile );
  small.Block::serialize( f , eBlockFile );
  NcMemFile::Image many;
  f.close( many );

  auto b0 = Block::deserialize( many , 0 );
  auto b1 = Block::deserialize( many , 1 );
  assert( b0 && same_block( *b0 , *block ) );
  assert( b1 && same( model_of( *b1 ) , model_of( small ) ) );
  assert( ! Block::deserialize( many , 2 ) );
  delete b0;
  delete b1;

  // serialize( image ) starts a new image, with the Block alone
  small.Block::serialize( many , eBlockFile );
  auto only = Block::deserialize( many , 0 );
  assert( only && same( model_of( *only ) , model_of( small ) ) );
  assert( ! Block::deserialize( many , 1 ) );
  delete only;
  }

 // what is not an image gives no Block, and a wrong type throws leaving
 // the image as it was
 {
  NcMemFile::Image junk( 100 , 7 );
  assert( ! Block::deserialize( junk ) );
  assert( ! Block::deserialize( NcMemFile::Image() ) );

  NcMemFile::Image kept( image0 );
  bool thrown = false;
  try { block->Block::serialize( kept , 7 ); }
  catch( std::invalid_argument & ) { thrown = true; }
  assert( thrown && ( kept == image0 ) );

  // a netCDF image with no SMS++_file_type is not an SMS++ file
  NcMemFile f;
  f.addGroup( "Block_0" );
  NcMemFile::Image plain;
  f.close( plain );
  assert( ! plain.empty() );
  assert( ! Block::deserialize( plain ) );
  thrown = false;
  try { f.close( plain ); }
  catch( std::logic_error & ) { thrown = true; }
  assert( thrown );

  thrown = false;
  try { NcMemFile g( junk ); }
  catch( netCDF::exceptions::NcException & ) { thrown = true; }
  assert( thrown );
  }

 // NcMemFile on its own: any netCDF content, read back from the image,
 // whose bytes are left as they were also when the file is closed with
 // close( Image & )
 {
  NcMemFile f;
  auto d = f.addDim( "n" , 3 );
  auto var = f.addVar( "v" , netCDF::ncDouble , d );
  const double v[] = { 1.5 , -2 , 3 };
  var.putVar( v );
  f.putAtt( "a" , netCDF::NcInt() , 4 );
  NcMemFile::Image img;
  f.close( img );
  const auto img0 = img;

  NcMemFile h( img );
  double rv[ 3 ];
  h.getVar( "v" ).getVar( rv );
  assert( ( rv[ 0 ] == 1.5 ) && ( rv[ 1 ] == -2 ) && ( rv[ 2 ] == 3 ) );
  int a = 0;
  h.getAtt( "a" ).getValues( & a );
  assert( a == 4 );
  NcMemFile::Image back;
  h.close( back );
  assert( ( img == img0 ) && ( back == img0 ) );
  }

 delete block;
 std::cout << "Block in memory: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* copy_group() and write_via_memory() give what writing in place gives: the
 * same groups, dimensions, variables and attributes, in the same order,
 * whatever their types (strings, scalars, a dimension of length 0, an
 * unlimited one) and however many the sub-groups (enough of them for an
 * order by name to differ from the order of creation); a dimension of an
 * ancestor of the copied group is found by name in the destination, and its
 * absence throws, as a variable of a user-defined type does; a copy that
 * would take a file over the groups netCDF-C can hold throws and copies
 * nothing. */

// writes in g a part with some of everything, i telling the parts apart:
// in the values, and in an attribute if index

static void write_part( netCDF::NcGroup & g , int i , bool index = true )
{
 g.putAtt( "type" , "Part" );
 if( index )
  g.putAtt( "index" , netCDF::NcInt() , i );

 auto t = g.addDim( "T" , 5 );
 auto z = g.addDim( "Z" , 0 );
 auto u = g.addDim( "U" );  // unlimited

 const double x[] = { 1.5 + i , -2 , 3 , 0 , 1e-300 };
 auto vx = g.addVar( "x" , netCDF::NcDouble() , t );
 vx.putVar( x );
 vx.putAtt( "unit" , "MW" );
 vx.putAtt( "scale" , netCDF::NcDouble() , 0.5 );

 g.addVar( "k" , netCDF::NcInt() ).putVar( & i );
 g.addVar( "z" , netCDF::NcDouble() , z );

 const char * s[] = { "a" , "" , "ccc" , "dd" , "e" };
 g.addVar( "s" , netCDF::NcString() , t ).putVar( s );

 const double m[] = { 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 };
 g.addVar( "m" , netCDF::NcDouble() , { t , g.addDim( "R" , 2 ) } )
  .putVar( { 0 , 0 } , { 5 , 2 } , m );

 const int w[] = { 7 , 8 , 9 };
 g.addVar( "w" , netCDF::NcInt() , u ).putVar( { 0 } , { 3 } , w );

 for( int j = 0 ; j < 12 ; ++j ) {
  auto sg = g.addGroup( "Sub_" + std::to_string( j ) );
  sg.putAtt( "type" , "Sub" );
  auto n = sg.addDim( "N" , 2 );
  const double v[] = { double( i ) , double( j ) };
  sg.addVar( "v" , netCDF::NcDouble() , n ).putVar( v );
  }
 }

/*--------------------------------------------------------------------------*/
// a text with everything in grp, in the order netCDF gives it

static void dump( int grp , std::ostream & out )
{
 char name[ NC_MAX_NAME + 1 ];
 nc_inq_grpname( grp , name );
 out << "group " << name << "\n";

 auto atts = [ & ]( int var ) {
  int na = 0;
  nc_inq_varnatts( grp , var , & na );
  for( int a = 0 ; a < na ; ++a ) {
   nc_type type;
   std::size_t len;
   nc_inq_attname( grp , var , a , name );
   nc_inq_att( grp , var , name , & type , & len );
   out << " att " << name << " " << type << " " << len << ":";
   if( type == NC_CHAR ) {
    std::string v( len , ' ' );
    nc_get_att_text( grp , var , name , v.data() );
    out << v;
    }
   else {
    std::vector< double > v( len );
    nc_get_att_double( grp , var , name , v.data() );
    for( auto d : v )
     out << " " << d;
    }
   out << "\n";
   }
  };
 atts( NC_GLOBAL );

 int nd = 0;
 nc_inq_dimids( grp , & nd , nullptr , 0 );
 std::vector< int > dims( nd );
 nc_inq_dimids( grp , & nd , dims.data() , 0 );
 int nu = 0;
 nc_inq_unlimdims( grp , & nu , nullptr );
 std::vector< int > unlim( nu );
 nc_inq_unlimdims( grp , & nu , unlim.data() );
 for( auto d : dims ) {
  std::size_t len;
  nc_inq_dim( grp , d , name , & len );
  out << " dim " << name << " " << len
      << ( std::ranges::find( unlim , d ) != unlim.end() ? " unlimited" : "" )
      << "\n";
  }

 int nv = 0;
 nc_inq_varids( grp , & nv , nullptr );
 std::vector< int > vars( nv );
 nc_inq_varids( grp , & nv , vars.data() );
 for( auto v : vars ) {
  nc_type type;
  int ndims;
  int vdims[ NC_MAX_VAR_DIMS ];
  nc_inq_var( grp , v , name , & type , & ndims , vdims , nullptr );
  out << " var " << name << " " << type;
  std::size_t n = 1;
  for( int k = 0 ; k < ndims ; ++k ) {
   char dn[ NC_MAX_NAME + 1 ];
   std::size_t len;
   nc_inq_dim( grp , vdims[ k ] , dn , & len );
   out << " " << dn;
   n *= len;
   }
  out << ":";
  if( type == NC_STRING ) {
   std::vector< char * > s( n );
   if( n )
    nc_get_var_string( grp , v , s.data() );
   for( auto c : s )
    out << " '" << c << "'";
   nc_free_string( n , s.data() );
   }
  else {
   std::vector< double > d( n );
   if( n )
    nc_get_var_double( grp , v , d.data() );
   for( auto e : d )
    out << " " << e;
   }
  out << "\n";
  atts( v );
  }

 int ng = 0;
 nc_inq_grps( grp , & ng , nullptr );
 std::vector< int > grps( ng );
 nc_inq_grps( grp , & ng , grps.data() );
 for( auto g : grps )
  dump( g , out );
 }

static std::string dump( const netCDF::NcGroup & g )
{
 std::ostringstream out;
 dump( g.getId() , out );
 return( out.str() );
 }

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

static void test_copy_group( void )
{
 const auto base = std::filesystem::temp_directory_path() /
  ( "smspp_NetCDF_test_copy_" + std::to_string(
    std::chrono::steady_clock::now().time_since_epoch().count() ) );
 const auto fa = base.string() + "_a.nc4";
 const auto fb = base.string() + "_b.nc4";

 // 15 parts, written in place and by way of memory: the same file, also
 // once it is read back from disk
 std::vector< std::string > names;
 for( int i = 0 ; i < 15 ; ++i )
  names.push_back( "Part_" + std::to_string( i ) );
 {
  netCDF::NcFile a( fa , netCDF::NcFile::replace );
  for( int i = 0 ; i < 15 ; ++i ) {
   auto g = a.addGroup( names[ i ] );
   write_part( g , i );
   }
  netCDF::NcFile b( fb , netCDF::NcFile::replace );
  write_via_memory( b , names , [ & ]( std::size_t i , netCDF::NcGroup & g ) {
		     write_part( g , int( i ) ); } );
  assert( dump( a ) == dump( b ) );
  }
 {
  netCDF::NcFile a( fa , netCDF::NcFile::read );
  netCDF::NcFile b( fb , netCDF::NcFile::read );
  assert( dump( a ) == dump( b ) );
  assert( dump( a ).find( "group Sub_11" ) > dump( a ).find( "group Sub_2" ) );
  }

 // no parts: nothing written
 {
  NcMemFile f;
  write_via_memory( f , {} , []( std::size_t , netCDF::NcGroup & ) {
		     assert( false ); } );
  assert( f.getGroupCount() == 0 );
  }

 // a dimension of an ancestor of the copied group is found by name in the
 // destination and in its ancestors, and its absence throws
 {
  NcMemFile from;
  auto n = from.addDim( "N" , 3 );
  auto g = from.addGroup( "G" );
  const double v[] = { 1 , 2 , 3 };
  g.addVar( "v" , netCDF::NcDouble() , n ).putVar( v );

  NcMemFile to;
  to.addDim( "N" , 3 );
  auto h = to.addGroup( "H" ).addGroup( "K" );
  copy_group( g , h );
  double r[ 3 ];
  h.getVar( "v" ).getVar( r );
  assert( ( r[ 0 ] == 1 ) && ( r[ 1 ] == 2 ) && ( r[ 2 ] == 3 ) );
  assert( h.getVar( "v" ).getDim( 0 ).getSize() == 3 );

  NcMemFile none;
  auto k = none.addGroup( "K" );
  bool thrown = false;
  try { copy_group( g , k ); }
  catch( netCDF::exceptions::NcException & ) { thrown = true; }
  assert( thrown );
  }

 // a variable of a user-defined type throws
 {
  // with the C API, netCDF::NcGroup::addVlenType() throwing on its own
  // type in netcdf-cxx4 4.3.1
  NcMemFile from;
  nc_type vl;
  int w;
  nc_def_vlen( from.getId() , "vl" , NC_INT , & vl );
  nc_def_var( from.getId() , "w" , vl , 0 , nullptr , & w );
  NcMemFile to;
  bool thrown = false;
  try { copy_group( from , to ); }
  catch( std::invalid_argument & ) { thrown = true; }
  assert( thrown );
  }

 // the groups of a file: a copy that would take it over 32767 throws, and
 // the file keeps what it had
 {
  NcMemFile to;
  for( int j = 1 ; j < 32760 ; ++j )
   to.addGroup( "G_" + std::to_string( j ) );
  NcMemFile from;
  for( int j = 0 ; j < 8 ; ++j )
   from.addGroup( "S_" + std::to_string( j ) );
  auto before = to.getGroupCount();
  bool thrown = false;
  try { copy_group( from , to ); }
  catch( std::length_error & ) { thrown = true; }
  assert( thrown && ( to.getGroupCount() == before ) );

  // one group fewer fits exactly, the root included
  NcMemFile fits;
  for( int j = 0 ; j < 7 ; ++j )
   fits.addGroup( "S_" + std::to_string( j ) );
  copy_group( fits , to );
  assert( to.getGroupCount() == before + 7 );
  }

 std::filesystem::remove( fa );
 std::filesystem::remove( fb );
 std::cout << "copy_group and write_via_memory: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/* write_stacked() writes identical parts (differing in the values of their
 * variables only) as one group with the dimension stacked_dimension in
 * front of every variable, and read_stacked() gives back each of them as it
 * was written, strings, scalars, an empty and an
 * unlimited dimension included, and so does it with a stack nested in the
 * parts, whose dimension has another name; parts of several trees are written in one stack per tree, with
 * the index of where each part is; parts whose root has a dimension
 * stacked_dimension are written as write_via_memory() writes them; parts that together take more groups than a file holds fit
 * in one when stacked; a part out of the stack throws. */

static void test_stacked( void )
{
 // the 15 parts of write_part(), each as it is written in place
 std::vector< std::string > names;
 std::vector< std::string > alone;
 for( int i = 0 ; i < 15 ; ++i ) {
  names.push_back( "Part_" + std::to_string( i ) );
  NcMemFile f;
  auto g = f.addGroup( "part" );
  write_part( g , i , false );
  alone.push_back( dump( g ) );
  }

 {
  NcMemFile f;
  assert( write_stacked( f , "Parts" , names ,
			 [ & ]( std::size_t i , netCDF::NcGroup & g ) {
			  write_part( g , int( i ) , false ); } ) );
  auto s = f.getGroup( "Parts" );
  assert( f.getGroupCount() == 1 );
  assert( stacked_size( s ) == 15 );
  assert( s.getVar( "x" ).getDim( 0 ).getName() == stacked_dimension );
  for( std::size_t i = 0 ; i < 15 ; ++i )
   read_stacked( s , i , [ & ]( const netCDF::NcGroup & g ) {
    assert( dump( g ) == alone[ i ] ); } );
  bool thrown = false;
  try { read_stacked( s , 15 , []( const netCDF::NcGroup & ) {} ); }
  catch( std::out_of_range & ) { thrown = true; }
  assert( thrown );
  assert( stacked_size( f ) == 0 );
  }

 // parts of three trees, mixed: one stack per tree, in the order of their
 // first part, and each part given back as it was written
 {
  auto write = [ & ]( std::size_t i , netCDF::NcGroup & g ) {
   write_part( g , int( i ) , false );
   if( i % 3 == 1 )
    g.addVar( "extra" , netCDF::NcInt() );
   if( i % 3 == 2 )
    g.addGroup( "More" ).addDim( "M" , i );
   };
  NcMemFile f;
  assert( write_stacked( f , "Parts" , names , write ) );
  auto s = f.getGroup( "Parts" );
  assert( ( f.getGroupCount() == 1 ) && ( stacked_size( s ) == 15 ) );
  assert( stacked_size( s.getGroup( "Stack_0" ) ) == 5 );
  assert( stacked_size( s.getGroup( "Stack_1" ) ) == 5 );
  assert( s.getGroup( "Stack_1" ).getVar( "extra" ).getDimCount() == 1 );
  // the dimension M differs in length, so every part with More is apart
  assert( s.getGroupCount() == 2 + 5 );
  for( std::size_t i = 0 ; i < 15 ; ++i )
   read_stacked( s , i , [ & ]( const netCDF::NcGroup & g ) {
    NcMemFile e;
    auto r = e.addGroup( "part" );
    write( i , r );
    assert( dump( g ) == dump( r ) ); } );
  }

 // no parts: nothing written
 {
  NcMemFile f;
  assert( ! write_stacked( f , "Parts" , {} ,
			   []( std::size_t , netCDF::NcGroup & ) {} ) );
  assert( f.getGroupCount() == 0 );
  }

 // parts whose root already has the dimension: the parts as they are
 {
  NcMemFile f;
  assert( ! write_stacked( f , "Parts" , { "A" , "B" } ,
			   [ & ]( std::size_t , netCDF::NcGroup & g ) {
			    g.addDim( stacked_dimension , 2 ); } ) );
  assert( f.getGroupCount() == 2 );
  }

 // parts with a stack nested in them: stacked, each part given back with
 // its nested stack, and each element of this as it was written
 {
  auto write = [ & ]( std::size_t i , netCDF::NcGroup & g ) {
   g.putAtt( "type" , "Outer" );
   write_stacked( g , "Inner" , { "a" , "b" , "c" } ,
		  [ & ]( std::size_t j , netCDF::NcGroup & h ) {
		   write_part( h , int( 10 * i + j ) , false ); } );
   };
  NcMemFile f;
  assert( write_stacked( f , "Parts" , names , write ) );
  auto s = f.getGroup( "Parts" );
  assert( s.getVar( "x" ).isNull() );
  // the two dimensions of the stacks have different names, the outer one
  // first
  auto x = s.getGroup( "Inner" ).getVar( "x" );
  assert( ( x.getDimCount() == 3 ) && ( x.getDim( 0 ).getSize() == 15 ) &&
	  ( x.getDim( 1 ).getSize() == 3 ) );
  assert( ( x.getDim( 0 ).getName() == std::string( stacked_dimension ) +
	    "_1" ) && ( x.getDim( 1 ).getName() == stacked_dimension ) );
  assert( stacked_size( s ) == 15 );
  for( std::size_t i = 0 ; i < 15 ; ++i )
   read_stacked( s , i , [ & ]( const netCDF::NcGroup & g ) {
    auto inner = g.getGroup( "Inner" );
    assert( stacked_size( inner ) == 3 );
    for( std::size_t j = 0 ; j < 3 ; ++j )
     read_stacked( inner , j , [ & ]( const netCDF::NcGroup & h ) {
      NcMemFile e;
      auto r = e.addGroup( "part" );
      write_part( r , int( 10 * i + j ) , false );
      assert( dump( h ) == dump( r ) ); } );
    } );
  }

 // 40 parts of 1000 groups each: 40040 groups written apart, which no file
 // holds, and 1001 stacked
 {
  std::vector< std::string > many;
  for( int i = 0 ; i < 40 ; ++i )
   many.push_back( "P_" + std::to_string( i ) );
  auto write = [ & ]( std::size_t i , netCDF::NcGroup & g ) {
   for( int j = 0 ; j < 1000 ; ++j ) {
    auto sg = g.addGroup( "U_" + std::to_string( j ) );
    const double v = double( i ) + j;
    sg.addVar( "v" , netCDF::NcDouble() ).putVar( & v );
    }
   };
  NcMemFile f;
  assert( write_stacked( f , "Parts" , many , write ) );
  auto s = f.getGroup( "Parts" );
  assert( s.getGroupCount() == 1000 );
  read_stacked( s , 39 , [ & ]( const netCDF::NcGroup & g ) {
   double v;
   g.getGroup( "U_999" ).getVar( "v" ).getVar( & v );
   assert( v == 39 + 999 ); } );
  }

 std::cout << "write_stacked and read_stacked: OK" << std::endl;
 }

/*--------------------------------------------------------------------------*/
/*----------------------------------- MAIN ---------------------------------*/
/*--------------------------------------------------------------------------*/

int main( void )
{
 test_abstract_block();
 test_empty_abstract_block();
 test_polyhedral_function_block();
 test_polyhedral_function_state();
 test_benders_function();
 test_benders_state();
 test_lagrangian_function();
 test_constants();
 test_append_to_file();
 test_in_memory();
 test_copy_group();
 test_stacked();

 while( ! given.empty() )
  given.pop_back();

 std::filesystem::remove( nc_file() );

 std::cout << "all tests passed" << std::endl;

 return( 0 );
 }

/*--------------------------------------------------------------------------*/
/*------------------------ End File tests_NetCDF.cpp -----------------------*/
/*--------------------------------------------------------------------------*/