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