NLPInterfacePack_NLPSerialPreprocess.hpp

// @HEADER
// ***********************************************************************
// 
// Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
//                  Copyright (2003) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library 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
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
// 
// ***********************************************************************
// @HEADER

#ifndef NLP_FULL_TO_REDUCED_H
#define NLP_FULL_TO_REDUCED_H

#include <valarray>

#include "NLPInterfacePack_NLPFirstOrder.hpp"
#include "NLPInterfacePack_NLPVarReductPerm.hpp"
#include "AbstractLinAlgPack_SpVectorClass.hpp"
#include "AbstractLinAlgPack_VectorMutableDense.hpp"
#include "AbstractLinAlgPack_PermutationSerial.hpp"
#include "DenseLinAlgPack_DVectorClass.hpp"
#include "DenseLinAlgPack_IVector.hpp"

namespace NLPInterfacePack {

class NLPSerialPreprocess
  : virtual public NLPObjGrad
  , virtual public NLPVarReductPerm
{
public:


  class InconsistantBounds : public std::logic_error
  {public: InconsistantBounds(const std::string& what_arg) : std::logic_error(what_arg) {}};


  NLPSerialPreprocess();

  static value_type fixed_var_mult(); 
  

  void force_xinit_in_bounds(bool force_xinit_in_bounds);
  bool force_xinit_in_bounds() const;
  void initialize(bool test_setup); 
  bool is_initialized() const;
  size_type n() const;
  size_type m() const;
  vec_space_ptr_t space_x() const;
  vec_space_ptr_t space_c() const;
  size_type num_bounded_x() const;
  const Vector& xl() const;
  const Vector& xu() const;
  const Vector& xinit() const;
  void get_init_lagrange_mult(
    VectorMutable*   lambda
    ,VectorMutable*  nu
    ) const;
  void scale_f( value_type scale_f );
  value_type scale_f() const;
  void report_final_solution(
    const Vector&    x
    ,const Vector*   lambda
    ,const Vector*   nu
    ,bool            is_optimal
    );
  virtual size_type ns() const;
  vec_space_ptr_t space_c_breve() const;
  vec_space_ptr_t space_h_breve() const;
  const Vector& hl_breve() const;
  const Vector& hu_breve() const;
  const Permutation& P_var() const;
  const Permutation& P_equ() const;



  const perm_fcty_ptr_t factory_P_var() const;
  const perm_fcty_ptr_t factory_P_equ() const;
  Range1D var_dep() const;
  Range1D var_indep() const;
  Range1D equ_decomp() const;
  Range1D equ_undecomp() const;
    bool nlp_selects_basis() const;
  bool get_next_basis(
    Permutation*  P_var,   Range1D* var_dep
    ,Permutation* P_equ,   Range1D* equ_decomp
    );
  void set_basis(
    const Permutation   &P_var,   const Range1D  &var_dep
    ,const Permutation  *P_equ,   const Range1D  *equ_decomp
    );
  void get_basis(
    Permutation*  P_var,   Range1D* var_dep
    ,Permutation* P_equ,   Range1D* equ_decomp
    ) const;


protected:


  void imp_calc_f(
    const Vector            &x
    ,bool                   newx
    ,const ZeroOrderInfo    &zero_order_info
    ) const;
  void imp_calc_c(
    const Vector            &x
    ,bool                   newx
    ,const ZeroOrderInfo    &zero_order_info
    ) const;
  void imp_calc_c_breve(
    const Vector            &x
    ,bool                   newx
    ,const ZeroOrderInfo    &zero_order_info_breve
    ) const;
  void imp_calc_h_breve(
    const Vector            &x
    ,bool                   newx
    ,const ZeroOrderInfo    &zero_order_info_breve
    ) const;



  void imp_calc_Gf(
    const Vector            &x
    ,bool                   newx
    ,const ObjGradInfo      &obj_grad_info
    ) const;



  struct ZeroOrderInfoSerial {
  public:
        ZeroOrderInfoSerial() : f(NULL)
    {}
    ZeroOrderInfoSerial( value_type* f_in, DVector* c_in, DVector* h_in )
      : f(f_in), c(c_in), h(h_in)
    {}
    value_type*    f;
    DVector*        c;
    DVector*        h;
  }; // end struct ZeroOrderInfoSerial

  struct ObjGradInfoSerial {
  public:
    ObjGradInfoSerial() : f(NULL)
    {}
    ObjGradInfoSerial( DVector* Gf_in, const ZeroOrderInfoSerial& first_order_info_in )
      : Gf(Gf_in), f(first_order_info_in.f), c(first_order_info_in.c), h(first_order_info_in.h)
    {}
    DVector*        Gf;
    value_type*    f;
    DVector*        c;
    DVector*        h;
  }; // end struct ObjGradInfoSerial



  virtual bool imp_nlp_has_changed() const { return true; }
  virtual size_type imp_n_orig() const = 0;
  virtual size_type imp_m_orig() const = 0;
  virtual size_type imp_mI_orig() const = 0;
  virtual const DVectorSlice imp_xinit_orig() const = 0;
  virtual bool imp_has_var_bounds() const = 0;
  virtual const DVectorSlice imp_xl_orig() const = 0;
  virtual const DVectorSlice imp_xu_orig() const = 0;
  virtual const DVectorSlice imp_hl_orig() const = 0;
  virtual const DVectorSlice imp_hu_orig() const = 0;
  virtual void imp_calc_f_orig(
    const DVectorSlice           &x_full
    ,bool                        newx
    ,const ZeroOrderInfoSerial   &zero_order_info
    ) const = 0;
  virtual void imp_calc_c_orig(
    const DVectorSlice           &x_full
    ,bool                        newx
    ,const ZeroOrderInfoSerial   &zero_order_info
    ) const = 0;
  virtual void imp_calc_h_orig(
    const DVectorSlice           &x_full
    ,bool                        newx
    ,const ZeroOrderInfoSerial   &zero_order_info
    ) const = 0;
  virtual void imp_calc_Gf_orig(
    const DVectorSlice           &x_full
    ,bool                        newx
    ,const ObjGradInfoSerial     &obj_grad_info
    ) const = 0;
  virtual bool imp_get_next_basis(
    IVector      *var_perm_full
    ,IVector     *equ_perm_full
    ,size_type   *rank_full
    ,size_type   *rank
    );
  virtual void imp_report_orig_final_solution(
    const DVectorSlice      &x_full
    ,const DVectorSlice     *lambda_orig
    ,const DVectorSlice     *lambdaI_orig
    ,const DVectorSlice     *nu_orig
    ,bool                   optimal
    )
  {}



  void set_not_initialized();

  void assert_initialized() const;

  void set_x_full(const DVectorSlice& x, bool newx, DVectorSlice* x_full) const;

  DVectorSlice x_full() const;

  const ZeroOrderInfoSerial zero_order_orig_info() const;

  const ObjGradInfoSerial obj_grad_orig_info() const;
  
  const IVector& var_remove_fixed_to_full() const;

  const IVector& var_full_to_remove_fixed() const;

  const IVector& var_perm() const;

  const IVector& equ_perm() const;

  const IVector& inv_equ_perm() const;

  // Perform the mapping from a full variable vector to the reduced permuted variable vector
  void var_from_full( DVectorSlice::const_iterator vec_full, DVectorSlice::iterator vec ) const;

  // Perform the mapping from a reduced permuted variable vector the full variable vector
  void var_to_full(DVectorSlice::const_iterator vec, DVectorSlice::iterator vec_full) const;

  // Perform the mapping from c_orig, h_orig, s_orig to the permuted constraint vector c
  void equ_from_full(
    const DVectorSlice   &c_orig
    ,const DVectorSlice  &h_orig
    ,const DVectorSlice  &s_orig
    ,DVectorSlice        *c_full
    ) const;


private:

  // ///////////////////////////
  // Private data members

  mutable value_type          f_orig_;    // Computed by subclasses
  mutable DVector           c_orig_;    // ...
  mutable DVector           h_orig_;    // ...
  mutable DVector           Gf_full_;   // ...

  bool initialized_;
  // Flag for if the NLP has has been properly initialized

  bool force_xinit_in_bounds_;
  // Determine if the initial point will be adjusted between bounds

  value_type scale_f_;
  // Set the scaling of the objective function used.

  IVector   var_full_to_fixed_;
  // Holds the indices of variables that are fixed by bounds and those
  // that are not (Length = n_full_).  These partitions are not
  // necessarly sorted in assending order as var_perm and con_perm are.
  //
  //  var_full_to_fixed_ =  [ not fixed by bounds | fixed by bounds  ]
  //              [1 ..         n_|n_ + 1 ..    n_full_]
  //

  IVector   inv_var_full_to_fixed_;
  // Inverse of var_full_to_fixed_.  If inv_var_full_to_fixed_(i) > n_ then this variable
  // is fixed between bounds, else inv_var_full_to_fixed_(i) is the indice of the 
  // variable in the unsorted x (not permuted to the current basis).

  IVector   var_perm_;
  // Variable permutations (length = n_) from the vector of unstorted variables not fixed
  // by bounds as defined by var_full_to_fixed_
  //
  // var_perm_  = [ dependent variables | independent variables ]
  //          [1..                r_|r_+1...              n_]
  //

  IVector   equ_perm_;
  // Equality Constraint permutations (length = m_)
  //
  // equ_perm_  = [ decomposed equalities | undecomposed equalities ]
  //          [1..                  r_|r_+1...                m_]
  //

  IVector   inv_equ_perm_;
  // Inverse of equ_perm
  //

  mutable DVector x_full_;
  DVector         xinit_full_;
  DVector         xl_full_;
  DVector         xu_full_;
  // The full vector (length = n_full_).  This vector may include
  // slack variables if mI_orig > 0:
  //
  //    [ x_orig; s_orig ]
  //

  perm_fcty_ptr_t            factory_P_var_;
  perm_fcty_ptr_t            factory_P_equ_;
  VectorSpaceSerial          space_x_;
  VectorSpaceSerial          space_c_;
  VectorSpaceSerial          space_c_breve_;
  VectorSpaceSerial          space_h_breve_;
  size_type                  num_bounded_x_;
  VectorMutableDense         xinit_; // Initial point of the shrunken NLP
  VectorMutableDense         xl_;    // Lower bounds of transformed NLP
  VectorMutableDense         xu_;    // Uppers bounds of transformed NLP
  VectorMutableDense         hl_breve_;// Lower bounds for general inequalities of transformed NLP
  VectorMutableDense         hu_breve_;// Uppers bounds for general inequalitiess of transformed NLP
  PermutationSerial          P_var_;
  PermutationSerial          P_equ_;
  size_type              n_orig_;  // Number of variables in the original NLP
  size_type              m_orig_;  // Number of general equality constraints in the original NLP
  size_type              mI_orig_; // Number of general inequality constraints in the original NLP
  size_type              n_full_;  // Number of variables in the transformed NLP (before fixed variables are removed)
  size_type              m_full_;  // Number of general equality constraints in the transformed NLP
  size_type              n_;       // Number of variables in the transformed NLP (with slacks and not fixed by bounds)
  size_type              r_;       // Number of independent equations in the transformed NLP

  int                    basis_selection_num_; // Number of the basis to select next

  // ///////////////////////////
  // Private member functions

  // Get the next basis (or first basis) from the NLP subclass and remove the
  // fixed variables.  Note that this function does not modify var_perm_, equ_perm_
  // or r_.  You must do that yourself by calling assert_and_set_basis.
  bool get_next_basis_remove_fixed(
    IVector* var_perm, IVector* equ_perm, size_type* rank );
  
  // Assert (throw std::length_error, NLPVarReductPerm::InvalidBasis) and set a basis selection
  // If &var_perm == &var_perm_ and/or &equ_perm == &equ_perm_ then the unneeded copy
  // is avoided.
  void assert_and_set_basis(
    const IVector& var_perm, const IVector& equ_perm, size_type rank );

  // Assert that there are bounds on the variables (throw NLP::NoBoundsOnVariables)
  void assert_bounds_on_variables() const;

  // Adjust initial point this->xinit_ to be within bound
  void do_force_xinit_in_bounds();

};  // end class NLPSerialPreprocess

// //////////////////////////////////////////////////
// Inline member functions

// protected

inline
void NLPSerialPreprocess::set_not_initialized()
{
  initialized_ = false;
}

inline
DVectorSlice NLPSerialPreprocess::x_full() const
{
  return x_full_();
}

inline
const NLPSerialPreprocess::ZeroOrderInfoSerial
NLPSerialPreprocess::zero_order_orig_info() const
{
  return ZeroOrderInfoSerial( &f_orig_, &c_orig_, &h_orig_ );
}

inline
const NLPSerialPreprocess::ObjGradInfoSerial
NLPSerialPreprocess::obj_grad_orig_info() const
{
  return ObjGradInfoSerial( &Gf_full_, zero_order_orig_info() );
}

inline
const IVector& NLPSerialPreprocess::var_remove_fixed_to_full() const
{
  return var_full_to_fixed_;
}

inline
const IVector& NLPSerialPreprocess::var_full_to_remove_fixed() const
{
  return inv_var_full_to_fixed_;
}

inline
const IVector& NLPSerialPreprocess::var_perm() const
{
  return var_perm_;
}

inline
const IVector& NLPSerialPreprocess::equ_perm() const
{
  return equ_perm_;
}

inline
const IVector& NLPSerialPreprocess::inv_equ_perm() const
{
  return inv_equ_perm_;
}

} // end namespace NLPInterfacePack 

#endif // NLP_FULL_TO_REDUCED_H