source: branches/devel/Bonmin/src/Interfaces/BonCurvatureEstimator.cpp @ 104

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

Add option to choose branching strategy,
include cmath in BonCurvatureEstimator?

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 16.9 KB
Line 
1// Copyright (C) 2006 International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: BonCurvatureEstimator.cpp 99 2006-10-31 04:55:45Z pbonami $
6//
7// Author:   Andreas Waechter                 IBM    2006-10-11
8
9#include "BonCurvatureEstimator.hpp"
10#include "IpGenTMatrix.hpp"
11#include "IpIdentityMatrix.hpp"
12#include "IpZeroMatrix.hpp"
13#include "IpDenseVector.hpp"
14#include "IpBlas.hpp"
15#include <cmath>
16#ifdef HAVE_MA27
17# include "IpMa27TSolverInterface.hpp"
18#endif
19#ifdef HAVE_MA57
20# include "IpMa57TSolverInterface.hpp"
21#endif
22#ifdef HAVE_MC19
23# include "IpMc19TSymScalingMethod.hpp"
24#endif
25#ifdef HAVE_PARDISO
26# include "IpPardisoSolverInterface.hpp"
27#endif
28#ifdef HAVE_TAUCS
29# include "IpTAUCSSolverInterface.hpp"
30#endif
31#ifdef HAVE_WSMP
32# include "IpWsmpSolverInterface.hpp"
33#endif
34#ifdef HAVE_MUMPS
35# include "IpMumpsSolverInterface.hpp"
36#endif
37
38namespace Bonmin
39{
40  using namespace Ipopt;
41
42  // ToDo: Consider NLP scaling?
43
44  CurvatureEstimator::CurvatureEstimator
45    (SmartPtr<Journalist> jnlst,
46     SmartPtr<OptionsList> options,
47     SmartPtr<TNLP> tnlp)
48      :
49      jnlst_(jnlst),
50      options_(options),
51      prefix_(""),
52      tnlp_(tnlp),
53      grad_f_(NULL),
54      irows_jac_(NULL),
55      jcols_jac_(NULL),
56      jac_vals_(NULL),
57      irows_hess_(NULL),
58      jcols_hess_(NULL),
59      hess_vals_(NULL),
60      x_l_(NULL),
61      x_u_(NULL),
62      g_l_(NULL),
63      g_u_(NULL),
64      x_free_map_(NULL),
65      g_fixed_map_(NULL),
66      lambda_(NULL),
67      initialized_(false)
68  {
69    DBG_ASSERT(IsValid(jnlst));
70    DBG_ASSERT(IsValid(options));
71    DBG_ASSERT(IsValid(tnlp));
72
73    ////////////////////////////////////////////////////////////////////
74    // Create a strategy object for solving the linear system for the //
75    // projection matrix                                              //
76    ////////////////////////////////////////////////////////////////////
77
78    // The following linear are taken from AlgBuilder in Ipopt
79    SmartPtr<SparseSymLinearSolverInterface> SolverInterface;
80    std::string linear_solver;
81    options->GetStringValue("linear_solver", linear_solver, prefix_);
82    if (linear_solver=="ma27") {
83#ifdef HAVE_MA27
84      SolverInterface = new Ma27TSolverInterface();
85#else
86
87      THROW_EXCEPTION(OPTION_INVALID,
88                      "Selected linear solver MA27 not available.");
89#endif
90
91    }
92    else if (linear_solver=="ma57") {
93#ifdef HAVE_MA57
94      SolverInterface = new Ma57TSolverInterface();
95#else
96
97      THROW_EXCEPTION(OPTION_INVALID,
98                      "Selected linear solver MA57 not available.");
99#endif
100
101    }
102    else if (linear_solver=="pardiso") {
103#ifdef HAVE_PARDISO
104      SolverInterface = new PardisoSolverInterface();
105#else
106
107      THROW_EXCEPTION(OPTION_INVALID,
108                      "Selected linear solver Pardiso not available.");
109#endif
110
111    }
112    else if (linear_solver=="taucs") {
113#ifdef HAVE_TAUCS
114      SolverInterface = new TAUCSSolverInterface();
115#else
116
117      THROW_EXCEPTION(OPTION_INVALID,
118                      "Selected linear solver TAUCS not available.");
119#endif
120
121    }
122    else if (linear_solver=="wsmp") {
123#ifdef HAVE_WSMP
124      SolverInterface = new WsmpSolverInterface();
125#else
126
127      THROW_EXCEPTION(OPTION_INVALID,
128                      "Selected linear solver WSMP not available.");
129#endif
130
131    }
132    else if (linear_solver=="mumps") {
133#ifdef HAVE_MUMPS
134      SolverInterface = new MumpsSolverInterface();
135#else
136
137      THROW_EXCEPTION(OPTION_INVALID,
138                      "Selected linear solver MUMPS not available.");
139#endif
140
141    }
142
143    SmartPtr<TSymScalingMethod> ScalingMethod;
144    std::string linear_system_scaling;
145    if (!options->GetStringValue("linear_system_scaling",
146                                 linear_system_scaling, prefix_)) {
147      // By default, don't use mc19 for non-HSL solvers
148      if (linear_solver!="ma27" && linear_solver!="ma57") {
149        linear_system_scaling="none";
150      }
151    }
152    if (linear_system_scaling=="mc19") {
153#ifdef HAVE_MC19
154      ScalingMethod = new Mc19TSymScalingMethod();
155#else
156
157      THROW_EXCEPTION(OPTION_INVALID,
158                      "Selected linear system scaling method MC19 not available.");
159#endif
160
161    }
162
163    tsymlinearsolver_ = new TSymLinearSolver(SolverInterface, ScalingMethod);
164    // End of lines from AlgBuilder
165  }
166
167  CurvatureEstimator::~CurvatureEstimator()
168  {
169    if (initialized_) {
170      delete [] grad_f_;
171      delete [] irows_jac_;
172      delete [] jcols_jac_;
173      delete [] jac_vals_;
174      delete [] irows_hess_;
175      delete [] jcols_hess_;
176      delete [] hess_vals_;
177      delete [] x_l_;
178      delete [] x_u_;
179      delete [] g_l_;
180      delete [] g_u_;
181      delete [] x_free_map_;
182      delete [] g_fixed_map_;
183      delete [] lambda_;
184    }
185  }
186
187  bool CurvatureEstimator::Initialize()
188  {
189    DBG_ASSERT(!initialized_);
190
191    //////////////////////////////////////
192    // Prepare internal data structures //
193    //////////////////////////////////////
194
195    // Get sizes
196    TNLP::IndexStyleEnum index_style;
197    if (!tnlp_->get_nlp_info(n_, m_, nnz_jac_, nnz_hess_, index_style)) {
198      return false;
199    }
200
201    // Space for gradient
202    grad_f_ = new Number[n_];
203
204    // Get nonzero entries in the matrices
205    irows_jac_ = new Index[nnz_jac_];
206    jcols_jac_ = new Index[nnz_jac_];
207    if (!tnlp_->eval_jac_g(n_, NULL, false, m_, nnz_jac_,
208                           irows_jac_, jcols_jac_, NULL)) {
209      return false;
210    }
211    if (index_style == TNLP::FORTRAN_STYLE) {
212      for (Index i=0; i<nnz_jac_; i++) {
213        irows_jac_[i]--;
214        jcols_jac_[i]--;
215      }
216    }
217    jac_vals_ = new Number[nnz_jac_];
218
219    irows_hess_ = new Index[nnz_hess_];
220    jcols_hess_ = new Index[nnz_hess_];
221    if (!tnlp_->eval_h(n_, NULL, false, 1., m_, NULL, false, nnz_hess_,
222                       irows_hess_, jcols_hess_, NULL)) {
223      return false;
224    }
225    if (index_style == TNLP::FORTRAN_STYLE) {
226      for (Index i=0; i<nnz_hess_; i++) {
227        irows_hess_[i]--;
228        jcols_hess_[i]--;
229      }
230    }
231    hess_vals_ = NULL; // We set it to NULL, so that we know later
232    // that we still need to compute the values
233
234    // Space for bounds
235    x_l_ = new Number[n_];
236    x_u_ = new Number[n_];
237    g_l_ = new Number[m_];
238    g_u_ = new Number[m_];
239
240    // Get space for the activities maps
241    x_free_map_ = new Index[n_];
242    g_fixed_map_ = new Index[m_];
243
244    // Get space for the multipliers
245    lambda_ = new Number[m_];
246
247    initialized_ = true;
248    return true;
249  }
250
251  bool
252  CurvatureEstimator::ComputeNullSpaceCurvature(
253    bool new_bounds,
254    int n,
255    const Number* x,
256    bool new_x,
257    const Number* z_L,
258    const Number* z_U,
259    int m,
260    const Number* lam,
261    bool new_mults,
262    const Number* orig_d,
263    Number* projected_d,
264    Number& gradLagTd,
265    Number& dTHLagd)
266  {
267    if (!initialized_) {
268      Initialize();
269      new_bounds = true;
270      new_mults = true;
271    }
272
273    DBG_ASSERT(n == n_);
274    DBG_ASSERT(m == m_);
275
276    // If necessary, get new Jacobian values (for the original matrix)
277    if (new_x) {
278      if (!tnlp_->eval_jac_g(n_, x, new_x, m_, nnz_jac_,
279                             NULL, NULL, jac_vals_)) {
280        return false;
281      }
282    }
283
284    if (new_bounds) {
285      // Get bounds
286      if (!tnlp_->get_bounds_info(n_, x_l_, x_u_, m_, g_l_, g_u_)) {
287        return false;
288      }
289    }
290
291    // If necessary, determine new activities
292    bool new_activities = false;
293    if (new_bounds || new_mults) {
294      new_activities = true;
295      active_x_.clear();
296      const Number zTol = 1e-4; // ToDo: make this more elaborate
297      jnlst_->Printf(J_MOREDETAILED, J_NLP,
298                     "List of variables considered fixed (with z_L and z_U)\n");
299      for (Index i=0; i<n; i++) {
300        if(z_L[i] > zTol || z_U[i] > zTol) {
301          active_x_.push_back(i+1);
302          jnlst_->Printf(J_MOREDETAILED, J_NLP,
303                         "x[%5d] (%e,%e)\n", i, z_L[i], z_U[i]);
304        }
305      }
306
307      active_g_.clear();
308      // Compute the product of the direction with the constraint Jacobian
309      // This could be done more efficient if we have d in sparse format
310      Number* jacTd = new Number[m];
311      const Number zero = 0.;
312      IpBlasDcopy(m, &zero, 0, jacTd, 1);
313      for (Index i=0; i<nnz_jac_; i++) {
314        const Index &irow = irows_jac_[i];
315        const Index &jcol = jcols_jac_[i];
316        jacTd[irow] += jac_vals_[i]*orig_d[jcol];
317      }
318
319      const Number lamTol = 1e-4;
320      jnlst_->Printf(J_MOREDETAILED, J_NLP,
321                     "List of constraints considered fixed (with lam and jacTd)\n");
322      for (Index j=0; j<m; j++) {
323        if (g_l_[j] < g_u_[j] && fabs(lam[j]) > lamTol) {
324          if (lam[j]*jacTd[j] > 0) {
325            active_g_.push_back(j+1);
326            jnlst_->Printf(J_MOREDETAILED, J_NLP,
327                           "g[%5d] (%e,%e)\n", j, lam[j], jacTd[j]);
328          }
329        }
330      }
331      delete [] jacTd;
332    }
333
334    // Check if the structure of the matrix has changed
335    if (new_activities) {
336      if (!PrepareNewMatrixStructure(new_activities)) {
337        return false;
338      }
339    }
340
341    bool new_lambda = false;
342    if (new_x || new_activities) {
343      if (!PrepareNewMatrixValues(new_activities, x, new_x)) {
344        return false;
345      }
346
347      // Compute least square multipliers for the given activities
348      if (!tnlp_->eval_grad_f(n_, x, new_x, grad_f_)) {
349        return false;
350      }
351      if (!SolveSystem(grad_f_, NULL, NULL, lambda_)) {
352        return false;
353      }
354      IpBlasDscal(m_, -1., lambda_, 1);
355      new_lambda = true;
356    }
357
358    // Compute the projection of the direction
359    if (!SolveSystem(orig_d, NULL, projected_d, NULL)) {
360      return false;
361    }
362
363    // Compute the product with the gradient of the Lagrangian
364    gradLagTd = IpBlasDdot(n, projected_d, 1, grad_f_, 1);
365    for (Index i=0; i<nnz_jac_; i++) {
366      const Index &irow = irows_jac_[i];
367      const Index &jcol = jcols_jac_[i];
368      gradLagTd += lambda_[irow]*jac_vals_[i]*projected_d[jcol];
369    }
370
371    // Compute the product with the Hessian of the Lagrangian
372    if (!Compute_dTHLagd(projected_d, x, new_x, lambda_, new_lambda, dTHLagd)) {
373      return false;
374    }
375
376    return true;
377  }
378
379  bool CurvatureEstimator::PrepareNewMatrixStructure(
380    bool new_activities)
381  {
382    if (new_activities) {
383      // Deterimine which variables are free
384      for (Index i=0; i<n_; i++) {
385        x_free_map_[i] = 0;
386      }
387      // fixed by activities
388      for (std::vector<int>::iterator i=active_x_.begin();
389           i != active_x_.end(); i++) {
390        DBG_ASSERT(*i != 0 && *i<=n_ && *i>=-n_);
391        if (*i<0) {
392          x_free_map_[(-*i)-1] = -1;
393          DBG_ASSERT(x_l_[(-*i)-1] > -1e19);
394        }
395        else {
396          x_free_map_[(*i)-1] = -1;
397          DBG_ASSERT(x_u_[(-*i)-1] < 1e19);
398        }
399      }
400      // fixed by bounds
401      for (Index i=0; i<n_; i++) {
402        if (x_l_[i] == x_u_[i]) {
403          x_free_map_[i] = -1;
404        }
405      }
406      // Correct the numbering in the x map and determine number of
407      // free variables
408      Index nx_free_ = 0;
409      for (Index i=0; i<n_; i++) {
410        if (x_free_map_[i] == 0) {
411          x_free_map_[i] = nx_free_;
412          nx_free_++;
413        }
414      }
415
416      // Determine which constraints are fixed
417      for (Index j=0; j<m_; j++) {
418        g_fixed_map_[j] = -1;
419      }
420      // fixed by activities
421      for (std::vector<int>::iterator i=active_g_.begin();
422           i != active_g_.end(); i++) {
423        DBG_ASSERT(*i != 0 && *i<=m_ && *i>=-m_);
424        if (*i<0) {
425          g_fixed_map_[(-*i)-1] = 0;
426          DBG_ASSERT(g_l_[(-*i)-1] > -1e19); //ToDo look at option?
427        }
428        else {
429          g_fixed_map_[(*i)-1] = 0;
430          DBG_ASSERT(g_u_[(*i)-1] < 1e19);
431        }
432      }
433      // fixed by bounds
434      for (Index j=0; j<m_; j++) {
435        if (g_l_[j] == g_u_[j]) {
436          g_fixed_map_[j] = 0;
437        }
438      }
439      // Correct the numbering in the g map and determine number of
440      // fixed constraints
441      Index ng_fixed_ = 0;
442      for (Index j=0; j<m_; j++) {
443        if (g_fixed_map_[j] == 0) {
444          g_fixed_map_[j] = ng_fixed_;
445          ng_fixed_++;
446        }
447      }
448
449      // Determine the row and column indices for the Jacobian of the fixed
450      // constraints in the space of the free variables
451      Index* iRows = new Index[nnz_jac_];
452      Index* jCols = new Index[nnz_jac_];
453      Index nnz_proj_jac = 0;
454      for (Index i=0; i<nnz_jac_; i++) {
455        const Index &irow = irows_jac_[i];
456        const Index &jcol = jcols_jac_[i];
457        if (x_free_map_[jcol] >= 0 && g_fixed_map_[irow] >= 0) {
458          iRows[nnz_proj_jac] = g_fixed_map_[irow] + 1;
459          jCols[nnz_proj_jac] = x_free_map_[jcol] + 1;
460          nnz_proj_jac++;
461        }
462      }
463
464      // Create the matrix space for the Jacobian matrices
465      SmartPtr<GenTMatrixSpace> proj_jac_space =
466        new GenTMatrixSpace(ng_fixed_, nx_free_, nnz_proj_jac, iRows, jCols);
467      delete [] iRows;
468      delete [] jCols;
469
470      // Create Matrix space for the projection matrix
471      comp_proj_matrix_space_ =
472        new CompoundSymMatrixSpace(2, nx_free_+ng_fixed_);
473      comp_proj_matrix_space_->SetBlockDim(0, nx_free_);
474      comp_proj_matrix_space_->SetBlockDim(1, ng_fixed_);
475      SmartPtr<SymMatrixSpace> identity_space =
476        new IdentityMatrixSpace(nx_free_);
477      comp_proj_matrix_space_->SetCompSpace(0, 0, *identity_space, true);
478      comp_proj_matrix_space_->SetCompSpace(1, 0, *proj_jac_space, true);
479      //      SmartPtr<MatrixSpace> zero_space =
480      //        new ZeroMatrixSpace(ng_fixed_, ng_fixed_);
481      //      comp_proj_matrix_space_->SetCompSpace(1, 1, *zero_space, true);
482
483      // Create a Vector space for the rhs and sol
484      comp_vec_space_ = new CompoundVectorSpace(2, nx_free_+ng_fixed_);
485      SmartPtr<DenseVectorSpace> x_space = new DenseVectorSpace(nx_free_);
486      comp_vec_space_->SetCompSpace(0, *x_space);
487      SmartPtr<DenseVectorSpace> g_space = new DenseVectorSpace(ng_fixed_);
488      comp_vec_space_->SetCompSpace(1, *g_space);
489    } 
490
491    return true;
492  }
493
494  bool CurvatureEstimator::PrepareNewMatrixValues(
495    bool new_activities,
496    const Number* x,
497    bool new_x)
498  {
499    if (new_x || new_activities) {
500      comp_proj_matrix_ = comp_proj_matrix_space_->MakeNewCompoundSymMatrix();
501      SmartPtr<Matrix> jac = comp_proj_matrix_->GetCompNonConst(1, 0);
502      SmartPtr<GenTMatrix> tjac = dynamic_cast<GenTMatrix*> (GetRawPtr(jac));
503      Number* vals = tjac->Values();
504      Index inz=0;
505      for (Index i=0; i<nnz_jac_; i++) {
506        Index irow = irows_jac_[i];
507        Index jcol = jcols_jac_[i];
508        if (x_free_map_[jcol] >= 0 && g_fixed_map_[irow] >= 0) {
509          vals[inz++] = jac_vals_[i];
510        }
511      }
512      comp_proj_matrix_->Print(*jnlst_,J_WARNING,J_MAIN,"compmat");
513      DBG_ASSERT(inz == tjac->Nonzeros());
514
515      // We need to reset the linear solver object, so that it knows
516      // that the structure of the linear system has changed
517      tsymlinearsolver_->ReducedInitialize(*jnlst_, *options_, prefix_);
518    }
519
520    return true;
521  }
522
523  bool CurvatureEstimator::SolveSystem(const Number* rhs_x,
524                                       const Number* rhs_g,
525                                       Number* sol_x, Number* sol_g)
526  {
527    // Create a vector for the RHS
528    SmartPtr<CompoundVector> rhs = comp_vec_space_->MakeNewCompoundVector();
529    SmartPtr<Vector> vrhs_x = rhs->GetCompNonConst(0);
530    SmartPtr<Vector> vrhs_g = rhs->GetCompNonConst(1);
531    // Now fill this vector with the values, extracting the relevant entries
532    if (rhs_x) {
533      SmartPtr<DenseVector> drhs_x =
534        dynamic_cast<DenseVector*> (GetRawPtr(vrhs_x));
535      Number* xvals = drhs_x->Values();
536      for (Index i=0; i<n_; i++) {
537        if (x_free_map_[i]>=0) {
538          xvals[x_free_map_[i]] = rhs_x[i];
539        }
540      }
541    }
542    else {
543      vrhs_x->Set(0.);
544    }
545    if (rhs_g) {
546      SmartPtr<DenseVector> drhs_g =
547        dynamic_cast<DenseVector*> (GetRawPtr(vrhs_g));
548      Number* gvals = drhs_g->Values();
549      for (Index j=0; j<m_; j++) {
550        if (g_fixed_map_[j]>=0) {
551          gvals[g_fixed_map_[j]] = rhs_g[j];
552        }
553      }
554    }
555    else {
556      vrhs_g->Set(0.);
557    }
558
559    // Solve the linear system
560    SmartPtr<CompoundVector> sol = comp_vec_space_->MakeNewCompoundVector();
561    ESymSolverStatus solver_retval =
562      tsymlinearsolver_->Solve(*comp_proj_matrix_, *rhs, *sol, false, 0);
563    if (solver_retval != SYMSOLVER_SUCCESS) {
564      return false;
565    }
566
567    // and get the solution out of it
568    if (sol_x) {
569      SmartPtr<Vector> vsol_x = sol->GetCompNonConst(0);
570      SmartPtr<const DenseVector> dsol_x =
571        dynamic_cast<const DenseVector*> (GetRawPtr(vsol_x));
572      const Number* xvals = dsol_x->Values();
573      for (Index i=0; i<n_; i++) {
574        if (x_free_map_[i]>=0) {
575          sol_x[i] = xvals[x_free_map_[i]];
576        }
577        else {
578          sol_x[i] = 0.;
579        }
580      }
581    }
582    if (sol_g) {
583      SmartPtr<Vector> vsol_g = sol->GetCompNonConst(1);
584      SmartPtr<const DenseVector> dsol_g =
585        dynamic_cast<const DenseVector*> (GetRawPtr(vsol_g));
586      const Number* gvals = dsol_g->Values();
587      for (Index j=0; j<m_; j++) {
588        if (g_fixed_map_[j]>=0) {
589          sol_g[j] = gvals[g_fixed_map_[j]];
590        }
591        else {
592          sol_g[j] = 0.;
593        }
594      }
595    }
596
597    return true;
598  }
599
600  bool CurvatureEstimator::Compute_dTHLagd(
601    const Number* d, const Number* x, bool new_x, const Number* lambda,
602    bool new_lambda,  Number& dTHLagd)
603  {
604    if (new_x || new_lambda || !hess_vals_) {
605      hess_vals_ = new Number[nnz_hess_];
606      if (!tnlp_->eval_h(n_, x, new_x, 1., m_, lambda, new_lambda, nnz_hess_,
607                         NULL, NULL, hess_vals_)) {
608        return false;
609      }
610    }
611    dTHLagd = 0.;
612    for (Index i=0; i<nnz_hess_; i++) {
613      Index irow = irows_hess_[i];
614      Index jcol = jcols_hess_[i];
615      if (irow == jcol) {
616        dTHLagd += d[irow]*d[irow]*hess_vals_[i];
617      }
618      else {
619        dTHLagd += 2.*d[irow]*d[jcol]*hess_vals_[i];
620      }
621    }
622    return true;
623  }
624
625
626} // namespace Bonmin
Note: See TracBrowser for help on using the repository browser.