NLPInterfacePack_ExampleNLPBanded.cpp

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 <math.h>
00030 
00031 #include "NLPInterfacePack_ExampleNLPBanded.hpp"
00032 #include "DenseLinAlgPack_PermVecMat.hpp"
00033 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
00034 #include "Teuchos_TestForException.hpp"
00035 
00036 namespace NLPInterfacePack {
00037 
00038 // Constructors / initializers
00039 
00040 ExampleNLPBanded::ExampleNLPBanded(
00041   size_type     nD
00042   ,size_type    nI
00043   ,size_type    bw
00044   ,size_type    mU
00045   ,size_type    mI
00046   ,value_type   xo
00047   ,value_type   xDl
00048   ,value_type   xDu
00049   ,value_type   xIl
00050   ,value_type   xIu
00051   ,value_type   hl
00052   ,value_type   hu
00053   ,bool         nlp_selects_basis
00054   ,value_type   diag_scal
00055   ,value_type   diag_vary
00056   ,bool         sym_basis
00057   ,value_type   f_offset
00058   ,value_type   co
00059   ,bool         ignore_constraints
00060   )
00061   :is_initialized_(false)
00062   ,nlp_selects_basis_(nlp_selects_basis)
00063   ,basis_selection_was_given_(false)
00064   ,has_var_bounds_(false)
00065   ,f_offset_(f_offset)
00066   ,nD_(nD)
00067   ,nI_(nI)
00068   ,bw_(bw)
00069   ,mU_(mU)
00070   ,mI_(mI)
00071   ,ignore_constraints_(ignore_constraints)
00072   ,diag_scal_(diag_scal)
00073   ,diag_vary_(diag_vary)
00074   ,fu_( sym_basis ? 3 : 6 )
00075 {
00076 #ifdef TEUCHOS_DEBUG  
00077   const char msg_err_head[] = "ExampleNLPBanded::ExampleNLPBanded(...) : Error";
00078   TEST_FOR_EXCEPTION(
00079     nD <= 0, std::invalid_argument
00080     ,msg_err_head<<"!" );
00081   TEST_FOR_EXCEPTION(
00082     nI <= 0 || nD < nI, std::invalid_argument
00083     ,msg_err_head<<"!" );
00084   TEST_FOR_EXCEPTION(
00085     bw < 1 || nD < bw, std::invalid_argument
00086     ,msg_err_head<<"!" );
00087   TEST_FOR_EXCEPTION(
00088     mU < 0, std::invalid_argument
00089     ,msg_err_head<<"!" );
00090   TEST_FOR_EXCEPTION(
00091     mI < 0, std::invalid_argument
00092     ,msg_err_head<<"!" );
00093   TEST_FOR_EXCEPTION(
00094     mU != 0, std::invalid_argument
00095     ,msg_err_head<<", can't handle undecomposed equalities yet!" );
00096 #endif
00097   Gc_orig_nz_ = nD_ * ( 1 + 2*(bw_-1) + nI_ ); // Overestimate, ToDo: compute the exact value!
00098   Gh_orig_nz_ = 2*mI_;
00099   //
00100   xinit_orig_.resize(nD_ + nI_);
00101   xl_orig_.resize(xinit_orig_.dim());
00102   xu_orig_.resize(xinit_orig_.dim());
00103   hl_orig_.resize(mI_);
00104   hu_orig_.resize(mI_);
00105   co_orig_.resize(nD_ + mU_);
00106   //
00107   xinit_orig_             = xo;
00108   xl_orig_(1,nD_)         = xDl;
00109   xu_orig_(1,nD_)         = xDu;
00110   xl_orig_(nD_+1,nD_+nI_) = xIl;
00111   xu_orig_(nD_+1,nD_+nI_) = xIu;
00112   if( mI_ ) {
00113     hl_orig_    = hl;
00114     hu_orig_    = hu;
00115   }
00116   co_orig_ = co;
00117   //
00118   const value_type inf = NLP::infinite_bound();
00119   if( xDl > -inf || xDu < +inf || xIl > -inf || xIu < +inf )
00120     has_var_bounds_ = true;
00121   else
00122     has_var_bounds_ = false;
00123 }
00124 
00125 // Overridden public members from NLP
00126 
00127 void ExampleNLPBanded::initialize(bool test_setup)
00128 {
00129   if(is_initialized_) {
00130     NLPSerialPreprocessExplJac::initialize(test_setup);
00131     return;
00132   }
00133   // Nothing to initialize?
00134   NLPSerialPreprocessExplJac::initialize(test_setup);
00135   is_initialized_ = true;
00136 }
00137 
00138 bool ExampleNLPBanded::is_initialized() const
00139 {
00140   return is_initialized_;
00141 }
00142 
00143 value_type ExampleNLPBanded::max_var_bounds_viol() const
00144 {
00145   return +1e+20; // Functions defined everywhere!
00146 }
00147 
00148 // Overridden protected methods from NLPSerialPreprocess
00149 
00150 bool ExampleNLPBanded::imp_nlp_has_changed() const
00151 {
00152   return !is_initialized_;
00153 }
00154 
00155 size_type ExampleNLPBanded::imp_n_orig() const
00156 {
00157   return nD_ + nI_;
00158 }
00159 
00160 size_type ExampleNLPBanded::imp_m_orig() const
00161 {
00162   if(ignore_constraints_) return 0;
00163   return nD_ + mU_;
00164 }
00165 
00166 size_type ExampleNLPBanded::imp_mI_orig() const
00167 {
00168   if(ignore_constraints_) return 0;
00169   return mI_;
00170 }
00171 
00172 const DVectorSlice ExampleNLPBanded::imp_xinit_orig() const
00173 {
00174   return xinit_orig_();
00175 }
00176 
00177 bool ExampleNLPBanded::imp_has_var_bounds() const
00178 {
00179   return has_var_bounds_;
00180 }
00181 
00182 const DVectorSlice ExampleNLPBanded::imp_xl_orig() const
00183 {
00184   return xl_orig_();
00185 }
00186 
00187 const DVectorSlice ExampleNLPBanded::imp_xu_orig() const
00188 {
00189   return xu_orig_();
00190 }
00191 
00192 const DVectorSlice ExampleNLPBanded::imp_hl_orig() const
00193 {
00194   return hl_orig_();
00195 }
00196 
00197 const DVectorSlice ExampleNLPBanded::imp_hu_orig() const
00198 {
00199   return hu_orig_();
00200 }
00201 
00202 void ExampleNLPBanded::imp_calc_f_orig(
00203   const DVectorSlice            &x_full
00204   ,bool                        newx
00205   ,const ZeroOrderInfoSerial   &zero_order_info
00206   ) const
00207 {
00208   inform_new_point(newx);
00209   const DVectorSlice x_orig = x_full(1,imp_n_orig());
00210   *zero_order_info.f = f_offset_ + ( 1.0 / 2.0 ) * DenseLinAlgPack::dot( x_orig, x_orig );
00211 }
00212 
00213 void ExampleNLPBanded::imp_calc_c_orig(
00214   const DVectorSlice            &x_full
00215   ,bool                        newx
00216   ,const ZeroOrderInfoSerial   &zero_order_info
00217   ) const
00218 {
00219   inform_new_point(newx);
00220   if(c_orig_updated_)
00221     return; // c(x) is already computed in *zero_order_info.c
00222   TEST_FOR_EXCEPT( !( zero_order_info.c ) );
00223   DVector
00224     &c  = *zero_order_info.c;
00225   const size_type
00226     num_I_per_D = nD_ / nI_,  // Integer division (rounds down)
00227     I_remainder = nD_ % nI_;
00228   size_type j = 0;
00229   const value_type
00230     ds_alpha = nD_ > 1 ? diag_scal_ * (diag_vary_ - 1.0)/(nD_ - 1.0) : 0.0,
00231     ds_beta = diag_scal_;
00232   for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
00233     const size_type num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
00234     for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
00235       ++j;
00236       const size_type
00237         klu = ( j - bw_     >= 0   ? bw_-1 : j-1   ),
00238         kuu = ( j + bw_ - 1 <= nD_ ? bw_-1 : nD_-j );
00239       const value_type
00240         ds_j = ds_alpha * (j-1) + ds_beta;
00241       value_type
00242         &c_j = (c(j) = ds_j * x_full(j));
00243       {for( size_type k = 1; k <= klu; ++k ) {
00244         c_j -= (3.0 / k) * x_full(j-k);
00245       }}
00246       {for( size_type k = 1; k <= kuu; ++k ) {
00247         c_j -= (fu_ / k) * x_full(j+k);
00248       }}
00249       const value_type
00250         term = x_full(nD_ + q_i) + 1;
00251       c_j *= (term * term);
00252       c_j += co_orig_(j);
00253     }
00254   }
00255   c_orig_updated_ = true;
00256 }
00257 
00258 void ExampleNLPBanded::imp_calc_h_orig(
00259   const DVectorSlice            &x_full
00260   ,bool                        newx
00261   ,const ZeroOrderInfoSerial   &zero_order_info
00262   ) const
00263 {
00264   inform_new_point(newx);
00265   DVector
00266     &h  = *zero_order_info.h;
00267   const size_type
00268     num_I_per_D = nD_ / nI_,  // Integer division (rounds down)
00269     I_remainder = nD_ % nI_;
00270   size_type jI = 0;
00271   for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
00272     const value_type  x_q = x_full(nD_ + q_i);
00273     const size_type   num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
00274     for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
00275       ++jI;
00276       if( jI > mI_ ) return;
00277       h(jI) = x_full(jI) - x_q;
00278     }
00279   }
00280  }
00281 
00282 void ExampleNLPBanded::imp_calc_Gf_orig(
00283   const DVectorSlice            &x_full
00284   ,bool                        newx
00285   ,const ObjGradInfoSerial     &obj_grad_info
00286   ) const
00287 {
00288   inform_new_point(newx);
00289   const Range1D var_orig(1,imp_n_orig());
00290   (*obj_grad_info.Gf)(var_orig) = x_full(var_orig);
00291 }
00292 
00293 bool ExampleNLPBanded::imp_get_next_basis(
00294   IVector      *var_perm_full
00295   ,IVector     *equ_perm_full
00296   ,size_type   *rank_full
00297   ,size_type   *rank
00298   )
00299 {
00300   if(basis_selection_was_given_)
00301     return false; // Already gave this basis selection.
00302   // Select the first nD variables as basis variables which gives
00303   // a nice banded matrix for the basis matrix C.
00304   // Also, if the general inequality constraints are begin
00305   // converted to equalities with slacks, make the slack variables
00306   // basic variables also (after the nD variables).
00307 #ifdef TEUCHOS_DEBUG
00308   TEST_FOR_EXCEPT( !( var_perm_full ) );
00309   TEST_FOR_EXCEPT( !( equ_perm_full ) );
00310   TEST_FOR_EXCEPT( !( rank ) );
00311 #endif
00312   const size_type    n_orig = nD_    + nI_;
00313   const size_type    n_full = n_orig + mI_;
00314   const size_type    m_full = nD_    + mI_;
00315   var_perm_full->resize(n_full);
00316   equ_perm_full->resize(m_full);
00317   if( mI_ ) {
00318     index_type k, i_perm = 1;;
00319     // basic variables
00320     for( k = 1; k <= nD_; ++k, ++i_perm )  // Put xD variables first
00321       (*var_perm_full)(i_perm) = k;
00322     for( k = 1; k <= mI_; ++k, ++i_perm )  // Put slacks after xD
00323       (*var_perm_full)(i_perm) = n_orig + k;
00324     // non-basic variables
00325     for( k = 1; k <= nI_; ++k, ++i_perm )  // Followed by nI
00326       (*var_perm_full)(i_perm) = nD_ + k;
00327   }
00328   else {
00329     DenseLinAlgPack::identity_perm( var_perm_full );
00330   }
00331   DenseLinAlgPack::identity_perm( equ_perm_full );
00332   // Count the number of fixed basic variables to reduce
00333   // the rank of the basis.
00334   index_type num_fixed = 0;
00335   for( index_type k = 1; k <= nD_; ++k ) {
00336     if( xl_orig_(k) == xu_orig_(k) )
00337       ++num_fixed;
00338   }
00339   if( mI_ ) {
00340     for( index_type k = 1; k <= mI_; ++k ) {
00341       if( hl_orig_(k) == hu_orig_(k) )
00342         ++num_fixed;
00343     }
00344   }
00345   // Set the rank
00346   *rank_full = m_full;
00347   *rank      = m_full - num_fixed;
00348   basis_selection_was_given_ = true;
00349   return true;
00350 }
00351 
00352 void ExampleNLPBanded::imp_report_orig_final_solution(
00353   const DVectorSlice      &x_orig
00354   ,const DVectorSlice     *lambda_orig
00355   ,const DVectorSlice     *lambdaI_orig
00356   ,const DVectorSlice     *nu_orig
00357   ,bool                   is_optimal
00358   ) const
00359 {
00360   // ToDo: Do something with the final soltuion?
00361 }
00362 
00363 bool ExampleNLPBanded::nlp_selects_basis() const
00364 {
00365   return nlp_selects_basis_;
00366 }
00367 
00368 // Overridden protected methods from NLPSerialPreprocessExplJac
00369 
00370 size_type ExampleNLPBanded::imp_Gc_nz_orig() const
00371 {
00372   return Gc_orig_nz_;
00373 }
00374 
00375 size_type ExampleNLPBanded::imp_Gh_nz_orig() const
00376 {
00377   return Gh_orig_nz_;
00378 }
00379 
00380 void ExampleNLPBanded::imp_calc_Gc_orig(
00381   const DVectorSlice& x_full, bool newx
00382   , const FirstOrderExplInfo& first_order_expl_info
00383   ) const
00384 {
00385   inform_new_point(newx);
00386   // Compute c(x) if not already (will compute in temp if needed)
00387   this->imp_calc_c_orig( x_full, newx, zero_order_orig_info() );
00388   TEST_FOR_EXCEPT( !( first_order_expl_info.c ) );
00389   DVector
00390     &c = *first_order_expl_info.c;
00391   // Get references/pointers to data for Gc to be computed/updated.
00392   index_type
00393     &Gc_nz = *first_order_expl_info.Gc_nz;
00394   value_type
00395     *Gc_val = &(*first_order_expl_info.Gc_val)[0];
00396   index_type
00397     *Gc_ivect = ( first_order_expl_info.Gc_ivect
00398             ? &(*first_order_expl_info.Gc_ivect)[0] : NULL ),
00399     *Gc_jvect = ( first_order_expl_info.Gc_jvect
00400             ? &(*first_order_expl_info.Gc_jvect)[0] : NULL );
00401   TEST_FOR_EXCEPT( !(  (Gc_ivect != NULL) == (Gc_jvect != NULL)  ) );
00402   // Set nonzeros for Gc (in sorted compressed column format w.r.t., i.e. grouped by constraints)
00403   const size_type
00404     num_I_per_D = nD_ / nI_,  // Integer division (rounds down)
00405     I_remainder = nD_ % nI_;
00406   Gc_nz = 0;
00407   size_type j = 0;
00408   const value_type
00409     ds_alpha = nD_ > 1 ? diag_scal_ * (diag_vary_ - 1.0)/(nD_ - 1.0) : 0.0,
00410     ds_beta = diag_scal_;
00411   for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
00412     const value_type
00413       x_q = x_full(nD_ + q_i),
00414       x_q_term = (x_q + 1) * (x_q + 1);
00415     const size_type num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
00416     for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
00417       ++j;
00418       const size_type
00419         klu = ( j - bw_     >= 0   ? bw_-1 : j-1   ),
00420         kuu = ( j + bw_ - 1 <= nD_ ? bw_-1 : nD_-j );
00421       const value_type
00422         ds_j = ds_alpha * (j-1) + ds_beta;
00423       //
00424       {for( index_type k = klu; k >= 1; --k ) {
00425         ++Gc_nz;
00426         *Gc_val++ = -3.0 / k * x_q_term;
00427         if(Gc_ivect) {
00428           *Gc_ivect++ = j - k;
00429           *Gc_jvect++ = j;
00430         }
00431       }}
00432       //
00433       ++Gc_nz;
00434       *Gc_val++ = ds_j * x_q_term;
00435       if(Gc_ivect) {
00436         *Gc_ivect++ = j;
00437         *Gc_jvect++ = j;
00438       }
00439       //
00440       {for( index_type k = 1; k <= kuu; ++k ) {
00441         ++Gc_nz;
00442         *Gc_val++ = -fu_ / k * x_q_term;
00443         if(Gc_ivect) {
00444           *Gc_ivect++ = j + k;
00445           *Gc_jvect++ = j;
00446         }
00447       }}
00448       //
00449       ++Gc_nz;
00450       *Gc_val++ = 2.0 * (c(j) - co_orig_(j)) / (x_q + 1);
00451       if(Gc_ivect) {
00452         *Gc_ivect++ = nD_ + q_i;
00453         *Gc_jvect++ = j;
00454       }
00455     }
00456   }
00457 }
00458 
00459 void ExampleNLPBanded::imp_calc_Gh_orig(
00460   const DVectorSlice& x_full, bool newx
00461   , const FirstOrderExplInfo& first_order_expl_info
00462   ) const
00463 {
00464   inform_new_point(newx);
00465   // Get references/pointers to data for Gh to be computed/updated.
00466   index_type
00467     &Gh_nz = *first_order_expl_info.Gh_nz;
00468   value_type
00469     *Gh_val = &(*first_order_expl_info.Gh_val)[0];
00470   index_type
00471     *Gh_ivect = ( first_order_expl_info.Gh_ivect
00472             ? &(*first_order_expl_info.Gh_ivect)[0] : NULL ),
00473     *Gh_jvect = ( first_order_expl_info.Gh_jvect
00474             ? &(*first_order_expl_info.Gh_jvect)[0] : NULL );
00475   TEST_FOR_EXCEPT( !(  (Gh_ivect != NULL) == (Gh_jvect != NULL)  ) );
00476   // Set nonzeros for Gh (in sorted compressed column format w.r.t., i.e. grouped by constraints)
00477   const size_type
00478     num_I_per_D = nD_ / nI_,  // Integer division (rounds down)
00479     I_remainder = nD_ % nI_;
00480   Gh_nz = 0;
00481   size_type jI = 0;
00482   for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
00483     const size_type   nD_q_i = nD_ + q_i;
00484     const size_type   num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
00485     for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
00486       ++jI;
00487       if( jI > mI_ ) return;
00488       // w.r.t. x(jI)
00489       ++Gh_nz;
00490       *Gh_val++ = 1.0;
00491       if(Gh_ivect) {
00492         *Gh_ivect++ = jI;
00493         *Gh_jvect++ = jI;
00494       }
00495       // w.r.t. x(nD+q(jI))
00496       ++Gh_nz;
00497       *Gh_val++ = -1.0;
00498       if(Gh_ivect) {
00499         *Gh_ivect++ = nD_q_i;
00500         *Gh_jvect++ = jI;
00501       }
00502     }
00503   }
00504 }
00505 
00506 // private
00507 
00508 void ExampleNLPBanded::assert_is_initialized() const
00509 {
00510   TEST_FOR_EXCEPT(true); //  ToDo: Implemenet!
00511 }
00512 
00513 void ExampleNLPBanded::inform_new_point(bool newx) const
00514 {
00515   if(newx) {
00516     c_orig_updated_ = false;
00517   }
00518 }
00519 
00520 } // end namespace NLPInterfacePack

Generated on Wed May 12 21:57:50 2010 for MOOCHO by  doxygen 1.4.7