MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NLPInterfacePack_NLPSerialPreprocess.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
44 #include <iostream> // Debug only
46 
47 #include <algorithm>
48 #include <sstream>
49 #include <limits>
50 #include <stdio.h>
51 #include <fstream>
52 
63 #include "Teuchos_Assert.hpp"
65 #include "Teuchos_dyn_cast.hpp"
66 
67 namespace LinAlgOpPack {
69 }
70 
71 namespace NLPInterfacePack {
72 
73 // NLPSerialPreprocess
74 
75 // Static public members
76 
79 {
80  return std::numeric_limits<DenseLinAlgPack::DVector::value_type>::max()-100; // Don't know what to use?
81 }
82 
83 // Constructors / nitializers
84 
86  )
87  :initialized_(false)
88  ,force_xinit_in_bounds_(true)
89  ,scale_f_(1.0)
90  ,basis_selection_num_(0)
91 
92 {}
93 
94 // Overridden public members from NLP
95 
96 void NLPSerialPreprocess::force_xinit_in_bounds(bool force_xinit_in_bounds)
97 {
99 }
100 
102 {
103  return force_xinit_in_bounds_;
104 }
105 
106 void NLPSerialPreprocess::initialize(bool test_setup)
107 {
108  namespace mmp = MemMngPack;
109 
110  const value_type inf_bnd = NLP::infinite_bound();
111 
113 
114  if( initialized_ && !imp_nlp_has_changed() ) {
115  // The subclass NLP has not changed so we can just
116  // slip this preprocessing.
117  NLPObjGrad::initialize(test_setup);
118  return;
119  }
120 
121  // Get the dimensions of the original problem
122 
123  n_orig_ = imp_n_orig();
124  m_orig_ = imp_m_orig(); // This may be zero!
125  mI_orig_ = imp_mI_orig(); // This may be zero!
126 
127  // Get the dimensions of the full problem
128 
131 
132  // Initialize the storage for the intermediate quanities
133 
139  h_orig_.resize(mI_orig_);
141  var_full_to_fixed_.resize(n_full_);
142  equ_perm_.resize(m_full_);
143  inv_equ_perm_.resize(m_full_);
144  space_c_.initialize(m_full_);
145  space_c_breve_.initialize(m_orig_);
146  space_h_breve_.initialize(mI_orig_);
149 
150  // Intialize xinit_full_, xl_full_ and xu_full_ for the initial point which will set the
151  // fixed elements which will not change during the optimization.
155  if( n_full_ > n_orig_ ) { // Include slack varaibles
156  xinit_full_(n_orig_+1,n_full_) = 0.0;
159  }
160 
161  const bool has_var_bounds = imp_has_var_bounds() || n_full_ > n_orig_;
162 
163  // Force the initial point in bounds if it is not.
164  if( force_xinit_in_bounds() && has_var_bounds ) {
166  VectorMutableDense( xl_full_(), Teuchos::null )
167  ,VectorMutableDense( xu_full_(), Teuchos::null )
168  ,&VectorMutableDense( x_full_(), Teuchos::null )
169  );
170  }
171 
172  // Determine which variables are fixed by bounds!
173  size_type
174  xl_nz = 0,
175  xu_nz = 0,
176  num_bnd_x = 0;
177  if( has_var_bounds ) {
178  // Determine which variables are fixed by bounds and
179  // adjust the bounds if needed.
180  DVector::iterator
181  xl_full = xl_full_.begin(),
182  xu_full = xu_full_.begin();
183  n_ = 0;
184  size_type num_fixed = 0;
185  for(int i = 1; i <= n_full_; ++i, ++xl_full, ++xu_full) {
187  *xl_full > *xu_full, InconsistantBounds
188  ,"NLPSerialPreprocess::initialize() : Error, Inconsistant bounds: xl_full("
189  << i << ") > xu_full(" << i << ")" );
190  if(*xl_full == *xu_full) {
191  //
192  // Fixed between bounds
193  //
194  var_full_to_fixed_(n_full_ - num_fixed) = i;
195  num_fixed++;
196  }
197  else {
198  //
199  // Not Fixed between bounds
200  //
201  // Adjust the bounds if needed
202  *xl_full = *xl_full < -inf_bnd ? -inf_bnd : *xl_full;
203  *xu_full = *xu_full > +inf_bnd ? +inf_bnd : *xu_full;
204  //
205  n_++;
206  var_full_to_fixed_(n_) = i;
207  // Check if xl is bounded
208  if( *xl_full != -inf_bnd )
209  ++xl_nz;
210  // Check if xu is bounded
211  if( *xu_full != inf_bnd )
212  ++xu_nz;
213  if( *xl_full != -inf_bnd || *xu_full != inf_bnd )
214  ++num_bnd_x;
215  }
216  }
217  }
218  else {
219  // None of the variables are fixed by bounds because there are no bounds
220  n_ = n_full_;
222  }
223 
224 // std::cerr << "n_ =" << n_ << std::endl;
225 // std::cerr << "var_full_to_fixed_ =\n" << var_full_to_fixed_;
226 
227  num_bounded_x_ = num_bnd_x;
228 
229  // Validate that we still have a solvable problem
232  ,"NLPSerialPreprocess::initialize() : Error, after removing fixed "
233  << "variables, n = " << n_ << " < m = " << m_full_
234  << ", and the NLP is over determined and can not be solved!" );
235 
236  // Initialize inverse of var_full_to_fixed_
238 
239 // std::cerr << "inv_var_full_to_fixed_ =\n" << inv_var_full_to_fixed_;
240 
241  var_perm_.resize(n_);
242  space_x_.initialize(n_);
243 
244  // Resize xinit, xl, xu, hl and hu
245  xinit_.initialize(n_);
246  xl_.initialize(n_);
247  xu_.initialize(n_);
248  if(mI_orig_) {
249  hl_breve_.initialize(mI_orig_);
250  hu_breve_.initialize(mI_orig_);
251  }
252 
253  if( m_full_ ) {
254  // Get the first basis
255  if( !nlp_selects_basis() ) {
256  // The NLP is not selecting the first basis so set to the initial basis to
257  // the indentity permutations and assume full column rank for Gc.
259 // std::cerr << "var_perm_ =\n" << var_perm_;
261 // std::cerr << "equ_perm_ =\n" << equ_perm_;
263 // std::cerr << "inv_equ_perm_ =\n" << inv_equ_perm_;
264  r_ = m_full_;
265  var_from_full( xinit_full_().begin(), xinit_.set_vec().begin() );
266  if(has_var_bounds) {
267  var_from_full( xl_full_().begin(), xl_.set_vec().begin() );
268  var_from_full( xu_full_().begin(), xu_.set_vec().begin() );
270  }
271  else {
272  xl_ = -inf_bnd;
273  xu_ = +inf_bnd;
274  }
275  }
276  else {
277  // The nlp subclass is selecting the first basis.
278 
279  // make intialized_ true temporaraly so you can call get_next_basis()
280  // and assert_initialized() called in it will not throw an exception.
281  initialized_ = true;
282 
283  try {
284  size_type rank;
285  const bool
286  get_next_basis_return = get_next_basis_remove_fixed( &var_perm_, &equ_perm_, &rank );
288  !get_next_basis_return, std::logic_error
289  ,"NLPSerialPreprocess::initialize(): "
290  " If nlp_selects_basis() is true then imp_get_next_basis() "
291  " must return true for the first call" );
293 // std::cerr << "var_perm_ =\n" << var_perm_;
294 // std::cerr << "equ_perm_ =\n" << equ_perm_;
295  }
296  catch(...) {
297  // In case an exception was thrown I don't want to leave #this#
298  // in an inconsistant state.
299  initialized_ = false;
300  throw;
301  }
302 
303  initialized_ = false; // resize to false to continue initialization
304  }
305  }
306  else {
308  r_ = 0;
309  var_from_full( xinit_full_().begin(), xinit_.set_vec().begin() );
310  if(has_var_bounds) {
311  var_from_full( xl_full_().begin(), xl_.set_vec().begin() );
312  var_from_full( xu_full_().begin(), xu_.set_vec().begin() );
314  }
315  else {
316  xl_ = -inf_bnd;
317  xu_ = +inf_bnd;
318  }
319  }
320 
321 // std::cerr << "n_full_ = " << n_full_ << std::endl;
322 // std::cerr << "n_ = " << n_ << std::endl;
323 // std::cerr << "var_full_to_fixed_ =\n" << var_full_to_fixed_;
324 // std::cerr << "inv_var_full_to_fixed_ =\n" << inv_var_full_to_fixed_;
325 // std::cerr << "var_perm_ =\n" << var_perm_;
326 // std::cerr << "equ_perm_ =\n" << equ_perm_;
327 
328  // If you get here then the initialization went Ok.
329  NLPObjGrad::initialize(test_setup);
330  initialized_ = true;
331 }
332 
334 {
335  return initialized_;
336 }
337 
339 {
341  return n_;
342 }
343 
345 {
347  return m_full_;
348 }
349 
351 {
352  namespace mmp = MemMngPack;
353  return Teuchos::rcp(&space_x_,false);
354 }
355 
357 {
358  namespace mmp = MemMngPack;
359  return ( m_full_ ? Teuchos::rcp(&space_c_,false) : Teuchos::null );
360 }
361 
363 {
364  return num_bounded_x_;
365 }
366 
368 {
370  return xl_;
371 }
372 
374 {
376  return xu_;
377 }
378 
380 {
382  return xinit_;
383 }
384 
386  VectorMutable* lambda
387  ,VectorMutable* nu
388  ) const
389 {
391  // ToDo: Get subclass to define what these are!
392  if(lambda)
393  *lambda = 0.0;
394  if(nu)
395  *nu = 0.0;
396 }
397 
399 {
401  scale_f_ = scale_f;
402 }
403 
405 {
407  return scale_f_;
408 }
409 
411  const Vector& x
412  ,const Vector* lambda
413  ,const Vector* nu
414  ,bool is_optimal
415  )
416 {
418  // set x_full
419  VectorDenseEncap x_d(x);
421  x_full = x_full_; // set any fixed components (as well as the others at first)
422  var_to_full( x_d().begin(), x_full().begin() ); // set the nonfixed components
423  // set lambda_full
424  DVector lambda_full;
425  if( lambda ) {
426  VectorDenseEncap lambda_d(*lambda);
427  DVectorSlice lambda = lambda_d();
428  lambda_full.resize(m_full_);
429  for(size_type j = 1; j <= m_full_; j++)
430  lambda_full(equ_perm_(j)) = lambda(j);
431  }
432  // set nu_full
433  DVector nu_full(n_full_);
434  if(nu) {
435  nu_full = 0.0; // We don't give lagrange multipliers for fixed varaibles!
436  // ToDo: Define a special constrant for multiplier values for fixed variables
437  VectorDenseEncap nu_d(*nu);
438  var_to_full( nu_d().begin(), nu_full().begin() ); // set the nonfixed components
439  }
440  // Report the final solution
442  lambda_orig = lambda && m_orig_ ? lambda_full(1,m_orig_) : DVectorSlice(),
443  lambdaI_orig = ( lambda && m_full_ > m_orig_
444  ? lambda_full(m_orig_+1,m_full_)
445  : DVectorSlice() ),
446  nu_orig = nu ? nu_full(1,n_orig_) : DVectorSlice();
448  x_full()
449  ,lambda_orig.dim() ? &lambda_orig : NULL
450  ,lambdaI_orig.dim() ? &lambdaI_orig : NULL
451  ,nu_orig.dim() ? &nu_orig : NULL
452  ,is_optimal
453  );
454 }
455 
457 {
459  return mI_orig_;
460 }
461 
464 {
465  namespace mmp = MemMngPack;
467  return ( m_orig_ ? Teuchos::rcp(&space_c_breve_,false) : Teuchos::null );
468 }
471 {
472  namespace mmp = MemMngPack;
474  return ( mI_orig_ ? Teuchos::rcp(&space_h_breve_,false) : Teuchos::null );
475 }
476 
478 {
480  return hl_breve_;
481 }
482 
484 {
486  return hu_breve_;
487 }
488 
489 const Permutation& NLPSerialPreprocess::P_var() const
490 {
492  return P_var_;
493 }
494 
495 const Permutation& NLPSerialPreprocess::P_equ() const
496 {
498  return P_equ_;
499 }
500 
501 // Overridden public members from NLPVarReductPerm
502 
505 {
506  return factory_P_var_;
507 }
508 
511 {
512  return factory_P_equ_;
513 }
514 
516 {
518  return r_ ? Range1D(1,r_) : Range1D::Invalid;
519 }
520 
522 {
524  return Range1D(r_+1,n_);
525 }
526 
528 {
530  return r_ ? Range1D(1,r_) : Range1D::Invalid;
531 }
532 
534 {
536  return r_ < m_full_ ? Range1D(r_+1,m_full_) : Range1D::Invalid;
537 }
538 
540  {
541  // Check if the user has supplied a basis from a file
542  char ind[17];
543  sprintf(ind, "%d", basis_selection_num_);
544  std::string fname = "basis_";
545  fname += ind;
546  fname += ".sel";
547 
548  std::ifstream basis_file(fname.c_str());
549  if (basis_file)
550  {
551  return true;
552  }
553 
554  return false;
555  }
556 
558  Permutation* P_var, Range1D* var_dep
559  ,Permutation* P_equ, Range1D* equ_decomp
560  )
561 {
563  size_type rank = 0;
564  const bool
565  get_next_basis_return = get_next_basis_remove_fixed( &var_perm_, &equ_perm_, &rank );
566  if(get_next_basis_return)
568  else
569  return false; // The NLP subclass did not have a new basis to give us!
570  this->get_basis(P_var,var_dep,P_equ,equ_decomp);
571  return true;
572 }
573 
575  const Permutation &P_var, const Range1D &var_dep
576  ,const Permutation *P_equ, const Range1D *equ_decomp
577  )
578 {
579  namespace mmp = MemMngPack;
580  using Teuchos::dyn_cast;
582  (m_full_ > 0 && (P_equ == NULL || equ_decomp == NULL))
583  ,std::invalid_argument
584  ,"NLPSerialPreprocess::set_basis(...) : Error!" );
586  m_full_ > 0 && var_dep.size() != equ_decomp->size()
587  ,InvalidBasis
588  ,"NLPSerialPreprocess::set_basis(...) : Error!" );
589  // Get the concrete types
590  const PermutationSerial
591  &P_var_s = dyn_cast<const PermutationSerial>(P_var),
592  *P_equ_s = m_full_ ? &dyn_cast<const PermutationSerial>(*P_equ) : NULL;
593  // Get the underlying permutation vectors
595  var_perm = Teuchos::rcp_const_cast<IVector>(P_var_s.perm()),
596  equ_perm = ( m_full_
597  ? Teuchos::rcp_const_cast<IVector>(P_equ_s->perm())
598  : Teuchos::null );
600  (m_full_ > 0 && equ_perm.get() == NULL)
601  ,std::invalid_argument
602  ,"NLPSerialPreprocess::set_basis(...) : Error, P_equ is not initialized properly!" );
603  // Set the basis
604  assert_and_set_basis( *var_perm, *equ_perm, var_dep.size() );
605 }
606 
608  Permutation* P_var, Range1D* var_dep
609  ,Permutation* P_equ, Range1D* equ_decomp
610  ) const
611 {
612  namespace mmp = MemMngPack;
613  using Teuchos::dyn_cast;
616  P_var == NULL || var_dep == NULL
617  || (m_full_ > 0 && (P_equ == NULL || equ_decomp == NULL))
618  ,std::invalid_argument
619  ,"NLPSerialPreprocess::get_basis(...) : Error!" );
620  // Get the concrete types
621  PermutationSerial
622  &P_var_s = dyn_cast<PermutationSerial>(*P_var),
623  *P_equ_s = m_full_ ? &dyn_cast<PermutationSerial>(*P_equ) : NULL;
624  // Get the underlying permutation vectors
626  var_perm = Teuchos::rcp_const_cast<IVector>(P_var_s.perm()),
627  equ_perm = ( m_full_
628  ? Teuchos::rcp_const_cast<IVector>(P_equ_s->perm())
629  : Teuchos::null );
630  // Allocate permutation vectors if none allocated yet or someone else has reference to them
631  if( var_perm.get() == NULL || var_perm.total_count() > 2 ) // P_var reference and my reference
632  var_perm = Teuchos::rcp( new IVector(n_) );
633  if( m_full_ && ( equ_perm.get() == NULL || equ_perm.total_count() > 2 ) ) // P_equ reference and my reference
634  equ_perm = Teuchos::rcp( new IVector(m_full_) );
635  // Copy the basis selection
636  (*var_perm) = var_perm_;
637  (*var_dep) = Range1D(1,r_);
638  if(m_full_) {
639  (*equ_perm) = equ_perm_;
640  (*equ_decomp) = Range1D(1,r_);
641  }
642  // Reinitialize the Permutation arguments.
643  P_var_s.initialize( var_perm, Teuchos::null, true ); // Allocate the inverse permuation as well!
644  if(m_full_)
645  P_equ_s->initialize( equ_perm, Teuchos::null, true );
646 }
647 
648 // Overridden protected members from NLP
649 
651  const Vector &x
652  ,bool newx
653  ,const ZeroOrderInfo &zero_order_info
654  ) const
655 {
657  VectorDenseEncap x_d(x);
658  set_x_full( x_d(), newx, &x_full_() );
660  *zero_order_info.f = scale_f_ * f_orig_;
661 }
662 
664  const Vector &x
665  ,bool newx
666  ,const ZeroOrderInfo &zero_order_info
667  ) const
668 {
670  VectorDenseEncap x_d(x);
671  set_x_full( x_d(), newx, &x_full_() );
672  if( m_orig_ )
674  if( mI_orig_ )
676  VectorDenseMutableEncap c_d(*zero_order_info.c);
679  ,mI_orig_ ? h_orig_() : DVectorSlice()
680  ,mI_orig_ ? x_full()(n_orig_+1,n_full_) : DVectorSlice() // s_orig
681  ,&c_d()
682  );
683 }
684 
686  const Vector &x
687  ,bool newx
688  ,const ZeroOrderInfo &zero_order_info_breve
689  ) const
690 {
692  VectorDenseEncap x_d(x);
693  set_x_full( x_d(), newx, &x_full_() );
694  if( m_orig_ )
696  if( mI_orig_ )
698  VectorDenseMutableEncap c_breve_d(*zero_order_info_breve.c);
699  c_breve_d() = c_orig_();
700 }
701 
703  const Vector &x
704  ,bool newx
705  ,const ZeroOrderInfo &zero_order_info_breve
706  ) const
707 {
708  // If this function gets called then this->mI() > 0 must be true
709  // which means that convert_inequ_to_equ must be false!
711  VectorDenseEncap x_d(x);
712  set_x_full( x_d(), newx, &x_full_() );
714  VectorDenseMutableEncap h_breve_d(*zero_order_info_breve.h);
715  h_breve_d() = h_orig_(); // Nothing fancy right now
716 }
717 
718 // Overridden protected members from NLPObjGrad
719 
721  const Vector &x
722  ,bool newx
723  ,const ObjGradInfo &obj_grad_info
724  ) const
725 {
726  using DenseLinAlgPack::Vt_S;
728  VectorDenseEncap x_d(x);
729  set_x_full( x_d(), newx, &x_full_() );
730  if( n_full_ > n_orig_ ) Gf_full_(n_orig_+1,n_full_) = 0.0; // Initialize terms for slacks to zero!
732  VectorDenseMutableEncap Gf_d(*obj_grad_info.Gf);
733  var_from_full( Gf_full_.begin(), Gf_d().begin() );
734  Vt_S( &Gf_d(), scale_f_ );
735 }
736 
737 // protected members
738 
740  IVector *var_perm_full
741  ,IVector *equ_perm_full
742  ,size_type *rank_full
743  ,size_type *rank
744  )
745 {
746  return false; // default is that the subclass does not select the basis
747 }
748 
750 {
753  ,"NLPSerialPreprocess : The nlp has not been initialized yet" );
754 }
755 
757  const DVectorSlice& x, bool newx
758  ,DVectorSlice* x_full
759  ) const
760 {
762  if(newx)
763  var_to_full(x.begin(), x_full->begin());
764 }
765 
767  DVectorSlice::const_iterator vec_full
768  ,DVectorSlice::iterator vec
769  ) const
770 {
771 // std::cout << "\nvar_from_full(...) : ";
772  for(size_type i = 1; i <= n_; i++) {
773  *vec++ = vec_full[var_full_to_fixed_(var_perm_(i)) - 1];
774 // std::cout
775 // << "\ni = " << i
776 // << "\nvar_perm_(i) = " << var_perm_(i)
777 // << "\nvar_full_to_fixed_(var_perm_(i)) = " << var_full_to_fixed_(var_perm_(i))
778 // << "\nvec_full[var_full_to_fixed_(var_perm_(i)) - 1] = " << vec_full[var_full_to_fixed_(var_perm_(i)) - 1]
779 // << "\nvec[i] = " << *(vec-1) << "\n\n";
780  }
781 }
782 
783 void NLPSerialPreprocess::var_to_full( DVectorSlice::const_iterator vec, DVectorSlice::iterator vec_full ) const
784 {
785  for(size_type i = 1; i <= n_; i++)
786  vec_full[var_full_to_fixed_(var_perm_(i)) - 1] = *vec++;
787 }
788 
790  const DVectorSlice &c_orig
791  ,const DVectorSlice &h_orig
792  ,const DVectorSlice &s_orig
793  ,DVectorSlice *c_full
794  ) const
795 {
796  size_type i;
797  // c_full = [ c_orig; h_orig - s_orig ]
798  for(i = 1; i <= m_orig_; ++i)
799  (*c_full)(inv_equ_perm_(i)) = c_orig(i);
800  for(i = 1; i <= mI_orig_; ++i)
801  (*c_full)(inv_equ_perm_(m_orig_+i)) = h_orig(i) - s_orig(i);
802 }
803 
804 // private members
805 
807  IVector* var_perm, IVector* equ_perm, size_type* rank
808  )
809 {
810  IVector var_perm_full(n_full_);
811  equ_perm->resize(m_full_);
812  size_type rank_full, rank_fixed_removed;
813  if( imp_get_next_basis( &var_perm_full, equ_perm, &rank_full, &rank_fixed_removed ) ) {
814 // std::cerr << "var_perm_full =\n" << var_perm_full;
815 // std::cerr << "equ_perm =\n" << *equ_perm;
816 // std::cerr << "rank_full = " << rank_full << std::endl;
817  //
818  // The NLP subclass has another basis to select
819  //
820  // Translate the basis by removing variables fixed by bounds.
821  // This is where it is important that var_perm_full is
822  // sorted in assending order for basis and non-basis variables
823  //
824  // This is a little bit of a tricky algorithm. We have to
825  // essentially loop through the set of basic and non-basic
826  // variables, remove fixed variables and adjust the indexes
827  // of the remaining variables. Since the set of indexes in
828  // the basic and non-basis sets are sorted, this will not
829  // be too bad of an algorithm.
830 
831  // Iterator for the fixed variables that we are to remove
832  IVector::const_iterator fixed_itr = var_full_to_fixed_.begin() + n_;
833  IVector::const_iterator fixed_end = var_full_to_fixed_.end();
834 
835  // Iterator for the basis variables
836  IVector::iterator basic_itr = var_perm_full.begin();
837  IVector::iterator basic_end = basic_itr + rank_full;
838 
839  // Iterator for the non-basis variables
840  IVector::iterator nonbasic_itr = basic_end;
841  IVector::iterator nonbasic_end = var_perm_full.end();
842 
843  // Count the number of fixed basic and non-basic variables
844  index_type
845  count_fixed = 0,
846  count_basic_fixed = 0,
847  count_nonbasic_fixed = 0;
848 
849  // Loop through all of the fixed variables and remove and compact
850  for( ; fixed_itr != fixed_end; ++fixed_itr ) {
851  const index_type
852  next_fixed = ( fixed_itr +1 != fixed_end ? *(fixed_itr+1) : n_full_+1);
853  // Bring the basic and nonbasic variables up to this fixed variable
854  for( ; *basic_itr < *fixed_itr; ++basic_itr )
855  *(basic_itr - count_basic_fixed) = *basic_itr - count_fixed;
856  for( ; *nonbasic_itr < *fixed_itr; ++nonbasic_itr )
857  *(nonbasic_itr - count_nonbasic_fixed) = *nonbasic_itr - count_fixed;
858  // Update the count of the fixed variables
859  if( *basic_itr == *fixed_itr ) {
860  ++count_basic_fixed;
861  ++basic_itr;
862  }
863  else {
864  TEUCHOS_TEST_FOR_EXCEPT( !( *nonbasic_itr == *fixed_itr ) ); // If basic was not fixed then nonbasic better be!
865  ++count_nonbasic_fixed;
866  ++nonbasic_itr;
867 
868  }
869  ++count_fixed;
870  // Now update the indexes until the next fixed variable
871  for( ; *basic_itr < next_fixed; ++basic_itr )
872  *(basic_itr - count_basic_fixed) = *basic_itr - count_fixed;
873  for( ; *nonbasic_itr < next_fixed; ++nonbasic_itr )
874  *(nonbasic_itr - count_nonbasic_fixed) = *nonbasic_itr - count_fixed;
875  }
876  TEUCHOS_TEST_FOR_EXCEPT( !( count_fixed == n_full_ - n_ ) ); // Basic check
877 
878  var_perm->resize(n_);
879  std::copy(
880  var_perm_full.begin()
881  ,var_perm_full.begin() + rank_fixed_removed
882  ,var_perm->begin()
883  );
884  std::copy(
885  var_perm_full.begin() + rank_full
886  ,var_perm_full.begin() + rank_full + (n_-rank_fixed_removed)
887  ,var_perm->begin() + rank_fixed_removed
888  );
889  *rank = rank_fixed_removed;
890  return true;
891  }
892  else {
893 
894  // try to find the file giving the basis...
895  char ind[17];
896  sprintf(ind, "%d", basis_selection_num_);
897  std::string fname = "basis_";
898  fname += ind;
899  fname += ".sel";
901 
902  std::ifstream basis_file(fname.c_str());
903  if (basis_file)
904  {
905  // try to read the basis file
906  std::string tags;
907 
908  int n;
909  basis_file >> tags;
911  tags != "n", std::logic_error
912  ,"Incorrect basis file format - \"n\" expected, \"" << tags << "\" found.");
913  basis_file >> n;
915  n <= 0, std::logic_error
916  , "Incorrect basis file format - n > 0 \"" << n << "\" found.");
917 
918  int m;
919  basis_file >> tags;
921  tags != "m", std::logic_error
922  ,"Incorrect basis file format - \"m\" expected, \"" << tags << "\" found.");
923  basis_file >> m;
925  m > n , std::logic_error
926  ,"Incorrect basis file format - 0 < m <= n expected, \"" << m << "\" found.");
927 
928  int r;
929  basis_file >> tags;
931  tags != "rank", std::logic_error
932  ,"Incorrect basis file format - \"rank\" expected, \"" << tags << "\" found.");
933  basis_file >> r;
935  r > m, std::logic_error
936  ,"Incorrect basis file format - 0 < rank <= m expected, \"" << r << "\" found.");
937  if (rank)
938  { *rank = r; }
939 
940  // var_permutation
941  basis_file >> tags;
943  tags != "var_perm", std::logic_error
944  ,"Incorrect basis file format -\"var_perm\" expected, \"" << tags << "\" found.");
945  var_perm->resize(n);
946  {for (int i=0; i < n; i++)
947  {
948  int var_index;
949  basis_file >> var_index;
951  var_index < 1 || var_index > n, std::logic_error
952  ,"Incorrect basis file format for var_perm: 1 <= indice <= n expected, \"" << n << "\" found.");
953  (*var_perm)[i] = var_index;
954  }}
955 
956  // eqn_permutation
957  basis_file >> tags;
959  tags != "equ_perm", std::logic_error
960  ,"Incorrect basis file format -\"equ_perm\" expected, \"" << tags << "\" found.");
961  equ_perm->resize(r);
962  {for (int i=0; i < r; i++)
963  {
964  int equ_index;
965  basis_file >> equ_index;
967  equ_index < 1 || equ_index > m, std::logic_error
968  ,"Incorrect basis file format for equ_perm: 1 <= indice <= m expected, \"" << m << "\" found.");
969  (*equ_perm)[i] = equ_index;
970  }}
971 
972  return true;
973  }
974  }
975 
976  return false;
977 }
978 
980  const IVector& var_perm, const IVector& equ_perm, size_type rank
981  )
982 {
983  namespace mmp = MemMngPack;
984 
985  // Assert that this is a valid basis and set the internal basis. Also repivot 'xinit',
986  // 'xl', and 'xu'.
987 
989 
990  // Assert the preconditions
992  var_perm.size() != n_ || equ_perm.size() != m_full_, std::length_error
993  ,"NLPSerialPreprocess::set_basis(): The sizes "
994  "of the permutation vectors does not match the size of the NLP" );
996  rank > m_full_, InvalidBasis
997  ,"NLPSerialPreprocess::set_basis(): The rank "
998  "of the basis can not be greater that the number of constraints" );
999 
1000  // Set the basis
1001  r_ = rank;
1002  if( &var_perm_ != &var_perm )
1003  var_perm_ = var_perm;
1004  if( &equ_perm_ != &equ_perm )
1005  equ_perm_ = equ_perm;
1007 
1008  var_from_full( xinit_full_().begin(), xinit_.set_vec().begin() );
1009  if(num_bounded_x_) {
1010  var_from_full( xl_full_().begin(), xl_.set_vec().begin() );
1011  var_from_full( xu_full_().begin(), xu_.set_vec().begin() );
1013  }
1014  else {
1015  xl_ = -NLP::infinite_bound();
1016  xu_ = +NLP::infinite_bound();
1017  }
1018  P_var_.initialize(Teuchos::rcp(new IVector(var_perm)),Teuchos::null);
1019  P_equ_.initialize(Teuchos::rcp(new IVector(equ_perm)),Teuchos::null);
1020 }
1021 
1023 {
1026  ,"There are no bounds on the variables for this NLP" );
1027 }
1028 
1030 {
1032 }
1033 
1034 } // end namespace NLPInterfacePack
void initialize(bool test_setup)
Initialize the NLP for its first use.
void set_basis(const Permutation &P_var, const Range1D &var_dep, const Permutation *P_equ, const Range1D *equ_decomp)
void var_to_full(DVectorSlice::const_iterator vec, DVectorSlice::iterator vec_full) const
bool get_next_basis_remove_fixed(IVector *var_perm, IVector *equ_perm, size_type *rank)
virtual bool imp_nlp_has_changed() const
Return if the definition of the NLP has changed since the last call to initialize() ...
AbstractLinAlgPack::size_type size_type
VectorMutable * c
Pointer to constraints residual c (Will be NULL if not set)
Struct for gradient (objective), objective and constriants (pointers)
void imp_calc_h_breve(const Vector &x, bool newx, const ZeroOrderInfo &zero_order_info_breve) const
virtual size_type imp_m_orig() const =0
Return the number of general equality constraints in the original problem.
Abstract interface for immutable, finite dimensional, coordinate vectors {abstract}.
const IVector & equ_perm() const
Permutes from the original constriant ordering to the current basis selection.
void var_from_full(DVectorSlice::const_iterator vec_full, DVectorSlice::iterator vec) const
void Vt_S(VectorMutable *v_lhs, const value_type &alpha)
v_lhs *= alpha
void inv_perm(const IVector &perm, IVector *inv_perm)
VectorMutable * Gf
Pointer to gradient of objective function Gf (may be NULL if not set)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void Vp_StV(VectorMutable *v_lhs, const value_type &alpha, const Vector &v_rhs)
v_lhs = alpha * v_rhs + v_lhs
virtual size_type imp_n_orig() const =0
Return the number of variables in the original problem (including those fixed by bounds) ...
int resize(OrdinalType length_in)
T * get() const
virtual void imp_calc_c_orig(const DVectorSlice &x_full, bool newx, const ZeroOrderInfoSerial &zero_order_info) const =0
Calculate the vector for all of the general equality constaints in the original NLP.
const IVector & var_perm() const
Permutes from the compated variable vector (removing fixed variables) to the current basis selection...
void assert_vs_sizes(size_type size1, size_type size2)
DVectorSlice x_full() const
Give reference to current x_full.
virtual void imp_calc_Gf_orig(const DVectorSlice &x_full, bool newx, const ObjGradInfoSerial &obj_grad_info) const =0
Calculate the vector for the gradient of the objective in the original NLP.
virtual void imp_calc_h_orig(const DVectorSlice &x_full, bool newx, const ZeroOrderInfoSerial &zero_order_info) const =0
Calculate the vector for all of the general inequality constaints in the original NLP...
void report_final_solution(const Vector &x, const Vector *lambda, const Vector *nu, bool is_optimal)
Overridden to permute the variables back into an order that is natural to the subclass.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void get_basis(Permutation *P_var, Range1D *var_dep, Permutation *P_equ, Range1D *equ_decomp) const
virtual const DVectorSlice imp_xinit_orig() const =0
Return the original initial point (size imp_n_orig()).
void copy(const f_int &N, const f_dbl_prec *X, const f_int &INCX, f_dbl_prec *Y, const f_int &INCY)
void imp_calc_f(const Vector &x, bool newx, const ZeroOrderInfo &zero_order_info) const
virtual const DVectorSlice imp_xu_orig() const =0
Return the original upper variable bounds (size imp_n_orig()).
T_To & dyn_cast(T_From &from)
void imp_calc_c(const Vector &x, bool newx, const ZeroOrderInfo &zero_order_info) const
virtual bool imp_has_var_bounds() const =0
Return if the NLP has bounds.
void assert_initialized() const
Assert if we have been initizlized (throws UnInitialized)
Struct for objective and constriants (pointer).
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)
To be overridden by subclasses to report the final solution in the original ordering natural to the s...
virtual const DVectorSlice imp_hu_orig() const =0
Return the original upper general inequality bounds (size imp_mI_orig()).
virtual void imp_calc_f_orig(const DVectorSlice &x_full, bool newx, const ZeroOrderInfoSerial &zero_order_info) const =0
Calculate the objective function for the original NLP.
Thrown if any member functions are called before initialize() has been called.
virtual size_type imp_mI_orig() const =0
Return the number of general inequality constraints in the original problem.
void identity_perm(IVector *perm)
void imp_calc_c_breve(const Vector &x, bool newx, const ZeroOrderInfo &zero_order_info_breve) const
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
void assert_and_set_basis(const IVector &var_perm, const IVector &equ_perm, size_type rank)
void Vt_S(DVectorSlice *vs_lhs, value_type alpha)
vs_lhs *= alpha (BLAS xSCAL) (*** Note that alpha == 0.0 is handeled as vs_lhs = 0.0)
AbstractLinAlgPack::value_type value_type
Abstract interface for mutable coordinate vectors {abstract}.
value_type * f
Pointer to objective function f (Will be NULL if not set)
void equ_from_full(const DVectorSlice &c_orig, const DVectorSlice &h_orig, const DVectorSlice &s_orig, DVectorSlice *c_full) const
void force_in_bounds(const Vector &xl, const Vector &xu, VectorMutable *x)
Force a vector in bounds.
virtual bool imp_get_next_basis(IVector *var_perm_full, IVector *equ_perm_full, size_type *rank_full, size_type *rank)
Return the next basis selection (default returns false).
void get_init_lagrange_mult(VectorMutable *lambda, VectorMutable *nu) const
virtual const DVectorSlice imp_xl_orig() const =0
Return the original lower variable bounds (size imp_n_orig()).
Thrown some bounds do not existe.
Thrown from initialize() if some logical error occured.
void set_x_full(const DVectorSlice &x, bool newx, DVectorSlice *x_full) const
Set the full x vector if newx == true
static value_type fixed_var_mult()
Gives the value of a Lagrange multipler for a fixed variable bound .that has been preprocessed out of...
static value_type infinite_bound()
Value for an infinite bound.
VectorMutable * h
Pointer to inequality constraints h (Will be NULL if not set)
int total_count() const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
virtual const DVectorSlice imp_hl_orig() const =0
Return the original lower general inequality bounds (size imp_mI_orig()).
void imp_calc_Gf(const Vector &x, bool newx, const ObjGradInfo &obj_grad_info) const
bool get_next_basis(Permutation *P_var, Range1D *var_dep, Permutation *P_equ, Range1D *equ_decomp)