PhaseIModel.cc

Go to the documentation of this file.

/*
 * $Id: PhaseIModel.cc 216 2014-11-04 09:00:49Z 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"PhaseIModel.hpp"

using namespace ipo::detail;

PhaseIModel::PhaseIModel( Model& model ) : model( model ){
  parameters = model.parameters;
  newtonDescentParameters = model.newtonDescentParameters;
  lineSearchParameters = model.lineSearchParameters;
  variablesChanged = true;
  setIndices(); // force model to add and recognise constraints, objective and variables
}

double
PhaseIModel::getInitial_s( gsl::vector const& vector ) const {
  // Make sure model is in sane state
  model.setIndices();

  if( not vector or vector.size() != model.indexMap.size() )
    throw IPOE( "PhaseIModel::getInitial_s(): vector size error." );

  double s { minusInfinity };
  
  // Check each constraint in base model
  for( size_t j { 0 }; j < model.constraints.size(); ++j ){
    auto& constraint = model.constraints[j];
    double const upperBound { constraint.getUpperBound() };
    double const lowerBound { constraint.getLowerBound() };
    if( upperBound < infinity or lowerBound > minusInfinity ){
      auto& constraintIndices_j = model.constraintIndices[j];
      auto subvector = constraintIndices_j->get( vector );
      double const functionValue { (*constraint.getFunction())( subvector )};
      if( upperBound < infinity )
    s = std::max( s, functionValue - upperBound );
      double const lowerBound { constraint.getLowerBound() };
      if( lowerBound > minusInfinity )
    s = std::max( s, lowerBound - functionValue );
    }
  }

  // Also each variable in base model
  for( auto& entry : model.indexMap ){
    auto& variable = entry.first;
    size_t index = entry.second;
    double const upperBound { variable.getUpperBound() };
    if( upperBound < infinity )
      s = std::max( s, vector[index] - upperBound );
    double const lowerBound { variable.getLowerBound() };
    if( lowerBound > minusInfinity )
      s = std::max( s, lowerBound - vector[index] );
  }
  
  // replace s with a strictly larger value.
  if( s < 1000 ) s += 0.1;
  else s *= 1.001;

  return s;
}

void 
PhaseIModel::setIndices(){
  if( model.variablesChanged == false and variablesChanged == false ) return; // nothing to do
  model.setIndices(); // Make sure model is up to date.

  // clear anything old out of this
  constraints.clear();
  equalityConstraints.clear();
  constraintIndices.clear();
  indexMap.clear();

  // We'll need to be able to look up Variable of this from Variable of model
  crossIndex.clear();

  // First create Variables for this
  for( auto& entry : model.indexMap ){
    auto& modelVariable = entry.first;
    size_t const INDEX { entry.second };
    // Create new unbounded variable
    Variable variable { *this, modelVariable.getName() };
    // Set to value of base variable
    variable.setValue( modelVariable.getValue() );
    // Add to indexMap and to crossIndex
    indexMap.insert( std::make_pair( variable, INDEX ) );
    crossIndex.insert( std::make_pair( modelVariable, variable ) );
  }
  // Add new variable for phase I
  Variable s { *this, "s: phase I" };
  size_t const INDEX = indexMap.size();
  indexMap.insert( std::make_pair( s, INDEX ) );

  // Now create an objective function
  ipo_function::Function const& objectiveFunction = *model.objective.getFunction();
  objective = Objective( *this, new ipo_function::detail::
             PhaseIObjectiveFunctionAndDerivatives( objectiveFunction ),
             "Phase I objective" );
  // Add variables to the new objective
  auto& modelObjVars = model.objective.getVariables();
  Array vars { *this, modelObjVars.getName() };
  for( auto& var : modelObjVars ) vars.push_back( crossIndex.find( var )->second );
  // Add new variable for phase I
  vars.push_back( s );
  objective.getVariables() = vars;

  // Create up to two constraints for each constraint in model
  for( auto& constraint : model.constraints ){
    ipo_function::Function& constraintFunction = *constraint.getFunction();
    double const upperBound = constraint.getUpperBound();
    if( upperBound < infinity ){
      // create a new constraint
      std::string name { constraint.getName() };
      name += " (upper bound)";
      Constraint c { *this, new ipo_function::detail::
	  PhaseIFunctionAndDerivatives { constraintFunction, upperBound, true }, name };
      // Add variables to the new constraint
      auto& cVars = constraint.getVariables();
      Array vars { *this, cVars.getName() };
      for( auto& var : cVars ) vars.push_back( crossIndex.find( var )->second );
      // Add new variable for phase I
      vars.push_back( s );
      c.getVariables() = vars;
      // Set upper bound to zero
      c.setUpperBound( 0.0 );
      // Add the constraint to the model
      constraints.push_back( c );
    }
    double const lowerBound { constraint.getLowerBound() };
    if( lowerBound > minusInfinity ){
      // create a new constraint
      std::string name { constraint.getName() };
      name += " (lower bound)";
      Constraint c { *this, new ipo_function::detail::
	  PhaseIFunctionAndDerivatives { constraintFunction, lowerBound, false }, name };
      // Add variables to the new constraint
      auto cVars = constraint.getVariables();
      Array vars { *this, cVars.getName() };
      for( auto& var : cVars ) vars.push_back( crossIndex.find( var )->second );
      // Add new variable for phase I
      vars.push_back( s );
      // And assign
      c.getVariables() = vars;
      // Set upper bound to zero
      c.setUpperBound( 0.0 );
      // Add the constraint to the model
      constraints.push_back( c );
    }
  }
  
  // Create a constraint for each equalityConstraint in model
  for( auto& equalityConstraint : model.equalityConstraints ){
    // create a new constraint
    LinearConstraint linearConstraint( *this, equalityConstraint.getName() );
    // Add variables to the new constraint
    auto cVars = equalityConstraint.getVariables();
    Array vars( *this, cVars.getName() );
    for( auto& var : cVars ) vars.push_back( crossIndex.find( var )->second );
    // Add new variable for phase I
    vars.push_back( s );
    // And assign
    linearConstraint.getVariables() = vars;
    // Set coefficients
    linearConstraint.setCoefficient( s, 0 ); // should be redundant
    for( auto& variable : cVars )
      linearConstraint.setCoefficient( crossIndex.find( variable )->second,
                       equalityConstraint.getCoefficient( variable ) );
    // Set bound
    linearConstraint.setValue( equalityConstraint.getUpperBound() );
    // Add the constraint to the model
    equalityConstraints.push_back( linearConstraint );
  }

  // Create up to two constraints for each variable constraint in model
  for( auto& entry : model.indexMap ){
    auto const& modelVar = entry.first;
    double const upperBound = modelVar.getUpperBound();
    if( upperBound < infinity ){
      // create a new constraint
      std::string name { modelVar.getName() };
      name += " (variable upper bound)";
      Constraint c { *this, new ipo_function::detail::
	  PhaseIBoundedVariableFunctionAndDerivatives( upperBound, true ), name };
      // Add variables to the new constraint
      Array vars { *this, name };
      // First the variable
      vars.push_back( crossIndex.find( modelVar )->second );
      // Then s
      vars.push_back( s );
      // And assign
      c.getVariables() = vars;
      // Set upper bound to zero
      c.setUpperBound( 0.0 );
      // Add the constraint to the model
      constraints.push_back( c );
    }
    double const lowerBound = modelVar.getLowerBound();
    if( lowerBound > minusInfinity ){
      // create a new constraint
      std::string name { modelVar.getName() };
      name += " (variable lower bound)";
      Constraint c { *this, new ipo_function::detail::
	  PhaseIBoundedVariableFunctionAndDerivatives { lowerBound, false }, name };
      // Add variables to the new constraint
      Array vars { *this, name };
      // First the variable
      vars.push_back( crossIndex.find( modelVar )->second );
      // Then s
      vars.push_back( s );
      // And assign
      c.getVariables() = vars;
      // Set upper bound to zero
      c.setUpperBound( 0.0 );
      // Add the constraint to the model
      constraints.push_back( c );
    }
  }
  // Next line need because objective/constraint.getVariables()
  // sets model.variablesChanged to true even though we have not
  // changed them.
  model.variablesChanged = false;

  // 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;
}

bool
PhaseIModel::findFeasibleSolution(){
  // Make sure phase I model is up to date.
  setIndices(); // This also sets the values of variables here

  if( not findEqualityConstraintFeasibleSolution() ) return false; // not even Phase I feasible
  if( isStrictlyFeasible() ) return true; // not a good idea to try to optimise.
  
  // tell phase I model to stop in Newton descent when a feasible solution is found
  getNewtonDescentParameters().setPhaseI( true );

  // get initial variable values.
  size_t const SIZE { indexMap.size() };
  gsl::vector vector( SIZE );
  for( auto& entry : indexMap ){
    auto& variable = entry.first;
    auto index = entry.second;
    vector[index] = variable.getValue();
  }
  // s would be set to zero or last-known value: correct this
  auto subvector = vector.subvector( 0, SIZE - 1 );
  vector[SIZE - 1] = getInitial_s( subvector );

  // set variables
  setVariablesFromVector( vector );
  
  // minimise and return
  if( vector[SIZE - 1] < 0 or minimise() ){ // minimise sets the values of variables in this
    // copy from variables of this to variables of model.
    for( auto& entry : model.indexMap ){
      auto& variable = entry.first;
      // copy variable values directly
      const_cast<Variable&>( variable )
    .setValue( crossIndex.find( variable )->second.getValue() );
    }
    return true;
  }
  // Just in case
  return false;
}

bool
PhaseIModel::findEqualityConstraintFeasibleSolution(){
  double constexpr ROOTEPS { std::sqrt( std::numeric_limits<double>::epsilon() ) };
  size_t const ROWS { equalityConstraints.size() };
  if( ROWS == 0 ) return true; // nothing to do
  setIndices(); // should return immediately: but allow possibility that indices not set
  size_t const COLS { model.indexMap.size() };
  if( COLS == 0 ) return true; // nothing to do

  // Create matrix A, vectors x and b from equality constraints
  gsl::matrix A { std::max( ROWS, COLS ), COLS };
  gsl::matrix V { COLS, COLS };
  A.set_all( 0 );
  gsl::vector x( COLS );
  gsl::vector S( COLS );
  gsl::vector b( std::max( ROWS, COLS ) );
  b.set_all( 0 );
  
  for( size_t i { 0 }; i < ROWS; ++i ){
    auto equalityConstraint = equalityConstraints[i];
    auto vars = equalityConstraint.getVariables();
    // set coefficients
    for( auto& var : vars ){
      size_t const index { indexMap[var] };
      if( index < COLS ) // ignore s
    A.set( i, indexMap[var], equalityConstraint.getCoefficient( var ) );
    }
    // set values
    b[i] = equalityConstraint.getUpperBound(); // == getLowerBound()
  }
  auto handler = gsl::exception::set_handler_gsl_exceptions();
  try {
    // Solve using SVD
    gsl::linalg::SV_decomp( A, V, S, x );
    for( auto& entry : S )
      if( std::fabs( entry ) < ROOTEPS ) entry = 0;
    gsl::linalg::SV_solve( A, V, S, b, x );
  } catch( gsl::exception e ){
    return false; // Singular value decomposition failed.
  } catch( ... ) {
    gsl::exception::set_handler( handler );
    throw IPOE( "ipo::detail::PhaseIModel::findEqualityConstraintFeasibleSolution(): "
        "unknown exception encountered." );
  }
  gsl::exception::set_handler( handler );
  
  // Copy values back to variables from x
  for( auto& p : model.indexMap ){
    auto var = p.first;
    auto q = crossIndex.find( var );
    if( crossIndex.end() == q )
      // This should not happen because PhaseIModel variables and Model variables should
      // not be modified independently
      throw IPOE( "ipo::detail::PhaseIModel::findEqualityConstraintFeasibleSolution(): "
          "variable in phase I model fails to match variable in model." );
    Variable& pIvar = const_cast<Variable&>( q->second );
    auto r = indexMap.find( pIvar );
    if( indexMap.end() == r )
      // This should not happen because PhaseIModel variables and Model variables should
      // not be modified independently
      throw IPOE( "ipo::detail::PhaseIModel::findEqualityConstraintFeasibleSolution(): "
          "variable in phase I model not contained in its index map." );
    size_t const index { r->second };
    if( index < COLS ){ // ignore s
      double const d { x[index] };
      // set value of variable: notice that s is not set
      pIvar.setValue( d );
      // also set in model
      var.setValue( d );
    }
  }
  
  return true;
}