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

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

Generated on Thu Aug 5 03:02:57 2010 by  doxygen 1.4.7