angel
mercurial changeset:
|
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 */