source: trunk/Bonmin/src/Interfaces/Ampl/BonAmplTMINLP.cpp @ 1536

Last change on this file since 1536 was 1536, checked in by pbonami, 10 years ago

Port Claudia's code to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to "Author Date Id Revision"
File size: 26.1 KB
Line 
1// (C) Copyright International Business Machines Corporation and
2// Carnegie Mellon University 2004, 2007
3//
4// All Rights Reserved.
5// This code is published under the Common Public License.
6//
7// Authors :
8// Carl D. Laird, Carnegie Mellon University,
9// Andreas Waechter, International Business Machines Corporation
10// Pierre Bonami, Carnegie Mellon University,
11//
12// Date : 12/01/2004
13#include "IpBlas.hpp"
14
15#include <list>
16
17#include "AmplTNLP.hpp"
18#include "BonAmplTMINLP.hpp"
19#include <iostream>
20#include <fstream>
21
22#include "CoinHelperFunctions.hpp"
23#include "BonExitCodes.hpp"
24namespace Bonmin{
25  void 
26  RegisteredOptions::fillAmplOptionList(ExtraCategoriesInfo which, Ipopt::AmplOptionsList * amplOptList){
27      std::list<int> test;
28      std::list< Ipopt::RegisteredOption * > options;
29      chooseOptions(which, options);
30      for(std::list< Ipopt::RegisteredOption * >::iterator i = options.begin();
31           i != options.end() ; i++)
32       {
33           std::string name = "bonmin.";
34           name += (*i)->Name();
35           Ipopt::RegisteredOptionType T = (*i)->Type();
36           Ipopt::AmplOptionsList::AmplOptionType  type;
37           switch(T){
38             case Ipopt::OT_Number: type = Ipopt::AmplOptionsList::Number_Option;
39                  break;
40             case Ipopt::OT_Integer: type = Ipopt::AmplOptionsList::Integer_Option;
41                  break;
42             case Ipopt::OT_String: type = Ipopt::AmplOptionsList::String_Option;
43                  break;
44             case Ipopt::OT_Unknown:
45             default:
46                throw CoinError("RegisteredOptions","fillAmplOptionList","Unknown option type");
47           }
48           amplOptList->AddAmplOption(name, name, type, (*i)->ShortDescription());
49       }
50}
51}
52#include "asl.h"
53#include "asl_pfgh.h"
54#include "getstub.h"
55
56namespace ampl_utils
57{
58  void sos_kludge(int nsos, int *sosbeg, double *sosref);
59}
60namespace Bonmin
61{
62
63  AmplTMINLP::AmplTMINLP()
64      :
65      TMINLP(),
66      appName_(),
67      upperBoundingObj_(-1),
68      ampl_tnlp_(NULL),
69      branch_(),
70      sos_(),
71      suffix_handler_(NULL),
72      constraintsConvexities_(NULL),
73      numberNonConvex_(0),
74      nonConvexConstraintsAndRelaxations_(NULL),
75      numberSimpleConcave_(0),
76      simpleConcaves_(NULL),
77      hasLinearObjective_(false)
78  {}
79
80
81  AmplTMINLP::AmplTMINLP(const SmartPtr<const Journalist>& jnlst,
82      const SmartPtr<Bonmin::RegisteredOptions> roptions,
83      const SmartPtr<OptionsList> options,
84      char**& argv,
85      AmplSuffixHandler* suffix_handler /*=NULL*/,
86      const std::string& appName,
87      std::string* nl_file_content /* = NULL */
88                        )
89      :
90      TMINLP(),
91      appName_(),
92      upperBoundingObj_(-1),
93      ampl_tnlp_(NULL),
94      branch_(),
95      sos_(),
96      suffix_handler_(NULL),
97      constraintsConvexities_(NULL),
98      numberNonConvex_(0),
99      nonConvexConstraintsAndRelaxations_(NULL),
100      numberSimpleConcave_(0),
101      simpleConcaves_(NULL),
102      hasLinearObjective_(false)
103  {
104    Initialize(jnlst, roptions, options, argv, suffix_handler, appName, nl_file_content);
105  }
106
107  void
108  AmplTMINLP::Initialize(const SmartPtr<const Journalist>& jnlst,
109      const SmartPtr<Bonmin::RegisteredOptions> roptions,
110      const SmartPtr<OptionsList> options,
111      char**& argv,
112      AmplSuffixHandler* suffix_handler /*=NULL*/,
113      const std::string& appName,
114      std::string* nl_file_content /* = NULL */
115                        )
116  {
117    appName_ = appName;
118    options->GetEnumValue("file_solution",writeAmplSolFile_,"bonmin.");
119    jnlst_ = jnlst;
120
121    if (suffix_handler==NULL)
122      suffix_handler_ = suffix_handler = new AmplSuffixHandler();
123
124    // Add the suffix handler for scaling
125    suffix_handler->AddAvailableSuffix("scaling_factor", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
126    suffix_handler->AddAvailableSuffix("scaling_factor", AmplSuffixHandler::Constraint_Source, AmplSuffixHandler::Number_Type);
127    suffix_handler->AddAvailableSuffix("scaling_factor", AmplSuffixHandler::Objective_Source, AmplSuffixHandler::Number_Type);
128
129    // Modified for warm-start from AMPL
130    suffix_handler->AddAvailableSuffix("ipopt_zL_out", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
131    suffix_handler->AddAvailableSuffix("ipopt_zU_out", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
132    suffix_handler->AddAvailableSuffix("ipopt_zL_in", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
133    suffix_handler->AddAvailableSuffix("ipopt_zU_in", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
134
135    // priority suffix
136    suffix_handler->AddAvailableSuffix("priority", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Index_Type);
137    suffix_handler->AddAvailableSuffix("direction", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
138    suffix_handler->AddAvailableSuffix("downPseudocost", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
139    suffix_handler->AddAvailableSuffix("upPseudocost", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
140
141
142
143    // sos suffixes
144    suffix_handler->AddAvailableSuffix("ref", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
145    suffix_handler->AddAvailableSuffix("sos",AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Index_Type);
146    suffix_handler->AddAvailableSuffix("sos",AmplSuffixHandler::Constraint_Source, AmplSuffixHandler::Index_Type);
147    suffix_handler->AddAvailableSuffix("sosno",AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
148    suffix_handler->AddAvailableSuffix("sosref",AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
149    suffix_handler->AddAvailableSuffix("sstatus", AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Index_Type);
150    suffix_handler->AddAvailableSuffix("sstatus", AmplSuffixHandler::Constraint_Source, AmplSuffixHandler::Index_Type);
151
152
153    // For marking convex/nonconvex constraints
154    suffix_handler->AddAvailableSuffix("id",AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Index_Type);
155    suffix_handler->AddAvailableSuffix("primary_var",AmplSuffixHandler::Constraint_Source, AmplSuffixHandler::Index_Type);
156
157    // For objectives
158    suffix_handler->AddAvailableSuffix("UBObj", AmplSuffixHandler::Objective_Source, AmplSuffixHandler::Index_Type);
159
160
161    // Perturbation radius
162    suffix_handler->AddAvailableSuffix("perturb_radius",AmplSuffixHandler::Variable_Source, AmplSuffixHandler::Number_Type);
163
164    SmartPtr<AmplOptionsList> ampl_options_list = new AmplOptionsList();
165    roptions->fillAmplOptionList(RegisteredOptions::BonminCategory, GetRawPtr(ampl_options_list));
166    roptions->fillAmplOptionList(RegisteredOptions::FilterCategory, GetRawPtr(ampl_options_list));
167    roptions->fillAmplOptionList(RegisteredOptions::BqpdCategory, GetRawPtr(ampl_options_list));
168    fillApplicationOptions(GetRawPtr(ampl_options_list) );
169    std::string options_id = appName + "_options";
170    ampl_tnlp_ = new AmplTNLP(jnlst, options, argv, suffix_handler, true,
171        ampl_options_list, options_id.c_str(),
172        appName.c_str(), appName.c_str(), nl_file_content);
173    /* Read suffixes */
174    read_obj_suffixes();
175    read_priorities();
176    read_convexities();
177    read_sos();
178
179
180    /* Determine if objective is linear.*/
181    Index n_non_linear_b= 0;
182    Index n_non_linear_bi= 0;
183    Index n_non_linear_c= 0;
184    Index n_non_linear_ci = 0;
185    Index n_non_linear_o= 0;
186    Index n_non_linear_oi = 0;
187    Index n_binaries = 0;
188    Index n_integers = 0;
189    ampl_tnlp_->get_discrete_info(n_non_linear_b, n_non_linear_bi, n_non_linear_c,
190        n_non_linear_ci, n_non_linear_o, n_non_linear_oi,
191        n_binaries, n_integers);
192    if (n_non_linear_b == 0 && n_non_linear_o == 0) {
193      hasLinearObjective_ = true;
194    }
195  }
196
197  AmplTMINLP::~AmplTMINLP()
198  {
199    delete [] constraintsConvexities_;
200    delete [] simpleConcaves_;
201    delete [] nonConvexConstraintsAndRelaxations_;
202    delete ampl_tnlp_;
203  }
204
205  void
206  AmplTMINLP::read_priorities()
207  {
208    int numcols, m, dummy1, dummy2;
209    TNLP::IndexStyleEnum index_style;
210    ampl_tnlp_->get_nlp_info(numcols, m, dummy1, dummy2, index_style);
211
212    const AmplSuffixHandler * suffix_handler = GetRawPtr(suffix_handler_);
213
214    const Index* pri = suffix_handler->GetIntegerSuffixValues("priority", AmplSuffixHandler::Variable_Source);
215    const Index* brac = suffix_handler->GetIntegerSuffixValues("direction", AmplSuffixHandler::Variable_Source);
216    const Number* upPs = suffix_handler->GetNumberSuffixValues("upPseudocost", AmplSuffixHandler::Variable_Source);
217    const Number* dwPs = suffix_handler->GetNumberSuffixValues("downPseudocost", AmplSuffixHandler::Variable_Source);
218
219
220    branch_.gutsOfDestructor();
221    branch_.size = numcols;
222    if (pri) {
223      branch_.priorities = new int[numcols];
224      for (int i = 0 ; i < numcols ; i++) {
225        branch_.priorities [i] = -pri[i] + 9999;
226      }
227    }
228    if (brac) {
229      branch_.branchingDirections = CoinCopyOfArray(brac,numcols);
230    }
231    if (upPs && !dwPs) dwPs = upPs;
232    else if (dwPs && !upPs) upPs = dwPs;
233
234    if (upPs) {
235      branch_.upPsCosts = CoinCopyOfArray(upPs,numcols);
236    }
237    if (dwPs) {
238      branch_.downPsCosts = CoinCopyOfArray(dwPs,numcols);
239    }
240
241    const double* perturb_radius =
242      suffix_handler->GetNumberSuffixValues("perturb_radius", AmplSuffixHandler::Variable_Source);
243    perturb_info_.SetPerturbationArray(numcols, perturb_radius);
244  }
245
246  void
247  AmplTMINLP::read_sos()
248  {
249    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
250
251    int i = 0;//ASL_suf_sos_explict_free;
252    int copri[2], **p_sospri;
253    copri[0] = 0;
254    copri[1] = 0;
255    int * starts = NULL;
256    int * indices = NULL;
257    char * types = NULL;
258    double * weights = NULL;
259    int * priorities = NULL;
260    p_sospri = &priorities;
261
262    sos_.gutsOfDestructor();
263
264    sos_.num = suf_sos(i, &sos_.numNz, &types, p_sospri, copri,
265        &starts, &indices, &weights);
266    if (sos_.num) {
267      //Copy sos information
268      sos_.priorities = CoinCopyOfArray(priorities,sos_.num);
269      sos_.starts = CoinCopyOfArray(starts, sos_.num + 1);
270      sos_.indices = CoinCopyOfArray(indices, sos_.numNz);
271      sos_.types = CoinCopyOfArray(types, sos_.num);
272      sos_.weights = CoinCopyOfArray(weights, sos_.numNz);
273
274      ampl_utils::sos_kludge(sos_.num, sos_.starts, sos_.weights);
275      for (int ii=0;ii<sos_.num;ii++) {
276        int ichar = sos_.types[ii];
277        if (ichar != '1' && ichar != '2') {
278          std::cerr<<"Unsuported type of sos constraint: "<<sos_.types[ii]<<std::endl;
279          throw;
280        }
281        sos_.types[ii]= 1;
282      }
283    }
284  }
285
286  void
287  AmplTMINLP::read_obj_suffixes()
288  {
289    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
290    DBG_ASSERT(asl);
291    //Get the values
292    if (n_obj < 2) return;
293
294    const AmplSuffixHandler * suffix_handler = GetRawPtr(suffix_handler_);
295
296    const Index* UBObj = suffix_handler->GetIntegerSuffixValues("UBObj", AmplSuffixHandler::Objective_Source);
297    if (UBObj) {
298      //get index of lower bounding objective
299      int lowerBoundingObj = (UBObj[0] == 1)? 1:0;
300      // Pass information to Ipopt
301      ampl_tnlp_->set_active_objective(lowerBoundingObj);
302
303      //get the index of upper bounding objective
304      for (int i = 0; i < n_obj; i++) {
305        if (UBObj[i]==1) {
306          if (upperBoundingObj_ != -1) {
307            jnlst_->Printf(J_ERROR, J_MAIN,
308                "Too many objectives for upper-bounding");
309          }
310          upperBoundingObj_ = i;
311        }
312      }
313    }
314    else {
315      ampl_tnlp_->set_active_objective(0);
316    }
317  }
318  /** To store all data stored in the nonconvex suffixes.*/
319  struct NonConvexSuff
320  {
321    /** Default constructor.*/
322    NonConvexSuff():
323        cIdx(-1),relIdx(-1),scXIdx(-1),scYIdx(-1)
324    {}
325    /** Copy constructor.*/
326    NonConvexSuff(const NonConvexSuff& other):
327        cIdx(other.cIdx), relIdx(other.relIdx),
328        scXIdx(other.scXIdx), scYIdx(other.scYIdx)
329    {}
330    /** Index of the nonconvex constraint.*/
331    int cIdx;
332    /** Index of its relaxation.*/
333    int relIdx;
334    /** Index of variable x in a simple concave constraint of type y >= F(x).*/
335    int scXIdx;
336    /** Index of variable y in a simple concave constraint of type y >= F(x).*/
337    int scYIdx;
338  };
339  void AmplTMINLP::read_convexities()
340  {
341    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
342    DBG_ASSERT(asl);
343
344    const AmplSuffixHandler * suffix_handler = GetRawPtr(suffix_handler_);
345    const Index * id = suffix_handler->GetIntegerSuffixValues("id", AmplSuffixHandler::Variable_Source);
346    const Index * primary_var = suffix_handler->GetIntegerSuffixValues("primary_var", AmplSuffixHandler::Constraint_Source);
347
348
349    if (primary_var!= NULL) {
350      if (constraintsConvexities_ != NULL) {
351        delete [] constraintsConvexities_;
352      }
353      constraintsConvexities_ = new TMINLP::Convexity[n_con];
354      if (id == NULL) {
355        std::cerr<<"Incorrect suffixes description in ampl model. id's are not declared "<<std::endl;
356        exit(ERROR_IN_AMPL_SUFFIXES);
357      }
358      int numberSimpleConcave = 0;
359      std::map<int, int> id_map;
360
361      for (int i = 0 ; i < n_var ; i++) {
362        id_map[id[i]] = i;
363      }
364
365
366      for (int i = 0 ; i < n_con ; i++) {
367        if (primary_var[i] != 0) {
368          constraintsConvexities_[i] = TMINLP::SimpleConcave;
369          numberSimpleConcave++;
370        }
371        else constraintsConvexities_[i] = TMINLP::Convex;
372      }
373      simpleConcaves_ = new SimpleConcaveConstraint[numberSimpleConcave];
374      nonConvexConstraintsAndRelaxations_ = new MarkedNonConvex[numberSimpleConcave];
375      numberSimpleConcave = 0;
376      int * jCol = new int[n_var];
377      for (int i = 0 ; i < n_con ; i++) {
378        if (primary_var[i] != 0) {
379          nonConvexConstraintsAndRelaxations_[numberSimpleConcave].cIdx = i;
380          nonConvexConstraintsAndRelaxations_[numberSimpleConcave].cRelaxIdx = -1;
381          simpleConcaves_[numberSimpleConcave].cIdx = i;
382          simpleConcaves_[numberSimpleConcave].yIdx = id_map[primary_var[i]];
383
384          //Now get gradient of i to get xIdx.
385          int nnz;
386          int & yIdx = simpleConcaves_[numberSimpleConcave].yIdx;
387          int & xIdx = simpleConcaves_[numberSimpleConcave].xIdx;
388          eval_grad_gi(n_var, NULL, false, i, nnz, jCol, NULL);
389          if (nnz != 2) {//Error in ampl model
390            std::cout<<"Incorrect suffixes description in ampl model. Constraint with id "
391            <<id<<" is simple concave and should have only two nonzero elements"<<std::endl;
392            exit(ERROR_IN_AMPL_SUFFIXES);
393          }
394          if (jCol[0] - 1== yIdx) {
395            xIdx = jCol[1] - 1;
396          }
397          else {
398            if (jCol[1] - 1!= yIdx) {//Error in ampl model
399              std::cout<<"Incorrect suffixes description in ampl model. Constraint with id "
400              <<id<<" : variable marked as y does not appear in the constraint."<<std::endl;
401              exit(ERROR_IN_AMPL_SUFFIXES);
402            }
403            xIdx = jCol[0] - 1;
404          }
405          numberSimpleConcave++;
406        }
407      }
408      delete [] jCol;
409      numberSimpleConcave_ = numberSimpleConcave;
410      numberNonConvex_ = numberSimpleConcave;
411    }
412
413  }
414
415  bool AmplTMINLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, Index& nnz_h_lag, TNLP::IndexStyleEnum& index_style)
416  {
417    return ampl_tnlp_->get_nlp_info(n, m, nnz_jac_g, nnz_h_lag, index_style);
418  }
419
420  bool AmplTMINLP::get_variables_types(Index n, VariableType* var_types)
421  {
422    // The variables are sorted by type in AMPL, so all we need to
423    // know are the counts of each type.
424
425
426    Index n_non_linear_b= 0;
427    Index n_non_linear_bi= 0;
428    Index n_non_linear_c= 0;
429    Index n_non_linear_ci = 0;
430    Index n_non_linear_o= 0;
431    Index n_non_linear_oi = 0;
432    Index n_binaries = 0;
433    Index n_integers = 0;
434    ampl_tnlp_->get_discrete_info(n_non_linear_b, n_non_linear_bi, n_non_linear_c,
435        n_non_linear_ci, n_non_linear_o, n_non_linear_oi,
436        n_binaries, n_integers);
437    int totalNumberOfNonContinuous = 0;
438
439
440    //begin with variables non-linear both in the objective function and in some constraints
441    // The first ones are continuous
442    Index start = 0;
443    Index end = n_non_linear_b - n_non_linear_bi;
444    for (int i=start; i<end; i++) {
445      var_types[i] = CONTINUOUS;
446    }
447
448    // The second ones are integers
449    start = end;
450    end = start + n_non_linear_bi;
451    for (int i=start; i < end; i++) {
452      var_types[i] = INTEGER;
453      totalNumberOfNonContinuous++;
454    }
455    //next variables non-linear in some constraints
456    // The first ones are continuous
457    start = end;
458    end = std::max(start,end + n_non_linear_c - n_non_linear_ci - n_non_linear_b);
459    for (int i=start; i<end; i++) {
460      var_types[i] = CONTINUOUS;
461    }
462
463    // The second ones are integers
464    start = end;
465    end = start + n_non_linear_ci;
466    for (int i=start; i < end; i++) {
467      var_types[i] = INTEGER;
468      totalNumberOfNonContinuous++;
469    }
470
471    //next variables non-linear in the objective function
472    // The first ones are continuous
473    start = end;
474    end = std::max(start,end + n_non_linear_o - std::max(n_non_linear_b, n_non_linear_c) - n_non_linear_oi);
475    for (int i=start; i<end; i++) {
476      var_types[i] = CONTINUOUS;
477    }
478
479    // The second ones are integers
480    start = end;
481    end = start + n_non_linear_oi;
482    for (int i=start; i < end; i++) {
483      var_types[i] = INTEGER;
484      totalNumberOfNonContinuous++;
485    }
486
487    //At last the linear variables
488    // The first ones are continuous
489    start = end;
490    end = n - n_binaries - n_integers;
491    for (int i=start; i<end; i++) {
492      var_types[i] = CONTINUOUS;
493    }
494
495    // The second ones are binaries
496    start = end;
497    end = start + n_binaries;
498    for (int i=start; i < end; i++) {
499      var_types[i] = BINARY;
500      totalNumberOfNonContinuous++;
501    }
502
503    // The third ones are integers
504    start = end;
505    end = start + n_integers;
506    for (int i=start; i < end; i++) {
507      var_types[i] = INTEGER;
508      totalNumberOfNonContinuous++;
509    }
510    return true;
511  }
512
513  bool AmplTMINLP::get_variables_linearity(Index n, Ipopt::TNLP::LinearityType* var_types)
514  {
515    // The variables are sorted by type in AMPL, so all we need to
516    // know are the counts of each type.
517
518
519    Index n_non_linear_b= 0;
520    Index n_non_linear_bi= 0;
521    Index n_non_linear_c= 0;
522    Index n_non_linear_ci = 0;
523    Index n_non_linear_o= 0;
524    Index n_non_linear_oi = 0;
525    Index n_binaries = 0;
526    Index n_integers = 0;
527    ampl_tnlp_->get_discrete_info(n_non_linear_b, n_non_linear_bi, n_non_linear_c,
528        n_non_linear_ci, n_non_linear_o, n_non_linear_oi,
529        n_binaries, n_integers);
530
531    //Compute the number of non linear variables:
532    int n_non_linear = n_non_linear_c + n_non_linear_o - n_non_linear_b;
533
534    int start = 0;
535    int end = n_non_linear;
536    for (int i=start; i<end; i++) {
537      var_types[i] = Ipopt::TNLP::NON_LINEAR;
538    }
539
540    //At last the linear variables
541    // The first ones are continuous
542    start = end;
543    end = n;
544    for (int i=start; i<end; i++) {
545      var_types[i] = Ipopt::TNLP::LINEAR;
546    }
547    return true;
548  }
549
550  /** Returns the constraint linearity.
551  * array should be alocated with length at least n.*/
552  bool
553  AmplTMINLP::get_constraints_linearity(Index m,
554      Ipopt::TNLP::LinearityType* const_types)
555  {
556    return ampl_tnlp_->get_constraints_linearity(m, const_types);
557  }
558
559  bool AmplTMINLP::get_bounds_info(Index n, Number* x_l, Number* x_u,
560      Index m, Number* g_l, Number* g_u)
561  {
562    return ampl_tnlp_->get_bounds_info(n, x_l, x_u, m, g_l, g_u);
563  }
564
565  bool AmplTMINLP::get_starting_point(Index n, bool init_x, Number* x,
566      bool init_z, Number* z_L, Number* z_U,
567      Index m, bool init_lambda, Number* lambda)
568  {
569    return ampl_tnlp_->get_starting_point(n, init_x, x,
570        init_z, z_L, z_U,
571        m, init_lambda, lambda);
572  }
573
574  bool AmplTMINLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
575  {
576    return ampl_tnlp_->eval_f(n, x, new_x, obj_value);
577  }
578
579  bool AmplTMINLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
580  {
581    return ampl_tnlp_->eval_grad_f(n, x, new_x, grad_f);
582  }
583
584  bool AmplTMINLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
585  {
586    return ampl_tnlp_->eval_g(n, x, new_x, m, g);
587  }
588
589  bool AmplTMINLP::eval_jac_g(Index n, const Number* x, bool new_x,
590      Index m, Index nele_jac, Index* iRow,
591      Index *jCol, Number* values)
592  {
593    return ampl_tnlp_->eval_jac_g(n, x, new_x,
594        m, nele_jac, iRow, jCol,
595        values);
596  }
597
598  bool AmplTMINLP::eval_h(Index n, const Number* x, bool new_x,
599      Number obj_factor, Index m, const Number* lambda,
600      bool new_lambda, Index nele_hess, Index* iRow,
601      Index* jCol, Number* values)
602  {
603    return ampl_tnlp_->eval_h(n, x, new_x, obj_factor,
604        m, lambda, new_lambda, nele_hess, iRow,
605        jCol, values);
606  }
607
608  bool AmplTMINLP::eval_gi(Index n, const Number* x, bool new_x,
609      Index i, Number& gi)
610  {
611    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
612
613    // ignore new_x for now
614    xunknown();
615
616    fint nerror = 0;
617    gi = conival(i, const_cast<real*>(x), &nerror);
618    if (nerror!=0) {
619      return false;
620    }
621    else {
622      return true;
623    }
624  }
625
626  bool AmplTMINLP::eval_grad_gi(Index n, const Number* x, bool new_x,
627      Index i, Index& nele_grad_gi, Index* jCol,
628      Number* values)
629  {
630    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
631
632    if (jCol) {
633      // Only compute the number of nonzeros and the indices
634      DBG_ASSERT(!values);
635      nele_grad_gi = 0;
636      for (cgrad* cg=Cgrad[i]; cg; cg = cg->next) {
637        jCol[nele_grad_gi++] = cg->varno + 1;
638      }
639      return true;
640    }
641    DBG_ASSERT(values);
642
643    // ignore new_x for now
644    xunknown();
645
646    asl->i.congrd_mode = 1;
647    fint nerror = 0;
648    congrd(i, const_cast<real*>(x), values, &nerror);
649    if (nerror!=0) {
650      return false;
651    }
652    else {
653      return true;
654    }
655  }
656
657  void AmplTMINLP::finalize_solution(TMINLP::SolverReturn status,
658      Index n, const Number* x, Number obj_value)
659  {
660    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
661    std::string message;
662    std::string status_str;
663    if (status == TMINLP::SUCCESS) {
664      status_str = "\t\"Finished\"";
665      message = "\n" + appName_ +": Optimal";
666      solve_result_num = 3;
667    }
668    else if (status == TMINLP::INFEASIBLE) {
669      status_str = "\t\"Finished\"";
670      message = "\n" + appName_ + ": Infeasible problem";
671      solve_result_num = 220;
672    }
673    else if (status == TMINLP::CONTINUOUS_UNBOUNDED) {
674      status_str = "\t\"Finished\"";
675      message = "\n" + appName_ +" Continuous relaxation is unbounded.";
676      solve_result_num = 300;
677    }
678    else if (status == TMINLP::LIMIT_EXCEEDED) {
679      status_str = "\t\"Not finished\"";
680      message = "\n" + appName_ + ": Optimization interupted on limit.";
681      if(x)
682        solve_result_num = 421; /* Limit reached with integer feasible solution.*/
683      else
684        solve_result_num = 410; /* Limit reached without solution.*/
685    }
686    else if (status == TMINLP::MINLP_ERROR) {
687      status_str = "\t\"Aborted\"";
688      message = "\n" + appName_ + ": Error encountered in optimization.";
689      solve_result_num = 500;
690    }
691    if (writeAmplSolFile_) {
692      write_solution(message, x);
693      std::cout<<"\n "<<status_str<<std::endl;
694    }
695    else {
696      std::cout<<status_str<<message<<std::endl;
697      std::string fName = appName_ + ".sol";
698      std::ofstream of(fName.c_str());
699      for (int i = 0 ; i < n ; i++) {
700        of<<i<<"\t"<<x[i]<<std::endl;
701      }
702      of<<"-1\n";
703    }
704  }
705
706  void AmplTMINLP::write_solution(const std::string & message, const Number *x_sol)
707  {
708    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
709
710    DBG_ASSERT(asl);
711    //    DBG_ASSERT(x_sol);
712
713    // We need to copy the message into a non-const char array to make
714    // it work with the AMPL C function.
715    char* cmessage = new char[message.length()+1];
716    strcpy(cmessage, message.c_str());
717
718    // In order to avoid casting into non-const, we copy the data here...
719    double* x_sol_copy = NULL;
720    if (x_sol) {
721      x_sol_copy = new double[n_var];
722      for (int i=0; i<n_var; i++) {
723        x_sol_copy[i] = x_sol[i];
724      }
725    }
726    write_sol(cmessage, x_sol_copy, NULL, NULL);
727
728    delete [] x_sol_copy;
729    delete [] cmessage;
730
731  }
732
733
734  /** This methods gives the linear part of the objective function */
735  void AmplTMINLP::getLinearPartOfObjective(double * obj)
736  {
737    Index n, m, nnz_jac_g, nnz_h_lag;
738    TNLP::IndexStyleEnum index_style = TNLP::FORTRAN_STYLE;
739    get_nlp_info( n, m, nnz_jac_g, nnz_h_lag, index_style);
740    eval_grad_f(n, NULL, 0,obj);
741    Index n_non_linear_b= 0;
742    Index n_non_linear_bi= 0;
743    Index n_non_linear_c= 0;
744    Index n_non_linear_ci = 0;
745    Index n_non_linear_o= 0;
746    Index n_non_linear_oi = 0;
747    Index n_binaries = 0;
748    Index n_integers = 0;
749    ampl_tnlp_->get_discrete_info(n_non_linear_b, n_non_linear_bi, n_non_linear_c,
750        n_non_linear_ci, n_non_linear_o, n_non_linear_oi,
751        n_binaries, n_integers);
752
753    // Now get the coefficients of variables wich are linear in the objective
754    // The first ones are not
755    Index start = 0;
756    Index end = n_non_linear_b;
757    for (int i=start; i<end; i++) {
758      obj[i] = 0.;
759    }
760
761    //next variables should be linear in the objective just skip them
762    // (from current end to (end + n_non_linear_c - n_non_linear_ci - n_non_linear_b;)
763
764
765    // Those are nonlinear in the objective
766    start = end + n_non_linear_c;
767    end = start + n_non_linear_o;
768    for (int i=start; i < end; i++) {
769      obj[i]=0.;
770    }
771    //The rest is linear keep the values of the gradient
772  }
773
774
775
776  /** This method to returns the value of an alternative objective function for
777  upper bounding (if one has been declared by using the prefix UBObj).*/
778  bool
779  AmplTMINLP::eval_upper_bound_f(Index n, const Number* x,
780      Number& obj_value)
781  {
782    ASL_pfgh* asl = ampl_tnlp_->AmplSolverObject();
783    //xknown(x);    // This tells ampl to use a new x
784    fint nerror = -1;
785    obj_value = objval(upperBoundingObj_, const_cast<double *>(x), &nerror);
786    if (nerror > 0) {
787      jnlst_->Printf(J_ERROR, J_MAIN,
788          "Error in evaluating upper bounding objecting");
789      throw -1;
790    }
791    return nerror;
792
793  }
794  /** Return the ampl solver object (ASL*) */
795  const ASL_pfgh*
796  AmplTMINLP::AmplSolverObject() const
797  {
798    return ampl_tnlp_->AmplSolverObject();
799  }
800
801} // namespace Bonmin
Note: See TracBrowser for help on using the repository browser.