source: trunk/Bonmin/experimental/Bcp/BM_tm.cpp @ 508

Last change on this file since 508 was 508, checked in by pbonami, 13 years ago

Make bcp-bonmin work with new setup

File size: 8.7 KB
Line 
1#include "OsiClpSolverInterface.hpp"
2// (C) Copyright International Business Machines Corporation 2006, 2007
3// All Rights Reserved.
4// This code is published under the Common Public License.
5//
6// Authors :
7// Laszlo Ladanyi, International Business Machines Corporation
8// Pierre Bonami, Carnegie Mellon University
9
10#include "BonIpoptSolver.hpp"
11#ifdef COIN_HAS_FILTERSQP
12# include "BonFilterSolver.hpp"
13#endif
14#include "BCP_problem_core.hpp"
15#include "BM.hpp"
16
17
18#include "BonOACutGenerator2.hpp"
19#include "BonEcpCuts.hpp"
20#include "BonOaNlpOptim.hpp"
21
22#include "BonAmplSetup.hpp"
23
24//#############################################################################
25
26void
27BM_tm::readIpopt()
28{
29    if (par.entry(BM_par::IpoptParamfile).length() != 0) {
30        FILE* ipopt = fopen(par.entry(BM_par::IpoptParamfile).c_str(), "r");
31        if (!ipopt) {
32            throw BCP_fatal_error("Non-existent IpoptParamfile\n");
33        }
34        fseek(ipopt, 0, SEEK_END);
35        int len = ftell(ipopt) + 1;
36        fseek(ipopt, 0, SEEK_SET);
37        char* ipopt_content = new char[len+1];
38        len = fread(ipopt_content, 1, len, ipopt);
39        ipopt_file_content.assign(ipopt_content, len);
40        delete[] ipopt_content;
41    }
42
43    FILE* nl = fopen(par.entry(BM_par::NL_filename).c_str(), "r");
44    if (!nl) {
45        throw BCP_fatal_error("Non-existent NL_filename\n");
46    }
47    fseek(nl, 0, SEEK_END);
48    int len = ftell(nl) + 1;
49    fseek(nl, 0, SEEK_SET);
50    char* nl_content = new char[len+1];
51    len = fread(nl_content, 1, len, nl);
52    nl_file_content.assign(nl_content, len);
53    delete[] nl_content;
54}
55
56/****************************************************************************/
57
58void
59BM_tm::initialize_core(BCP_vec<BCP_var_core*>& vars,
60                       BCP_vec<BCP_cut_core*>& cuts,
61                       BCP_lp_relax*& matrix)
62{
63    char* argv_[3];
64    char** argv = argv_;
65    argv[0] = NULL;
66    argv[1] = strdup(par.entry(BM_par::NL_filename).c_str());
67    argv[2] = NULL;
68   
69    /* Get the basic options. */
70    Bonmin::BonminAmplSetup bonmin;
71    bonmin.readOptionsFile(par.entry(BM_par::IpoptParamfile).c_str());
72    bonmin.initializeBonmin(argv);   
73   
74    Bonmin::OsiTMINLPInterface& nlpSolver = *bonmin.nonlinearSolver();
75   
76    free(argv[1]);
77   
78    OsiSolverInterface& clp  = *bonmin.continuousSolver();
79   
80    const int numCols = clp.getNumCols();
81    const int numRows = clp.getNumRows();
82
83    double* obj = new double[numCols];
84    if (bonmin.getAlgorithm() == Bonmin::B_BB /* pure B&B */) {
85      std::cout<<"Doing branch and bound"<<std::endl;
86      CoinFillN(obj, numCols, 0.0);
87    }
88    else {
89      std::cout<<"Doing hybrid"<<std::endl;
90      CoinDisjointCopyN(clp.getObjCoefficients(), numCols, obj);
91    }
92
93    vars.reserve(numCols);
94    const double* clb = clp.getColLower();
95    const double* cub = clp.getColUpper();
96    for (int i = 0; i < numCols; ++i)
97        {
98            BCP_var_t type = BCP_ContinuousVar;
99            if (clp.isBinary(i)) type = BCP_BinaryVar;
100            if (clp.isInteger(i)) type = BCP_IntegerVar;
101            vars.push_back(new BCP_var_core(type, obj[i], clb[i], cub[i]));
102        }
103
104    if (bonmin.getAlgorithm() == Bonmin::B_BB /* pure B&B */) {
105        // Just fake something into the core matrix. In this case: 0 <= 1
106        BCP_vec<double> OBJ(obj, numCols);
107        BCP_vec<double> CLB(clb, numCols);
108        BCP_vec<double> CUB(cub, numCols);
109        BCP_vec<double> RLB(1, 0.0);
110        BCP_vec<double> RUB(1, 1.0);
111        BCP_vec<int> VB;
112        BCP_vec<int> EI;
113        BCP_vec<double> EV;
114        VB.push_back(0);
115        VB.push_back(2);
116        EI.push_back(0);         EV.push_back(0.0);
117        EI.push_back(numCols-1); EV.push_back(0.0);
118        matrix = new BCP_lp_relax(false, VB, EI, EV, OBJ, CLB, CUB, RLB, RUB);
119        cuts.push_back(new BCP_cut_core(0.0, 1.0));
120    } else {
121        cuts.reserve(numRows);
122        const double* rlb = clp.getRowLower();
123        const double* rub = clp.getRowUpper();
124        for (int i = 0; i < numRows; ++i)
125            {
126                cuts.push_back(new BCP_cut_core(rlb[i], rub[i]));
127            }
128
129        matrix = new BCP_lp_relax(true /*column major ordered*/);
130        matrix->copyOf(*clp.getMatrixByCol(), obj, clb, cub, rlb, rub);
131    }
132    delete[] obj;
133}
134
135/****************************************************************************/
136void
137BM_tm::create_root(BCP_vec<BCP_var*>& added_vars,
138                   BCP_vec<BCP_cut*>& added_cuts,
139                   BCP_user_data*& user_data)
140{
141    BM_node* data = new BM_node;
142    data->numNlpFailed_ = 0;
143    user_data = data;
144}
145
146void
147BM_tm::write_AMPL_solution(const BCP_solution* sol,
148                           bool write_file, bool write_screen)
149{
150  const BM_solution* bs = dynamic_cast<const BM_solution*>(sol);
151  if (!bs) {
152    throw BCP_fatal_error("Trying to pack non-BM_solution.\n");
153  }
154  /* Parse again the input file so that we have a nice and clean ampl
155     setup */ 
156 
157  char* argv_[3];
158  char** argv = argv_;
159  argv[0] = NULL;
160  argv[1] = strdup(par.entry(BM_par::NL_filename).c_str());
161  argv[2] = NULL;
162  Bonmin::BonminAmplSetup bonmin;
163  bonmin.initializeBonmin(argv);   
164 
165  Bonmin::OsiTMINLPInterface& nlpSolver = *bonmin.nonlinearSolver();
166  minlpParams_.extractParams(&nlpSolver);
167  free(argv[1]);
168  OsiSolverInterface& clp = *bonmin.continuousSolver();
169
170  const int numCols = clp.getNumCols();
171
172  int i;
173 
174  /* This will give the vector of core variables we have created in
175     BM_tm::initialize_core */
176  const BCP_vec<BCP_var_core*>& vars = getTmProblemPointer()->core->vars;
177
178  /* Create a dense vector with the value of each variable. Round the
179     integer vars (they were tested to be near enough to integrality) so
180     in the printouts we won't have 9.99999991234, etc.) */
181  double* dsol = new double[numCols];
182  for (i = 0; i < numCols; ++i) {
183    dsol[i] = 0.0;
184  }
185  const int size = bs->_ind.size();
186  for (i = 0; i < size; ++i) {
187    const int ind = bs->_ind[i];
188    const double v = bs->_values[i];
189    const BCP_var_t type = vars[ind]->var_type();
190    dsol[ind] = (type == BCP_ContinuousVar) ? v : floor(v+0.5);
191  }
192
193  if (write_screen) {
194    /* Display the solution on the screen */
195    printf("bonmin: feasible solution found.  Objective value: %f\n",
196           bs->_objective);
197    for (i = 0; i < size; ++i) {
198      printf("    index: %5i   value: %f\n", bs->_ind[i], dsol[bs->_ind[i]]);
199    }
200    printf("\n");
201  }
202
203  if (write_file) {
204    /* create the AMPL solfile */
205    nlpSolver.model()->finalize_solution(Bonmin::TMINLP::SUCCESS, nlpSolver.getNumCols(), 
206                                         dsol, bs->_objective);
207  }
208  delete[] dsol;
209}
210
211//#############################################################################
212
213const BCP_proc_id*
214BM_tm::process_id() const
215{
216}
217
218void
219BM_tm::send_message(const BCP_proc_id* const target, const BCP_buffer& buf)
220{
221}
222
223void
224BM_tm::broadcast_message(const BCP_process_t proc_type, const BCP_buffer& buf)
225{
226}
227
228void
229BM_tm::process_message(BCP_buffer& buf)
230{
231}
232
233/****************************************************************************/
234
235void
236BM_tm::display_final_information(const BCP_lp_statistics& lp_stat)
237{
238  bool write_screen = false;
239  BCP_tm_prob *p = getTmProblemPointer();
240  if (p->param(BCP_tm_par::TmVerb_FinalStatistics)) {
241    printf("TM: Running time: %.3f\n", CoinCpuTime() - p->start_time);
242    printf("TM: search tree size: %i   ( processed %i )   max depth: %i\n",
243           int(p->search_tree.size()), int(p->search_tree.processed()),
244           p->search_tree.maxdepth());
245    lp_stat.display();
246
247    if (! p->feas_sol) {
248      printf("TM: No feasible solution is found\n");
249    } else {
250      printf("TM: The best solution found has value %f\n",
251             p->feas_sol->objective_value());
252      if (p->param(BCP_tm_par::TmVerb_BestFeasibleSolution)) {
253        write_screen = true;
254      }
255    }
256  }
257  if (p->feas_sol) {
258    write_AMPL_solution(p->feas_sol, true, write_screen);
259  }
260}
261
262/****************************************************************************/
263
264void
265BM_tm::display_feasible_solution(const BCP_solution* sol)
266{
267  // For now, we want to write also the AMPL solution file...(?)
268  // This needs to be changed though
269  write_AMPL_solution(sol, true, true);
270}
271
272struct BMSearchTreeCompareBest {
273  inline bool operator()(const CoinTreeSiblings* x,
274                         const CoinTreeSiblings* y) const {
275    double allowable_gap = 1e-8;
276    return ((x->currentNode()->quality_ <
277             y->currentNode()->quality_-allowable_gap) ||
278            ((x->currentNode()->quality_ <= y->currentNode()->quality_+allowable_gap) &&
279             x->currentNode()->depth_ > y->currentNode()->depth_));
280    }
281};
282
283void
284BM_tm::init_new_phase(int phase,
285                      BCP_column_generation& colgen,
286                      CoinSearchTreeBase*& candidates)
287{
288  BCP_tm_prob* p = getTmProblemPointer();
289  colgen = BCP_DoNotGenerateColumns_Fathom;
290  switch (p->param(BCP_tm_par::TreeSearchStrategy)) {
291  case BCP_BestFirstSearch:
292    candidates = new CoinSearchTree<BMSearchTreeCompareBest>;
293    break;
294  case BCP_BreadthFirstSearch:
295    candidates = new CoinSearchTree<CoinSearchTreeCompareBreadth>;
296    break;
297  case BCP_DepthFirstSearch:
298    candidates = new CoinSearchTree<CoinSearchTreeCompareDepth>;
299    break;
300  }
301}
302
Note: See TracBrowser for help on using the repository browser.