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