MoochoPack_PBFGS_helpers.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 <ostream>
00030 #include <iomanip>
00031 
00032 #include "MoochoPack_PBFGS_helpers.hpp"
00033 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
00034 #include "AbstractLinAlgPack/src/AbstractLinAlgPack_SpVectorClass.hpp"
00035 #include "AbstractLinAlgPack_SpVectorOp.hpp"
00036 #include "MiWorkspacePack.h"
00037 
00038 namespace MoochoPack {
00039 
00040 bool PBFGSPack::act_set_calmed_down( 
00041   const ActSetStats         &stats
00042   ,const value_type         act_set_frac_proj_start
00043   ,EJournalOutputLevel      olevel
00044   ,std::ostream             &out
00045   )
00046 {
00047   typedef ActSetStats ASS;
00048   const size_type
00049     num_active_indep = stats.num_active_indep(),
00050     num_adds_indep   = stats.num_adds_indep(),
00051     num_drops_indep  = stats.num_drops_indep();
00052   const value_type
00053     frac_same
00054     = ( num_adds_indep == ASS::NOT_KNOWN || num_drops_indep == ASS::NOT_KNOWN || num_active_indep == 0
00055       ? 0.0
00056       : std::_MAX(((double)(num_active_indep)-num_adds_indep-num_drops_indep) / num_active_indep, 0.0 ) );
00057   const bool act_set_calmed_down = ( num_active_indep > 0 && frac_same >= act_set_frac_proj_start );
00058   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00059     out << "\nnum_active_indep = " << num_active_indep;
00060     if( num_active_indep ) {
00061       out << "\nmax(num_active_indep-num_adds_indep-num_drops_indep,0)/(num_active_indep) = "
00062         << "max("<<num_active_indep<<"-"<<num_adds_indep<<"-"<<num_drops_indep<<",0)/("<<num_active_indep<<") = "
00063         << frac_same;
00064       if( act_set_calmed_down )
00065         out << " >= ";
00066       else
00067         out << " < ";
00068       out << "act_set_frac_proj_start = " << act_set_frac_proj_start;
00069       if( act_set_calmed_down )
00070         out << "\nThe active set has calmed down enough\n";
00071       else
00072         out << "\nThe active set has not calmed down enough\n";
00073     }
00074   }
00075   return act_set_calmed_down;
00076 }
00077 
00078 void PBFGSPack::init_i_x_free_sRTsR_sRTyR(
00079   const SpVectorSlice        &nu_indep
00080   ,const DVectorSlice         &s
00081   ,const DVectorSlice         &y
00082   ,size_type                 *n_pz_R
00083   ,size_type                 i_x_free[]
00084   ,value_type                *sRTsR
00085   ,value_type                *sRTyR
00086   )
00087 {
00088   const size_type
00089     n_pz = nu_indep.size();
00090   SpVectorSlice::const_iterator
00091     nu_indep_itr   = nu_indep.begin(),
00092     nu_indep_end   = nu_indep.end();
00093   const SpVectorSlice::difference_type
00094     o = nu_indep.offset();
00095   *n_pz_R = 0;
00096   *sRTsR  = 0.0;
00097   *sRTyR  = 0.0;
00098   for( size_type i = 1; i <= n_pz; ++i ) {
00099     if( nu_indep_itr != nu_indep_end && nu_indep_itr->indice() + o == i ) {
00100       ++nu_indep_itr;
00101     }
00102     else {
00103       ++(*n_pz_R);
00104       *i_x_free++ = i;
00105       const value_type s_i = s(i);
00106       (*sRTsR) += s_i*s_i;
00107       (*sRTyR) += s_i*y(i);
00108     }
00109   }
00110 }
00111 
00112 void PBFGSPack::sort_fixed_max_cond_viol(
00113   const SpVectorSlice        &nu_indep
00114   ,const DVectorSlice         &s
00115   ,const DVectorSlice         &y
00116   ,const DVectorSlice         &B_XX
00117   ,const value_type          sRTBRRsR
00118   ,const value_type          sRTyR
00119   ,value_type                *sXTBXXsX
00120   ,value_type                *sXTyX
00121   ,size_type                 l_x_fixed_sorted[]
00122   )
00123 {
00124   using Teuchos::Workspace;
00125   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00126 
00127   const size_type
00128     n_pz = nu_indep.size();
00129   
00130   *sXTBXXsX = 0.0;
00131   *sXTyX    = 0.0;
00132 
00133   // Initial spare vector so we can sort this stuff
00134   typedef SpVector::element_type ele_t;
00135   Workspace<ele_t> sort_array(wss,nu_indep.nz());
00136   {
00137     SpVectorSlice::const_iterator
00138       nu_indep_itr = nu_indep.begin();
00139     ele_t
00140       *itr         = &sort_array[0];
00141     for( size_type l = 1 ; l <= nu_indep.nz(); ++l, ++itr, ++nu_indep_itr ) {
00142       const size_type i = nu_indep_itr->indice() + nu_indep.offset();
00143       const value_type
00144         s_i          = s(i),
00145         y_i          = y(i),
00146         B_ii         = B_XX[l-1],
00147         s_i_B_ii_s_i = s_i*B_ii*s_i,
00148         s_i_y_i      = s_i*y_i;
00149       *sXTBXXsX += s_i_B_ii_s_i;
00150       *sXTyX    += s_i_y_i;
00151       itr->initialize( l, s_i_B_ii_s_i/sRTBRRsR + ::fabs(s_i_y_i)/sRTyR );
00152     }
00153   }
00154   // Sort this sparse vector in decending order
00155   std::sort(
00156     &sort_array[0], &sort_array[0] + sort_array.size()
00157     , AbstractLinAlgPack::SortByDescendingAbsValue()
00158     );
00159   // Extract this ordering
00160   {
00161     for( size_type l = 0; l < nu_indep.nz(); ++l )
00162       l_x_fixed_sorted[l] = sort_array[l].indice();
00163   }
00164 }
00165 
00166 void PBFGSPack::choose_fixed_free(
00167   const value_type                       project_error_tol
00168   ,const value_type                      super_basic_mult_drop_tol
00169   ,const SpVectorSlice                   &nu_indep
00170   ,const DVectorSlice                     &s
00171   ,const DVectorSlice                     &y
00172   ,const DVectorSlice                     &B_XX
00173   ,const size_type                       l_x_fixed_sorted[]
00174   ,EJournalOutputLevel                   olevel
00175   ,std::ostream                          &out
00176   ,value_type                            *sRTBRRsR
00177   ,value_type                            *sRTyR
00178   ,value_type                            *sXTBXXsX
00179   ,value_type                            *sXTyX
00180   ,size_type                             *n_pz_X
00181   ,size_type                             *n_pz_R
00182   ,size_type                             i_x_free[]
00183   ,size_type                             i_x_fixed[]
00184   ,ConstrainedOptPack::EBounds  bnd_fixed[]
00185   )
00186 {
00187   using std::setw;
00188   using std::endl;
00189   using std::right;
00190   using AbstractLinAlgPack::norm_inf;
00191   namespace COP = ConstrainedOptPack;
00192   using Teuchos::Workspace;
00193   Teuchos::WorkspaceStore* wss = Teuchos::get_default_workspace_store().get();
00194 
00195   const size_type
00196     n_pz = nu_indep.size();
00197 
00198   // Store the status of the variables so that we can put sorted i_x_free[]
00199   // and i_x_fixed[] together at the end
00200   Workspace<long int>  i_x_status(wss,n_pz);  // free if > 0 , fixed if < 0 , error if 0
00201   std::fill_n( &i_x_status[0], n_pz, 0 );
00202   {for( size_type l = 0; l < (*n_pz_R); ++l ) {
00203     i_x_status[i_x_free[l]-1] = +1;
00204   }
00205   // Adjust i_x_free[] and i_x_fixed to meat the projection conditions
00206   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00207     out << "\nDetermining which fixed variables to put in superbasis and which to leave out (must have at least one in superbasis)...\n";
00208   }
00209   const value_type
00210     max_nu_indep = norm_inf(nu_indep);
00211   const bool
00212     all_fixed = n_pz == nu_indep.nz();
00213   if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ACTIVE_SET) ) {
00214     out << "\nmax{|nu_k(indep)|,i=r+1...n} = " << max_nu_indep            << std::endl
00215       << "super_basic_mult_drop_tol    = " << super_basic_mult_drop_tol << std::endl
00216       << "project_error_tol            = " << project_error_tol         << std::endl;
00217   }
00218   if( super_basic_mult_drop_tol > 1.0 ) {
00219     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00220       out << "super_basic_mult_drop_tol = " << super_basic_mult_drop_tol << " > 1"
00221         << "\nNo variables will be removed from the super basis!  (You might consider decreasing super_basic_mult_drop_tol < 1)\n";
00222     }
00223   }
00224   else {
00225     const int prec = out.precision();
00226     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ACTIVE_SET) ) {
00227       out << endl
00228         << right << setw(10)      << "i"
00229         << right << setw(prec+12) << "nu_indep(i)"
00230         << right << setw(1)       << " "
00231         << right << setw(prec+12) << "s(i)*B(ii)*s(i)"
00232         << right << setw(prec+12) << "s_R'*B_RR*s_R"
00233         << right << setw(prec+12) << "s_X'*B_XX*s_X"
00234         << right << setw(1)       << " "
00235         << right << setw(prec+12) << "s(i)*y(i)"
00236         << right << setw(prec+12) << "s_R'*y_R"
00237         << right << setw(prec+12) << "s_X'*y_X"
00238         << right << setw(1)       << " "
00239         << right << setw(14)      << "status"
00240         << endl
00241         << right << setw(10)      << "--------"
00242         << right << setw(prec+12) << "---------------"
00243         << right << setw(1)       << " "
00244         << right << setw(prec+12) << "---------------"
00245         << right << setw(prec+12) << "---------------"
00246         << right << setw(prec+12) << "---------------"
00247         << right << setw(1)       << " "
00248         << right << setw(prec+12) << "---------------"
00249         << right << setw(prec+12) << "---------------"
00250         << right << setw(prec+12) << "---------------"
00251         << right << setw(1)       << " "
00252         << right << setw(14)      << "------------"
00253         << endl;
00254     }
00255     // Loop through the fixed variables in decending order of the violation.
00256     bool kept_one = false;
00257     for( size_type k = 0; k < nu_indep.nz(); ++k ) {
00258       const size_type
00259         l = l_x_fixed_sorted[k];
00260       const SpVectorSlice::element_type
00261         &nu_i = *(nu_indep.begin() + (l-1));
00262       const size_type
00263         i = nu_i.indice() + nu_indep.offset();
00264       const value_type
00265         abs_val_nu   = ::fabs(nu_i.value()),
00266         rel_val_nu   = abs_val_nu / max_nu_indep,
00267         s_i          = s(i),
00268         y_i          = y(i),
00269         B_ii         = B_XX[l-1],
00270         s_i_B_ii_s_i = s_i*B_ii*s_i,
00271         s_i_y_i      = s_i*y_i;
00272       const bool
00273         nu_cond       =  rel_val_nu < super_basic_mult_drop_tol,
00274         sXTBXXsX_cond = (*sXTBXXsX) / (*sRTBRRsR) > project_error_tol,
00275         sXTyX_cond    = ::fabs(*sXTyX) / ::fabs(*sRTyR) > project_error_tol,
00276         keep = ( (all_fixed && abs_val_nu == max_nu_indep && !kept_one)
00277              || nu_cond || sXTBXXsX_cond || nu_cond );
00278       if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ACTIVE_SET) ) {
00279         out << right << setw(10)      << i
00280           << right << setw(prec+12) << nu_i.value()
00281           << right << setw(1)       << (nu_cond ? "*" : " ")
00282           << right << setw(prec+12) << s_i_B_ii_s_i
00283           << right << setw(prec+12) << (*sRTBRRsR)
00284           << right << setw(prec+12) << (*sXTBXXsX)
00285           << right << setw(1)       << (sXTBXXsX_cond ? "*" : " ")
00286           << right << setw(prec+12) << s_i_y_i
00287           << right << setw(prec+12) << (*sRTyR)
00288           << right << setw(prec+12) << (*sXTyX)
00289           << right << setw(1)       << (sXTyX_cond ? "*" : " ")
00290           << right << setw(14)      << (keep ? "superbasic" : "nonbasic")
00291           << endl;
00292       }
00293       if(keep) {
00294         kept_one = true;
00295         *sRTBRRsR += s_i_B_ii_s_i;
00296         *sXTBXXsX -= s_i_B_ii_s_i;
00297         *sRTyR    += s_i_y_i;
00298         *sXTyX    -= s_i_y_i;
00299         i_x_status[i-1] = +1;
00300         ++(*n_pz_R);
00301       }
00302       else {
00303         i_x_status[i-1] = -1;
00304         ++(*n_pz_X);
00305       }
00306     }}
00307     if( static_cast<int>(olevel) >= static_cast<int>(PRINT_ALGORITHM_STEPS) ) {
00308       out << "\nFinal selection: n_pz_X = " << (*n_pz_X) << " nonbasic variables and n_pz_R = " << (*n_pz_R) << " superbasic variables\n";
00309     }
00310   }
00311   // Set the final set of i_x_free[] and i_x_fixed[]
00312   SpVector::const_iterator
00313     nu_itr = nu_indep.begin(),
00314     nu_end = nu_indep.end();
00315   SpVector::difference_type
00316     nu_o = nu_indep.offset();
00317   size_type
00318     *i_x_free_itr  = i_x_free,
00319     *i_x_fixed_itr = i_x_fixed;
00320   ConstrainedOptPack::EBounds
00321     *bnd_fixed_itr = bnd_fixed;
00322   {for( size_type i = 1; i <= n_pz; ++i ) {
00323     long int status = i_x_status[i-1];
00324     TEST_FOR_EXCEPT( !( status ) ); // should not be zero!
00325     if( status > 0 ) {
00326       // A superbasic variable
00327       *i_x_free_itr++ = i;
00328     }
00329     else {
00330       // A nonbasic variable
00331       for( ; nu_itr->indice() + nu_o < i; ++nu_itr ); // Find the multiplier
00332       TEST_FOR_EXCEPT( !( nu_itr != nu_end ) );
00333       TEST_FOR_EXCEPT( !( nu_itr->indice() + nu_o == i  ) );
00334       *i_x_fixed_itr++ = i;
00335       *bnd_fixed_itr++ = ( nu_itr->value() > 0.0 ? COP::UPPER : COP::LOWER );
00336     }
00337   }}
00338   TEST_FOR_EXCEPT( !(  i_x_free_itr  - i_x_free  == *n_pz_R  ) );
00339   TEST_FOR_EXCEPT( !(  i_x_fixed_itr - i_x_fixed == *n_pz_X  ) );
00340   TEST_FOR_EXCEPT( !(  bnd_fixed_itr - bnd_fixed == *n_pz_X  ) );
00341 }
00342 
00343 } // end namespace MoochoPack

Generated on Tue Oct 20 12:51:47 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7