Model.hpp

Go to the documentation of this file.

/*
 * $Id: Model.hpp 225 2014-11-17 17:52:40Z 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.
 */

#ifndef IPO_MODEL_HPP
#define IPO_MODEL_HPP

#include<map>
#include"Objective.hpp"
#include"Constraint.hpp"
#include"../ipo_function/DerivativesEstimates.hpp"
#include"detail/NewtonDescent.hpp"
#include"../ipo_function/concrete/Null.hpp"
#include"LinearConstraint.hpp"

namespace ipo {
  // Forward declarations
  class LinearConstraint;
  namespace detail {
    class BarrierFunction;
    class PhaseIModel;
  }
  class Model : public detail::ModelBase {
    // Classes used internally to get subvectors of a vector
  public:
    class Parameters {
    public:
      Parameters();
      double getMu() const { return mu; }
      double getEpsilon() const { return epsilon; }
      std::ostream* getOutputStream() const { return outputStream; }
      void setMu( double const mu );
      void setEpsilon( double const epsilon );
      void setOutputStream( std::ostream* outputStream );
    private:
      double mu;
      double epsilon;
      std::ostream* outputStream;
    };
  private:
    Parameters parameters;
#ifdef TESTSUBVECTOR
  public:
#endif

    class Subvector {
    public:
      Subvector(){}
      virtual gsl::vector const get( gsl::vector const& vector ) const = 0;
      virtual void set( gsl::vector& vector, gsl::vector const& subvector,
            double const scale = 1.0 ) const = 0;
      virtual void add( gsl::vector& vector, gsl::vector const& subvector,
            double const scale = 1.0 ) const = 0;
      virtual void set( gsl::matrix& matrix, gsl::matrix const& submatrix,
            double const scale = 1.0 ) const = 0;
      virtual void add( gsl::matrix& matrix, gsl::matrix const& submatrix,
            double const scale = 1.0 ) const = 0;
      virtual void set( gsl::matrix& matrix, gsl::vector const& vector,
            double const scale = 1.0 ) const = 0;
      virtual void add( gsl::matrix& matrix, gsl::vector const& vector,
            double const scale = 1.0 ) const = 0;
    };
    class DirectSubvector : public Subvector {
    public:
      DirectSubvector( size_t const offset, size_t const length );
      virtual gsl::vector const get( gsl::vector const& vector ) const {
    return vector.subvector( offset, length ); }
      virtual void set( gsl::vector& vector, gsl::vector const& subvector,
            double const scale = 1.0 ) const {
    for( size_t i { 0 }; i < length; ++i ) vector[i + offset] = scale * subvector[i]; }
      virtual void add( gsl::vector& vector, gsl::vector const& subvector, 
            double const scale = 1.0 ) const {
    for( size_t i { 0 }; i < length; ++i ) vector[i + offset] += scale * subvector[i]; }
      virtual void set( gsl::matrix& matrix, gsl::matrix const& submatrix,
            double const scale = 1.0 ) const {
    for( size_t i { 0 }; i < length; ++i )
      for( size_t j { i }; j < length; ++j ){
        double const d { scale * submatrix.get( i, j ) };
        matrix.set( i + offset, j + offset, d );
        if( j > i ) matrix.set( j + offset, i + offset, d );
      }
      }
      virtual void add( gsl::matrix& matrix, gsl::matrix const& submatrix,
            double const scale = 1 ) const {
    for( size_t i = 0; i < length; ++i )
      for( size_t j = i; j < length; ++j ){
        double const d = scale * submatrix.get( i, j )
          + matrix.get( i + offset, j + offset );
        matrix.set( i + offset, j + offset, d );
        if( j > i ) matrix.set( j + offset, i + offset, d );
      }
      }
      virtual void set( gsl::matrix& matrix, gsl::vector const& vector,
            double const scale = 1 ) const {
    for( size_t i = 0; i < length; ++i )
      for( size_t j = i; j < length; ++j ){
        double const d = scale * vector[i] * vector[j];
        matrix.set( i + offset, j + offset, d );
        if( j > i ) matrix.set( j + offset, i + offset, d );
      }
      }
      virtual void add( gsl::matrix& matrix, gsl::vector const& vector,
            double const scale = 1.0 ) const {
    for( size_t i { 0 }; i < length; ++i )
      for( size_t j { i }; j < length; ++j ){
        double const d { scale * vector[i] * vector[j]
        + matrix.get( i + offset, j + offset ) };
        matrix.set( i + offset, j + offset, d );
        if( j > i ) matrix.set( j + offset, i + offset, d );
      }
      }
    private:
      size_t const offset;
      size_t const length;
    };
    class IndirectSubvector : public Subvector, public std::vector<size_t> {
    public:
      IndirectSubvector( size_t const size );
      virtual gsl::vector const get( gsl::vector const& vector ) const {
    gsl::vector subvector( size() );
    for( size_t index { 0 }; index < size(); ++index )
      subvector[index] = vector[(*this)[index]];
    return subvector;
      }
      virtual void set( gsl::vector& vector, gsl::vector const& subvector,
            double const scale = 1.0 ) const {
    for( size_t index { 0 }; index < size(); ++index )
      vector[(*this)[index]] = scale * subvector[index]; }
      virtual void add( gsl::vector& vector, gsl::vector const& subvector,
            double const scale = 1.0 ) const {
    for( size_t index { 0 }; index < size(); ++index )
      vector[(*this)[index]] += scale * subvector[index]; }
      virtual void set( gsl::matrix& matrix, gsl::matrix const& submatrix,
            double const scale = 1.0 ) const {
    for( size_t i { 0 }; i < size(); ++i )
      for( size_t j = i; j < size(); ++j ){
        size_t const r { (*this)[i] };
        size_t const c { (*this)[j] };
        double const d { scale * submatrix.get( i , j ) };
        matrix.set( r, c, d );
        if( j > i ) matrix.set( c, r, d );
      }
      }
      virtual void add( gsl::matrix& matrix, gsl::matrix const& submatrix,
            double const scale = 1 ) const {
    for( size_t i = 0; i < size(); ++i )
      for( size_t j = i; j < size(); ++j ){
        size_t const r = (*this)[i];
        size_t const c = (*this)[j];
        double const d = scale * submatrix.get( i , j ) + matrix.get( r, c );
        matrix.set( r, c, d );
        if( j > i ) matrix.set( c, r, d );
      }
      }
      virtual void set( gsl::matrix& matrix, gsl::vector const& vector,
            double const scale = 1.0 ) const {
    for( size_t i { 0 }; i < size(); ++i )
      for( size_t j { i }; j < size(); ++j ){
        size_t const r { (*this)[i] };
        size_t const c { (*this)[j] };
        double const d { scale * vector[i] * vector[j] };
        matrix.set( r, c, d );
        if( j > i ) matrix.set( c, r, d );
      }
      }
      virtual void add( gsl::matrix& matrix, gsl::vector const& vector,
            double const scale = 1.0 ) const {
    for( size_t i { 0 }; i < size(); ++i )
      for( size_t j = i; j < size(); ++j ){
        size_t const r { (*this)[i] };
        size_t const c { (*this)[j] };
        double const d { scale * vector[i] * vector[j] + matrix.get( r, c ) };
        matrix.set( r, c, d );
        if( j > i ) matrix.set( c, r, d );
      }
      }
    };
    friend class detail::BarrierFunction;
    friend class detail::PhaseIModel;
  public:
    Model();
    virtual void notify();
    void setObjective( Objective& objective );
    Objective getObjective();
    void addConstraint( Constraint& constraint );
    bool removeConstraint( Constraint& constraint );
    std::vector<Constraint> const& getConstraints() const;
    std::vector<LinearConstraint> const& getEqualityConstraints() const;
    bool isFeasible();
    bool isStrictlyFeasible();
    bool minimise(){ return optimise( true ); }
    bool maximise(){ return optimise( false ); }
    virtual bool findFeasible();
    Parameters& getParameters();
    detail::NewtonDescent::Parameters& getNewtonDescentParameters();
    detail::LineSearch::Parameters& getLineSearchParameters();
    bool checkSize() const;
    void testSizes() const;
    void summary( std::ostream& ostream = std::cout ) const;
  protected:
    gsl::vector getVectorFromVariables() const;
    void setVariablesFromVector( gsl::vector const & vector );
    bool stoppingCriterion( double const t ) const;
    virtual bool optimise( bool const minimise );
    virtual void setIndices();
    inline void setObjectiveIndices(){
      size_t const OSIZE { objective.getVariables().size() }; /* identifies vars not in 
                                   * objective */
      bool const DIRECTSUBVECTORFOROBJECTIVEINDICES { OSIZE == indexMap.size() };
      // Set objectiveIndices
      if( DIRECTSUBVECTORFOROBJECTIVEINDICES )
    objectiveIndices = std::shared_ptr<Subvector>( new DirectSubvector( 0, OSIZE ) );
      else {
    IndirectSubvector* oIndices = new IndirectSubvector( OSIZE );
    objectiveIndices = std::shared_ptr<Subvector>( oIndices );
    size_t index { 0 };
    for( auto& variable : objective.getVariables() ){
      (*oIndices)[index] = indexMap[variable];
      ++index;
    }
      }
    }
     void setConstraintIndices();  // local function 
  protected:
    std::vector<Constraint> constraints;
    std::vector<ipo::LinearConstraint> equalityConstraints;
    Objective objective;
    std::shared_ptr<Subvector> objectiveIndices;
    std::vector<std::shared_ptr<Subvector> > constraintIndices;
    std::vector<std::shared_ptr<Subvector> > equalityConstraintIndices;
    std::map<Variable,size_t> indexMap;
    bool variablesChanged = false;
    detail::NewtonDescent::Parameters newtonDescentParameters;
    detail::LineSearch::Parameters lineSearchParameters;
  };
}

#endif