BarrierFunction.cc

Go to the documentation of this file.

/*
 * $Id: BarrierFunction.cc 222 2014-11-09 14:04:41Z jdl3 $
 * Copyright (C) 2013 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"BarrierFunction.hpp"

using namespace ipo::detail;

BarrierFunction::BarrierFunction( Model& model, Sign const sign, double const t )
  : model( model ), sign { sign }, t{ t }
{
  size_t const SIZE { model.indexMap.size() };
  functionGradient = gsl::vector( SIZE );
  functionHessian = gsl::matrix( SIZE, SIZE );
}

void
BarrierFunction::set_t( double const t ){
  this->t = t;
}

double
BarrierFunction::get_t() const {
  return t;
}

double
BarrierFunction::operator()( gsl::vector const& vector ){
  // First, the objective function
  gsl::vector subvector { model.objectiveIndices->get( vector ) };
  // Update
  initialiseResultFromObjective( subvector );

  // Then the explicit constraints
  for( size_t j { 0 }; j < model.constraints.size(); ++j ){
    auto& constraintIndices_j = model.constraintIndices[j];
    subvector = constraintIndices_j->get( vector );
    auto& constraint = model.constraints[j];
    // get bounds
    double const upperBound { constraint.getUpperBound() };
    double const lowerBound { constraint.getLowerBound() };
    if( upperBound >= infinity and lowerBound <= minusInfinity ) continue;
    double const f { (*constraint.getFunction())( subvector ) };
    // Update
    updateResultFromConstraint( constraint, subvector, lowerBound, f, upperBound );
  }

  // Then any variable constraints
  for( auto& entry : model.indexMap )
    // Update
    updateResultFromVariable( entry, vector );

  // Finally, return result
  return functionValue;
}

gsl::vector
BarrierFunction::gradient( gsl::vector const& vector ){
  // First, the objective function
  gsl::vector subvector { model.objectiveIndices->get( vector ) };
  // Update
  initialiseGradientFromObjective( subvector );

  // Then the explicit constraints
  for( size_t j { 0 }; j < model.constraints.size(); ++j ){
    auto& constraintIndices_j = model.constraintIndices[j];
    subvector = constraintIndices_j->get( vector );
    auto& constraint = model.constraints[j];
    // get bounds
    double const upperBound { constraint.getUpperBound() };
    double const lowerBound { constraint.getLowerBound() };
    if( upperBound >= infinity and lowerBound <= minusInfinity ) continue;
    double const f { (*constraint.getFunction())( subvector ) };
    // Update
    updateGradientFromConstraint( j, subvector, lowerBound, f, upperBound );
  }

  // Then any variable constraints
  for( auto& entry : model.indexMap )
    // Update
    updateGradientFromVariable( entry, vector );

  // Finally, return result
  return functionGradient;
}

gsl::matrix
BarrierFunction::hessian( gsl::vector const& vector ){
  // First, the objective function
  gsl::vector subvector { model.objectiveIndices->get( vector ) };
  // Update
  initialiseHessianFromObjective( subvector );

  // Then the explicit constraints
  for( size_t j = 0; j < model.constraints.size(); ++j ){
    auto& constraintIndices_j = model.constraintIndices[j];
    subvector = constraintIndices_j->get( vector );
    auto& constraint = model.constraints[j];
    // get bounds
    double const upperBound { constraint.getUpperBound() };
    double const lowerBound { constraint.getLowerBound() };
    if( upperBound >= infinity and lowerBound <= minusInfinity ) continue;
    double const f { (*constraint.getFunction())( subvector ) };
    // Update
    gsl::vector gradient { constraint.getFunction()->gradient( subvector ) };
    updateHessianFromConstraint( j, subvector, lowerBound, f, upperBound, gradient );
  }

  // Then any variable constraints
  for( auto& entry : model.indexMap )
    // Update
    updateHessianFromVariable( entry, vector );
  
  // Finally, return result
  return functionHessian;
}

void
BarrierFunction::initialiseResultFromObjective( gsl::vector& subvector, bool const simul ){
  functionValue = 0;
  DerivativesEstimates* de
    = dynamic_cast<DerivativesEstimates*>( model.objective.getFunction() );
  if( simul and nullptr != de ){
    de->setVector( subvector );
    if( Sign::positive == sign )
      functionValue += t * de->value();
    else
      functionValue -= t * de->value();
  } else {
    if( Sign::positive == sign )
      functionValue += t * (*model.objective.getFunction())( subvector );
    else
      functionValue -= t * (*model.objective.getFunction())( subvector );
  }
}

void
BarrierFunction::updateResultFromConstraint( Constraint& constraint, gsl::vector& subvector,
                         double const& lowerBound,
                         double const& functionValue,
                         double const& upperBound ){
  // Add constraints for upper and lower bounds.
  if( upperBound < infinity )
    this->functionValue -= std::log( upperBound - functionValue );
  if( lowerBound > minusInfinity )
    this->functionValue -= std::log( functionValue - lowerBound );
}

void
BarrierFunction::updateResultFromVariable( std::map<Variable,size_t>::value_type& entry,
                       gsl::vector const& vector ){
  Variable const& variable { entry.first };
  // get bounds
  double const upperBound { variable.getUpperBound() };
  double const lowerBound { variable.getLowerBound() };
  if( upperBound >= infinity and lowerBound <= minusInfinity ) return;
  double const f = vector[entry.second];
  // Add results for upper and lower bounds.
  if( upperBound < infinity )
    functionValue -= std::log( upperBound - f );
  if( lowerBound > minusInfinity )
    functionValue -= std::log( f - lowerBound );
}

gsl::vector
BarrierFunction::initialiseGradientFromObjective( gsl::vector& subvector, bool const simul ){
  DerivativesEstimates*
    de = dynamic_cast<DerivativesEstimates*>( model.objective.getFunction() );
  gsl::vector gradient { simul and nullptr != de ? de->gradient()
      : model.objective.getFunction()->gradient( subvector ) };
  // In any case initialise functionGradient
  if( functionGradient.size() != model.indexMap.size() )
    functionGradient = gsl::vector( model.indexMap.size() );
  functionGradient.set_all( 0 );
  // and copy to functionGradient;
  if( Sign::positive == sign ){
    model.objectiveIndices->set( functionGradient, gradient, t );
  } else {
    model.objectiveIndices->set( functionGradient, gradient, -t );
  }
  // Return objective gradient
  return gradient;
}

gsl::vector
BarrierFunction::updateGradientFromConstraint( size_t const constraintIndex,
                           gsl::vector& subvector,
                           double const& lowerBound,
                           double const& functionValue,
                           double const& upperBound, bool const simul ){
  auto& constraint = model.constraints[constraintIndex];
  DerivativesEstimates*
    de = dynamic_cast<DerivativesEstimates*>( constraint.getFunction() );
  gsl::vector gradient { simul and nullptr != de ? de->gradient() 
      : constraint.getFunction()->gradient( subvector ) };
  // Now gradient is set correctly
  if( upperBound < infinity ){
    double const scale_inv { upperBound - functionValue };
    double const scale { 1.0 / scale_inv };
    model.constraintIndices[constraintIndex]->add( functionGradient, gradient, scale );
  }
  if( lowerBound > minusInfinity ){
    double const scale_inv { lowerBound - functionValue };
    double const scale { 1.0 / scale_inv };
    model.constraintIndices[constraintIndex]->add( functionGradient, gradient, scale );
  }
  return gradient;
}

void
BarrierFunction::updateGradientFromVariable( std::map<Variable,size_t>::value_type& entry,
                         gsl::vector const& vector ){
  Variable const& variable { entry.first };
  // get bounds
  double const upperBound { variable.getUpperBound() };
  double const lowerBound { variable.getLowerBound() };
  if( upperBound >= infinity and lowerBound <= minusInfinity ) return;
  double const f { vector[entry.second] };
  // Add results for upper and lower bounds.
  if( upperBound < infinity ){
    double const scale_inv { upperBound - f };
    double const scale { 1.0 / scale_inv };
    functionGradient[entry.second] += scale;
  }
  if( lowerBound > minusInfinity ){
    double const scale_inv { lowerBound - f };
    double const scale { 1.0 / scale_inv };
    functionGradient[entry.second] += scale;
  }
}

void
BarrierFunction::initialiseHessianFromObjective( gsl::vector& subvector, bool const simul ){
  size_t const SIZE { model.indexMap.size() };
  // In any case initialise functionHessian
  if( functionHessian.size1() != SIZE or functionHessian.size2() != SIZE )
    functionHessian = gsl::matrix( SIZE, SIZE );
  functionHessian.set_all( 0 );
  // set hessian of objective
  DerivativesEstimates*
    de = dynamic_cast<DerivativesEstimates*>( model.objective.getFunction() );
  gsl::matrix hessian { simul and nullptr != de ? de->hessian()
      : model.objective.getFunction()->hessian( subvector ) };
  // Copy to functionHessian;
  if( Sign::positive == sign ){
    model.objectiveIndices->set( functionHessian, hessian, t );
  } else {
    model.objectiveIndices->set( functionHessian, hessian, -t );
  }
}

void
BarrierFunction::updateHessianFromConstraint( size_t const constraintIndex,
                          gsl::vector& subvector,
                          double const& lowerBound,
                          double const& functionValue,
                          double const& upperBound,
                          gsl::vector& gradient, bool const simul ){
  auto& constraint = model.constraints[constraintIndex];
  DerivativesEstimates*
    de = dynamic_cast<DerivativesEstimates*>( constraint.getFunction() );
  gsl::matrix hessian { simul and nullptr != de ? de->hessian()
      : constraint.getFunction()->hessian( subvector ) };
  // Now update from hessian
  if( upperBound < infinity ){
    double const scale_inv { upperBound - functionValue };
    double const scale { 1.0 / scale_inv };
    double const scale2 { scale * scale };
    model.constraintIndices[constraintIndex]->add( functionHessian, gradient, scale2 );
    model.constraintIndices[constraintIndex]->add( functionHessian, hessian, scale );
  }
  if( lowerBound > minusInfinity ){
    double const scale_inv { functionValue - lowerBound };
    double const scale { 1.0 / scale_inv };
    double const scale2 { scale * scale };
    model.constraintIndices[constraintIndex]->add( functionHessian, gradient, scale2 );
    model.constraintIndices[constraintIndex]->add( functionHessian, hessian, scale );
  }
}

void
BarrierFunction::updateHessianFromVariable( std::map<Variable,size_t>::value_type& entry,
                        gsl::vector const& vector ){
  Variable const& variable { entry.first };
  // get bounds
  double const upperBound { variable.getUpperBound() };
  double const lowerBound { variable.getLowerBound() };
  if( upperBound >= infinity and lowerBound <= minusInfinity ) return;
  double const f { vector[entry.second] };
  // Add results for upper and lower bounds.
  if( upperBound < infinity ){
    double const scale_inv { f - upperBound };
    double const scale { 1.0 / scale_inv };
    double const scale2 { scale * scale };
    double const d { functionHessian.get( entry.second, entry.second ) + scale2 };
    functionHessian.set( entry.second, entry.second, d );
  }
  if( lowerBound > minusInfinity ){
    double const scale_inv { lowerBound - f };
    double const scale { 1.0 / scale_inv };
    double const scale2 { scale * scale };
    double const d { functionHessian.get( entry.second, entry.second ) + scale2 };
    functionHessian.set( entry.second, entry.second, d );
  }
}

void
BarrierFunction::setVector( gsl::vector const& vector ){
  // vector is not actually stored anywhere
  
  // First, the objective function
  auto subvector = model.objectiveIndices->get( vector );
  // Update
  initialiseResultFromObjective( subvector, true );
  initialiseGradientFromObjective( subvector, true );
  initialiseHessianFromObjective( subvector, true );

  // Then the explicit constraints
  for( size_t j { 0 }; j < model.constraints.size(); ++j ){
    subvector = model.constraintIndices[j]->get( vector );
    auto& constraint = model.constraints[j];
    // get bounds
    double const upperBound { constraint.getUpperBound() };
    double const lowerBound { constraint.getLowerBound() };
    if( upperBound >= infinity and lowerBound <= minusInfinity ) continue;
    // Try cross-casting
    DerivativesEstimates* de = dynamic_cast<DerivativesEstimates*>( constraint.getFunction() );
    // Make sure constraint value and derivatives are calculated
    if( nullptr != de ) de->setVector( subvector );
    double const f {  nullptr == de ? (*constraint.getFunction())( subvector ) : de->value() };
    // Update
    updateResultFromConstraint( constraint, subvector, lowerBound, f, upperBound );
    auto gradient
      = updateGradientFromConstraint( j, subvector, lowerBound, f, upperBound, true );
    updateHessianFromConstraint( j, subvector, lowerBound, f, upperBound, gradient, true);
  }

  // Then any variable constraints
  for( auto& entry : model.indexMap ){
    // Update
    updateResultFromVariable( entry, vector );
    updateGradientFromVariable( entry, vector );
    updateHessianFromVariable( entry, vector );
  }
}