/home/coin/SVN-release/OS-2.4.0/Couenne/src/readnl/CouenneAmplInterface.cpp

Go to the documentation of this file.
00001 // $Id: CouenneAmplInterface.cpp 694 2011-06-18 20:13:17Z stefan $
00002 //
00003 // (C) Copyright XXX 2009
00004 // All Rights Reserved.
00005 // This code is published under the Eclipse Public License (EPL).
00006 //
00007 // Authors :
00008 // Pietro Belotti, Lehigh University
00009 // Stefan Vigerske, Humboldt University
00010 //
00011 // Date : 07/18/2009
00012 
00013 #include "CouenneConfig.h"
00014 #include "CouenneAmplInterface.hpp"
00015 #include "CoinPragma.hpp"
00016 
00017 #include <cstdlib>
00018 
00019 #if   defined HAVE_CSTDINT
00020 #include <cstdint>
00021 #elif defined HAVE_STDINT_H
00022 #include <stdint.h>
00023 #endif
00024 
00025 #include <string>
00026 
00027 #include "BonAmplTMINLP.hpp"
00028 #include "BonCbc.hpp"
00029 
00030 #include "CouenneProblem.hpp"
00031 #include "CouenneTypes.hpp"
00032 
00033 #include "CouenneExprClone.hpp"
00034 #include "CouenneExprGroup.hpp"
00035 #include "CouenneExprAbs.hpp"
00036 #include "CouenneExprSum.hpp"
00037 #include "CouenneExprSub.hpp"
00038 #include "CouenneExprMul.hpp"
00039 #include "CouenneExprDiv.hpp"
00040 #include "CouenneExprInv.hpp"
00041 #include "CouenneExprSin.hpp"
00042 #include "CouenneExprPow.hpp"
00043 #include "CouenneExprLog.hpp"
00044 #include "CouenneExprOpp.hpp"
00045 #include "CouenneExprCos.hpp"
00046 #include "CouenneExprExp.hpp"
00047 
00048 #include "asl.h"
00049 #include "nlp.h"
00050 #include "getstub.h"
00051 #include "opcode.hd"
00052 
00053 // get ASL op. code relative to function pointer passed as parameter 
00054 int getOperator (efunc *);
00055 
00056 #define OBJ_DE    ((const ASL_fg *) asl) -> I.obj_de_
00057 #define VAR_E     ((const ASL_fg *) asl) -> I.var_e_
00058 #define CON_DE    ((const ASL_fg *) asl) -> I.con_de_
00059 #define OBJ_sense ((const ASL_fg *) asl) -> i.objtype_
00060 
00061 #include "r_opn.hd" /* for N_OPS */
00062 
00063 static fint timing = 0;
00064 
00065 static
00066 keyword keywds[] = { // must be alphabetical
00067    KW(const_cast<char*>("timing"), L_val, &timing, const_cast<char*>("display timings for the run")),
00068 };
00069 
00070 static
00071 Option_Info Oinfo = { const_cast<char*>("testampl"), const_cast<char*>("ANALYSIS TEST"),
00072                       const_cast<char*>("concert_options"), keywds, nkeywds, 0, const_cast<char*>("ANALYSIS TEST") };
00073 
00074 
00075 // (C++) code starts here ///////////////////////////////////////////////////////////////////////////
00076 
00077 using namespace Couenne;
00078 
00079 void CouenneAmplInterface::registerOptions(Ipopt::SmartPtr<Bonmin::RegisteredOptions> roptions) {
00080         roptions->AddStringOption1("nlfile", "name of an ampl .nl file to get the problem from", "", "*", "name of .nl file");
00081 }
00082 
00083 CouenneAmplInterface::~CouenneAmplInterface() {
00084         delete problem;
00085         
00086         if (asl) {
00087           delete[] X0;
00088           delete[] havex0;
00089           delete[] pi0;
00090           delete[] havepi0;             
00091                 ASL_free(&asl);
00092         }
00093 }
00094 
00095 // create an AMPL problem by using ASL interface to the .nl file
00096 CouenneProblem* CouenneAmplInterface::getCouenneProblem() {
00097         if (problem)
00098                 return problem;
00099 
00100   if (!readASLfg())
00101         return NULL;
00102 
00103   problem = new CouenneProblem;
00104 
00105   if (!readnl()) {
00106         delete problem;
00107         problem = NULL;
00108         return NULL;
00109   }
00110   
00111   return problem;
00112 }
00113 
00114 Ipopt::SmartPtr<Bonmin::TMINLP> CouenneAmplInterface::getTMINLP() {
00115         if (IsValid(tminlp))
00116                 return tminlp;
00117         
00118         if (IsNull(roptions)) {
00119                 jnlst->Printf(Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "Error: Need registered options to create AmplTMINLP object!\n");
00120                 return NULL;
00121         }
00122         
00123         std::string nlfile;
00124         options->GetStringValue("nlfile", nlfile, "");
00125         char** argv = new char*[3];
00126         argv[0] = const_cast<char*>("dummy");
00127         argv[1] = strdup(nlfile.c_str());
00128         argv[2] = NULL;
00129         tminlp = new Bonmin::AmplTMINLP(GetRawPtr(jnlst), roptions, options, argv);
00130         
00131         free(argv[1]);
00132         delete[] argv;
00133         
00134         return tminlp;
00135 }
00136 
00137 bool CouenneAmplInterface::writeSolution(Bonmin::Bab& bab) {
00138         const char* message;
00139 
00140         //TODO setup a nicer message
00141         if (bab.bestSolution()) {
00142                 message = "Couenne found a solution.\n";
00143         } else {
00144                 message = "Couenne could not found a solution.\n";
00145         }
00146 
00147         write_sol(const_cast<char*>(message), const_cast<double*>(bab.bestSolution()), NULL, NULL);
00148         
00149         return true;
00150 }
00151 
00152 bool CouenneAmplInterface::readASLfg() {
00153         assert(asl == NULL);
00154         
00155   std::string nlfile;
00156   options->GetStringValue("nlfile", nlfile, "");
00157   
00158   if (nlfile == "")
00159         return false;
00160 
00161   char** argv = new char*[3];
00162         argv[0] = const_cast<char*>("dummy");
00163         argv[1] = strdup(nlfile.c_str());
00164         argv[2] = NULL;
00165         
00166   // Create the ASL structure
00167   asl = (ASL*) ASL_alloc (ASL_read_fg);
00168 
00169   char* stub = getstub (&argv, &Oinfo);
00170 
00171   // Although very intuitive, we shall explain why the second argument
00172   // is passed with a minus sign: it is to tell the ASL to retrieve
00173   // the nonlinear information too.
00174   FILE* nl = jac0dim (stub, - (fint) strlen (stub));
00175 
00176   // Set options in the asl structure
00177   want_xpi0 = 1 | 2;  // allocate initial values for primal and dual if available
00178   obj_no = 0;         // always want to work with the first (and only?) objective
00179 
00180   // allocate space for initial values
00181   X0      = new real [n_var];
00182   havex0  = new char [n_var];
00183   pi0     = new real [n_con];
00184   havepi0 = new char [n_con];
00185 
00186   // read the rest of the nl file
00187   fg_read (nl, ASL_return_read_err | ASL_findgroups);
00188   
00189   //FIXME freeing argv and argv[1] gives segfault !!!
00190 //  free(argv[1]);
00191 //  delete[] argv;
00192 
00193   return true;
00194 }
00195 
00196 // check if an expression is a null pointer or equals zero
00197 //static inline bool is_expr_zero (expr* e) {
00198 //      return ((e==NULL) || ((((Intcast (e->op)) == OPNUM) && 
00199 //                        (fabs (((expr_n *)e) -> v)  < COUENNE_EPS) 
00200 //                        //  && (fabs (e -> dL) < COUENNE_EPS)
00201 //                        // *** CHECK THIS! dL is the derivative
00202 //                        )));
00203 //} 
00204 
00205 
00206 // Reads a MINLP from an AMPL .nl file through the ASL methods
00207 bool CouenneAmplInterface::readnl() {
00208 
00209   std::string nlfile;
00210   options->GetStringValue("nlfile", nlfile, "");
00211   problem -> setProblemName (nlfile);
00212 
00213   // number of defined variables (aka common expressions)
00214   problem -> setNDefVars(como + comc + comb + como1 + comc1);
00215 
00216   // see "hooking your solver to AMPL", by David M. Gay, tables 3, 4, and 5
00217 
00218   // nonlinear in both objectives and constraints
00219   if (nlvb >= 0) {
00220     for (int i = 0; i < nlvb - nlvbi; i++) problem -> addVariable (false, problem -> domain ());
00221     for (int i = 0; i < nlvbi;        i++) problem -> addVariable (true,  problem -> domain ());
00222   }
00223 
00224   // nonlinear in either objectives or constraints
00225   if (nlvo > nlvc) {
00226     for (int i = 0; i < nlvc - (nlvb + nlvci); i++) problem -> addVariable (false, problem -> domain ());
00227     for (int i = 0; i < nlvci;                 i++) problem -> addVariable (true,  problem -> domain ());
00228     for (int i = 0; i < nlvo - (nlvc + nlvoi); i++) problem -> addVariable (false, problem -> domain ());
00229     for (int i = 0; i < nlvoi;                 i++) problem -> addVariable (true,  problem -> domain ());
00230   } else {
00231     for (int i = 0; i < nlvo - (nlvb + nlvoi); i++) problem -> addVariable (false, problem -> domain ());
00232     for (int i = 0; i < nlvoi;                 i++) problem -> addVariable (true,  problem -> domain ());
00233     for (int i = 0; i < nlvc - (nlvo + nlvci); i++) problem -> addVariable (false, problem -> domain ());
00234     for (int i = 0; i < nlvci;                 i++) problem -> addVariable (true,  problem -> domain ());
00235   }
00236 
00237   for (int i = 0; i < nwv; i++)                                  problem -> addVariable(false, problem -> domain ());//arc
00238   for (int i = n_var - (CoinMax (nlvc,nlvo) +niv+nbv+nwv); i--;) problem -> addVariable(false, problem -> domain ());//other
00239   for (int i = 0; i < nbv; i++)                                  problem -> addVariable(true,  problem -> domain ());//binary
00240   for (int i = 0; i < niv; i++)                                  problem -> addVariable(true,  problem -> domain ());//int.
00241 
00242   // add space for common expressions
00243   for (int i = problem->nDefVars(); i--;)  problem -> addVariable(false, problem -> domain ());
00244 
00245   // common expressions (or defined variables) ///////////////////////////////////////
00246 
00247 #ifdef DEBUG
00248   printf ("tot var = %d\n", variables_ . size ());
00249   printf ("c_vars_ = %d\n", ((const ASL_fg *) asl) -> i.c_vars_ );
00250   printf ("comb_ = %d\n",   ((const ASL_fg *) asl) -> i.comb_  );
00251   printf ("combc_ = %d\n",  ((const ASL_fg *) asl) -> i.combc_ );
00252   printf ("comc1_ = %d\n",  ((const ASL_fg *) asl) -> i.comc1_ );
00253   printf ("comc_ = %d\n",   ((const ASL_fg *) asl) -> i.comc_  );
00254   printf ("como1_ = %d\n",  ((const ASL_fg *) asl) -> i.como1_ );
00255   printf ("como_ = %d\n",   ((const ASL_fg *) asl) -> i.como_  );
00256 #endif
00257 
00258   // Each has a linear and a nonlinear part (thanks to Dominique
00259   // Orban: http://www.gerad.ca/~orban/drampl/def-vars.html)
00260 
00261   try {
00262   for (int i = 0; i < como + comc + comb; i++) {
00263 
00264     struct cexp *common = ((const ASL_fg *) asl) -> I.cexps_ + i;
00265     expression *nle = nl2e (common -> e);
00266 
00267 #ifdef DEBUG
00268     printf ("cexp  %d [%d]: ", i, problem -> nVars()); nle -> print ();  printf (" ||| ");
00269 #endif
00270 
00271     int nlin = common -> nlin;  // Number of linear terms
00272 
00273     if (nlin > 0) {
00274 
00275       int       *indexL = new int       [nlin+1];
00276       CouNumber *coeff  = new CouNumber [nlin];
00277 
00278       linpart *L = common -> L;
00279 
00280       for (int j = 0; j < nlin; j++) {
00281         //vp = (expr_v *)((char *)L->v.rp - ((char *)&ev.v - (char *)&ev));
00282         //Printf( " %-g x[%-d]", L->fac, (int)(vp - VAR_E) );   
00283         coeff [j] = L [j]. fac;
00284         indexL [j] = ((uintptr_t) (L [j].v.rp) - (uintptr_t) VAR_E) / sizeof (expr_v);
00285 #ifdef DEBUG
00286         Printf( " %+g x_%-3d", L [j]. fac, 
00287                 (expr_v *) (L [j].v.rp) - VAR_E //((const ASL_fg *) asl) -> I.cexps_
00288                 //L [j]. v.i
00289                 );
00290 #endif
00291       }
00292 
00293       indexL [nlin] = -1;
00294 
00295       expression **al = new expression * [1];
00296       *al = nle;
00297 
00298       std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
00299       problem -> indcoe2vector (indexL, coeff, lcoeff);
00300 
00301       expression *eg = exprGroup::genExprGroup (0, lcoeff, al, 1);
00302       problem -> commonExprs (). push_back (eg);
00303     } 
00304     else problem -> commonExprs () . push_back (nle);
00305 #ifdef DEBUG
00306     printf ("\n");
00307 #endif
00308   }
00309 
00310   for (int i = 0; i < como1 + comc1; i++) {
00311 
00312     struct cexp1 *common = ((const ASL_fg *) asl) -> I.cexps1_ + i;
00313     expression *nle = nl2e (common -> e);
00314 
00315 #ifdef DEBUG
00316     printf ("cexp1 %d [%d]: ", i, variables_ . size ()); nle -> print ();  printf (" ||| ");
00317 #endif
00318 
00319     int nlin = common -> nlin;  // Number of linear terms
00320 
00321     if (nlin > 0) {
00322 
00323       int       *indexL = new int       [nlin+1];
00324       CouNumber *coeff  = new CouNumber [nlin];
00325 
00326       linpart *L = common -> L;
00327 
00328       for (int j = 0; j < nlin; j++) {
00329         //vp = (expr_v *)((char *)L->v.rp - ((char *)&ev.v - (char *)&ev));
00330         coeff  [j] = L [j]. fac;
00331         indexL [j] = ((uintptr_t) (L [j].v.rp) - (uintptr_t) VAR_E) / sizeof (expr_v);
00332 #ifdef DEBUG
00333         Printf( " %+g x_%-3d", L [j]. fac, 
00334                 (expr_v *) (L [j].v.rp) - VAR_E //((const ASL_fg *) asl) -> I.cexps_
00335                 //L [j]. v.i
00336                 );
00337 #endif
00338       }
00339 
00340       indexL [nlin] = -1;
00341 
00342       expression **al = new expression * [1];
00343       *al = nle;
00344 
00345       std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
00346       problem -> indcoe2vector (indexL, coeff, lcoeff);
00347 
00348       expression *eg = exprGroup::genExprGroup (0, lcoeff, al, 1);
00349       problem -> commonExprs () . push_back (eg);
00350     } 
00351     else problem -> commonExprs () . push_back (nle);
00352 #ifdef DEBUG
00353     printf ("\n");
00354 #endif
00355     //    addAuxiliary (nl2e (((const ASL_fg *) asl) -> I.cexps1_ [i] . e));
00356   }
00357 
00358   // objective functions /////////////////////////////////////////////////////////////
00359 
00360   for (int i = 0; i < n_obj; i++) {
00361 
00363     int nterms = 0;
00364 
00365     // count nonzero terms in linear part
00366  
00367     for (ograd *objgrad = Ograd [i];
00368          objgrad;
00369          objgrad = objgrad -> next)
00370       if (fabs (objgrad -> coef) > COUENNE_EPS)
00371         nterms++;
00372 
00373     expression 
00374       *body,
00375       *nl = nl2e (OBJ_DE [i] . e);
00376 
00377     if (nterms) { // have linear terms
00378 
00379       int       *indexL = new int       [nterms+1];
00380       CouNumber *coeff  = new CouNumber [nterms];
00381 
00382       for (ograd *objgrad = Ograd [i]; objgrad; objgrad = objgrad -> next)
00383         if (fabs (objgrad -> coef) > COUENNE_EPS) {
00384 
00385           *indexL++ = objgrad -> varno;
00386           *coeff++  = objgrad -> coef;
00387         }
00388 
00389       *indexL = -1;
00390 
00391       indexL -= nterms;
00392       coeff  -= nterms;
00393 
00394       std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
00395       problem -> indcoe2vector (indexL, coeff, lcoeff);
00396 
00397       if (nl -> code () == COU_EXPRSUM) {
00398         body = exprGroup::genExprGroup (0., lcoeff, nl -> ArgList (), nl -> nArgs ());
00399         // delete node without deleting children (they are now in body)
00400         nl -> ArgList (NULL);
00401         delete nl;
00402       }
00403       else {
00404 
00405         expression **nll = new expression * [1];
00406 
00407         *nll = nl;
00408 
00409         // apparently, objconst (i) is included in the obj expression
00410         body = exprGroup::genExprGroup (0., lcoeff, nll, 1);
00411         //body = new exprGroup (objconst (i), indexL, coeff, nll, 1);
00412       }
00413 
00414       delete [] indexL;
00415       delete [] coeff;
00416 
00417     } else
00418       // apparently, objconst (i) is included in the obj expression
00419       body = nl;
00420       //if (fabs (objconst (i) > COUENNE_EPS))
00421       //body = new exprSum (nl, new exprConst (objconst (i)));
00422       //else 
00423 
00425 
00426     expression *subst = body -> simplify ();
00427 
00428     if (subst) {
00429       delete body; // VALGRIND
00430       body = subst;
00431     }
00432 
00433     // ThirdParty/ASL/solvers/asl.h, line 336: 0 is minimization, 1 is maximization
00434     problem -> addObjective (body, (OBJ_sense [i] == 0) ? "min" : "max");
00435   }
00436 
00437   // constraints ///////////////////////////////////////////////////////////////////
00438 
00439   int *nterms = new int [n_con];
00440 
00441   // allocate space for argument list of all constraints' summations
00442   // of linear and nonlinear terms
00443 
00444   // init array with # terms of each constraint
00445   for (int i = n_con; i--;) 
00446     *nterms++ = 0;
00447   nterms -= n_con;
00448 
00449   cgrad *congrad;
00450 
00451   // count all linear terms
00452   if (A_colstarts && A_vals)         // Constraints' linear info is stored in A_vals
00453     for (register int j = A_colstarts [n_var]; j--;) {
00454 
00455       real coeff = A_vals [j];
00456 
00457       if (fabs (coeff) > COUENNE_EPS)
00458         nterms [A_rownos [j]] ++;
00459     }
00460   else {                             // Constraints' linear info is stored in Cgrad
00461     for (register int i = 0; i < n_con; i++)
00462       for (congrad = Cgrad [i]; 
00463            congrad; 
00464            congrad = congrad -> next) 
00465         if (fabs (congrad -> coef) > COUENNE_EPS) 
00466           nterms [i] ++;
00467   }
00468 
00469 
00470   // vectors of the linear part
00471   CouNumber **coeff  = new CouNumber * [n_con];
00472   int       **indexL = new int       * [n_con];
00473 
00474   for (register int i = n_con; i--;) 
00475     *indexL++ = NULL;
00476 
00477   indexL -= n_con;
00478 
00479 
00480   // set linear terms
00481 
00482   if (A_colstarts && A_vals)         // Constraints' linear info is stored in A_vals
00483     for (int j = 0; j < n_var; j++)
00484       for (register int i = A_colstarts [j], k = A_colstarts [j+1] - i; k--; i++) {
00485 
00486         int rowno = A_rownos [i],
00487             nt    = nterms [rowno] --;
00488 
00489         CouNumber **cline = coeff  + rowno;
00490         int       **iline = indexL + rowno;
00491 
00492         if (*iline==NULL) {
00493           *cline = new CouNumber [nt];
00494           *iline = new int       [nt+1];
00495           (*iline) [nt] = -1;
00496         }
00497 
00498         (*cline) [--nt] = A_vals [i];
00499         (*iline)   [nt] = j;
00500 
00501       }
00502   else {                             // Constraints' linear info is stored in Cgrad
00503     for (int i=0; i < n_con; i++) {
00504 
00505       int nt = nterms [i];
00506 
00507       CouNumber **cline = coeff + i;
00508       int       **iline = indexL + i;
00509 
00510       *cline = new CouNumber [nt];
00511       *iline = new int       [nt+1];
00512       (*iline) [nt] = -1;
00513 
00514       for (congrad = Cgrad [i]; congrad; congrad = congrad -> next) 
00515         if (fabs (congrad -> coef) > COUENNE_EPS) {
00516           (*cline) [--nt] = congrad -> coef;
00517           (*iline)   [nt] = congrad -> varno;
00518         }
00519     }
00520   }
00521 
00522   // set constraints' bound and sign and store nonlinear part ///////////////////////////////
00523 
00524   for (int i = 0; i < n_con; i++) {
00525 
00526     enum con_sign sign;
00527     double lb, ub;
00528 
00529     if (Urhsx) {
00530       lb = LUrhs [i];
00531       ub = Urhsx [i];
00532     } else {
00533       int j = 2*i;
00534       lb = LUrhs [j];
00535       ub = LUrhs [j+1];
00536     }
00537 
00538     // set constraint sign
00539     if (lb > negInfinity)
00540       if (ub < Infinity) sign = COUENNE_RNG;
00541       else               sign = COUENNE_GE;
00542     else                 sign = COUENNE_LE;
00543 
00544     // this is an equality constraint  
00545     if (fabs (lb - ub) < COUENNE_EPS)
00546       sign = COUENNE_EQ;
00547 
00548     expression *body;
00549 
00550     expression **nll = new expression * [1];
00551     *nll = nl2e (CON_DE [i] . e);
00552 
00553     if (indexL [i] && (*(indexL [i]) >= 0)) {
00554 
00555       int code = (*nll) -> code ();
00556 
00557       std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
00558       problem -> indcoe2vector (indexL [i], coeff [i], lcoeff);
00559 
00560       /*std::vector <std::pair <exprVar *, CouNumber> > lcoeff;
00561       for (int i=0, *ind = indexL; *ind >= 0; *ind++, i++)
00562       lcoeff.push_back (std::pair <exprVar *, CouNumber> (Var (*ind), coeff [i]));*/
00563 
00564       if ((code == COU_EXPRSUM) || 
00565           (code == COU_EXPRGROUP)) {
00566 
00567         body    = exprGroup::genExprGroup (0., lcoeff, (*nll) -> ArgList (), (*nll) -> nArgs ());
00568         // delete node without deleting children (they are now in body)
00569         (*nll) -> ArgList (NULL);
00570         delete *nll;
00571         delete [] nll;
00572       }
00573       else body = exprGroup::genExprGroup (0., lcoeff, nll, 1);
00574     }
00575     else {
00576       body = *nll;
00577       delete [] nll;
00578     }
00579 
00580     expression *subst = body -> simplify ();
00581     if (subst) {
00582       delete body; // VALGRIND
00583       body = subst;
00584     }
00585 
00586     // add them (and set lower-upper bound)
00587     switch (sign) {
00588 
00589     case COUENNE_EQ:  problem -> addEQConstraint  (body, new exprConst (ub)); break;
00590     case COUENNE_LE:  problem -> addLEConstraint  (body, new exprConst (ub)); break;
00591     case COUENNE_GE:  problem -> addGEConstraint  (body, new exprConst (lb)); break;
00592     case COUENNE_RNG: problem -> addRNGConstraint (body, new exprConst (lb), 
00593                                                    new exprConst (ub)); break;
00594     default: jnlst->Printf(Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "Error: could not recognize constraint\n"); return false;
00595     }
00596 
00597     delete [] indexL [i];
00598     delete [] coeff  [i];
00599   }
00600 
00601   delete [] indexL;
00602   delete [] coeff;
00603   delete [] nterms;
00604   
00605   } catch (...) {
00606         return false;
00607   }
00608 
00609   // create room for problem's variables and bounds
00610   CouNumber 
00611     *x  = (CouNumber *) malloc ((n_var + problem -> nDefVars() ) * sizeof (CouNumber)),
00612     *lb = (CouNumber *) malloc ((n_var + problem -> nDefVars() ) * sizeof (CouNumber)),
00613     *ub = (CouNumber *) malloc ((n_var + problem -> nDefVars() ) * sizeof (CouNumber));
00614 
00615   for (int i = n_var + problem -> nDefVars(); i--;) {
00616     x  [i] =  0.;
00617     lb [i] = -COUENNE_INFINITY;
00618     ub [i] =  COUENNE_INFINITY;
00619   }
00620 
00621   problem -> domain () -> push (n_var + problem -> nDefVars(), x, lb, ub);
00622   free (x); free (lb); free (ub);
00623 
00624   // suggested:
00625   // problem -> domain () -> push (n_var + problem -> nDefVars(), x, lb, ub, false);
00626   // //free (x); free (lb); free (ub);
00627   // saves three allocations (default last parameter is true, which copies x,l,b)
00628 
00629   // lower and upper bounds ///////////////////////////////////////////////////////////////
00630 
00631   if (LUv) {
00632 
00633     real *Uvx_copy = Uvx;
00634 
00635     if (!Uvx_copy)
00636       for (register int i=0; i<n_var; i++) {
00637 
00638         register int j = 2*i;
00639 
00640         problem -> Lb (i) = LUv[j]   <= -COUENNE_INFINITY ? -COUENNE_INFINITY : LUv[j]  ;
00641         problem -> Ub (i) = LUv[j+1] >=  COUENNE_INFINITY ?  COUENNE_INFINITY : LUv[j+1];
00642       }
00643     else
00644       for (register int i=n_var; i--;) {
00645         problem -> Lb (i) = LUv [i]      <= -COUENNE_INFINITY ? -COUENNE_INFINITY : LUv[i];
00646         problem -> Ub (i) = Uvx_copy [i] >=  COUENNE_INFINITY ?  COUENNE_INFINITY : Uvx_copy[i];
00647       }
00648 
00649   } else
00650     for (register int i=n_var; i--;) {
00651         problem -> Lb (i) = - COUENNE_INFINITY;
00652         problem -> Ub (i) =   COUENNE_INFINITY;
00653     }
00654 
00655   // initial x ////////////////////////////////////////////////////////////////////
00656 
00657   for (register int i=n_var; i--;) 
00658 
00659     if (X0 && havex0 [i]) problem -> X (i) = X0 [i]; 
00660 
00661     else {
00662 
00663       CouNumber x, l = problem -> Lb (i), u = problem -> Ub (i);
00664 
00665       if      (l < - COUENNE_INFINITY)
00666         if    (u >   COUENNE_INFINITY)  x = 0.;
00667         else                            x = u;
00668       else if (u >   COUENNE_INFINITY)  x = l;
00669       else                              x = 0.5 * (l+u);
00670 
00671       problem -> X (i) = x;
00672     }
00673 
00674   for (register int i=n_var; i < problem -> nDefVars() ; i++) {  //FIXME: shouldn't this loop go until n_var + problem -> nDefVars() ?
00675 
00676         problem -> X  (i) =  0.;
00677         problem -> Lb (i) = -COUENNE_INFINITY;
00678         problem -> Ub (i) =  COUENNE_INFINITY;
00679   }
00680 
00681   return true;
00682 }
00683 
00684 
00685 // warning for non-implemented functions -- return 0 constant expression
00686 //expression *notimpl (const std::string &fname) {
00687 //static void notimpl (const std::string &fname) {
00688 //  std::cerr << "*** Error: " << fname << " not implemented" << std::endl;
00689 //  exit (-1);
00690 //}
00691 
00692 // converts an AMPL expression (sub)tree into an expression* (sub)tree
00693 expression *CouenneAmplInterface::nl2e(expr *e) {
00694 
00695   switch (getOperator (e -> op)) {
00696 
00697   case OPPLUS:  return new exprSum (nl2e (e -> L.e), nl2e (e -> R.e));
00698   case OPMINUS: return new exprSub (nl2e (e -> L.e), nl2e (e -> R.e));
00699   case OPMULT:  return new exprMul (nl2e (e -> L.e), nl2e (e -> R.e));
00700   case OPDIV:   return new exprDiv (nl2e (e -> L.e), nl2e (e -> R.e));
00701     //case OPREM:   notimpl ("remainder");
00702   case OPPOW:   return new exprPow (nl2e (e -> L.e), nl2e (e -> R.e));
00703     //case OPLESS:  notimpl ("less");
00704     //case MINLIST: notimpl ("min");
00705     //case MAXLIST: notimpl ("max");
00706     //case FLOOR:   notimpl ("floor");
00707     //case CEIL:    notimpl ("ceil");
00708   case ABS:     return new exprAbs (nl2e (e -> L.e));
00709   case OPUMINUS:return new exprOpp (nl2e (e -> L.e));
00710     //          return new exprOpp (nl2e (e -> L.e -> L.e));
00711     //case OPIFnl:  { notimpl ("ifnl");
00712 
00713     // see ASL/solvers/rops.c, IfNL
00714     //}
00715 
00716   case OP_tanh: return new exprDiv 
00717       (new exprSub (new exprExp (nl2e (e -> L.e)),
00718                     new exprExp (new exprOpp (nl2e (e->L.e)))),
00719        new exprSum (new exprExp (nl2e (e -> L.e)),
00720                     new exprExp (new exprOpp (nl2e (e->L.e)))));
00721 
00722   case OP_tan: {
00723     expression *arg;
00724     arg = nl2e (e -> L.e);
00725     return new exprDiv (new exprSin (arg), new exprCos (new exprClone (arg)));
00726   }
00727   case OP_sqrt:    return new exprPow (nl2e (e -> L.e), new exprConst (0.5));
00728   case OP_sinh:    return new exprMul (new exprConst (0.5),
00729                                        new exprSub (new exprExp (nl2e (e -> L.e)),
00730                                                     new exprExp (new exprOpp (nl2e (e->L.e)))));
00731   case OP_sin:     return new exprSin (nl2e (e -> L.e));
00732   case OP_log10:   return new exprMul (new exprConst (1.0 / log (10.0)), 
00733                                        new exprLog (nl2e (e -> L.e)));
00734   case OP_log:     return new exprLog (nl2e (e -> L.e));
00735   case OP_exp:     return new exprExp (nl2e (e -> L.e));
00736   case OP_cosh:    return new exprMul (new exprConst (0.5),
00737                                        new exprSum (new exprExp (nl2e (e -> L.e)),
00738                                                     new exprExp (new exprOpp (nl2e (e->L.e)))));
00739 
00740   case OP_cos:   return new exprCos (nl2e (e -> L.e));
00741     //case OP_atanh: notimpl ("atanh");
00742     //case OP_atan2: notimpl ("atan2");
00743     //case OP_atan:  notimpl ("atan");
00744     //case OP_asinh: notimpl ("asinh");
00745     //case OP_asin:  notimpl ("asin");
00746     //case OP_acosh: notimpl ("acosh");
00747     //case OP_acos:  notimpl ("acos");
00748 
00749   case OPSUMLIST: {
00750     int i=0;
00751     expression **al = new expression * [(e->R.ep - e->L.ep)];
00752     for (expr **ep = e->L.ep; ep < e->R.ep; ep++)
00753       al [i++] = nl2e (*ep);
00754     return new exprSum (al, i);
00755   }
00756     //case OPintDIV: notimpl ("intdiv");
00757     //case OPprecision: notimpl ("precision");
00758     //case OPround:  notimpl ("round");
00759     //case OPtrunc:  notimpl ("trunc");
00760 
00761   case OP1POW: return new exprPow (nl2e (e -> L.e), 
00762                                    new exprConst (((expr_n *)e->R.e)->v));
00763   case OP2POW: return new exprPow (nl2e (e -> L.e), 
00764                                    new exprConst (2.));
00765   case OPCPOW: return new exprPow (new exprConst (((expr_n *)e->L.e)->v),
00766                                    nl2e (e -> R.e));
00767     //case OPFUNCALL: notimpl ("function call");
00768   case OPNUM:     return new exprConst (((expr_n *)e)->v);
00769     //case OPPLTERM:  notimpl ("plterm");
00770     //case OPIFSYM:   notimpl ("ifsym");
00771     //case OPHOL:     notimpl ("hol");
00772   case OPVARVAL:  {
00773 
00774     int j = ((expr_v *) e) -> a;
00775 
00776     if (j >= problem -> nOrigVars()) // common expression
00777       // use base pointer otherwise the .a field returns an awkward, out-of-bound index
00778       j = ((expr_v *) e) - ((const ASL_fg *) asl) -> I.var_e_; 
00779 
00780     if (j >= problem -> nOrigVars() + problem -> nDefVars()) {
00781       jnlst -> Printf (Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "Error: unknown variable x_%d\n", j);
00782       throw -1;
00783     }
00784 
00785     return new exprClone (problem -> Variables() [j]);
00786   }
00787 
00788   default:
00789     jnlst -> Printf (Ipopt::J_ERROR, Ipopt::J_INITIALIZATION, "ERROR: unknown operator (address %p), aborting.\n", Intcast (e -> op));
00790     throw -2;
00791   }
00792 
00793   return new exprConst (0.);
00794 }

Generated on Thu Sep 22 03:05:59 2011 by  doxygen 1.4.7