source: trunk/Apps/AmplSolver/AmplTNLP.cpp @ 2

Last change on this file since 2 was 2, checked in by andreasw, 15 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.6 KB
Line 
1// Copyright (C) 2004, International Business Machines and others.
2// All Rights Reserved.
3// This code is published under the Common Public License.
4//
5// $Id: AmplTNLP.cpp 2 2004-10-21 01:03:09Z andreasw $
6//
7// Authors:  Carl Laird, Andreas Waechter     IBM    2004-08-13
8
9// TODO:
10// - clean up the boolean for initialization
11// - pass in a tag so we can know when x has changed
12// - look closer at the interface - it should pass in non-zeros as well as m in the jacobian stuff (maybe hessian stuff)
13
14/* AMPL includes */
15extern "C"
16{
17#include "asl.h"
18#include "asl_pfgh.h"
19}
20
21#include "AmplTNLP.hpp"
22#include "IpDenseVector.hpp"
23#include "IpGenTMatrix.hpp"
24#include "IpSymTMatrix.hpp"
25#include "getstub.h"
26
27namespace Ipopt
28{
29
30  AmplTNLP::AmplTNLP(const SmartPtr<const Journalist>& jnlst, char**& argv, SmartPtr<AmplSuffixHandler> suffix_handler)
31      :
32      TNLP(),
33      jnlst_(jnlst),
34      asl_(NULL),
35      obj_sign_(1),
36      nz_h_full_(-1),
37      non_const_x_(NULL),
38      objval_called_with_current_x_(false),
39      conval_called_with_current_x_(false)
40  {
41
42    // The ASL include files #define certain
43    // variables that they expect you to work with.
44    // These variables then appear as though they are
45    // global variables when, in fact, they are not
46    // Most of them are data members of an asl object
47
48    // Create the ASL structure
49    ASL_pfgh* asl = (ASL_pfgh*)ASL_alloc(ASL_read_pfgh);
50    DBG_ASSERT(asl);
51    asl_ = asl; // keep the pointer for ourselves to use later...
52
53    // Read the options and stub
54    // ToDo: Figure out the options stuff
55    char* stub = getstub(&argv, NULL); // need to deal with options here
56    if (!stub) {
57      printf("No .nl file given!\n");
58      exit(1);
59    }
60
61    jnlst_->Printf(J_SUMMARY, J_MAIN, "Ampl Model: %s\n", stub);
62
63    // Parse the first part of the nl file
64    //char* stub = argv[0];
65    FILE* nl = jac0dim(stub, (fint)strlen(stub));
66    DBG_ASSERT(nl);
67    // check the problem statistics (see Table 1 in AMPL doc)
68    DBG_ASSERT(n_var > 0); // need some continuous variables
69    DBG_ASSERT(nbv == 0); // Cannot handle binary variables
70    DBG_ASSERT(niv == 0); // Cannot handle integer variables
71    // n_con can be >= 0
72    DBG_ASSERT(n_obj == 1); // Currently can handle only 1 objective
73    DBG_ASSERT(nlo == 0 || nlo == 1); // Can handle nonlinear obj.
74    DBG_ASSERT(nwv == 0); // Don't know what "linear arc" variables are
75    DBG_ASSERT(nlnc == 0); // Don't know what "nonlinear network"constraints are
76    DBG_ASSERT(lnc == 0); // Don't know what "linear network" constraints are
77
78    // Set options in the asl structure
79    want_xpi0 = 1 | 2;  // allocate initial values for primal and dual if available
80    DBG_ASSERT((want_xpi0 & 1) == 1 && (want_xpi0 & 2) == 2);
81    obj_no = 0; // always want to work with the first (and only?) objective
82
83    // allocate space for initial values
84    X0 = new real[n_var];
85    havex0 = new char[n_var];
86    pi0 = new real[n_con];
87    havepi0 = new char[n_con];
88
89    // prepare for suffixes
90    if (IsValid(suffix_handler)) {
91      suffix_handler->PrepareAmplForSuffixes(asl_);
92    }
93
94    // read the rest of the nl file
95    int retcode = pfgh_read(nl, ASL_return_read_err | ASL_findgroups);
96
97    // see "changes" in solvers directory of ampl code...
98    hesset(1,0,1,0,nlc);
99
100    switch (retcode) {
101      case ASL_readerr_none : {}
102      break;
103      case ASL_readerr_nofile : {
104        printf("Cannot open .nl file\n");
105        exit(1);
106      }
107      break;
108      case ASL_readerr_nonlin : {
109        DBG_ASSERT(false); // this better not be an error!
110        printf("model involves nonlinearities (ed0read)\n");
111        exit(1);
112      }
113      break;
114      case  ASL_readerr_argerr : {
115        printf("user-defined function with bad args\n");
116        exit(1);
117      }
118      break;
119      case ASL_readerr_unavail : {
120        printf("user-defined function not available\n");
121        exit(1);
122      }
123      break;
124      case ASL_readerr_corrupt : {
125        printf("corrupt .nl file\n");
126        exit(1);
127      }
128      break;
129      case ASL_readerr_bug : {
130        printf("bug in .nl reader\n");
131        exit(1);
132      }
133      break;
134      default: {
135        printf("Unknown error in stub file read\n");
136        exit(1);
137      }
138      break;
139    }
140
141    obj_sign_ = 1; // minimization
142    if (objtype[obj_no] != 0) {
143      obj_sign_ = -1;
144    }
145
146    // find the nonzero structure for the hessian
147    // parameters to sphsetup:
148    int coeff_obj = 1; // coefficient of the objective fn ???
149    int mult_supplied = 1; // multipliers will be supplied
150    int uptri = 1; // only need the upper triangular part
151    nz_h_full_ = sphsetup(-1, coeff_obj, mult_supplied, uptri);
152  }
153
154  AmplTNLP::~AmplTNLP()
155  {
156    ASL_pfgh* asl = asl_;
157
158    if (asl) {
159      if (X0) {
160        delete [] X0;
161        X0 = NULL;
162      }
163      if (havex0) {
164        delete [] havex0;
165        havex0 = NULL;
166      }
167      if (pi0) {
168        delete [] pi0;
169        pi0 = NULL;
170      }
171      if (havepi0) {
172        delete [] havepi0;
173        havepi0 = NULL;
174      }
175      ASL_free((ASL**)&asl_);
176      asl_ = NULL;
177    }
178
179    delete [] non_const_x_;
180  }
181
182  bool AmplTNLP::get_nlp_info(Index& n, Index& m, Index& nnz_jac_g, Index& nnz_h_lag)
183  {
184    ASL_pfgh* asl = asl_;
185    DBG_ASSERT(asl_);
186
187    n = n_var; // # of variables (variable types have been asserted in the constructor
188    m = n_con; // # of constraints
189    nnz_jac_g = nzc; // # of non-zeros in the jacobian
190    nnz_h_lag = nz_h_full_; // # of non-zeros in the hessian
191
192    return true;
193  }
194
195  bool AmplTNLP::get_bounds_info(Index n, Number* x_l, Number* x_u, Index m, Number* d_l, Number* d_u)
196  {
197    ASL_pfgh* asl = asl_;
198    DBG_ASSERT(asl_);
199
200    DBG_ASSERT(n == n_var);
201    DBG_ASSERT(m == n_con);
202
203    for (Index i=0; i<n; i++) {
204      x_l[i] = LUv[2*i];
205      x_u[i] = LUv[2*i+1];
206    }
207
208    for (Index i=0; i<m; i++) {
209      d_l[i] = LUrhs[2*i];
210      d_u[i] = LUrhs[2*i+1];
211    }
212
213    return true;
214  }
215
216  bool AmplTNLP::get_starting_point(Index n, bool init_x, Number* x, bool init_z, Number* z_L, Number* z_U, Index m, bool init_lambda, Number* lambda)
217  {
218    ASL_pfgh* asl = asl_;
219    DBG_ASSERT(asl_);
220    DBG_ASSERT(n == n_var);
221    DBG_ASSERT(m == n_con);
222
223    if (init_x) {
224      for (Index i=0; i<n; i++) {
225        if (havex0[i]) {
226          x[i] = X0[i];
227        }
228        else {
229          x[i] = 0.0;
230        }
231      }
232    }
233
234    if (init_z) {
235      for (Index i=0; i<n; i++) {
236        z_L[i] = z_U[i] = 1.0;
237      }
238    }
239
240    if (init_lambda) {
241      for (Index i=0; i<m; i++) {
242        if (havepi0[i]) {
243          lambda[i] = pi0[i];
244        }
245        else {
246          lambda[i] = 0.0;
247        }
248      }
249    }
250
251    return true;
252  }
253
254  bool AmplTNLP::eval_f(Index n, const Number* x, bool new_x, Number& obj_value)
255  {
256    apply_new_x(new_x, n, x);
257
258    return internal_objval(obj_value);
259  }
260
261  bool AmplTNLP::eval_grad_f(Index n, const Number* x, bool new_x, Number* grad_f)
262  {
263    ASL_pfgh* asl = asl_;
264    DBG_ASSERT(asl_);
265
266    apply_new_x(new_x, n, x);
267
268    fint nerror = 0;
269    objgrd(0, non_const_x_, grad_f, &nerror);
270    if (nerror == 0) {
271      return true;
272    }
273    return false;
274  }
275
276  bool AmplTNLP::eval_g(Index n, const Number* x, bool new_x, Index m, Number* g)
277  {
278    ASL_pfgh* asl = asl_;
279    DBG_ASSERT(asl_);
280    DBG_ASSERT(n == n_var);
281    DBG_ASSERT(m == n_con);
282
283    apply_new_x(new_x, n, x);
284
285    fint nerror = 0;
286    return internal_conval(m, g);
287  }
288
289  bool AmplTNLP::eval_jac_g(Index n, const Number* x, bool new_x,
290                            Index nele_jac, Index* iRow, Index *jCol,
291                            Number* values)
292  {
293    ASL_pfgh* asl = asl_;
294    DBG_ASSERT(asl_);
295    DBG_ASSERT(n == n_var);
296
297    if (iRow && jCol && !values) {
298      // setup the structure
299      Index current_nz = 0;
300      for (Index i=0; i<n_con; i++) {
301        for (cgrad* cg=Cgrad[i]; cg; cg = cg->next) {
302          iRow[cg->goff] = i + 1;
303          jCol[cg->goff] = cg->varno + 1;
304          //                            iRow[current_nz] = i + 1;
305          //                            jCol[current_nz] = cg->varno+1;
306          current_nz++;
307        }
308      }
309      DBG_ASSERT(current_nz == nele_jac);
310      return true;
311    }
312    else if (!iRow && !jCol && values) {
313      apply_new_x(new_x, n, x);
314
315      fint nerror = 0;
316      jacval(non_const_x_, values, &nerror);
317
318      if (nerror == 0) {
319        return true;
320      }
321    }
322    else {
323      DBG_ASSERT(false && "Invalid combination of iRow, jCol, and values pointers");
324    }
325
326    return false;
327  }
328
329  bool AmplTNLP::eval_h(Index n, const Number* x, bool new_x,
330                        Number obj_factor, Index m, const Number* lambda,
331                        bool new_lambda, Index nele_hess, Index* iRow,
332                        Index* jCol, Number* values)
333  {
334    ASL_pfgh* asl = asl_;
335    DBG_ASSERT(asl_);
336    DBG_ASSERT(n == n_var);
337    DBG_ASSERT(m == n_con);
338
339    if (iRow && jCol && !values) {
340      // setup the structure
341      int k=0;
342      for (int i=0; i<n; i++) {
343        for (int j=sputinfo->hcolstarts[i]; j<sputinfo->hcolstarts[i+1]; j++) {
344          iRow[k] = i + 1;
345          jCol[k] = sputinfo->hrownos[j]+1;
346          k++;
347        }
348      }
349      DBG_ASSERT(k==nele_hess);
350      return true;
351    }
352    else if (!iRow & !jCol && values) {
353      apply_new_x(new_x, n, x);
354      if (!objval_called_with_current_x_) {
355        Number dummy;
356        internal_objval(dummy);
357      }
358      if (!conval_called_with_current_x_) {
359        internal_conval(m);
360      }
361      // copy lambda to non_const_lambda - note, we do not store a copy like
362      // we do with x since lambda is only used here and not in other calls
363      Number* non_const_lambda = new Number[m];
364      for (Index i=0; i<m; i++) {
365        non_const_lambda[i] = lambda[i];
366      }
367
368      fint nerror = 0;
369      real ow=obj_sign_*obj_factor;
370      sphes(values, -1, &ow, non_const_lambda);
371
372      delete [] non_const_lambda;
373      return true;
374    }
375    else {
376      DBG_ASSERT(false && "Invalid combination of iRow, jCol, and values pointers");
377    }
378
379    return false;
380  }
381
382  bool AmplTNLP::internal_objval(Number& obj_val)
383  {
384    ASL_pfgh* asl = asl_;
385    DBG_ASSERT(asl_);
386    objval_called_with_current_x_ = false; // in case the call below fails
387
388    fint nerror = 0;
389    Number retval = objval(0, non_const_x_, &nerror);
390    if (nerror == 0) {
391      obj_val = obj_sign_*retval;
392      objval_called_with_current_x_ = true;
393      return true;
394    }
395
396    //DBG_ASSERT(false && "Error evaluating AMPL objective.\n");
397    return false;
398  }
399
400  bool AmplTNLP::internal_conval(Index m, Number* g)
401  {
402    DBG_START_METH("AmplTNLP::internal_conval", 0);
403    ASL_pfgh* asl = asl_;
404    DBG_ASSERT(asl_);
405    DBG_ASSERT(m == n_con);
406    conval_called_with_current_x_ = false; // in case the call below fails
407
408    bool allocated = false;
409    if (!g) {
410      g = new double[m];
411      allocated = true;
412    }
413
414    fint nerror = 0;
415    conval(non_const_x_, g, &nerror);
416
417    if (allocated) {
418      delete [] g;
419      g = NULL;
420    }
421
422    if (nerror == 0) {
423      conval_called_with_current_x_ = true;
424      return true;
425    }
426    return false;
427  }
428
429
430  void AmplTNLP::apply_new_x(bool new_x, Index n, const Number* x)
431  {
432    ASL_pfgh* asl = asl_;
433    DBG_ASSERT(asl_);
434
435    if (new_x) {
436      // update the flags so these methods are called
437      // before evaluating the hessian
438      conval_called_with_current_x_ = false;
439      objval_called_with_current_x_ = false;
440
441      //copy the data to the non_const_x_
442      if (!non_const_x_) {
443        non_const_x_ = new Number[n];
444      }
445
446      for (Index i=0; i<n; i++) {
447        non_const_x_[i] = x[i];
448      }
449
450      // tell ampl that we have a new x
451      xknown(non_const_x_);
452    }
453  }
454
455  void AmplTNLP::write_solution_file(const std::string& message,
456                                     Number* x, Number* y) const
457  {
458    ASL_pfgh* asl = asl_;
459    DBG_ASSERT(asl);
460
461    // We need to copy the message into a non-const char array to make
462    // it work with the AMPL C function.
463    char* cmessage = new char[message.length()+1];
464    strcpy(cmessage, message.c_str());
465
466    write_sol(cmessage, x, y, NULL);
467
468    delete [] cmessage;
469  }
470
471  AmplSuffixHandler::AmplSuffixHandler()
472      :
473      asl_(NULL),
474      suftab_ (NULL)
475  {}
476
477  AmplSuffixHandler::~AmplSuffixHandler()
478  {
479    if (suftab_) {
480      Index n = suffix_ids_.size();
481      for (Index i=0; i<n; i++) {
482        delete [] suftab_[i].name;
483        suftab_[i].name = NULL;
484      }
485    }
486    delete [] suftab_;
487    suftab_ = NULL;
488  }
489
490  void AmplSuffixHandler::PrepareAmplForSuffixes(ASL_pfgh* asl)
491  {
492    DBG_ASSERT(asl);
493    asl_ = asl;
494
495    Index n = suffix_ids_.size();
496    suftab_ = new SufDecl[n];
497    for (Index i=0; i<n; i++) {
498      Index id_len = strlen(suffix_ids_[i].c_str());
499      suftab_[i].name = new char[id_len + 1];
500      strcpy(suftab_[i].name, suffix_ids_[i].c_str());
501
502      suftab_[i].table = 0;
503
504      if (suffix_sources_[i] == Variable_Source) {
505        suftab_[i].kind = ASL_Sufkind_var;
506      }
507      else if (suffix_sources_[i]  == Constraint_Source) {
508        suftab_[i].kind = ASL_Sufkind_con;
509      }
510      else {
511        DBG_ASSERT(false && "Unknown suffix source in PrepareAmplForSuffixes");
512      }
513
514      suftab_[i].nextra = 0;
515    }
516
517    suf_declare(suftab_, n);
518  }
519
520  const Index* AmplSuffixHandler::GetIntegerSuffixValues(std::string suffix_string, Suffix_Source source) const
521  {
522    ASL_pfgh* asl = asl_;
523    DBG_ASSERT(asl);
524
525    int kind = (source == Variable_Source) ? ASL_Sufkind_var : ASL_Sufkind_con;
526    SufDesc* dp = suf_get(suffix_string.c_str(), kind);
527    return dp->u.i;
528  }
529
530  const Number* AmplSuffixHandler::GetNumberSuffixValues(std::string suffix_string, Suffix_Source source) const
531  {
532    ASL_pfgh* asl = asl_;
533    DBG_ASSERT(asl);
534
535    int kind = (source == Variable_Source) ? ASL_Sufkind_var : ASL_Sufkind_con;
536    SufDesc* dp = suf_get(suffix_string.c_str(), kind);
537    return dp->u.r;
538  }
539
540
541} // namespace Ipopt
542
543
544
Note: See TracBrowser for help on using the repository browser.