NLPInterfacePack_NLPSerialPreprocess.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <assert.h>
00030 
00031 #include <iostream> // Debug only
00032 #include "DenseLinAlgPack_PermOut.hpp"
00033 
00034 #include <algorithm>
00035 #include <sstream>
00036 #include <limits>
00037 #include <stdio.h>
00038 #include <fstream>
00039 
00040 #include "NLPInterfacePack_NLPSerialPreprocess.hpp"
00041 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00042 #include "AbstractLinAlgPack_PermutationSerial.hpp"
00043 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00044 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00045 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00046 #include "DenseLinAlgPack_DVectorOp.hpp"
00047 #include "DenseLinAlgPack_IVector.hpp"
00048 #include "DenseLinAlgPack_PermVecMat.hpp"
00049 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00050 #include "Teuchos_TestForException.hpp"
00051 #include "Teuchos_AbstractFactoryStd.hpp"
00052 #include "Teuchos_dyn_cast.hpp"
00053 
00054 namespace LinAlgOpPack {
00055   using AbstractLinAlgPack::Vp_StV;
00056 }
00057 
00058 namespace NLPInterfacePack {
00059 
00060 // NLPSerialPreprocess
00061 
00062 // Static public members
00063 
00064 value_type
00065 NLPSerialPreprocess::fixed_var_mult()
00066 {
00067   return std::numeric_limits<DenseLinAlgPack::DVector::value_type>::max()-100; // Don't know what to use?
00068 }
00069 
00070 // Constructors / nitializers
00071 
00072 NLPSerialPreprocess::NLPSerialPreprocess(
00073   )
00074   :initialized_(false)
00075   ,force_xinit_in_bounds_(true)
00076     ,scale_f_(1.0)
00077   ,basis_selection_num_(0)
00078 
00079 {}
00080 
00081 // Overridden public members from NLP
00082 
00083 void NLPSerialPreprocess::force_xinit_in_bounds(bool force_xinit_in_bounds)
00084 {
00085   force_xinit_in_bounds_ = force_xinit_in_bounds;
00086 }
00087 
00088 bool NLPSerialPreprocess::force_xinit_in_bounds() const
00089 {
00090   return force_xinit_in_bounds_;
00091 }
00092 
00093 void NLPSerialPreprocess::initialize(bool test_setup)
00094 {
00095   namespace mmp = MemMngPack;
00096 
00097   const value_type inf_bnd = NLP::infinite_bound();
00098 
00099   basis_selection_num_ = 0;
00100 
00101   if( initialized_  && !imp_nlp_has_changed() ) {
00102     // The subclass NLP has not changed so we can just
00103     // slip this preprocessing.
00104     NLPObjGrad::initialize(test_setup);
00105     return;
00106   }
00107 
00108   // Get the dimensions of the original problem
00109 
00110   n_orig_  = imp_n_orig();
00111   m_orig_  = imp_m_orig();   // This may be zero!
00112   mI_orig_ = imp_mI_orig();  // This may be zero!
00113   
00114   // Get the dimensions of the full problem
00115 
00116   n_full_  = n_orig_ + mI_orig_;
00117   m_full_  = m_orig_ + mI_orig_;
00118 
00119   // Initialize the storage for the intermediate quanities
00120   
00121   xinit_full_.resize(n_full_);
00122   xl_full_.resize(n_full_);
00123   xu_full_.resize(n_full_);
00124   x_full_.resize(n_full_);
00125   c_orig_.resize(m_orig_);
00126   h_orig_.resize(mI_orig_);
00127   Gf_full_.resize(n_full_);
00128   var_full_to_fixed_.resize(n_full_);
00129   equ_perm_.resize(m_full_);
00130   inv_equ_perm_.resize(m_full_);
00131   space_c_.initialize(m_full_);
00132   space_c_breve_.initialize(m_orig_);
00133   space_h_breve_.initialize(mI_orig_);
00134   factory_P_var_   = Teuchos::rcp( new Teuchos::AbstractFactoryStd<Permutation,PermutationSerial>() );
00135   factory_P_equ_   = Teuchos::rcp( new Teuchos::AbstractFactoryStd<Permutation,PermutationSerial>() );
00136 
00137   // Intialize xinit_full_, xl_full_ and xu_full_ for the initial point which will set the
00138   // fixed elements which will not change during the optimization.
00139   xinit_full_(1,n_orig_)  = imp_xinit_orig();
00140   xl_full_(1,n_orig_)     = imp_xl_orig();
00141   xu_full_(1,n_orig_)     = imp_xu_orig();
00142   if( n_full_ > n_orig_ ) { // Include slack varaibles
00143     xinit_full_(n_orig_+1,n_full_)  = 0.0;
00144     xl_full_(n_orig_+1,n_full_)     = imp_hl_orig();
00145     xu_full_(n_orig_+1,n_full_)     = imp_hu_orig();
00146   }
00147 
00148   const bool has_var_bounds = imp_has_var_bounds() || n_full_ > n_orig_;
00149 
00150   // Force the initial point in bounds if it is not.
00151   if( force_xinit_in_bounds() && has_var_bounds ) {
00152     AbstractLinAlgPack::force_in_bounds(
00153       VectorMutableDense( xl_full_(), Teuchos::null )
00154       ,VectorMutableDense( xu_full_(), Teuchos::null )
00155       ,&VectorMutableDense( x_full_(), Teuchos::null )
00156       );
00157   }
00158   
00159   // Determine which variables are fixed by bounds!
00160   size_type
00161     xl_nz     = 0,
00162     xu_nz     = 0,
00163     num_bnd_x = 0;
00164   if( has_var_bounds ) {
00165     // Determine which variables are fixed by bounds and
00166     // adjust the bounds if needed.
00167     DVector::iterator
00168       xl_full   = xl_full_.begin(),
00169       xu_full   = xu_full_.begin();
00170     n_ = 0;
00171     size_type num_fixed = 0;
00172     for(int i = 1; i <= n_full_; ++i, ++xl_full, ++xu_full) {
00173       TEST_FOR_EXCEPTION(
00174         *xl_full > *xu_full, InconsistantBounds
00175         ,"NLPSerialPreprocess::initialize() : Error, Inconsistant bounds: xl_full("
00176         << i << ") > xu_full(" << i << ")" ); 
00177       if(*xl_full == *xu_full) {
00178         //
00179         // Fixed between bounds
00180         //
00181         var_full_to_fixed_(n_full_ - num_fixed) = i;
00182         num_fixed++;
00183       }
00184       else {
00185         //
00186         // Not Fixed between bounds
00187         //
00188         // Adjust the bounds if needed
00189         *xl_full = *xl_full < -inf_bnd ? -inf_bnd : *xl_full;
00190         *xu_full = *xu_full > +inf_bnd ? +inf_bnd : *xu_full;
00191         //
00192         n_++;
00193         var_full_to_fixed_(n_) = i;
00194         // Check if xl is bounded
00195         if( *xl_full != -inf_bnd )
00196           ++xl_nz;
00197         // Check if xu is bounded
00198         if( *xu_full != inf_bnd )
00199           ++xu_nz;
00200         if( *xl_full != -inf_bnd || *xu_full != inf_bnd )
00201           ++num_bnd_x;
00202       }
00203     }
00204   }
00205   else {
00206     // None of the variables are fixed by bounds because there are no bounds
00207     n_ = n_full_;
00208     DenseLinAlgPack::identity_perm( &var_full_to_fixed_ );
00209   }
00210 
00211 //  std::cerr << "n_ =" << n_ << std::endl;
00212 //  std::cerr << "var_full_to_fixed_ =\n" << var_full_to_fixed_;
00213   
00214   num_bounded_x_ = num_bnd_x;
00215 
00216   // Validate that we still have a solvable problem
00217   TEST_FOR_EXCEPTION(
00218     n_ < m_full_, InvalidInitialization
00219     ,"NLPSerialPreprocess::initialize() : Error, after removing fixed "
00220     << "variables, n = " << n_ << " < m = " << m_full_
00221     << ", and the NLP is over determined and can not be solved!" );
00222 
00223   // Initialize inverse of var_full_to_fixed_
00224   DenseLinAlgPack::inv_perm( var_full_to_fixed_, &inv_var_full_to_fixed_ );
00225 
00226 //  std::cerr << "inv_var_full_to_fixed_ =\n"  << inv_var_full_to_fixed_;
00227 
00228   var_perm_.resize(n_);
00229   space_x_.initialize(n_);
00230   
00231   // Resize xinit, xl, xu, hl and hu
00232   xinit_.initialize(n_);
00233   xl_.initialize(n_);
00234   xu_.initialize(n_);
00235   if(mI_orig_) {
00236     hl_breve_.initialize(mI_orig_);
00237     hu_breve_.initialize(mI_orig_);
00238   }
00239 
00240   if( m_full_ ) {
00241     // Get the first basis
00242     if( !nlp_selects_basis() ) {
00243       // The NLP is not selecting the first basis so set to the initial basis to
00244       // the indentity permutations and assume full column rank for Gc.
00245       DenseLinAlgPack::identity_perm(&var_perm_);
00246 //      std::cerr << "var_perm_ =\n" << var_perm_;
00247       DenseLinAlgPack::identity_perm(&equ_perm_);
00248 //      std::cerr << "equ_perm_ =\n" << equ_perm_;
00249       DenseLinAlgPack::identity_perm(&inv_equ_perm_);
00250 //      std::cerr << "inv_equ_perm_ =\n" << inv_equ_perm_;
00251       r_ = m_full_;
00252       var_from_full( xinit_full_().begin(), xinit_.set_vec().begin() );
00253       if(has_var_bounds) {
00254         var_from_full( xl_full_().begin(), xl_.set_vec().begin() );
00255         var_from_full( xu_full_().begin(), xu_.set_vec().begin() );
00256         do_force_xinit_in_bounds();
00257       }
00258       else {
00259         xl_ = -inf_bnd;
00260         xu_ = +inf_bnd;
00261       }
00262     }
00263     else {
00264       // The nlp subclass is selecting the first basis.
00265       
00266       // make intialized_ true temporaraly so you can call get_next_basis()
00267       // and assert_initialized() called in it will not throw an exception.
00268       initialized_ = true;
00269       
00270       try {
00271         size_type rank;
00272         const bool 
00273            get_next_basis_return = get_next_basis_remove_fixed( &var_perm_, &equ_perm_, &rank );
00274         TEST_FOR_EXCEPTION(
00275           !get_next_basis_return, std::logic_error
00276           ,"NLPSerialPreprocess::initialize():  "
00277           " If nlp_selects_basis() is true then imp_get_next_basis() "
00278           " must return true for the first call" );
00279         assert_and_set_basis( var_perm_, equ_perm_, rank );
00280 //        std::cerr << "var_perm_ =\n" << var_perm_;
00281 //        std::cerr << "equ_perm_ =\n" << equ_perm_;
00282       }
00283       catch(...) {
00284         // In case an exception was thrown I don't want to leave #this#
00285         // in an inconsistant state.
00286         initialized_ = false;
00287         throw;
00288       }
00289       
00290       initialized_ = false; // resize to false to continue initialization
00291     }
00292   }
00293   else {
00294     DenseLinAlgPack::identity_perm(&var_perm_);
00295     r_ = 0;
00296     var_from_full( xinit_full_().begin(), xinit_.set_vec().begin() );
00297     if(has_var_bounds) {
00298       var_from_full( xl_full_().begin(), xl_.set_vec().begin() );
00299       var_from_full( xu_full_().begin(), xu_.set_vec().begin() );
00300       do_force_xinit_in_bounds();
00301     }
00302     else {
00303       xl_ = -inf_bnd;
00304       xu_ = +inf_bnd;
00305     }
00306   }
00307   
00308 //  std::cerr << "n_full_ = " << n_full_ << std::endl;
00309 //  std::cerr << "n_ = " << n_ << std::endl;
00310 //  std::cerr << "var_full_to_fixed_ =\n" << var_full_to_fixed_;
00311 //  std::cerr << "inv_var_full_to_fixed_ =\n"  << inv_var_full_to_fixed_;
00312 //  std::cerr << "var_perm_ =\n" << var_perm_;
00313 //  std::cerr << "equ_perm_ =\n" << equ_perm_;
00314 
00315   // If you get here then the initialization went Ok.
00316   NLPObjGrad::initialize(test_setup);
00317   initialized_ = true;
00318 }
00319 
00320 bool NLPSerialPreprocess::is_initialized() const
00321 {
00322   return initialized_;
00323 }
00324 
00325 size_type NLPSerialPreprocess::n() const 
00326 {
00327   assert_initialized();
00328   return n_;
00329 }
00330 
00331 size_type NLPSerialPreprocess::m() const 
00332 { 
00333   assert_initialized(); 
00334   return m_full_;
00335 }
00336 
00337 NLP::vec_space_ptr_t NLPSerialPreprocess::space_x() const
00338 {
00339   namespace mmp = MemMngPack;
00340   return Teuchos::rcp(&space_x_,false);
00341 }
00342 
00343 NLP::vec_space_ptr_t NLPSerialPreprocess::space_c() const
00344 {
00345   namespace mmp = MemMngPack;
00346   return ( m_full_ ? Teuchos::rcp(&space_c_,false) : Teuchos::null );
00347 }
00348 
00349 size_type NLPSerialPreprocess::num_bounded_x() const 
00350 {
00351   return num_bounded_x_;
00352 }
00353 
00354 const Vector& NLPSerialPreprocess::xl() const 
00355 {
00356   assert_initialized();
00357   return xl_;
00358 }
00359 
00360 const Vector& NLPSerialPreprocess::xu() const 
00361 {
00362   assert_initialized();
00363   return xu_;
00364 }
00365 
00366 const Vector& NLPSerialPreprocess::xinit() const 
00367 {
00368   assert_initialized();
00369   return xinit_;
00370 }
00371 
00372 void NLPSerialPreprocess::get_init_lagrange_mult(
00373   VectorMutable*   lambda
00374   ,VectorMutable*  nu
00375   ) const
00376 {
00377   assert_initialized();
00378   // ToDo: Get subclass to define what these are!
00379   if(lambda)
00380     *lambda   = 0.0;
00381   if(nu)
00382     *nu      = 0.0;
00383 }
00384 
00385 void NLPSerialPreprocess::scale_f( value_type scale_f )
00386 {
00387   assert_initialized();
00388   scale_f_ = scale_f;
00389 }
00390 
00391 value_type NLPSerialPreprocess::scale_f() const
00392 {
00393   assert_initialized();
00394   return scale_f_;
00395 }
00396 
00397 void NLPSerialPreprocess::report_final_solution(
00398   const Vector&    x
00399   ,const Vector*   lambda
00400   ,const Vector*   nu
00401   ,bool            is_optimal
00402   )
00403 {
00404   assert_initialized();
00405   // set x_full
00406   VectorDenseEncap  x_d(x);
00407   DVector x_full(n_full_);
00408   x_full = x_full_; // set any fixed components (as well as the others at first)
00409   var_to_full( x_d().begin(), x_full().begin() ); // set the nonfixed components
00410   // set lambda_full
00411   DVector lambda_full;
00412   if( lambda ) {
00413     VectorDenseEncap lambda_d(*lambda);
00414     DVectorSlice      lambda = lambda_d();
00415     lambda_full.resize(m_full_);
00416     for(size_type j = 1; j <= m_full_; j++)
00417       lambda_full(equ_perm_(j)) = lambda(j);
00418   }
00419   // set nu_full
00420   DVector nu_full(n_full_);
00421   if(nu) {
00422     nu_full = 0.0; // We don't give lagrange multipliers for fixed varaibles!
00423     // ToDo: Define a special constrant for multiplier values for fixed variables 
00424     VectorDenseEncap nu_d(*nu);
00425     var_to_full( nu_d().begin(), nu_full().begin() ); // set the nonfixed components
00426   }
00427   // Report the final solution
00428   DVectorSlice
00429     lambda_orig   = lambda && m_orig_ ? lambda_full(1,m_orig_) : DVectorSlice(),
00430     lambdaI_orig  = ( lambda && m_full_ > m_orig_ 
00431               ? lambda_full(m_orig_+1,m_full_)
00432               : DVectorSlice() ),
00433     nu_orig       = nu ? nu_full(1,n_orig_) : DVectorSlice();
00434   imp_report_orig_final_solution(
00435     x_full()
00436     ,lambda_orig.dim()  ? &lambda_orig  : NULL
00437     ,lambdaI_orig.dim() ? &lambdaI_orig : NULL
00438     ,nu_orig.dim()      ? &nu_orig      : NULL
00439     ,is_optimal
00440     );
00441 }
00442 
00443 size_type NLPSerialPreprocess::ns() const
00444 {
00445   assert_initialized();
00446   return mI_orig_;
00447 }
00448 
00449 NLP::vec_space_ptr_t
00450 NLPSerialPreprocess::space_c_breve() const
00451 {
00452   namespace mmp = MemMngPack;
00453   assert_initialized();
00454   return ( m_orig_ ? Teuchos::rcp(&space_c_breve_,false) : Teuchos::null );
00455 } 
00456 NLP::vec_space_ptr_t
00457 NLPSerialPreprocess::space_h_breve() const
00458 {
00459   namespace mmp = MemMngPack;
00460   assert_initialized();
00461   return ( mI_orig_ ? Teuchos::rcp(&space_h_breve_,false) : Teuchos::null );
00462 }
00463 
00464 const Vector& NLPSerialPreprocess::hl_breve() const
00465 {
00466   assert_initialized();
00467   return hl_breve_;
00468 }
00469 
00470 const Vector& NLPSerialPreprocess::hu_breve() const
00471 {
00472   assert_initialized();
00473   return hu_breve_;
00474 }
00475 
00476 const Permutation& NLPSerialPreprocess::P_var() const
00477 {
00478   assert_initialized();
00479   return P_var_;
00480 }
00481 
00482 const Permutation& NLPSerialPreprocess::P_equ() const
00483 {
00484   assert_initialized();
00485   return P_equ_;
00486 }
00487 
00488 // Overridden public members from NLPVarReductPerm
00489 
00490 const NLPVarReductPerm::perm_fcty_ptr_t
00491 NLPSerialPreprocess::factory_P_var() const
00492 {
00493   return factory_P_var_;
00494 }
00495 
00496 const NLPVarReductPerm::perm_fcty_ptr_t
00497 NLPSerialPreprocess::factory_P_equ() const
00498 {
00499   return factory_P_equ_;
00500 }
00501 
00502 Range1D NLPSerialPreprocess::var_dep() const
00503 {
00504   assert_initialized();
00505   return r_ ? Range1D(1,r_) : Range1D::Invalid;
00506 }
00507 
00508 Range1D NLPSerialPreprocess::var_indep() const
00509 {
00510   assert_initialized();
00511   return Range1D(r_+1,n_);
00512 }
00513 
00514 Range1D NLPSerialPreprocess::equ_decomp() const
00515 {
00516   assert_initialized();
00517   return r_ ? Range1D(1,r_) : Range1D::Invalid;
00518 }
00519 
00520 Range1D NLPSerialPreprocess::equ_undecomp() const
00521 {
00522   assert_initialized();
00523   return r_ < m_full_ ? Range1D(r_+1,m_full_) : Range1D::Invalid;
00524 }
00525 
00526 bool NLPSerialPreprocess::nlp_selects_basis() const
00527   {
00528   // Check if the user has supplied a basis from a file
00529   char ind[17];
00530   sprintf(ind, "%d", basis_selection_num_);
00531   std::string fname = "basis_";
00532   fname += ind;
00533   fname += ".sel";
00534 
00535   std::ifstream basis_file(fname.c_str());
00536   if (basis_file)
00537     {
00538     return true;
00539     }
00540 
00541   return false;
00542   }
00543 
00544 bool NLPSerialPreprocess::get_next_basis(
00545   Permutation*  P_var,   Range1D* var_dep
00546   ,Permutation* P_equ,   Range1D* equ_decomp
00547   )
00548 {
00549   assert_initialized();
00550   size_type rank = 0;
00551   const bool 
00552     get_next_basis_return = get_next_basis_remove_fixed( &var_perm_, &equ_perm_, &rank );
00553   if(get_next_basis_return)
00554     assert_and_set_basis( var_perm_, equ_perm_, rank );
00555   else
00556     return false; // The NLP subclass did not have a new basis to give us!
00557   this->get_basis(P_var,var_dep,P_equ,equ_decomp);
00558   return true;
00559 }
00560 
00561 void NLPSerialPreprocess::set_basis(
00562   const Permutation   &P_var,   const Range1D  &var_dep
00563   ,const Permutation  *P_equ,   const Range1D  *equ_decomp
00564   )
00565 {
00566   namespace mmp = MemMngPack;
00567   using Teuchos::dyn_cast;
00568   TEST_FOR_EXCEPTION(
00569     (m_full_ > 0 && (P_equ == NULL || equ_decomp == NULL))
00570     ,std::invalid_argument
00571     ,"NLPSerialPreprocess::set_basis(...) : Error!" );
00572   TEST_FOR_EXCEPTION(
00573     m_full_ > 0 && var_dep.size() != equ_decomp->size()
00574     ,InvalidBasis
00575     ,"NLPSerialPreprocess::set_basis(...) : Error!" );
00576   // Get the concrete types
00577   const PermutationSerial
00578     &P_var_s   = dyn_cast<const PermutationSerial>(P_var),
00579     *P_equ_s   = m_full_  ? &dyn_cast<const PermutationSerial>(*P_equ)   : NULL;
00580   // Get the underlying permutation vectors
00581   Teuchos::RCP<IVector>
00582     var_perm   = Teuchos::rcp_const_cast<IVector>(P_var_s.perm()),
00583     equ_perm   = ( m_full_
00584              ? Teuchos::rcp_const_cast<IVector>(P_equ_s->perm())
00585              : Teuchos::null );
00586   TEST_FOR_EXCEPTION(
00587     (m_full_ > 0 && equ_perm.get() == NULL)
00588     ,std::invalid_argument
00589     ,"NLPSerialPreprocess::set_basis(...) : Error, P_equ is not initialized properly!" );
00590   // Set the basis
00591   assert_and_set_basis( *var_perm, *equ_perm, var_dep.size() );
00592 }
00593 
00594 void NLPSerialPreprocess::get_basis(
00595   Permutation*  P_var,   Range1D* var_dep
00596   ,Permutation* P_equ,   Range1D* equ_decomp
00597   ) const
00598 {
00599   namespace mmp = MemMngPack;
00600   using Teuchos::dyn_cast;
00601   assert_initialized();
00602   TEST_FOR_EXCEPTION(
00603     P_var == NULL || var_dep == NULL
00604     || (m_full_ > 0 && (P_equ == NULL || equ_decomp == NULL))
00605     ,std::invalid_argument
00606     ,"NLPSerialPreprocess::get_basis(...) : Error!" );
00607   // Get the concrete types
00608   PermutationSerial
00609     &P_var_s   = dyn_cast<PermutationSerial>(*P_var),
00610     *P_equ_s   = m_full_  ? &dyn_cast<PermutationSerial>(*P_equ)   : NULL;
00611   // Get the underlying permutation vectors
00612   Teuchos::RCP<IVector>
00613     var_perm   = Teuchos::rcp_const_cast<IVector>(P_var_s.perm()),
00614     equ_perm   = ( m_full_
00615              ? Teuchos::rcp_const_cast<IVector>(P_equ_s->perm())
00616              : Teuchos::null );
00617   // Allocate permutation vectors if none allocated yet or someone else has reference to them
00618   if( var_perm.get() == NULL || var_perm.count() > 2 ) // P_var reference and my reference
00619     var_perm = Teuchos::rcp( new IVector(n_) );
00620   if( m_full_ && ( equ_perm.get() == NULL || equ_perm.count() > 2 ) ) // P_equ reference and my reference
00621     equ_perm = Teuchos::rcp( new IVector(m_full_) );
00622   // Copy the basis selection
00623   (*var_perm)   = var_perm_;
00624   (*var_dep)    = Range1D(1,r_);
00625   if(m_full_) {
00626     (*equ_perm)   = equ_perm_;
00627     (*equ_decomp) = Range1D(1,r_);
00628   }
00629   // Reinitialize the Permutation arguments.
00630   P_var_s.initialize( var_perm, Teuchos::null, true );  // Allocate the inverse permuation as well!
00631   if(m_full_)
00632     P_equ_s->initialize( equ_perm, Teuchos::null, true );
00633 }
00634 
00635 // Overridden protected members from NLP
00636 
00637 void NLPSerialPreprocess::imp_calc_f(
00638   const Vector            &x
00639   ,bool                   newx
00640   ,const ZeroOrderInfo    &zero_order_info
00641   ) const
00642 {
00643   assert_initialized();
00644   VectorDenseEncap  x_d(x);
00645   set_x_full( x_d(), newx, &x_full_() );
00646   imp_calc_f_orig( x_full(), newx, zero_order_orig_info() );
00647   *zero_order_info.f = scale_f_ * f_orig_;
00648 }
00649 
00650 void NLPSerialPreprocess::imp_calc_c(
00651   const Vector            &x
00652   ,bool                   newx
00653   ,const ZeroOrderInfo    &zero_order_info
00654   ) const
00655 {
00656   assert_initialized();
00657   VectorDenseEncap  x_d(x);
00658   set_x_full( x_d(), newx, &x_full_() );
00659   if( m_orig_ )
00660     imp_calc_c_orig( x_full(), newx, zero_order_orig_info() );
00661   if( mI_orig_ )
00662     imp_calc_h_orig( x_full(), newx, zero_order_orig_info() );
00663   VectorDenseMutableEncap  c_d(*zero_order_info.c);
00664   equ_from_full(
00665     m_orig_   ? c_orig_()                     : DVectorSlice()
00666     ,mI_orig_ ? h_orig_()                     : DVectorSlice()
00667     ,mI_orig_ ? x_full()(n_orig_+1,n_full_)   : DVectorSlice() // s_orig
00668     ,&c_d()
00669     );
00670 }
00671 
00672 void NLPSerialPreprocess::imp_calc_c_breve(
00673   const Vector            &x
00674   ,bool                   newx
00675   ,const ZeroOrderInfo    &zero_order_info_breve
00676   ) const
00677 {
00678   assert_initialized();
00679   VectorDenseEncap x_d(x);
00680   set_x_full( x_d(), newx, &x_full_() );
00681   if( m_orig_ )
00682     imp_calc_c_orig( x_full(), newx, zero_order_orig_info() );
00683   if( mI_orig_ )
00684     imp_calc_h_orig( x_full(), newx, zero_order_orig_info() );
00685   VectorDenseMutableEncap  c_breve_d(*zero_order_info_breve.c);
00686   c_breve_d() = c_orig_();
00687 }
00688 
00689 void NLPSerialPreprocess::imp_calc_h_breve(
00690   const Vector            &x
00691   ,bool                   newx
00692   ,const ZeroOrderInfo    &zero_order_info_breve
00693   ) const
00694 {
00695   // If this function gets called then this->mI() > 0 must be true
00696   // which means that convert_inequ_to_equ must be false!
00697   assert_initialized();
00698   VectorDenseEncap  x_d(x);
00699   set_x_full( x_d(), newx, &x_full_() );
00700   imp_calc_h_orig( x_full(), newx, zero_order_orig_info() );
00701   VectorDenseMutableEncap  h_breve_d(*zero_order_info_breve.h);
00702   h_breve_d() = h_orig_(); // Nothing fancy right now
00703 }
00704 
00705 // Overridden protected members from NLPObjGrad
00706 
00707 void NLPSerialPreprocess::imp_calc_Gf(
00708   const Vector            &x
00709   ,bool                   newx
00710   ,const ObjGradInfo      &obj_grad_info
00711   ) const
00712 {
00713   using DenseLinAlgPack::Vt_S;
00714   assert_initialized();
00715   VectorDenseEncap  x_d(x);
00716   set_x_full( x_d(), newx, &x_full_() );
00717   if( n_full_ > n_orig_ ) Gf_full_(n_orig_+1,n_full_) = 0.0; // Initialize terms for slacks to zero!
00718   imp_calc_Gf_orig( x_full(), newx, obj_grad_orig_info() );
00719   VectorDenseMutableEncap  Gf_d(*obj_grad_info.Gf);
00720   var_from_full( Gf_full_.begin(), Gf_d().begin() );
00721   Vt_S( &Gf_d(), scale_f_ );
00722 }
00723 
00724 // protected members
00725 
00726 bool NLPSerialPreprocess::imp_get_next_basis(
00727   IVector      *var_perm_full
00728   ,IVector     *equ_perm_full
00729   ,size_type   *rank_full
00730   ,size_type   *rank
00731   )
00732 {
00733   return false; // default is that the subclass does not select the basis
00734 }
00735 
00736 void NLPSerialPreprocess::assert_initialized() const
00737 {
00738   TEST_FOR_EXCEPTION(
00739     !initialized_, UnInitialized
00740     ,"NLPSerialPreprocess : The nlp has not been initialized yet" );
00741 }
00742 
00743 void NLPSerialPreprocess::set_x_full(
00744   const DVectorSlice& x, bool newx
00745   ,DVectorSlice* x_full
00746   ) const
00747 {
00748   DenseLinAlgPack::assert_vs_sizes(x.dim(),n_);
00749   if(newx)
00750     var_to_full(x.begin(), x_full->begin());
00751 }
00752 
00753 void NLPSerialPreprocess::var_from_full(
00754   DVectorSlice::const_iterator   vec_full
00755   ,DVectorSlice::iterator        vec
00756   ) const
00757 {
00758 //  std::cout << "\nvar_from_full(...) : ";
00759   for(size_type i = 1; i <= n_; i++) {
00760     *vec++ = vec_full[var_full_to_fixed_(var_perm_(i)) - 1];
00761 //    std::cout
00762 //      << "\ni = " << i
00763 //      << "\nvar_perm_(i) = " << var_perm_(i)
00764 //      << "\nvar_full_to_fixed_(var_perm_(i)) = " << var_full_to_fixed_(var_perm_(i))
00765 //      << "\nvec_full[var_full_to_fixed_(var_perm_(i)) - 1] = " << vec_full[var_full_to_fixed_(var_perm_(i)) - 1]
00766 //      << "\nvec[i] = " << *(vec-1) << "\n\n";
00767   }   
00768 }
00769 
00770 void NLPSerialPreprocess::var_to_full( DVectorSlice::const_iterator vec, DVectorSlice::iterator vec_full ) const
00771 {
00772   for(size_type i = 1; i <= n_; i++)
00773     vec_full[var_full_to_fixed_(var_perm_(i)) - 1] = *vec++;
00774 }
00775 
00776 void NLPSerialPreprocess::equ_from_full(
00777   const DVectorSlice   &c_orig
00778   ,const DVectorSlice  &h_orig
00779   ,const DVectorSlice  &s_orig
00780   ,DVectorSlice        *c_full
00781   ) const
00782 {
00783   size_type i;
00784   // c_full = [ c_orig; h_orig - s_orig ]
00785    for(i = 1; i <= m_orig_; ++i)
00786     (*c_full)(inv_equ_perm_(i)) = c_orig(i);
00787    for(i = 1; i <= mI_orig_; ++i)
00788     (*c_full)(inv_equ_perm_(m_orig_+i)) = h_orig(i) - s_orig(i);
00789 }
00790 
00791 // private members
00792 
00793 bool NLPSerialPreprocess::get_next_basis_remove_fixed(
00794   IVector* var_perm, IVector* equ_perm, size_type* rank
00795   )
00796 {
00797   IVector var_perm_full(n_full_);
00798   equ_perm->resize(m_full_);
00799   size_type rank_full, rank_fixed_removed;
00800   if( imp_get_next_basis( &var_perm_full, equ_perm, &rank_full, &rank_fixed_removed ) ) {
00801 //    std::cerr << "var_perm_full =\n"  << var_perm_full;
00802 //    std::cerr << "equ_perm =\n"  << *equ_perm;
00803 //    std::cerr << "rank_full = "  << rank_full << std::endl;
00804     //
00805     // The NLP subclass has another basis to select
00806     //
00807     // Translate the basis by removing variables fixed by bounds.
00808     // This is where it is important that var_perm_full is
00809     // sorted in assending order for basis and non-basis variables
00810     //
00811     // This is a little bit of a tricky algorithm.  We have to
00812     // essentially loop through the set of basic and non-basic
00813     // variables, remove fixed variables and adjust the indexes
00814     // of the remaining variables.  Since the set of indexes in
00815     // the basic and non-basis sets are sorted, this will not
00816     // be too bad of an algorithm.
00817 
00818     // Iterator for the fixed variables that we are to remove
00819     IVector::const_iterator     fixed_itr     = var_full_to_fixed_.begin() + n_;
00820     IVector::const_iterator     fixed_end     = var_full_to_fixed_.end();
00821 
00822     // Iterator for the basis variables
00823     IVector::iterator           basic_itr  = var_perm_full.begin();
00824     IVector::iterator           basic_end  = basic_itr + rank_full;
00825 
00826     // Iterator for the non-basis variables
00827     IVector::iterator           nonbasic_itr  = basic_end;
00828     IVector::iterator           nonbasic_end  = var_perm_full.end();
00829     
00830     // Count the number of fixed basic and non-basic variables
00831     index_type
00832       count_fixed          = 0,
00833       count_basic_fixed    = 0,
00834       count_nonbasic_fixed = 0;
00835 
00836     // Loop through all of the fixed variables and remove and compact
00837     for( ; fixed_itr != fixed_end; ++fixed_itr ) {
00838       const index_type
00839         next_fixed = ( fixed_itr +1 != fixed_end ? *(fixed_itr+1) : n_full_+1);
00840       // Bring the basic and nonbasic variables up to this fixed variable
00841       for( ; *basic_itr < *fixed_itr; ++basic_itr )
00842         *(basic_itr - count_basic_fixed) = *basic_itr - count_fixed;
00843       for( ; *nonbasic_itr < *fixed_itr; ++nonbasic_itr )
00844         *(nonbasic_itr - count_nonbasic_fixed) = *nonbasic_itr - count_fixed;
00845       // Update the count of the fixed variables
00846       if( *basic_itr == *fixed_itr ) {
00847         ++count_basic_fixed;
00848         ++basic_itr;
00849       }
00850       else {
00851         TEST_FOR_EXCEPT( !( *nonbasic_itr == *fixed_itr ) ); // If basic was not fixed then nonbasic better be!
00852         ++count_nonbasic_fixed;
00853         ++nonbasic_itr;
00854 
00855       }
00856       ++count_fixed;
00857       // Now update the indexes until the next fixed variable
00858       for( ; *basic_itr < next_fixed; ++basic_itr )
00859         *(basic_itr - count_basic_fixed) = *basic_itr - count_fixed;
00860       for( ; *nonbasic_itr < next_fixed; ++nonbasic_itr )
00861         *(nonbasic_itr - count_nonbasic_fixed) = *nonbasic_itr - count_fixed;
00862     }
00863     TEST_FOR_EXCEPT( !( count_fixed == n_full_ - n_ ) ); // Basic check
00864 
00865     var_perm->resize(n_);
00866     std::copy(
00867       var_perm_full.begin()
00868       ,var_perm_full.begin() + rank_fixed_removed
00869       ,var_perm->begin()
00870       );
00871     std::copy(
00872       var_perm_full.begin() + rank_full
00873       ,var_perm_full.begin() + rank_full + (n_-rank_fixed_removed)
00874       ,var_perm->begin() + rank_fixed_removed
00875       );
00876     *rank = rank_fixed_removed;
00877     return true;
00878   }
00879   else {
00880   
00881   // try to find the file giving the basis...
00882   char ind[17];
00883   sprintf(ind, "%d", basis_selection_num_);
00884   std::string fname = "basis_";
00885   fname += ind;
00886   fname += ".sel";
00887   basis_selection_num_++;
00888 
00889   std::ifstream basis_file(fname.c_str());
00890   if (basis_file)
00891     {
00892     // try to read the basis file
00893     std::string tags;
00894 
00895     int n;
00896     basis_file >> tags;
00897     TEST_FOR_EXCEPTION(
00898       tags != "n", std::logic_error
00899       ,"Incorrect basis file format - \"n\" expected, \"" << tags << "\" found.");
00900     basis_file >> n;
00901     TEST_FOR_EXCEPTION(
00902       n <= 0, std::logic_error
00903       , "Incorrect basis file format - n > 0 \"" << n << "\" found.");
00904 
00905     int m;
00906     basis_file >> tags;
00907     TEST_FOR_EXCEPTION(
00908       tags != "m", std::logic_error
00909       ,"Incorrect basis file format - \"m\" expected, \"" << tags << "\" found.");
00910     basis_file >> m;
00911     TEST_FOR_EXCEPTION(
00912       m > n , std::logic_error
00913       ,"Incorrect basis file format - 0 < m <= n expected, \"" << m << "\" found.");
00914     
00915     int r;
00916     basis_file >> tags;
00917     TEST_FOR_EXCEPTION(
00918       tags != "rank", std::logic_error
00919       ,"Incorrect basis file format - \"rank\" expected, \"" << tags << "\" found.");
00920     basis_file >> r;
00921     TEST_FOR_EXCEPTION(
00922       r > m, std::logic_error
00923       ,"Incorrect basis file format - 0 < rank <= m expected, \"" << r << "\" found.");   
00924     if (rank)
00925       { *rank = r; }
00926 
00927     // var_permutation
00928     basis_file >> tags;
00929     TEST_FOR_EXCEPTION(
00930       tags != "var_perm", std::logic_error
00931       ,"Incorrect basis file format -\"var_perm\" expected, \"" << tags << "\" found.");
00932     var_perm->resize(n);
00933     {for (int i=0; i < n; i++)
00934       {
00935       int var_index;
00936       basis_file >> var_index;
00937       TEST_FOR_EXCEPTION(
00938         var_index < 1 || var_index > n, std::logic_error
00939         ,"Incorrect basis file format for var_perm: 1 <= indice <= n expected, \"" << n << "\" found.");
00940       (*var_perm)[i] = var_index;
00941       }}
00942 
00943     // eqn_permutation
00944     basis_file >> tags;
00945     TEST_FOR_EXCEPTION(
00946       tags != "equ_perm", std::logic_error
00947       ,"Incorrect basis file format -\"equ_perm\" expected, \"" << tags << "\" found.");
00948     equ_perm->resize(r);
00949     {for (int i=0; i < r; i++)
00950       {
00951       int equ_index;
00952       basis_file >> equ_index;
00953       TEST_FOR_EXCEPTION(
00954         equ_index < 1 || equ_index > m, std::logic_error
00955         ,"Incorrect basis file format for equ_perm: 1 <= indice <= m expected, \"" << m << "\" found.");
00956       (*equ_perm)[i] = equ_index;
00957       }}
00958 
00959     return true;
00960     }
00961   }
00962 
00963   return false;
00964 }
00965 
00966 void NLPSerialPreprocess::assert_and_set_basis(
00967   const IVector& var_perm, const IVector& equ_perm, size_type rank
00968   )
00969 {
00970   namespace mmp = MemMngPack;
00971 
00972   // Assert that this is a valid basis and set the internal basis.  Also repivot 'xinit', 
00973   // 'xl', and 'xu'.
00974 
00975   const value_type inf_bnd = NLPSerialPreprocess::infinite_bound();
00976 
00977   // Assert the preconditions
00978   TEST_FOR_EXCEPTION(
00979     var_perm.size() != n_ || equ_perm.size() != m_full_, std::length_error
00980     ,"NLPSerialPreprocess::set_basis():  The sizes "
00981     "of the permutation vectors does not match the size of the NLP" );
00982   TEST_FOR_EXCEPTION(
00983     rank > m_full_, InvalidBasis
00984     ,"NLPSerialPreprocess::set_basis():  The rank "
00985     "of the basis can not be greater that the number of constraints" );
00986     
00987   // Set the basis
00988   r_ = rank;
00989   if( &var_perm_ != &var_perm )
00990     var_perm_ = var_perm;
00991   if( &equ_perm_ != &equ_perm )
00992     equ_perm_ = equ_perm;
00993   DenseLinAlgPack::inv_perm( equ_perm_, &inv_equ_perm_ );
00994 
00995   var_from_full( xinit_full_().begin(), xinit_.set_vec().begin() );
00996   if(num_bounded_x_) {
00997     var_from_full( xl_full_().begin(), xl_.set_vec().begin() );
00998     var_from_full( xu_full_().begin(), xu_.set_vec().begin() );
00999     do_force_xinit_in_bounds();
01000   }
01001   else {
01002     xl_ = -NLP::infinite_bound();
01003     xu_ = +NLP::infinite_bound();
01004   }
01005   P_var_.initialize(Teuchos::rcp(new IVector(var_perm)),Teuchos::null);
01006   P_equ_.initialize(Teuchos::rcp(new IVector(equ_perm)),Teuchos::null);
01007 }
01008 
01009 void NLPSerialPreprocess::assert_bounds_on_variables() const
01010 {
01011   TEST_FOR_EXCEPTION(
01012     !(imp_has_var_bounds() || n_full_ > n_orig_), NLP::NoBounds
01013     ,"There are no bounds on the variables for this NLP" );
01014 }
01015 
01016 void NLPSerialPreprocess::do_force_xinit_in_bounds()
01017 {
01018   AbstractLinAlgPack::force_in_bounds( xl_, xu_, &xinit_ );
01019 }
01020 
01021 } // end namespace NLPInterfacePack

Generated on Tue Jul 13 09:30:54 2010 for MOOCHO (Single Doxygen Collection) by  doxygen 1.4.7