MOOCHO (Single Doxygen Collection) Version of the Day
MoochoPack_PBFGS_helpers.cpp
Go to the documentation of this file.
00001 #if 0
00002 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00007 //                  Copyright (2003) Sandia Corporation
00008 // 
00009 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 // license for use of this work by or on behalf of the U.S. Government.
00011 // 
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission.
00026 //
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //
00039 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00040 // 
00041 // ***********************************************************************
00042 // @HEADER
00043 
00044 #include <ostream>
00045 #include <iomanip>
00046 
00047 #include "MoochoPack_PBFGS_helpers.hpp"
00048 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
00049 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00050 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00051 #include "MiWorkspacePack.h"
00052 
00053 namespace MoochoPack {
00054 
00055 bool PBFGSPack::act_set_calmed_down( 
00056   const ActSetStats         &stats
00057   ,const value_type         act_set_frac_proj_start
00058   ,EJournalOutputLevel      olevel
00059   ,std::ostream             &out
00060   )
00061 {
00062   typedef ActSetStats ASS;
00063   const size_type
00064     num_active_indep = stats.num_active_indep(),
00065     num_adds_indep   = stats.num_adds_indep(),
00066     num_drops_indep  = stats.num_drops_indep();
00067   const value_type
00068     frac_same
00069     = ( num_adds_indep == ASS::NOT_KNOWN || num_drops_indep == ASS::NOT_KNOWN || num_active_indep == 0
00070       ? 0.0
00071       : std::_MAX(((double)(num_active_indep)-num_adds_indep-num_drops_indep) / num_active_indep, 0.0 ) );
00072   const bool act_set_calmed_down = ( num_active_indep > 0 && frac_same >= act_set_frac_proj_start );
00073   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00074     out << "\nnum_active_indep = " << num_active_indep;
00075     if( num_active_indep ) {
00076       out << "\nmax(num_active_indep-num_adds_indep-num_drops_indep,0)/(num_active_indep) = "
00077         << "max("<<num_active_indep<<"-"<<num_adds_indep<<"-"<<num_drops_indep<<",0)/("<<num_active_indep<<") = "
00078         << frac_same;
00079       if( act_set_calmed_down )
00080         out << " >= ";
00081       else
00082         out << " < ";
00083       out << "act_set_frac_proj_start = " << act_set_frac_proj_start;
00084       if( act_set_calmed_down )
00085         out << "\nThe active set has calmed down enough\n";
00086       else
00087         out << "\nThe active set has not calmed down enough\n";
00088     }
00089   }
00090   return act_set_calmed_down;
00091 }
00092 
00093 void PBFGSPack::init_i_x_free_sRTsR_sRTyR(
00094   const SpVectorSlice        &nu_indep
00095   ,const DVectorSlice         &s
00096   ,const DVectorSlice         &y
00097   ,size_type                 *n_pz_R
00098   ,size_type                 i_x_free[]
00099   ,value_type                *sRTsR
00100   ,value_type                *sRTyR
00101   )
00102 {
00103   const size_type
00104     n_pz = nu_indep.size();
00105   SpVectorSlice::const_iterator
00106     nu_indep_itr   = nu_indep.begin(),
00107     nu_indep_end   = nu_indep.end();
00108   const SpVectorSlice::difference_type
00109     o = nu_indep.offset();
00110   *n_pz_R = 0;
00111   *sRTsR  = 0.0;
00112   *sRTyR  = 0.0;
00113   for( size_type i = 1; i <= n_pz; ++i ) {
00114     if( nu_indep_itr != nu_indep_end && nu_indep_itr->indice() + o == i ) {
00115       ++nu_indep_itr;
00116     }
00117     else {
00118       ++(*n_pz_R);
00119       *i_x_free++ = i;
00120       const value_type s_i = s(i);
00121       (*sRTsR) += s_i*s_i;
00122       (*sRTyR) += s_i*y(i);
00123     }
00124   }
00125 }
00126 
00127 void PBFGSPack::sort_fixed_max_cond_viol(
00128   const SpVectorSlice        &nu_indep
00129   ,const DVectorSlice         &s
00130   ,const DVectorSlice         &y
00131   ,const DVectorSlice         &B_XX
00132   ,const value_type          sRTBRRsR
00133   ,const value_type          sRTyR
00134   ,value_type                *sXTBXXsX
00135   ,value_type                *sXTyX
00136   ,size_type                 l_x_fixed_sorted[]
00137   )
00138 {
00139   using Teuchos::Workspace;
00140   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00141 
00142   const size_type
00143     n_pz = nu_indep.size();
00144   
00145   *sXTBXXsX = 0.0;
00146   *sXTyX    = 0.0;
00147 
00148   // Initial spare vector so we can sort this stuff
00149   typedef SpVector::element_type ele_t;
00150   Workspace<ele_t> sort_array(wss,nu_indep.nz());
00151   {
00152     SpVectorSlice::const_iterator
00153       nu_indep_itr = nu_indep.begin();
00154     ele_t
00155       *itr         = &sort_array[0];
00156     for( size_type l = 1 ; l <= nu_indep.nz(); ++l, ++itr, ++nu_indep_itr ) {
00157       const size_type i = nu_indep_itr->indice() + nu_indep.offset();
00158       const value_type
00159         s_i          = s(i),
00160         y_i          = y(i),
00161         B_ii         = B_XX[l-1],
00162         s_i_B_ii_s_i = s_i*B_ii*s_i,
00163         s_i_y_i      = s_i*y_i;
00164       *sXTBXXsX += s_i_B_ii_s_i;
00165       *sXTyX    += s_i_y_i;
00166       itr->initialize( l, s_i_B_ii_s_i/sRTBRRsR + ::fabs(s_i_y_i)/sRTyR );
00167     }
00168   }
00169   // Sort this sparse vector in decending order
00170   std::sort(
00171     &sort_array[0], &sort_array[0] + sort_array.size()
00172     , AbstractLinAlgPack::SortByDescendingAbsValue()
00173     );
00174   // Extract this ordering
00175   {
00176     for( size_type l = 0; l < nu_indep.nz(); ++l )
00177       l_x_fixed_sorted[l] = sort_array[l].indice();
00178   }
00179 }
00180 
00181 void PBFGSPack::choose_fixed_free(
00182   const value_type                       project_error_tol
00183   ,const value_type                      super_basic_mult_drop_tol
00184   ,const SpVectorSlice                   &nu_indep
00185   ,const DVectorSlice                     &s
00186   ,const DVectorSlice                     &y
00187   ,const DVectorSlice                     &B_XX
00188   ,const size_type                       l_x_fixed_sorted[]
00189   ,EJournalOutputLevel                   olevel
00190   ,std::ostream                          &out
00191   ,value_type                            *sRTBRRsR
00192   ,value_type                            *sRTyR
00193   ,value_type                            *sXTBXXsX
00194   ,value_type                            *sXTyX
00195   ,size_type                             *n_pz_X
00196   ,size_type                             *n_pz_R
00197   ,size_type                             i_x_free[]
00198   ,size_type                             i_x_fixed[]
00199   ,ConstrainedOptPack::EBounds  bnd_fixed[]
00200   )
00201 {
00202   using std::setw;
00203   using std::endl;
00204   using std::right;
00205   using AbstractLinAlgPack::norm_inf;
00206   namespace COP = ConstrainedOptPack;
00207   using Teuchos::Workspace;
00208   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00209 
00210   const size_type
00211     n_pz = nu_indep.size();
00212 
00213   // Store the status of the variables so that we can put sorted i_x_free[]
00214   // and i_x_fixed[] together at the end
00215   Workspace<long int>  i_x_status(wss,n_pz);  // free if > 0 , fixed if < 0 , error if 0
00216   std::fill_n( &i_x_status[0], n_pz, 0 );
00217   {for( size_type l = 0; l < (*n_pz_R); ++l ) {
00218     i_x_status[i_x_free[l]-1] = +1;
00219   }
00220   // Adjust i_x_free[] and i_x_fixed to meat the projection conditions
00221   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00222     out << "\nDetermining which fixed variables to put in superbasis and which to leave out (must have at least one in superbasis)...\n";
00223   }
00224   const value_type
00225     max_nu_indep = norm_inf(nu_indep);
00226   const bool
00227     all_fixed = n_pz == nu_indep.nz();
00228   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ACTIVE_SET) ) {
00229     out << "\nmax{|nu_k(indep)|,i=r+1...n} = " << max_nu_indep            << std::endl
00230       << "super_basic_mult_drop_tol    = " << super_basic_mult_drop_tol << std::endl
00231       << "project_error_tol            = " << project_error_tol         << std::endl;
00232   }
00233   if( super_basic_mult_drop_tol > 1.0 ) {
00234     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00235       out << "super_basic_mult_drop_tol = " << super_basic_mult_drop_tol << " > 1"
00236         << "\nNo variables will be removed from the super basis!  (You might consider decreasing super_basic_mult_drop_tol < 1)\n";
00237     }
00238   }
00239   else {
00240     const int prec = out.precision();
00241     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ACTIVE_SET) ) {
00242       out << endl
00243         << right << setw(10)      << "i"
00244         << right << setw(prec+12) << "nu_indep(i)"
00245         << right << setw(1)       << " "
00246         << right << setw(prec+12) << "s(i)*B(ii)*s(i)"
00247         << right << setw(prec+12) << "s_R'*B_RR*s_R"
00248         << right << setw(prec+12) << "s_X'*B_XX*s_X"
00249         << right << setw(1)       << " "
00250         << right << setw(prec+12) << "s(i)*y(i)"
00251         << right << setw(prec+12) << "s_R'*y_R"
00252         << right << setw(prec+12) << "s_X'*y_X"
00253         << right << setw(1)       << " "
00254         << right << setw(14)      << "status"
00255         << endl
00256         << right << setw(10)      << "--------"
00257         << right << setw(prec+12) << "---------------"
00258         << right << setw(1)       << " "
00259         << right << setw(prec+12) << "---------------"
00260         << right << setw(prec+12) << "---------------"
00261         << right << setw(prec+12) << "---------------"
00262         << right << setw(1)       << " "
00263         << right << setw(prec+12) << "---------------"
00264         << right << setw(prec+12) << "---------------"
00265         << right << setw(prec+12) << "---------------"
00266         << right << setw(1)       << " "
00267         << right << setw(14)      << "------------"
00268         << endl;
00269     }
00270     // Loop through the fixed variables in decending order of the violation.
00271     bool kept_one = false;
00272     for( size_type k = 0; k < nu_indep.nz(); ++k ) {
00273       const size_type
00274         l = l_x_fixed_sorted[k];
00275       const SpVectorSlice::element_type
00276         &nu_i = *(nu_indep.begin() + (l-1));
00277       const size_type
00278         i = nu_i.indice() + nu_indep.offset();
00279       const value_type
00280         abs_val_nu   = ::fabs(nu_i.value()),
00281         rel_val_nu   = abs_val_nu / max_nu_indep,
00282         s_i          = s(i),
00283         y_i          = y(i),
00284         B_ii         = B_XX[l-1],
00285         s_i_B_ii_s_i = s_i*B_ii*s_i,
00286         s_i_y_i      = s_i*y_i;
00287       const bool
00288         nu_cond       =  rel_val_nu < super_basic_mult_drop_tol,
00289         sXTBXXsX_cond = (*sXTBXXsX) / (*sRTBRRsR) > project_error_tol,
00290         sXTyX_cond    = ::fabs(*sXTyX) / ::fabs(*sRTyR) > project_error_tol,
00291         keep = ( (all_fixed && abs_val_nu == max_nu_indep && !kept_one)
00292              || nu_cond || sXTBXXsX_cond || nu_cond );
00293       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ACTIVE_SET) ) {
00294         out << right << setw(10)      << i
00295           << right << setw(prec+12) << nu_i.value()
00296           << right << setw(1)       << (nu_cond ? "*" : " ")
00297           << right << setw(prec+12) << s_i_B_ii_s_i
00298           << right << setw(prec+12) << (*sRTBRRsR)
00299           << right << setw(prec+12) << (*sXTBXXsX)
00300           << right << setw(1)       << (sXTBXXsX_cond ? "*" : " ")
00301           << right << setw(prec+12) << s_i_y_i
00302           << right << setw(prec+12) << (*sRTyR)
00303           << right << setw(prec+12) << (*sXTyX)
00304           << right << setw(1)       << (sXTyX_cond ? "*" : " ")
00305           << right << setw(14)      << (keep ? "superbasic" : "nonbasic")
00306           << endl;
00307       }
00308       if(keep) {
00309         kept_one = true;
00310         *sRTBRRsR += s_i_B_ii_s_i;
00311         *sXTBXXsX -= s_i_B_ii_s_i;
00312         *sRTyR    += s_i_y_i;
00313         *sXTyX    -= s_i_y_i;
00314         i_x_status[i-1] = +1;
00315         ++(*n_pz_R);
00316       }
00317       else {
00318         i_x_status[i-1] = -1;
00319         ++(*n_pz_X);
00320       }
00321     }}
00322     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00323       out << "\nFinal selection: n_pz_X = " << (*n_pz_X) << " nonbasic variables and n_pz_R = " << (*n_pz_R) << " superbasic variables\n";
00324     }
00325   }
00326   // Set the final set of i_x_free[] and i_x_fixed[]
00327   SpVector::const_iterator
00328     nu_itr = nu_indep.begin(),
00329     nu_end = nu_indep.end();
00330   SpVector::difference_type
00331     nu_o = nu_indep.offset();
00332   size_type
00333     *i_x_free_itr  = i_x_free,
00334     *i_x_fixed_itr = i_x_fixed;
00335   ConstrainedOptPack::EBounds
00336     *bnd_fixed_itr = bnd_fixed;
00337   {for( size_type i = 1; i <= n_pz; ++i ) {
00338     long int status = i_x_status[i-1];
00339     TEUCHOS_TEST_FOR_EXCEPT( !( status ) ); // should not be zero!
00340     if( status > 0 ) {
00341       // A superbasic variable
00342       *i_x_free_itr++ = i;
00343     }
00344     else {
00345       // A nonbasic variable
00346       for( ; nu_itr->indice() + nu_o < i; ++nu_itr ); // Find the multiplier
00347       TEUCHOS_TEST_FOR_EXCEPT( !( nu_itr != nu_end ) );
00348       TEUCHOS_TEST_FOR_EXCEPT( !( nu_itr->indice() + nu_o == i  ) );
00349       *i_x_fixed_itr++ = i;
00350       *bnd_fixed_itr++ = ( nu_itr->value() > 0.0 ? COP::UPPER : COP::LOWER );
00351     }
00352   }}
00353   TEUCHOS_TEST_FOR_EXCEPT( !(  i_x_free_itr  - i_x_free  == *n_pz_R  ) );
00354   TEUCHOS_TEST_FOR_EXCEPT( !(  i_x_fixed_itr - i_x_fixed == *n_pz_X  ) );
00355   TEUCHOS_TEST_FOR_EXCEPT( !(  bnd_fixed_itr - bnd_fixed == *n_pz_X  ) );
00356 }
00357 
00358 } // end namespace MoochoPack
00359 
00360 #endif // 0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines