angel  mercurial changeset:
include/sa_impl.hpp
Go to the documentation of this file.
00001 // $Id: sa_impl.hpp,v 1.4 2004/04/23 12:59:11 gottschling Exp $
00002 /*
00003 #############################################################
00004 # This file is part of angel released under the BSD license #
00005 # The full COPYRIGHT notice can be found in the top         #
00006 # level directory of the angel distribution                 #
00007 #############################################################
00008 */
00009 
00010 #ifndef SA_IMPL_HPP
00011 #define SA_IMPL_HPP
00012 
00013 #ifdef USE_MPI
00014 #include "angel_comm.hpp"
00015 #endif // USE_MPI
00016 
00017 namespace angel {
00018   using std::vector;
00019 
00020 template <class Object_t, class Neighbor_t, class Cost_t, class Temp_t>
00021 int SA (Object_t& object, Neighbor_t neighbor, Cost_t costs, Temp_t temp, int max_iter) {
00022   int min_costs= costs (object), last_costs= min_costs; // initial costs
00023   Object_t min_object (object), last_object (object); // corresponding objects
00024 
00025   for (int i= 0; i < max_iter; i++) {
00026     Object_t new_object (last_object);
00027     if (!neighbor (new_object)) {
00028       cout << "No neighbor found!"; return -1;}
00029     int new_costs= costs (new_object);
00030     // cout << "LSA costs: " << new_costs;
00031 #ifdef GNUPLOT
00032     cout << i << "\t" << new_costs << '\n'; // for gnuplot
00033 #endif
00034     if (new_costs < last_costs) {
00035       if (new_costs < min_costs) {
00036         min_costs= new_costs; min_object= new_object; }
00037       last_costs= new_costs; last_object= new_object; 
00038       // cout << " accepted due to improvement";
00039     } else {
00040       double acc= SA_acceptance (last_costs-new_costs, i, temp);
00041       double ra= angel::random (1.0);
00042       if (ra < acc) {
00043         last_costs= new_costs; last_object= new_object; 
00044         // cout << " accepted due to Metropolis";
00045       } } 
00046     // cout << '\n'; 
00047   }
00048   object= min_object;
00049   return min_costs;
00050 }
00051 
00052 template <class Object_t, class Neighbor_t, class Cost_t, class Adaption_t>
00053 int ALSA (Object_t& object, Neighbor_t neighbor, Cost_t costs, Adaption_t adaption,
00054           double& gamma, int max_iter) {
00055 
00056   int min_costs= costs (object), last_costs= min_costs; // initial costs
00057   Object_t min_object (object), last_object (object); // corresponding objects
00058 
00059   for (int i= 0; i < max_iter; i++) {
00060     Object_t new_object (last_object);
00061     if (!neighbor (new_object)) {
00062       cout << "No neighbor found!"; return -1;}
00063     int new_costs= costs (new_object);
00064     //cout << "LSA costs: " << new_costs << '\n';
00065 #ifdef GNUPLOT
00066     cout << i << "\t" << new_costs << '\n'; // for gnuplot
00067 #endif
00068     adaption (new_costs, gamma);
00069     if (new_costs < last_costs) {
00070       if (new_costs < min_costs) {
00071         min_costs= new_costs; min_object= new_object; }
00072       last_costs= new_costs; last_object= new_object; 
00073       // cout << "new configuration accepted due to improvement\n";
00074     } else {
00075       double acc= exp ((last_costs-new_costs) / gamma * log (double (i+2)));
00076       double ra= angel::random (1.0);
00077       if (ra < acc) {
00078         last_costs= new_costs; last_object= new_object; 
00079         // cout << "new configuration accepted due to Metropolis\n";
00080       } } }
00081 
00082   object= min_object;
00083   return min_costs;
00084 }
00085 
00086 // =====================================================
00087 // neighbourhoods
00088 // =====================================================
00089 
00090   template <class Ad_graph_t, class El_spec_t>
00091   bool neighbor_last_removable_t::operator() (elimination_history_t<Ad_graph_t, El_spec_t>& eh) {
00092     bool eliminate= eh.seq.size() == 0 || angel::random (1.0) < .5;
00093     if (eliminate) {
00094       vector<El_spec_t>   el;
00095       eliminatable_objects (eh.cg, el);
00096       if (el.empty() && eh.seq.empty()) return false; // graph was always bipartite
00097       if (el.size() > 0) {  // if nothing to eliminate -> re-insert
00098         int nextn= angel::random (int (el.size()));
00099         El_spec_t next= el[nextn];
00100         return eh.elimination (next) > 0; } } // elimination successful ?
00101     eh.seq.resize (eh.seq.size()-1);
00102     return eh.rebuild_graph (); 
00103   }
00104 
00105 // -------------------------------------------------------------------------
00106 
00107   template <class Ad_graph_t, class El_spec_t>
00108   bool neighbor_multi_step_t::operator() (elimination_history_t<Ad_graph_t, El_spec_t>& eh) {
00109     bool eliminate= eh.seq.size() == 0 || angel::random (1.0) < .5;
00110     int steps= angel::random (1, max_steps);
00111     if (eliminate) {
00112       for (int i= 0; i < steps; i++) {
00113         vector<El_spec_t>   el;
00114         eliminatable_objects (eh.cg, el);
00115         if (el.empty() && eh.seq.empty()) return false; // graph was always bipartite
00116         if (el.size() == 0) return true;
00117         int nextn= angel::random (int (el.size()));
00118         El_spec_t next= el[nextn];
00119         bool okay= eh.elimination (next) > 0; // elimination successful ?
00120         if (!okay) return false; }
00121       return true; // enough eliminations
00122     } else {
00123       eh.seq.resize (std::max (int (eh.seq.size())-steps, 0));
00124       return eh.rebuild_graph (); } 
00125   }
00126 
00127 // -------------------------------------------------------------------------
00128 
00129   template <class Ad_graph_t, class El_spec_t>
00130   bool neighbor_sequence_check_t::operator() (elimination_history_t<Ad_graph_t, El_spec_t>& eh) {
00131     bool eliminate= eh.seq.size() == 0 || angel::random (1.0) < .5;
00132     if (eliminate) {
00133       vector<El_spec_t>   el;
00134       eliminatable_objects (eh.cg, el);
00135       if (el.empty() && eh.seq.empty()) return false; // graph was always bipartite
00136       if (el.size() > 0) {  // if nothing to eliminate -> re-insert
00137         int nextn= angel::random (int (el.size()));
00138         El_spec_t next= el[nextn];
00139         return eh.elimination (next) > 0; } } // elimination successful ?
00140 
00141     int next_re_ins;
00142     bool reinsertable= false;
00143     for (int c= 1; !reinsertable; c++) {
00144       next_re_ins= angel::random_high (int (eh.seq.size()), c);
00145       vector<El_spec_t> seq_copy (eh.seq);
00146       typename vector<El_spec_t>::iterator it= eh.seq.begin() + next_re_ins;
00147       eh.seq.erase (it);
00148       reinsertable= eh.rebuild_graph ();
00149       // if reinsertion failed then restore old seq, cg is unchanged then
00150       if (!reinsertable) eh.seq.swap (seq_copy); }
00151     return true; }
00152 
00153 // -------------------------------------------------------------------------
00154 
00155   template <class Ad_graph_t, class El_spec_t>
00156   bool neighbor_check_meta_t::operator() (elimination_history_t<Ad_graph_t, El_spec_t>& eh) {
00157     bool eliminate= eh.seq.size() == 0 || angel::random (1.0) < .5;
00158     if (eliminate) {
00159       vector<El_spec_t>   el;
00160       eliminatable_objects (eh.cg, el);
00161       if (el.empty() && eh.seq.empty()) return false; // graph was always bipartite
00162       if (el.size() > 0) {  // if nothing to eliminate -> re-insert
00163         int nextn= angel::random (int (el.size()));
00164         El_spec_t next= el[nextn];
00165         return eh.elimination (next) > 0; } } // elimination successful ?
00166 
00167     int next_re_ins;
00168     bool reinsertable= false;
00169     for (int c= 1; !reinsertable; c++) {
00170       next_re_ins= angel::random_high (int (eh.seq.size()), c);
00171       elimination_history_t<Ad_graph_t, El_spec_t> eh_copy (eh);
00172       typename vector<El_spec_t>::iterator it= eh_copy.seq.begin() + next_re_ins;
00173       El_spec_t re= *it;  // removed elimination
00174       eh_copy.seq.erase (it);
00175       reinsertable= eh_copy.rebuild_graph ();
00176       // if successful check if new graph (eh_copy.cg) is predecessor of eh.cg
00177       // eh_copy.cg --re--> eh.cg
00178       if (reinsertable) {
00179         Ad_graph_t copy_copy (eh_copy.cg);
00180         angel::eliminate (re, copy_copy);
00181         reinsertable= eh.cg == copy_copy; }
00182       if (reinsertable) eh= eh_copy; }
00183     return true; }
00184 
00185 
00186 // =====================================================
00187 // Gamma adaption operators
00188 // =====================================================
00189 
00190 /*  void gamma_adaption_max_t::operator() (int costs, double& gamma) {
00191     if (imp >= D) return;
00192     if (last_min == 0) { // very first iteration must be distinguished
00193       last_max= last_min= costs; return; }
00194     if (costs < last_min) {
00195       diff= last_max - costs; last_max= last_min= costs;
00196       if (diff > max_diff) max_diff= diff;
00197       if (++imp >= D) gamma= scaling * double (max_diff);
00198     } else if (costs > last_max) last_max= costs;
00199   }
00200 
00201   void gamma_adaption_average_t::operator() (int costs, double& gamma) {
00202     if (imp >= D) return;
00203     if (last_min == 0) { // very first iteration must be distinguished
00204       last_max= last_min= costs; return; }
00205     if (costs < last_min) {
00206       sum_diff+= last_max - costs; last_max= last_min= costs;
00207       if (++imp >= D) {gamma= scaling * double (sum_diff) / double (imp); }
00208     } else if (costs > last_max) last_max= costs;
00209   }
00210 */
00211 // ============ Parallel computations only if USE_MPI is defined ====================
00212 
00213 #ifdef USE_MPI
00214 
00215 template <class Temp_t>
00216 int pick_processor_sa (int my_costs, int it, Temp_t temp, const GMPI::Intracomm& comm) {
00217   int          me= comm.Get_rank(), size= comm.Get_size(),  picked;
00218   vector<int>  all_costs (size);
00219   comm.Gather (&my_costs, 1, MPI::INT, &all_costs[0], 1, MPI::INT, 0);
00220   if (me == 0) {
00221     int min_costs= *min_element (all_costs.begin(), all_costs.end());
00222     // compute acceptance probs and their partial sums
00223     vector<double>  acc_probs (size), psum_probs (size);
00224     for (size_t c= 0; c < size; c++)
00225       acc_probs[c]= SA_acceptance (min_costs - all_costs[c], it, temp);
00226     partial_sum (acc_probs.begin(), acc_probs.end(), psum_probs.begin());
00227     // pick one according to probabilities
00228     double sum_probs= psum_probs[size-1], ra= angel::random (sum_probs);
00229     for (picked= 0; picked < size; picked++)
00230       if (ra < psum_probs[picked]) break;
00231     THROW_DEBUG_EXCEPT_MACRO (picked == size, consistency_exception, "No processor picked");
00232 #ifdef WRITE_LOG_FILE
00233     THROW_DEBUG_EXCEPT_MACRO (!log_file.is_open(), io_exception, "Log file not opened");    
00234     log_file << "pick_processor_sa at iteration " << it << endl;
00235     write_vector (log_file, "Processor's costs: ", all_costs);
00236     write_vector (log_file, "Accumulated probabilities: ", psum_probs);
00237     log_file << "Random value: " << ra << " --> picked processor " << picked << endl;
00238 #endif
00239   }
00240   comm.Bcast (&picked, 1, MPI::INT, 0);
00241   return picked;
00242 }
00243 
00244 template <class Object_t, class Neighbor_t, class Cost_t, class Inner_temp_t, class Outer_temp_t,
00245           class Pre_comm_t, class Post_comm_t>
00246 int parallel_SA (Object_t& object, Neighbor_t neighbor, Cost_t costs, 
00247                  Inner_temp_t inner_temp, Outer_temp_t outer_temp, int inner_iter, int max_iter,
00248                  const GMPI::Intracomm& comm, Pre_comm_t pre_comm, Post_comm_t post_comm) {
00249 
00250   int min_costs= costs (object), last_costs= min_costs; // initial costs
00251   Object_t min_object (object), last_object (object); // corresponding objects
00252   GMPI::comm_ref_t<int, Object_t> comm_ref (last_object); // reference used in communication
00253 
00254   for (int i= 0; i < max_iter; i++) {
00255 
00256     int ii; bool inner_completion;
00257     for (ii= 0, inner_completion= false; !inner_completion; ) {
00258 
00259       Object_t new_object (last_object);
00260       if (!neighbor (new_object)) {
00261         cout << "No neighbor found!"; return -1;}
00262       int new_costs= costs (new_object);
00263       // log_file << "LSA costs: " << new_costs;
00264 #ifdef GNUPLOT
00265       log_file << i+ii << "\t" << new_costs << "   # gnuplot \n"; // different files per proc with MPI
00266 #endif
00267       if (new_costs < last_costs) {
00268         if (new_costs < min_costs) {
00269           min_costs= new_costs; min_object= new_object; }
00270         last_costs= new_costs; last_object= new_object; 
00271         // log_file << " accepted due to improvement";
00272       } else {
00273         double acc= SA_acceptance (last_costs-new_costs, i, inner_temp);
00274         double ra= angel::random (1.0);
00275         if (ra < acc) {
00276           last_costs= new_costs; last_object= new_object; 
00277           // log_file << " accepted due to Metropolis";
00278         } } 
00279       // log_file << '\n'; 
00280       if (++ii == inner_iter) { 
00281         inner_completion= true;
00282         mark_completion (comm.mpi_comm_ref());
00283       } else inner_completion= test_completion (comm.mpi_comm_ref());
00284     } // inner for loop
00285 
00286     i+= inner_iter; // all proc are handled as having finished inner loop here
00287     if (i < max_iter) {
00288       clean_completion_messages (comm.mpi_comm_ref());
00289       pre_comm (last_object);
00290       int sa_root= pick_processor_sa (last_costs, i, outer_temp, comm);
00291       comm.Bcast (comm_ref, sa_root); 
00292       post_comm (last_object); }
00293   } // outer for loop
00294 
00295   // now look for global minimum
00296   int me= comm.Get_rank(); 
00297   std::pair<int,int> my_min_costs_rank (min_costs, me), min_costs_rank;
00298   comm.Allreduce (&my_min_costs_rank, &min_costs_rank, 1, MPI::TWOINT, MPI::MINLOC); 
00299   int min_root= min_costs_rank.second;
00300   GMPI::comm_ref_t<int, Object_t> min_comm_ref (min_object); // reference to min
00301   comm.Bcast (min_comm_ref, min_root);
00302 
00303   object= min_object;
00304   return min_root;
00305 }
00306 
00307 #endif // USE_MPI
00308 
00309 } // namespace angel
00310 
00311 #endif /* SA_IMPL_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines