Ignore:
Timestamp:
Oct 30, 2006 4:29:42 PM (13 years ago)
Author:
andreasw
Message:

first version of BonChooseVariable? for Hessian-based strong branching; minor corrections

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/devel/Bonmin/src/Interfaces/BonCurvatureEstimator.cpp

    r66 r98  
    5050      options_(options),
    5151      prefix_(""),
    52       tnlp_(tnlp)
     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)
    5368  {
    5469    DBG_ASSERT(IsValid(jnlst));
     
    153168  {
    154169    if (initialized_) {
     170      delete [] grad_f_;
    155171      delete [] irows_jac_;
    156172      delete [] jcols_jac_;
     
    172188  {
    173189    DBG_ASSERT(!initialized_);
     190
    174191    //////////////////////////////////////
    175192    // Prepare internal data structures //
     
    181198      return false;
    182199    }
     200
     201    // Space for gradient
     202    grad_f_ = new Number[n_];
    183203
    184204    // Get nonzero entries in the matrices
     
    232252  CurvatureEstimator::ComputeNullSpaceCurvature(
    233253    bool new_bounds,
    234     std::vector<int>& active_x,
    235     std::vector<int>& active_g,
    236     bool new_activities,
    237254    int n,
    238255    const Number* x,
    239256    bool new_x,
     257    const Number* z_L,
     258    const Number* z_U,
     259    int m,
     260    const Number* lam,
     261    bool new_mults,
    240262    const Number* orig_d,
    241263    Number* projected_d,
    242     Number& dTHd)
    243   {
    244     DBG_ASSERT(n == n_);
    245 
     264    Number& gradLagTd,
     265    Number& dTHLagd)
     266  {
    246267    if (!initialized_) {
    247268      Initialize();
    248     }
    249 
    250     if (new_bounds) new_activities = true;
     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    }
    251333
    252334    // Check if the structure of the matrix has changed
    253335    if (new_activities) {
    254       if (!PrepareNewMatrixStructure(new_bounds, active_x, active_g,
    255                                      new_activities)) {
     336      if (!PrepareNewMatrixStructure(new_activities)) {
    256337        return false;
    257338      }
     
    265346
    266347      // Compute least square multipliers for the given activities
    267       Number* grad_f = new Number[n_];
    268       if (!tnlp_->eval_grad_f(n_, x, new_x, grad_f)) {
     348      if (!tnlp_->eval_grad_f(n_, x, new_x, grad_f_)) {
    269349        return false;
    270350      }
    271       if (!SolveSystem(grad_f, NULL, NULL, lambda_)) {
     351      if (!SolveSystem(grad_f_, NULL, NULL, lambda_)) {
    272352        return false;
    273353      }
    274       delete [] grad_f;
    275       IpBlasDscal(n_, -1., lambda_, 1);
     354      IpBlasDscal(m_, -1., lambda_, 1);
    276355      new_lambda = true;
    277356    }
     
    281360      return false;
    282361    }
    283     if (!Compute_dTHd(projected_d, x, new_x, lambda_, new_lambda, dTHd)) {
     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)) {
    284373      return false;
    285374    }
     
    289378
    290379  bool CurvatureEstimator::PrepareNewMatrixStructure(
    291     bool new_bounds,
    292     std::vector<int>& active_x,
    293     std::vector<int>& active_g,
    294380    bool new_activities)
    295381  {
    296     if (new_bounds) {
    297       // Get bounds
    298       if (!tnlp_->get_bounds_info(n_, x_l_, x_u_, m_, g_l_, g_u_)) {
    299         return false;
    300       }
    301     }
    302 
    303382    if (new_activities) {
    304383      // Deterimine which variables are free
     
    307386      }
    308387      // fixed by activities
    309       for (std::vector<int>::iterator i=active_x.begin();
    310            i != active_x.end(); i++) {
     388      for (std::vector<int>::iterator i=active_x_.begin();
     389           i != active_x_.end(); i++) {
    311390        DBG_ASSERT(*i != 0 && *i<=n_ && *i>=-n_);
    312391        if (*i<0) {
     
    340419      }
    341420      // fixed by activities
    342       for (std::vector<int>::iterator i=active_g.begin();
    343            i != active_g.end(); i++) {
     421      for (std::vector<int>::iterator i=active_g_.begin();
     422           i != active_g_.end(); i++) {
    344423        DBG_ASSERT(*i != 0 && *i<=m_ && *i>=-m_);
    345424        if (*i<0) {
     
    349428        else {
    350429          g_fixed_map_[(*i)-1] = 0;
    351           DBG_ASSERT(g_u_[(-*i)-1] < 1e19);
     430          DBG_ASSERT(g_u_[(*i)-1] < 1e19);
    352431        }
    353432      }
     
    361440      // fixed constraints
    362441      Index ng_fixed_ = 0;
    363       for (Index j=0; j<n_; j++) {
     442      for (Index j=0; j<m_; j++) {
    364443        if (g_fixed_map_[j] == 0) {
    365444          g_fixed_map_[j] = ng_fixed_;
     
    374453      Index nnz_proj_jac = 0;
    375454      for (Index i=0; i<nnz_jac_; i++) {
    376         Index irow = irows_jac_[i];
    377         Index jcol = jcols_jac_[i];
     455        const Index &irow = irows_jac_[i];
     456        const Index &jcol = jcols_jac_[i];
    378457        if (x_free_map_[jcol] >= 0 && g_fixed_map_[irow] >= 0) {
    379           iRows[nnz_proj_jac] = g_fixed_map_[irow];
    380           jCols[nnz_proj_jac] = x_free_map_[jcol];
     458          iRows[nnz_proj_jac] = g_fixed_map_[irow] + 1;
     459          jCols[nnz_proj_jac] = x_free_map_[jcol] + 1;
    381460          nnz_proj_jac++;
    382461        }
     
    392471      comp_proj_matrix_space_ =
    393472        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_);
    394475      SmartPtr<SymMatrixSpace> identity_space =
    395476        new IdentityMatrixSpace(nx_free_);
    396477      comp_proj_matrix_space_->SetCompSpace(0, 0, *identity_space, true);
    397478      comp_proj_matrix_space_->SetCompSpace(1, 0, *proj_jac_space, true);
    398       SmartPtr<MatrixSpace> zero_space =
    399         new ZeroMatrixSpace(ng_fixed_, ng_fixed_);
    400       comp_proj_matrix_space_->SetCompSpace(1, 1, *zero_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);
    401482
    402483      // Create a Vector space for the rhs and sol
     
    416497    bool new_x)
    417498  {
    418     // If necessary, get new Jacobian values (for the original matrix)
    419     if (new_x) {
    420       if (!tnlp_->eval_jac_g(n_, x, new_x, m_, nnz_jac_,
    421                              NULL, NULL, jac_vals_)) {
    422         return false;
    423       }
    424     }
    425 
    426499    if (new_x || new_activities) {
    427500      comp_proj_matrix_ = comp_proj_matrix_space_->MakeNewCompoundSymMatrix();
     
    437510        }
    438511      }
     512      comp_proj_matrix_->Print(*jnlst_,J_WARNING,J_MAIN,"compmat");
    439513      DBG_ASSERT(inz == tjac->Nonzeros());
    440514
     
    524598  }
    525599
    526   bool CurvatureEstimator::Compute_dTHd(
     600  bool CurvatureEstimator::Compute_dTHLagd(
    527601    const Number* d, const Number* x, bool new_x, const Number* lambda,
    528     bool new_lambda,  Number& dTHd)
     602    bool new_lambda,  Number& dTHLagd)
    529603  {
    530604    if (new_x || new_lambda || !hess_vals_) {
     
    535609      }
    536610    }
    537     dTHd = 0.;
     611    dTHLagd = 0.;
    538612    for (Index i=0; i<nnz_hess_; i++) {
    539613      Index irow = irows_hess_[i];
    540614      Index jcol = jcols_hess_[i];
    541615      if (irow == jcol) {
    542         dTHd += d[irow]*d[irow]*hess_vals_[i];
     616        dTHLagd += d[irow]*d[irow]*hess_vals_[i];
    543617      }
    544618      else {
    545         dTHd += 2.*d[irow]*d[jcol]*hess_vals_[i];
     619        dTHLagd += 2.*d[irow]*d[jcol]*hess_vals_[i];
    546620      }
    547621    }
Note: See TracChangeset for help on using the changeset viewer.