Model.cc

Go to the documentation of this file.

/*
 * $Id: Model.cc 238 2014-11-30 14:11:20Z jdl3 $
 * Copyright (C) 2011–2013, 2014 John D Lamb
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

#ifdef HAVE_CONFIG_H
#  include<config.h>
#endif

#include"Model.hpp"
#include"detail/BarrierFunction.hpp"
#include"detail/PhaseIModel.hpp"
#include"LinearConstraint.hpp"

using namespace ipo;

Model::DirectSubvector::DirectSubvector( size_t const offset, size_t const length ) 
  : offset { offset }, length { length }{}

Model::IndirectSubvector::IndirectSubvector( size_t const size )
  : std::vector<size_t>( size ){}

void
Model::setObjective( Objective& objective ){
  this->objective = objective;
}

Objective
Model::getObjective(){
  return objective;
}

void
Model::addConstraint( Constraint& constraint ){
  auto lc = dynamic_cast<LinearConstraint*>( &constraint );
  if( nullptr != lc and lc->isEqualityConstraint() )
    equalityConstraints.push_back( *lc );
  else
    constraints.push_back( constraint );
}

bool
Model::removeConstraint( Constraint& constraint ){
  /* Note the importance here of ensuring that a LinearCombination is allowed in
   * equalityConstraints if and only if it is an equality constaint and allowed in
   * constraints if and only if it is not. To ensure consistency, LinearCombination must
   * call this function before making a change.
   */
  auto lc = dynamic_cast<LinearConstraint*>( &constraint );
  if( nullptr != lc and lc->isEqualityConstraint() ){
    auto iter = std::find( equalityConstraints.begin(), equalityConstraints.end(), *lc );
    if( iter != equalityConstraints.end() ){
      equalityConstraints.erase( iter );
      return true;
    }
  } else {
    auto iter = std::find( constraints.begin(), constraints.end(), constraint );
    if( iter != constraints.end() ){
      constraints.erase( iter );
      return true;
    }
  }
  // constraint not in either list
  return false;
}

std::vector<Constraint> const&
Model::getConstraints() const {
  return constraints;
}

std::vector<LinearConstraint> const&
Model::getEqualityConstraints() const {
  return equalityConstraints;
}

void
Model::setIndices(){
  if( not variablesChanged ) return;
  // Check that numbers of variables supplied make sense
  try {
    testSizes();
  } catch( IPOException& ){
    RETHROW( "ipo::Model::setIndices()" );
  }
  // Clear indexMap so that any removed variables are no longer indexed.
  indexMap.clear();
  // Add objective variables
  size_t const OSIZE { objective.getVariables().size() }; // identifies vars not in objective
  for( auto& variable : objective.getVariables() ) indexMap[variable] = OSIZE;
  bool const DIRECTSUBVECTORFOROBJECTIVEINDICES { OSIZE == indexMap.size() };
  // Add constraint variables
  for( auto& constraint : constraints )
    for( auto& variable : constraint.getVariables() ) indexMap[variable] = OSIZE;
  for( auto& constraint : equalityConstraints )
    for( auto& variable : constraint.getVariables() ) indexMap[variable] = OSIZE;
  // Map keys to unique index value
  size_t index { 0 };
  if( DIRECTSUBVECTORFOROBJECTIVEINDICES )
    for( auto& variable : objective.getVariables() )
      indexMap[variable] = index++;
  for( auto& entry : indexMap ) if( OSIZE == entry.second ) entry.second = index++;
  // indexMap complete

  // Set objectiveIndices
  setObjectiveIndices();
  
  // Set constraintIndices (including constraints with no upper or lower bound)
  setConstraintIndices();

  // Now indices are up to date; no need to check again
  variablesChanged = false;

}

void Model::notify(){ variablesChanged = true; }

Model::Model()
  : objective( *this, new ipo_function::concrete::NullFunction() ){} 

bool
Model::isFeasible(){
  setIndices();
  // create a vector of suitable size: essentially reserves capacity
  gsl::vector superVector( indexMap.size() );
  // Copy values of variables to superVector
  for( auto& p : indexMap )
    superVector[p.second] = p.first.getValue();

  // Check variable bounds and copy values of variables to superVector
  for( auto& entry : indexMap ){
    Variable const& variable { entry.first };
    double value { variable.getValue() };
    if( value <= variable.getLowerBound() or value >= variable.getUpperBound() )
      return false;
    superVector[entry.second] = value;
  }

  // Check each constraint in turn
  for( size_t index { 0 }; index < constraints.size(); ++index ){
    auto& constraint = constraints[index];
    // Create subvector
    gsl::vector vector { constraintIndices[index]->get( superVector ) };
    double value { (*constraint.getFunction())( vector ) };
    if( value < constraint.getLowerBound() or value > constraint.getUpperBound() )
      return false;
  }
  // Repeat for equality constraints
  for( size_t index { 0 }; index < equalityConstraints.size(); ++index ){
    auto& equalityConstraint = equalityConstraints[index];
    // Create subvector
    gsl::vector vector { constraintIndices[index]->get( superVector ) };
    double value { (*equalityConstraint.getFunction())( vector ) };
    if( value != equalityConstraint.getUpperBound() ) // == constraint.getLowerBound()
      return false;
  }
  return true;
}

bool
Model::isStrictlyFeasible(){
  setIndices();
  // create a vector of suitable size: essentially reserves capacity
  gsl::vector superVector( indexMap.size() );
  // Copy values of variables to superVector
  for( auto& p : indexMap )
    superVector[p.second] = p.first.getValue();

  // Check variable bounds and copy values of variables to superVector
  for( auto& entry : indexMap ){
    Variable const& variable { entry.first };
    double value { variable.getValue() };
    if( value <= variable.getLowerBound() or value >= variable.getUpperBound() )
      return false;
    superVector[entry.second] = value;
  }

  // Check each constraint in turn
  for( size_t index { 0 }; index < constraints.size(); ++index ){
    auto& constraint = constraints[index];
    // Create subvector
    gsl::vector vector { constraintIndices[index]->get( superVector ) };
    double value { (*constraint.getFunction())( vector ) };
    if( value <= constraint.getLowerBound() or value >= constraint.getUpperBound() )
      return false;
  }
  // Repeat for equality constraints
  for( size_t index { 0 }; index < equalityConstraints.size(); ++index ){
    auto& equalityConstraint = equalityConstraints[index];
    // Create subvector
    gsl::vector vector { equalityConstraintIndices[index]->get( superVector ) };
    double value { (*equalityConstraint.getFunction())( vector ) };
    if( value != equalityConstraint.getUpperBound() ) // == constraint.getLowerBound()
      return false;
  }
  return true;
}

gsl::vector
Model::getVectorFromVariables() const {
  // Should usually be called after setIndices()
  size_t const SIZE = indexMap.size();
  gsl::vector vector( SIZE );
  for( auto& entry : indexMap )
    vector[entry.second] = entry.first.getValue();
  return vector;
}

void
Model::setVariablesFromVector( gsl::vector const & vector ){
  size_t const SIZE { indexMap.size() };
  if( vector.size() != SIZE ){
    std::ostringstream error;
    error << "ipo::Model::setVariablesFromVector(): ";
    error << "vector size is " << vector.size();
    error << ", but there are " << SIZE << " variables.";
    throw IPOE( error.str() );
  }
  for( auto& entry : indexMap )
    const_cast<Variable&>( entry.first ).setValue( vector[entry.second] );
}

bool
Model::stoppingCriterion( double const t ) const {
  // First identify number of constraints
  double numConstraints { static_cast<double>( constraints.size()
                           + equalityConstraints.size() ) };
  numConstraints = std::max( parameters.getEpsilon(), numConstraints );
  return numConstraints < t * parameters.getEpsilon();
}

bool
Model::optimise( bool const minimise ){
  //First find feasible solution
  if( not findFeasible() ) return false;
  
  // Set up
  auto const outputStream = parameters.getOutputStream();
  // Try to minimise barrier function
  try {
    // Create a BarrierFunction
    detail::BarrierFunction barrierFunction { *this, minimise ?
    detail::BarrierFunction::Sign::positive : detail::BarrierFunction::Sign::negative };
    // Get variables
    gsl::vector vector { getVectorFromVariables() };

    // Create a matrix and vector for equality constraints
    size_t const EQUALITYCONSTRAINTSSIZE { equalityConstraints.size() };
    gsl::matrix A { EQUALITYCONSTRAINTSSIZE, vector.size() };
    gsl::vector bounds( EQUALITYCONSTRAINTSSIZE );
    for( size_t i { 0 }; i < EQUALITYCONSTRAINTSSIZE; ++i ){
      // Allow for possibility that equality constraint does not contain all variables
      auto row = A.row( i );
      auto& vars = equalityConstraints[i].getVariables();
      for( auto& var : vars )
    row[indexMap[var]] = equalityConstraints[i].getCoefficient( var );
      bounds[i] = equalityConstraints[i].getUpperBound(); // also = lowerBound
    }
    // Create NewtonDescent object
    detail::NewtonDescent newtonDescent { barrierFunction, barrierFunction, A, bounds };

    // Set parameters
    newtonDescent.getParameters() = getNewtonDescentParameters();
    newtonDescent.getLineSearchParameters() = getLineSearchParameters();

    for( double t { 0.1 }; not stoppingCriterion( t ); t *= parameters.getMu() ){
      barrierFunction.set_t( t );
      
      if( nullptr != outputStream ) *outputStream << " t = " << t << std::endl;
      // Scales NewtonDescent and LineSearch output to approximate function values
      // Affects any output shown rather than optimisation
      newtonDescent.setFunctionScale( (minimise ? 1 : -1) / t );
      double const objectiveValue { barrierFunction( vector ) };
      newtonDescent( vector, objectiveValue ); // updates vector
      // early exit possible in phase I
      // std::cout << "* t = " << t << ", s = " << *vector.rbegin();
      // if( nullptr != dynamic_cast<detail::PhaseIModel*>( this ) )
      //    std::cout << " in phase 1" << std::endl;
      // else
      //    std::cout << " in phase 2" << std::endl;
      if( nullptr != dynamic_cast<detail::PhaseIModel*>( this ) and *vector.rbegin() < 0 )
    break;
    }
    
    // Set variables
    setVariablesFromVector( vector );
    // And return
    return true;

  } catch( ... ){
    RETHROW( "ipo::Model::optimise()" );
  }
}



Model::Parameters::Parameters(){
  mu = 10;
  epsilon = std::sqrt( std::numeric_limits<double>::epsilon() );
  outputStream = nullptr;
}

void
Model::Parameters::setMu( double const mu ){ if( mu > 1 ) this->mu = mu; }

void
Model::Parameters::setEpsilon( double const epsilon ){ if( epsilon > 0 )
    this->epsilon = epsilon; }

void
Model::Parameters::setOutputStream( std::ostream* outputStream ){
  this->outputStream = outputStream; }

detail::NewtonDescent::Parameters&
Model::getNewtonDescentParameters(){
  return newtonDescentParameters;
}

detail::LineSearch::Parameters&
Model::getLineSearchParameters(){
  return lineSearchParameters;
}

Model::Parameters&
Model::getParameters(){ return parameters; }

bool Model::checkSize() const {
  if( objective.checkSize() == false ) return false;
  for( auto& constraint : constraints )
    if( constraint.checkSize() == false ) return false;
  return true;
}

void
Model::testSizes() const {
  std::ostringstream ost;

  auto sizes = objective.getSizes();
  size_t const sizeF { std::get<0>( sizes ) };
  size_t const sizeV { std::get<1>( sizes ) };
  if( sizeF == 0 ){
    if( sizeV == 0 )
      ost << "Objective function " << objective.getName() << " has no variables." << std::endl;
  } else {
    if( sizeV != sizeF )
      ost << "Objective function " << objective.getName()
      << " needs " << sizeF << " variables, "
      << "but " << sizeV << " supplied." << std::endl;
  }
  for( auto& constraint : constraints ){
    auto sizes = constraint.getSizes();
    size_t const sizeF { std::get<0>( sizes ) };
    size_t const sizeV { std::get<1>( sizes ) };
    if( sizeF == 0 ){
      if( sizeV == 0 )
    ost << "Constraint " << constraint.getName() << " has no variables." << std::endl;
    } else {
      if( sizeV != sizeF )
    ost << "Constraint " << constraint.getName()
        << " needs " << sizeF << " variables, "
        << " but " << sizeV << " supplied." << std::endl;
    }
  }

  std::string const str { ost.str() };
  if( str.size() > 0 ){
    std::ostringstream ost { "ipo::Model::testSizes(): " };
    ost << "errors in objective and constraints:";
    ost << std::endl << str;
    throw IPOException { ost.str(), __FILE__, __LINE__ };
  }
}

void
Model::summary( std::ostream& ostream ) const {
  // Prefix
  std::string const prefix { " " };
  ostream << "Optimise" << std::endl;
  const_cast<Model&>( *this ).getObjective().summary( ostream, prefix );
  ostream << "subject to" << std::endl;
  
  // For collecting variables: not particularly efficient, but this should not matter here
  Array variables { const_cast<Model&>( *this ), 0, "variables" };

  // Add the model variables
  for( auto const& variable : const_cast<Model&>( *this ).getObjective().getVariables() ){
    if( not variables.contains( variable ) ) variables.push_back( variable );
  }
  
  // Summarise inequality constraints
  for( auto const& constraint : const_cast<Model&>( *this ).getConstraints() ){
    constraint.summary( ostream, prefix );
    for( auto const& variable : const_cast<Constraint&>( constraint ).getVariables() ){
      if( not variables.contains( variable ) ) variables.push_back( variable );
    }
  }
  
  // Summarise equality constraints
  for( auto const& constraint : const_cast<Model&>( *this ).getEqualityConstraints() ){
    constraint.summary( ostream, prefix );
    for( auto const& variable : const_cast<LinearConstraint&>( constraint ).getVariables() ){
      if( not variables.contains( variable ) ) variables.push_back( variable );
    }
  }
  
  // list variables
  ostream << "over";
  variables.summary( ostream, prefix );
  ostream << std::endl;
}

void
Model::setConstraintIndices(){  // local function
  auto indirectSubvector = [this]( Constraint& constraint, size_t const index,
                   bool const eqConstraint )->void {
    size_t const CSIZE { constraint.getVariables().size() };
    IndirectSubvector* cIndices = new IndirectSubvector( CSIZE );
    // Now fill with values
    size_t i = 0;
    for( auto& variable : constraint.getVariables() ){
      (*cIndices)[i] = this->indexMap[variable];
      ++i;
    }
    if( eqConstraint )
      equalityConstraintIndices[index] = std::shared_ptr<Subvector>( cIndices );
    else
      constraintIndices[index] = std::shared_ptr<Subvector>( cIndices );
  };

  size_t const OSIZE { objective.getVariables().size() }; /* identifies vars not in
                               * objective */
  bool const DIRECTSUBVECTORFOROBJECTIVEINDICES { OSIZE == indexMap.size() };
  // Resize
  constraintIndices.resize( constraints.size() );
  size_t index { 0 };
  // Set indices
  for( auto& constraint : constraints ){
    if( not DIRECTSUBVECTORFOROBJECTIVEINDICES )
      indirectSubvector( constraint, index, false );
    else if( constraint.getVariables().size() == 0 )
      // unlikely case; but this produces empty vector on demand
      indirectSubvector( constraint, index, false );
    else if( constraint.getVariables().size() > 1
         and constraint.getVariables().size() < OSIZE )
      indirectSubvector( constraint, index, false );
    else {
      if( constraint.getVariables().size() == 1 ){
    // constraint variable must exist in objective variables
    Variable& variable = constraint.getVariables().front();
    constraintIndices[index]
      = std::shared_ptr<Subvector>( new DirectSubvector( indexMap[variable], 1 ) );
      } else { // constraint.size() == OSIZE
    // Check variables are in the right order
    bool correctOrder { true };
    auto& variables = constraint.getVariables();
    for( size_t j { 0 }; j < variables.size(); ++j ){
      auto& variable = variables[j];
      if( indexMap[variable] != j ){
        correctOrder = false;
        break;
      }
    }
    if( correctOrder )
      constraintIndices[index]
        = std::shared_ptr<Subvector>( new DirectSubvector( 0, OSIZE ) );
    else
      indirectSubvector( constraint, index, false );
      }
    }
    ++index;
  } 
  // Repeat for equality constraints
  // Resize
  equalityConstraintIndices.resize( equalityConstraints.size() );
  index = 0;
  // Set indices
  for( auto& constraint : equalityConstraints ){
    if( not DIRECTSUBVECTORFOROBJECTIVEINDICES )
      indirectSubvector( constraint, index, true );
    else if( constraint.getVariables().size() == 0 )
      // unlikely case; but this produces empty vector on demand
      indirectSubvector( constraint, index, true );
    else if( constraint.getVariables().size() > 1
         and constraint.getVariables().size() < OSIZE )
      indirectSubvector( constraint, index, true );
    else {
      if( constraint.getVariables().size() == 1 ){
    // constraint variable must exist in objective variables
    Variable& variable = constraint.getVariables().front();
    equalityConstraintIndices[index]
      = std::shared_ptr<Subvector>( new DirectSubvector( indexMap[variable], 1 ) );
      } else { // constraint.size() == OSIZE
    // Check variables are in the right order
    bool correctOrder { true };
    auto& variables = constraint.getVariables();
    for( size_t j { 0 }; j < variables.size(); ++j ){
      auto& variable = variables[j];
      if( indexMap[variable] != j ){
        correctOrder = false;
        break;
      }
    }
    if( correctOrder )
      equalityConstraintIndices[index]
        = std::shared_ptr<Subvector>( new DirectSubvector( 0, OSIZE ) );
    else
      indirectSubvector( constraint, index, true );
      }
    }
    ++index;
  } 
}

bool
Model::findFeasible(){
  // First check for feasibility (but only if not PhaseIModel)
  if( nullptr == dynamic_cast<detail::PhaseIModel*>( this ) and not isStrictlyFeasible() ){
    detail::PhaseIModel phaseIModel( *this );
    // First find feasible solution for equality constraints
    if( not phaseIModel.findFeasibleSolution() ) return false;
    // Model now has a feasible solution.
  }
  return true;
}