source: trunk/Cbc/src/ClpAmplStuff.cpp @ 789

Last change on this file since 789 was 789, checked in by forrest, 12 years ago

trying to move ampl stuff away

File size: 37.9 KB
Line 
1// Copyright (C) 2007, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "ClpConfig.h"
5#include "CbcConfig.h"
6#ifdef COIN_HAS_ASL
7#include "CoinPragma.hpp"
8#include "CoinHelperFunctions.hpp"
9#include "CoinIndexedVector.hpp"
10#include "ClpFactorization.hpp"
11#include "ClpSimplex.hpp"
12#include "ClpAmplObjective.hpp"
13#include "ClpConstraintAmpl.hpp"
14#include "ClpMessage.hpp"
15#include "CoinUtilsConfig.h"
16#include "CoinHelperFunctions.hpp"
17#include "CoinWarmStartBasis.hpp"
18#include "OsiSolverInterface.hpp"
19#include "CbcSolver.hpp"
20#include "Cbc_ampl.h"
21#include "CoinTime.hpp"
22#include "CglStored.hpp"
23#include "CoinModel.hpp"
24#include "CbcLinked.hpp"
25class CbcAmpl  : public CbcUser {
26 
27public:
28  ///@name usage methods
29  //@{
30  /// Solve (whatever that means)
31  virtual void solve(CbcSolver * model, const char * options);
32  /// Returns true if function knows about option
33  virtual bool canDo(const char * options) ;
34  /** Import - gets full command arguments
35      Returns -1 - no action
36               0 - data read in without error
37               1 - errors
38  */
39  virtual int importData(CbcSolver * model, int & argc, char * argv[]);
40  /// Export 1 OsiClpSolver, 2 CbcModel - add 10 if infeasible from odd situation
41  virtual void exportSolution(CbcSolver * model, int mode,const char * message=NULL) ;
42  /// Export Data (i.e. at very end)
43  virtual void exportData(CbcSolver * model);
44  /// Get useful stuff
45  virtual void fillInformation(CbcSolver * model,
46                               CbcSolverUsefulData & info);
47  //@}
48  ///@name Constructors and destructors etc
49  //@{
50  /// Default Constructor
51  CbcAmpl(); 
52 
53  /** Copy constructor .
54   */ 
55  CbcAmpl(const CbcAmpl & rhs);
56 
57  /// Assignment operator
58  CbcAmpl & operator=(const CbcAmpl& rhs);
59
60  /// Clone
61  virtual CbcUser * clone() const;
62 
63  /// Destructor
64  virtual ~CbcAmpl ();
65  //@}
66private:
67  ///@name Private member data
68  //@{
69  /// AMPL info
70  ampl_info info_;
71  //@}
72};
73// Mechanicsburg stuff
74CbcAmpl::CbcAmpl()
75  : CbcUser()
76{
77  userName_ = "mech";
78  memset(&info_,0,sizeof(info_));
79}
80CbcAmpl::~CbcAmpl()
81{
82}
83// Copy constructor
84CbcAmpl::CbcAmpl ( const CbcAmpl & rhs)
85  : CbcUser(rhs)
86{
87  info_ = rhs.info_;
88}
89// Assignment operator
90CbcAmpl &
91CbcAmpl::operator=(const CbcAmpl & rhs)
92{
93  if (this != &rhs) {
94    CbcUser::operator=(rhs);
95    info_ = rhs.info_;
96  }
97  return *this;
98}
99// Clone
100CbcUser *
101CbcAmpl::clone() const
102{
103  return new CbcAmpl(*this);
104}
105// Solve (whatever that means)
106void 
107CbcAmpl::solve(CbcSolver * controlModel, const char * options)
108{
109  CbcModel * model = controlModel->model();
110  assert (model);
111  //OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model->solver());
112  //ClpSimplex * lpSolver = clpSolver->getModelPtr();
113  if (!strcmp(options,"cbc_load")) {
114  } else if (!strcmp(options,"cbc_quit")) {
115  } else {
116    printf("unknown option for CbcAmpl is %s\n",options);
117    abort();
118  }
119}
120// Returns true if function knows about option
121bool 
122CbcAmpl::canDo(const char * options) 
123{
124  return (!strcmp(options,"cbc_load")||!strcmp(options,"cbc_quit"));
125}
126/* Import - gets full command arguments
127   Returns -1 - no action
128            0 - data read in without error
129            1 - errors
130*/
131int 
132CbcAmpl::importData(CbcSolver * control, int &argc, char * argv[])
133{
134  CbcModel * babModel = control->model();
135  assert (babModel);
136  CoinMessageHandler * generalMessageHandler = babModel->messageHandler();
137  OsiClpSolverInterface * solver = dynamic_cast< OsiClpSolverInterface*> (control->model()->solver());
138  assert (solver);
139  CoinMessages generalMessages = solver->getModelPtr()->messages();
140  char generalPrint[10000];
141  OsiSolverLink * si = NULL;
142  ClpSimplex * lpSolver = solver->getModelPtr();
143  if (argc>2&&!strcmp(argv[2],"-AMPL")) {
144    // see if log in list
145    bool printing=false;
146    for (int i=1;i<argc;i++) {
147      if (!strncmp(argv[i],"log",3)) {
148        const char * equals = strchr(argv[i],'=');
149        if (equals&&atoi(equals+1)>0) {
150          printing=true;
151          info_.logLevel=atoi(equals+1);
152          control->setIntValue(LOGLEVEL,info_.logLevel);
153          // mark so won't be overWritten
154          info_.numberRows=-1234567;
155          break;
156        }
157      }
158    }
159    CoinModel * coinModel=NULL;
160    int returnCode = readAmpl(&info_,argc, argv,(void **) (& coinModel));
161    if (returnCode)
162      return returnCode;
163    control->setReadMode(3); // so will start with parameters
164    // see if log in list (including environment)
165    for (int i=1;i<info_.numberArguments;i++) {
166      if (!strcmp(info_.arguments[i],"log")) {
167        if (i<info_.numberArguments-1&&atoi(info_.arguments[i+1])>0)
168          printing=true;
169        break;
170      }
171    }
172    control->setPrinting(printing);
173    if (printing)
174      printf("%d rows, %d columns and %d elements\n",
175             info_.numberRows,info_.numberColumns,info_.numberElements);
176    if (!coinModel) {
177      solver->loadProblem(info_.numberColumns,info_.numberRows,info_.starts,
178                          info_.rows,info_.elements,
179                          info_.columnLower,info_.columnUpper,info_.objective,
180                          info_.rowLower,info_.rowUpper);
181      if (info_.numberSos) {
182        // SOS
183        solver->setSOSData(info_.numberSos,info_.sosType,info_.sosStart,
184                           info_.sosIndices,info_.sosReference);
185      }
186    } else {
187      // save
188      control->setOriginalCoinModel(coinModel);
189      // load from coin model
190      OsiSolverLink solver1;
191      OsiSolverInterface * solver2 = solver1.clone();
192      babModel->assignSolver(solver2,false);
193      si = dynamic_cast<OsiSolverLink *>(babModel->solver()) ;
194      assert (si != NULL);
195      si->setDefaultMeshSize(0.001);
196      // need some relative granularity
197      si->setDefaultBound(100.0);
198      double dextra3 = control->doubleValue(DEXTRA3);
199      if (dextra3)
200        si->setDefaultMeshSize(dextra3);
201      si->setDefaultBound(100000.0);
202      si->setIntegerPriority(1000);
203      si->setBiLinearPriority(10000);
204      CoinModel * model2 = (CoinModel *) coinModel;
205      int logLevel = control->intValue(LOGLEVEL);
206      si->load(*model2,true,logLevel);
207      // redo
208      solver = dynamic_cast< OsiClpSolverInterface*> (control->model()->solver());
209      lpSolver = solver->getModelPtr();
210      solver->messageHandler()->setLogLevel(0) ;
211      control->setIntValue(TESTOSI,0);
212      if (info_.cut) {
213        printf("Sorry - can't do cuts with LOS as ruins delicate row order\n");
214        abort();
215      }
216    }
217    if (info_.cut) {
218      int numberRows = info_.numberRows;
219      int * whichRow = new int [numberRows];
220      // Row copy
221      const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
222      const double * elementByRow = matrixByRow->getElements();
223      const int * column = matrixByRow->getIndices();
224      const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
225      const int * rowLength = matrixByRow->getVectorLengths();
226     
227      const double * rowLower = solver->getRowLower();
228      const double * rowUpper = solver->getRowUpper();
229      int nDelete=0;
230      CglStored storedAmpl;
231      for (int iRow=0;iRow<numberRows;iRow++) {
232        if (info_.cut[iRow]) {
233          whichRow[nDelete++]=iRow;
234          int start = rowStart[iRow];
235          storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
236                            rowLength[iRow],column+start,elementByRow+start);
237        }
238      }
239      control->addCutGenerator(&storedAmpl);
240      solver->deleteRows(nDelete,whichRow);
241      // and special matrix
242      si->cleanMatrix()->deleteRows(nDelete,whichRow);
243      delete [] whichRow;
244    }
245    // If we had a solution use it
246    if (info_.primalSolution) {
247      solver->setColSolution(info_.primalSolution);
248    }
249    // status
250    if (info_.rowStatus) {
251      unsigned char * statusArray = lpSolver->statusArray();
252      memset(statusArray,0,lpSolver->numberColumns()+lpSolver->numberRows());
253      int i;
254      for (i=0;i<info_.numberColumns;i++)
255        statusArray[i]=(char)info_.columnStatus[i];
256      statusArray+=info_.numberColumns;
257      for (i=0;i<info_.numberRows;i++)
258        statusArray[i]=(char)info_.rowStatus[i];
259      CoinWarmStartBasis * basis = lpSolver->getBasis();
260      solver->setWarmStart(basis);
261      delete basis;
262    }
263    freeArrays1(&info_);
264    // modify objective if necessary
265    solver->setObjSense(info_.direction);
266    solver->setDblParam(OsiObjOffset,info_.offset);
267    if (info_.offset) {
268      sprintf(generalPrint,"Ampl objective offset is %g",
269              info_.offset);
270      generalMessageHandler->message(CLP_GENERAL,generalMessages)
271        << generalPrint
272        <<CoinMessageEol;
273    }
274    // Set integer variables (unless nonlinear when set)
275    if (!info_.nonLinear) {
276      for (int i=info_.numberColumns-info_.numberIntegers;
277           i<info_.numberColumns;i++)
278        solver->setInteger(i);
279    }
280    // change argc etc
281    argc = info_.numberArguments;
282    argv = info_.arguments;
283    return 0;
284  } else {
285    return -1;
286  }
287  abort();
288  return -1;
289}
290// Export 1 OsiClpSolver, 2 CbcModel - add 10 if infeasible from odd situation
291void 
292CbcAmpl::exportSolution(CbcSolver * model, int mode,const char * message)
293{
294  OsiClpSolverInterface * solver = model->originalSolver();
295  if (!solver) {
296    solver = dynamic_cast< OsiClpSolverInterface*> (model->model()->solver());
297    assert (solver);
298  }
299  ClpSimplex * lpSolver = solver->getModelPtr();
300  int numberColumns=lpSolver->numberColumns();
301  int numberRows = lpSolver->numberRows();
302  double totalTime = CoinCpuTime()-model->startTime();
303  if (mode==1) {
304    double value = lpSolver->getObjValue()*lpSolver->getObjSense();
305    char buf[300];
306    int pos=0;
307    int iStat = lpSolver->status();
308    if (iStat==0) {
309      pos += sprintf(buf+pos,"optimal," );
310    } else if (iStat==1) {
311      // infeasible
312      pos += sprintf(buf+pos,"infeasible,");
313    } else if (iStat==2) {
314      // unbounded
315      pos += sprintf(buf+pos,"unbounded,");
316    } else if (iStat==3) {
317      pos += sprintf(buf+pos,"stopped on iterations or time,");
318    } else if (iStat==4) {
319      iStat = 7;
320      pos += sprintf(buf+pos,"stopped on difficulties,");
321    } else if (iStat==5) {
322      iStat = 3;
323      pos += sprintf(buf+pos,"stopped on ctrl-c,");
324    } else {
325      pos += sprintf(buf+pos,"status unknown,");
326      iStat=6;
327    }
328    info_.problemStatus=iStat;
329    info_.objValue = value;
330    pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
331                   value);
332    sprintf(buf+pos,"\n%d iterations",
333            lpSolver->getIterationCount());
334    free(info_.primalSolution);
335    info_.primalSolution = (double *) malloc(numberColumns*sizeof(double));
336    CoinCopyN(lpSolver->primalColumnSolution(),numberColumns,info_.primalSolution);
337    free(info_.dualSolution);
338    info_.dualSolution = (double *) malloc(numberRows*sizeof(double));
339    CoinCopyN(lpSolver->dualRowSolution(),numberRows,info_.dualSolution);
340    CoinWarmStartBasis * basis = lpSolver->getBasis();
341    free(info_.rowStatus);
342    info_.rowStatus = (int *) malloc(numberRows*sizeof(int));
343    free(info_.columnStatus);
344    info_.columnStatus = (int *) malloc(numberColumns*sizeof(int));
345    // Put basis in
346    int i;
347    // free,basic,ub,lb are 0,1,2,3
348    for (i=0;i<numberRows;i++) {
349      CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
350      info_.rowStatus[i]=status;
351    }
352    for (i=0;i<numberColumns;i++) {
353      CoinWarmStartBasis::Status status = basis->getStructStatus(i);
354      info_.columnStatus[i]=status;
355    }
356    // put buffer into info_
357    strcpy(info_.buffer,buf);
358    delete basis;
359  } else if (mode==2) {
360    CbcModel * babModel = model->model();
361    int iStat = babModel->status();
362    int iStat2 = babModel->secondaryStatus();
363    double value = babModel->getObjValue()*lpSolver->getObjSense();
364    char buf[300];
365    int pos=0;
366    if (iStat==0) {
367      if (babModel->getObjValue()<1.0e40) {
368        pos += sprintf(buf+pos,"optimal," );
369      } else {
370        // infeasible
371        iStat=1;
372        pos += sprintf(buf+pos,"infeasible,");
373      }
374    } else if (iStat==1) {
375      if (iStat2!=6)
376        iStat=3;
377      else
378        iStat=4;
379      std::string minor[]={"","","gap","nodes","time","","solutions","user ctrl-c"};
380      pos += sprintf(buf+pos,"stopped on %s,",minor[iStat2].c_str());
381    } else if (iStat==2) {
382      iStat = 7;
383      pos += sprintf(buf+pos,"stopped on difficulties,");
384    } else if (iStat==5) {
385      iStat = 3;
386      pos += sprintf(buf+pos,"stopped on ctrl-c,");
387    } else {
388      pos += sprintf(buf+pos,"status unknown,");
389      iStat=6;
390    }
391    info_.problemStatus=iStat;
392    info_.objValue = value;
393    if (babModel->getObjValue()<1.0e40) {
394      int precision = ampl_obj_prec();
395      if (precision>0)
396        pos += sprintf(buf+pos," objective %.*g",precision,
397                       value);
398      else
399        pos += sprintf(buf+pos," objective %g",value);
400    }
401    sprintf(buf+pos,"\n%d nodes, %d iterations, %g seconds",
402            babModel->getNodeCount(),
403            babModel->getIterationCount(),
404            totalTime);
405    if (babModel->bestSolution()) {
406      free(info_.primalSolution);
407      info_.primalSolution = (double *) malloc(numberColumns*sizeof(double));
408      CoinCopyN(lpSolver->primalColumnSolution(),numberColumns,info_.primalSolution);
409      free(info_.dualSolution);
410      info_.dualSolution = (double *) malloc(numberRows*sizeof(double));
411      CoinCopyN(lpSolver->dualRowSolution(),numberRows,info_.dualSolution);
412    } else {
413      info_.primalSolution=NULL;
414      info_.dualSolution=NULL;
415    }
416    // put buffer into info
417    strcpy(info_.buffer,buf);
418  } else if (mode==11||mode==12) {
419    // infeasible
420    info_.problemStatus=1;
421    info_.objValue = 1.0e100;
422    sprintf(info_.buffer,"%s",message);
423    info_.primalSolution=NULL;
424    info_.dualSolution=NULL;
425  }
426}
427// Export Data (i.e. at very end)
428void 
429CbcAmpl::exportData(CbcSolver * model)
430{
431  writeAmpl(&info_);
432  freeArrays2(&info_);
433  freeArgs(&info_);
434}
435// Get useful stuff
436void 
437CbcAmpl::fillInformation(CbcSolver * model,
438                         CbcSolverUsefulData & info)
439{
440  memset(&info,0,sizeof(info));
441  info.priorities_ = info_.priorities;
442  info.sosPriority_ = info_.sosPriority;
443  info.branchDirection_ = info_.branchDirection;
444  info.primalSolution_ = info_.primalSolution;
445  info.pseudoDown_ = info_.pseudoDown;
446  info.pseudoUp_ = info_.pseudoUp;
447}
448void addAmplToCbc(CbcSolver * control)
449{
450  CbcAmpl ampl;
451  control->addUserFunction(&ampl);
452}
453extern "C" {
454  //# include "getstub.h"
455# include "asl_pfgh.h"
456}
457//#############################################################################
458// Constructors / Destructor / Assignment
459//#############################################################################
460
461//-------------------------------------------------------------------
462// Default Constructor
463//-------------------------------------------------------------------
464ClpAmplObjective::ClpAmplObjective () 
465: ClpObjective()
466{
467  type_=12;
468  objective_=NULL;
469  amplObjective_=NULL;
470  gradient_ = NULL;
471  offset_ = 0.0;
472}
473// stolen from IPopt with changes
474typedef struct {
475  double obj_sign_;
476  ASL_pfgh * asl_;
477  double * non_const_x_;
478  int * column_; // for jacobian
479  int * rowStart_;
480  double * gradient_;
481  double * constraintValues_;
482  int nz_h_full_; // number of nonzeros in hessian
483  int nerror_;
484  bool objval_called_with_current_x_;
485  bool conval_called_with_current_x_;
486  bool jacval_called_with_current_x_;
487} CbcAmplInfo;
488#if 0
489static bool get_nlp_info(void * amplInfo,int & n, int & m, int & nnz_jac_g,
490                              int & nnz_h_lag)
491{
492  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
493  ASL_pfgh* asl = info->asl_;
494 
495  n = n_var; // # of variables
496  m = n_con; // # of constraints
497  nnz_jac_g = nzc; // # of non-zeros in the jacobian
498  nnz_h_lag = info->nz_h_full_; // # of non-zeros in the hessian
499 
500  return true;
501}
502
503static bool get_bounds_info(void * amplInfo,int  n, double * x_l,
504                            double * x_u, int  m, double * g_l, double * g_u)
505{
506  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
507  ASL_pfgh* asl = info->asl_;
508  assert(n == n_var);
509  assert(m == n_con);
510  int i;
511  for (i=0; i<n; i++) {
512    x_l[i] = LUv[2*i];
513    x_u[i] = LUv[2*i+1];
514  }
515 
516  for (i=0; i<m; i++) {
517    g_l[i] = LUrhs[2*i];
518    g_u[i] = LUrhs[2*i+1];
519  }
520  return true;
521}
522
523#endif
524bool get_constraints_linearity(void * amplInfo,int  n,
525      int * const_types)
526{
527  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
528  ASL_pfgh* asl = info->asl_;
529  //check that n is good
530  assert(n == n_con);
531  // check that there are no network constraints
532  assert(nlnc == 0 && lnc == 0);
533  //the first nlc constraints are non linear the rest is linear
534  int i;
535  for (i=0; i<nlc; i++) {
536    const_types[i]=1;
537  }
538  // the rest is linear
539  for (i=nlc; i<n_con; i++)
540    const_types[i]=0;
541  return true;
542}
543#if 0
544bool get_starting_point(int  n, bool init_x, double * x, bool init_z,
545                        double * z_L, double * z_U, int  m, bool init_lambda, double * lambda)
546{
547  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
548  ASL_pfgh* asl = info->asl_;
549  assert(n == n_var);
550  assert(m == n_con);
551  int i;
552 
553  if (init_x) {
554    for (i=0; i<n; i++) {
555      if (havex0[i]) {
556        x[i] = X0[i];
557      }
558      else {
559        x[i] = 0.0;
560      }
561    }
562  }
563 
564  if (init_z) {
565    for (i=0; i<n; i++) {
566      z_L[i] = z_U[i] = 1.0;
567    }
568  }
569 
570  if (init_lambda) {
571    for (i=0; i<m; i++) {
572      if (havepi0[i]) {
573        lambda[i] = pi0[i];
574      }
575      else {
576        lambda[i] = 0.0;
577      }
578    }
579  }
580 
581  return true;
582}
583#endif
584static bool internal_objval(CbcAmplInfo * info ,double & obj_val)
585{
586  ASL_pfgh* asl = info->asl_;
587  info->objval_called_with_current_x_ = false; // in case the call below fails
588
589  if (n_obj==0) {
590    obj_val = 0;
591    info->objval_called_with_current_x_ = true;
592    return true;
593  }  else {
594    double  retval = objval(0, info->non_const_x_, (fint*)info->nerror_);
595    if (!info->nerror_) {
596      obj_val = info->obj_sign_*retval;
597      info->objval_called_with_current_x_ = true;
598      return true;
599    } else {
600      abort();
601    }
602  }
603 
604  return false;
605}
606static bool internal_conval(CbcAmplInfo * info ,double * g)
607{
608  ASL_pfgh* asl = info->asl_;
609  info->conval_called_with_current_x_ = false; // in case the call below fails
610  assert (g);
611 
612  conval(info->non_const_x_, g, (fint*)info->nerror_);
613
614  if (!info->nerror_) {
615    info->conval_called_with_current_x_ = true;
616    return true;
617  } else {
618    abort();
619  }
620  return false;
621}
622
623static bool apply_new_x(CbcAmplInfo * info  ,bool new_x, int  n, const double * x)
624{
625  ASL_pfgh* asl = info->asl_;
626 
627  if (new_x) {
628    // update the flags so these methods are called
629    // before evaluating the hessian
630    info->conval_called_with_current_x_ = false;
631    info->objval_called_with_current_x_ = false;
632    info->jacval_called_with_current_x_ = false;
633
634    //copy the data to the non_const_x_
635    if (!info->non_const_x_) {
636      info->non_const_x_ = new double [n];
637    }
638
639    for (int  i=0; i<n; i++) {
640      info->non_const_x_[i] = x[i];
641    }
642   
643    // tell ampl that we have a new x
644    xknowne(info->non_const_x_, (fint*)info->nerror_);
645    return info->nerror_ ? false : true;
646  }
647 
648  return true;
649}
650static bool eval_f(void * amplInfo,int  n, const double * x, bool new_x, double & obj_value)
651{
652  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
653  if (!apply_new_x(info,new_x, n, x)) {
654    return false;
655  }
656 
657  return internal_objval(info,obj_value);
658}
659
660static bool eval_grad_f(void * amplInfo,int  n, const double * x, bool new_x, double * grad_f)
661{
662  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
663  ASL_pfgh* asl = info->asl_;
664  if (!apply_new_x(info,new_x, n, x)) {
665    return false;
666  }
667  int i;
668 
669  if (n_obj==0) {
670    for (i=0; i<n; i++) {
671      grad_f[i] = 0.;
672    }
673  }
674  else {
675    objgrd(0, info->non_const_x_, grad_f, (fint*)info->nerror_);
676    if (info->nerror_) {
677      return false;
678    }
679   
680    if (info->obj_sign_==-1) {
681      for (i=0; i<n; i++) {
682        grad_f[i] = -grad_f[i];
683      }
684    }
685  }
686  return true;
687}
688static bool eval_g(void * amplInfo,int  n, const double * x, bool new_x, double * g)
689{
690  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
691  ASL_pfgh* asl = info->asl_;
692  assert(n == n_var);
693 
694  if (!apply_new_x(info,new_x, n, x)) {
695    return false;
696  }
697 
698  return internal_conval(info, g);
699}
700
701static bool eval_jac_g(void * amplInfo,int  n, const double * x, bool new_x,
702                       double * values)
703{
704  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
705  ASL_pfgh* asl = info->asl_;
706  assert(n == n_var);
707 
708  assert (values);
709  if (!apply_new_x(info,new_x, n, x)) {
710    return false;
711  }
712   
713  jacval(info->non_const_x_, values, (fint*)info->nerror_);
714  if (!info->nerror_) {
715    return true;
716  } else {
717    abort();
718  }
719  return false;
720}
721#if 0
722
723static bool eval_h(void * amplInfo,int  n, const double * x, bool new_x,
724            double  obj_factor, int  m, const double * lambda,
725            bool new_lambda, int  nele_hess, int * iRow,
726            int * jCol, double * values)
727{
728  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
729  ASL_pfgh* asl = info->asl_;
730  assert(n == n_var);
731  assert(m == n_con);
732  int i;
733 
734  if (iRow && jCol && !values) {
735    // setup the structure
736    int k=0;
737    for (int i=0; i<n; i++) {
738      for (int j=sputinfo->hcolstarts[i]; j<sputinfo->hcolstarts[i+1]; j++) {
739        //iRow[k] = i + 1;
740        iRow[k] = i;
741        jCol[k] = sputinfo->hrownos[j]+1;
742        k++;
743      }
744    }
745    assert(k==nele_hess);
746    return true;
747  }
748  else if (!iRow & !jCol && values) {
749    if (!apply_new_x(info,new_x, n, x)) {
750      return false;
751    }
752    if (!info->objval_called_with_current_x_) {
753      double  dummy;
754      internal_objval(info,dummy);
755      internal_conval(info,m,NULL);
756    }
757    if (!info->conval_called_with_current_x_) {
758      internal_conval(info,m,NULL);
759    }
760    // copy lambda to non_const_lambda - note, we do not store a copy like
761    // we do with x since lambda is only used here and not in other calls
762    double * non_const_lambda = new double [m];
763    for (i=0; i<m; i++) {
764      non_const_lambda[i] = lambda[i];
765    }
766   
767    real ow=info->obj_sign_*obj_factor;
768    sphes(values, -1, &ow, non_const_lambda);
769   
770    delete [] non_const_lambda;
771    return true;
772  }
773  else {
774    assert(false && "Invalid combination of iRow, jCol, and values pointers");
775  }
776 
777  return false;
778}
779#endif
780//-------------------------------------------------------------------
781// Useful Constructor
782//-------------------------------------------------------------------
783ClpAmplObjective::ClpAmplObjective (void * amplInfo)
784  : ClpObjective()
785{
786  type_=12;
787  activated_=1;
788  gradient_ = NULL;
789  objective_ = NULL;
790  offset_ = 0.0;
791  amplObjective_ = amplInfo;
792}
793
794//-------------------------------------------------------------------
795// Copy constructor
796//-------------------------------------------------------------------
797ClpAmplObjective::ClpAmplObjective (const ClpAmplObjective & rhs) 
798: ClpObjective(rhs)
799{ 
800  amplObjective_ = rhs.amplObjective_;
801  offset_ = rhs.offset_;
802  type_=rhs.type_;
803  if (!amplObjective_) {
804    objective_=NULL;
805    gradient_ = NULL;
806  } else {
807    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
808    ASL_pfgh* asl = info->asl_;
809   
810    int numberColumns = n_var;;
811    if (rhs.objective_) {
812      objective_ = new double [numberColumns];
813      memcpy(objective_,rhs.objective_,numberColumns*sizeof(double));
814    } else {
815      objective_=NULL;
816    }
817    if (rhs.gradient_) {
818      gradient_ = new double [numberColumns];
819      memcpy(gradient_,rhs.gradient_,numberColumns*sizeof(double));
820    } else {
821      gradient_=NULL;
822    }
823  }
824}
825 
826
827//-------------------------------------------------------------------
828// Destructor
829//-------------------------------------------------------------------
830ClpAmplObjective::~ClpAmplObjective ()
831{
832  delete [] objective_;
833  delete [] gradient_;
834}
835
836//----------------------------------------------------------------
837// Assignment operator
838//-------------------------------------------------------------------
839ClpAmplObjective &
840ClpAmplObjective::operator=(const ClpAmplObjective& rhs)
841{
842  if (this != &rhs) {
843    delete [] objective_;
844    delete [] gradient_;
845    amplObjective_ = rhs.amplObjective_;
846    offset_ = rhs.offset_;
847    type_=rhs.type_;
848    if (!amplObjective_) {
849      objective_=NULL;
850      gradient_ = NULL;
851    } else {
852      CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
853      ASL_pfgh* asl = info->asl_;
854     
855      int numberColumns = n_var;;
856      if (rhs.objective_) {
857        objective_ = new double [numberColumns];
858        memcpy(objective_,rhs.objective_,numberColumns*sizeof(double));
859      } else {
860        objective_=NULL;
861      }
862      if (rhs.gradient_) {
863        gradient_ = new double [numberColumns];
864        memcpy(gradient_,rhs.gradient_,numberColumns*sizeof(double));
865      } else {
866        gradient_=NULL;
867      }
868    }
869  }
870  return *this;
871}
872
873// Returns gradient
874double * 
875ClpAmplObjective::gradient(const ClpSimplex * model,
876                                const double * solution, double & offset,bool refresh,
877                                int includeLinear)
878{
879  if (model)
880    assert (model->optimizationDirection()==1.0);
881  bool scaling=false;
882  if (model&&(model->rowScale()||
883              model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0))
884    scaling=true;
885  const double * cost = NULL;
886  if (model)
887    cost = model->costRegion();
888  if (!cost) {
889    // not in solve
890    cost = objective_;
891    scaling=false;
892  }
893  assert (!scaling);
894  if (!amplObjective_||!solution||!activated_) {
895    offset=offset_;
896    return objective_;
897  } else {
898    if (refresh||!gradient_) {
899      CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
900      ASL_pfgh* asl = info->asl_;
901      int numberColumns = n_var;;
902     
903      if (!gradient_) 
904        gradient_ = new double[numberColumns];
905      assert (solution);
906      eval_grad_f(amplObjective_,numberColumns,solution,true,gradient_);
907      // Is this best way?
908      double objValue=0.0;
909      eval_f(amplObjective_,numberColumns,solution,false,objValue);
910      double objValue2=0.0;
911      for (int i=0;i<numberColumns;i++)
912        objValue2 += gradient_[i]*solution[i];
913      offset_ = objValue2 - objValue; // or other way???
914      if (model&&model->optimizationDirection()!=1.0) {
915        offset *= model->optimizationDirection();
916        for (int i=0;i<numberColumns;i++)
917          gradient_[i] *= -1.0;
918      }
919    }
920    offset=offset_;
921    return gradient_;
922  }
923}
924 
925//-------------------------------------------------------------------
926// Clone
927//-------------------------------------------------------------------
928ClpObjective * ClpAmplObjective::clone() const
929{
930  return new ClpAmplObjective(*this);
931}
932// Resize objective
933void 
934ClpAmplObjective::resize(int newNumberColumns)
935{
936  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
937  ASL_pfgh* asl = info->asl_;
938  int numberColumns = n_var;;
939  if (numberColumns!=newNumberColumns) {
940    abort();
941  } 
942 
943}
944// Delete columns in  objective
945void 
946ClpAmplObjective::deleteSome(int numberToDelete, const int * which) 
947{
948  if (numberToDelete)
949    abort();
950}
951/* Returns reduced gradient.Returns an offset (to be added to current one).
952 */
953double 
954ClpAmplObjective::reducedGradient(ClpSimplex * model, double * region,
955                                       bool useFeasibleCosts)
956{
957  int numberRows = model->numberRows();
958  int numberColumns=model->numberColumns();
959 
960  //work space
961  CoinIndexedVector  * workSpace = model->rowArray(0);
962 
963  CoinIndexedVector arrayVector;
964  arrayVector.reserve(numberRows+1);
965 
966  int iRow;
967#ifdef CLP_DEBUG
968  workSpace->checkClear();
969#endif
970  double * array = arrayVector.denseVector();
971  int * index = arrayVector.getIndices();
972  int number=0;
973  const double * costNow = gradient(model,model->solutionRegion(),offset_,
974                                    true,useFeasibleCosts ? 2 : 1);
975  double * cost = model->costRegion();
976  const int * pivotVariable = model->pivotVariable();
977  for (iRow=0;iRow<numberRows;iRow++) {
978    int iPivot=pivotVariable[iRow];
979    double value;
980    if (iPivot<numberColumns)
981      value = costNow[iPivot];
982    else if (!useFeasibleCosts) 
983      value = cost[iPivot];
984    else 
985      value=0.0;
986    if (value) {
987      array[iRow]=value;
988      index[number++]=iRow;
989    }
990  }
991  arrayVector.setNumElements(number);
992
993  // Btran basic costs
994  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
995  double * work = workSpace->denseVector();
996  ClpFillN(work,numberRows,0.0);
997  // now look at dual solution
998  double * rowReducedCost = region+numberColumns;
999  double * dual = rowReducedCost;
1000  const double * rowCost = cost+numberColumns;
1001  for (iRow=0;iRow<numberRows;iRow++) {
1002    dual[iRow]=array[iRow];
1003  }
1004  double * dj = region;
1005  ClpDisjointCopyN(costNow,numberColumns,dj);
1006 
1007  model->transposeTimes(-1.0,dual,dj);
1008  for (iRow=0;iRow<numberRows;iRow++) {
1009    // slack
1010    double value = dual[iRow];
1011    value += rowCost[iRow];
1012    rowReducedCost[iRow]=value;
1013  }
1014  return offset_;
1015}
1016/* Returns step length which gives minimum of objective for
1017   solution + theta * change vector up to maximum theta.
1018   
1019   arrays are numberColumns+numberRows
1020*/
1021double 
1022ClpAmplObjective::stepLength(ClpSimplex * model,
1023                                  const double * solution,
1024                                  const double * change,
1025                                  double maximumTheta,
1026                                  double & currentObj,
1027                                  double & predictedObj,
1028                                  double & thetaObj)
1029{
1030  // Assume convex
1031  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1032  ASL_pfgh* asl = info->asl_;
1033 
1034  int numberColumns = n_var;;
1035  double * tempSolution = new double [numberColumns];
1036  double * tempGradient = new double [numberColumns];
1037  // current
1038  eval_f(amplObjective_,numberColumns,solution,true,currentObj);
1039  double objA = currentObj;
1040  double thetaA = 0.0;
1041  // at maximum
1042  int i;
1043  for (i=0;i<numberColumns;i++)
1044    tempSolution[i] = solution[i] + maximumTheta*change[i];
1045  eval_f(amplObjective_,numberColumns,tempSolution,true,thetaObj);
1046  double objC = thetaObj;
1047  double thetaC = maximumTheta;
1048  double objB=0.5*(objA+objC);
1049  double thetaB=0.5*maximumTheta;
1050  double gradientNorm=1.0e6;
1051  while (gradientNorm>1.0e-6&&thetaC-thetaA>1.0e-8) {
1052    for (i=0;i<numberColumns;i++)
1053      tempSolution[i] = solution[i] + thetaB*change[i];
1054    eval_grad_f(amplObjective_,numberColumns,tempSolution,true,tempGradient);
1055    eval_f(amplObjective_,numberColumns,tempSolution,false,objB);
1056    double changeObj=0.0;
1057    gradientNorm=0.0;
1058    for (i=0;i<numberColumns;i++) {
1059      changeObj += tempGradient[i]*change[i];
1060      gradientNorm += tempGradient[i]*tempGradient[i];
1061    }
1062    gradientNorm = fabs(changeObj)/sqrt(gradientNorm);
1063    // Should try and get quadratic convergence by interpolation
1064    if (changeObj<0.0) {
1065      // increasing is good
1066      thetaA = thetaB;
1067    } else {
1068      // decreasing is good
1069      thetaC = thetaB;
1070    }
1071    thetaB = 0.5*(thetaA+thetaC);
1072  }
1073  delete [] tempSolution;
1074  delete [] tempGradient;
1075  predictedObj = objB;
1076  return thetaB;
1077}
1078// Return objective value (without any ClpModel offset) (model may be NULL)
1079double 
1080ClpAmplObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
1081{
1082  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1083  ASL_pfgh* asl = info->asl_;
1084 
1085  int numberColumns = n_var;;
1086  // current
1087  double currentObj=0.0;
1088  eval_f(amplObjective_,numberColumns,solution,true,currentObj);
1089  return currentObj;
1090}
1091// Scale objective
1092void 
1093ClpAmplObjective::reallyScale(const double * columnScale) 
1094{
1095  abort();
1096}
1097/* Given a zeroed array sets nonlinear columns to 1.
1098   Returns number of nonlinear columns
1099*/
1100int 
1101ClpAmplObjective::markNonlinear(char * which)
1102{
1103  int iColumn;
1104  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1105  ASL_pfgh* asl = info->asl_;
1106  int nonLinear = CoinMax(nlvc,nlvo);
1107  for (iColumn=0;iColumn<nonLinear;iColumn++) {
1108    which[iColumn]=1;
1109  }
1110  int numberNonLinearColumns = 0;
1111  int numberColumns = n_var;;
1112  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1113    if(which[iColumn])
1114      numberNonLinearColumns++;
1115  }
1116  return numberNonLinearColumns;
1117}
1118// Say we have new primal solution - so may need to recompute
1119void 
1120ClpAmplObjective::newXValues() 
1121{
1122  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1123  info->conval_called_with_current_x_ = false;
1124  info->objval_called_with_current_x_ = false;
1125  info->jacval_called_with_current_x_ = false;
1126}
1127//#############################################################################
1128// Constructors / Destructor / Assignment
1129//#############################################################################
1130//-------------------------------------------------------------------
1131// Default Constructor
1132//-------------------------------------------------------------------
1133ClpConstraintAmpl::ClpConstraintAmpl () 
1134: ClpConstraint()
1135{
1136  type_=3;
1137  column_=NULL;
1138  coefficient_ = NULL;
1139  numberCoefficients_=0;
1140  amplInfo_ = NULL;
1141}
1142
1143//-------------------------------------------------------------------
1144// Useful Constructor
1145//-------------------------------------------------------------------
1146ClpConstraintAmpl::ClpConstraintAmpl (int row, void * amplInfo)
1147  : ClpConstraint()
1148{
1149  type_=3;
1150  rowNumber_=row;
1151  amplInfo_ = amplInfo;
1152  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1153  ASL_pfgh* asl = info->asl_;
1154  assert (rowNumber_<nlc);
1155  numberCoefficients_ = info->rowStart_[rowNumber_+1]-info->rowStart_[rowNumber_];
1156  column_ = CoinCopyOfArray(info->column_+info->rowStart_[rowNumber_],numberCoefficients_);
1157  coefficient_ = new double [numberCoefficients_];;
1158}
1159
1160//-------------------------------------------------------------------
1161// Copy constructor
1162//-------------------------------------------------------------------
1163ClpConstraintAmpl::ClpConstraintAmpl (const ClpConstraintAmpl & rhs)
1164: ClpConstraint(rhs)
1165{ 
1166  numberCoefficients_=rhs.numberCoefficients_;
1167  column_ = CoinCopyOfArray(rhs.column_,numberCoefficients_);
1168  coefficient_ = CoinCopyOfArray(rhs.coefficient_,numberCoefficients_);
1169}
1170 
1171
1172//-------------------------------------------------------------------
1173// Destructor
1174//-------------------------------------------------------------------
1175ClpConstraintAmpl::~ClpConstraintAmpl ()
1176{
1177  delete [] column_;
1178  delete [] coefficient_;
1179}
1180
1181//----------------------------------------------------------------
1182// Assignment operator
1183//-------------------------------------------------------------------
1184ClpConstraintAmpl &
1185ClpConstraintAmpl::operator=(const ClpConstraintAmpl& rhs)
1186{
1187  if (this != &rhs) {
1188    delete [] column_;
1189    delete [] coefficient_;
1190    numberCoefficients_=rhs.numberCoefficients_;
1191    column_ = CoinCopyOfArray(rhs.column_,numberCoefficients_);
1192    coefficient_ = CoinCopyOfArray(rhs.coefficient_,numberCoefficients_);
1193  }
1194  return *this;
1195}
1196//-------------------------------------------------------------------
1197// Clone
1198//-------------------------------------------------------------------
1199ClpConstraint * ClpConstraintAmpl::clone() const
1200{
1201  return new ClpConstraintAmpl(*this);
1202}
1203
1204// Returns gradient
1205int 
1206ClpConstraintAmpl::gradient(const ClpSimplex * model,
1207                              const double * solution, 
1208                              double * gradient,
1209                              double & functionValue, 
1210                              double & offset,
1211                              bool useScaling,
1212                              bool refresh) const
1213{
1214  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1215  ASL_pfgh* asl = info->asl_;
1216  int numberColumns = n_var;;
1217  // If not done then do all
1218  if (!info->jacval_called_with_current_x_) {
1219    bool getStuff = eval_g(amplInfo_,numberColumns,solution,true,info->constraintValues_);
1220    assert (getStuff);
1221    getStuff = eval_jac_g(amplInfo_,numberColumns,solution,false,info->gradient_);
1222    assert (getStuff);
1223    info->jacval_called_with_current_x_=getStuff;
1224  }
1225  if (refresh||!lastGradient_) {
1226    functionValue_ = info->constraintValues_[rowNumber_]; 
1227    offset_=functionValue_; // sign??
1228    if (!lastGradient_)
1229      lastGradient_ = new double[numberColumns];
1230    CoinZeroN(lastGradient_,numberColumns);
1231    bool scaling=(model&&model->rowScale()&&useScaling);
1232    assert (!scaling);
1233    int i;
1234    int start = info->rowStart_[rowNumber_];
1235    assert (numberCoefficients_==info->rowStart_[rowNumber_+1]-start);
1236    for (i=0;i<numberCoefficients_;i++) {
1237      int iColumn= column_[i];
1238      double valueS = solution[iColumn];
1239      double valueG = info->gradient_[start+i];
1240      lastGradient_[iColumn]=valueG;
1241      offset_ -= valueS*valueG;
1242    }
1243  }
1244  functionValue = functionValue_;
1245  offset = offset_;
1246  memcpy(gradient,lastGradient_,numberColumns*sizeof(double));
1247  return 0;
1248}
1249// Resize constraint
1250void 
1251ClpConstraintAmpl::resize(int newNumberColumns)
1252{
1253  abort();
1254}
1255// Delete columns in  constraint
1256void 
1257ClpConstraintAmpl::deleteSome(int numberToDelete, const int * which) 
1258{
1259  if (numberToDelete) {
1260    abort();
1261  }
1262}
1263// Scale constraint
1264void 
1265ClpConstraintAmpl::reallyScale(const double * columnScale) 
1266{
1267  abort();
1268}
1269/* Given a zeroed array sets nonlinear columns to 1.
1270   Returns number of nonlinear columns
1271*/
1272int 
1273ClpConstraintAmpl::markNonlinear(char * which) const
1274{
1275  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1276  ASL_pfgh* asl = info->asl_;
1277  int iColumn;
1278  int numberNon=0;
1279  int nonLinear = CoinMax(nlvc,nlvo);
1280  for (iColumn=0;iColumn<numberCoefficients_;iColumn++) {
1281    int jColumn = column_[iColumn];
1282    if (jColumn<nonLinear) {
1283      which[jColumn]=1;
1284      numberNon++;
1285    }
1286  }
1287  return numberNon;
1288}
1289/* Given a zeroed array sets possible nonzero coefficients to 1.
1290   Returns number of nonzeros
1291*/
1292int 
1293ClpConstraintAmpl::markNonzero(char * which) const
1294{
1295  int iColumn;
1296  for (iColumn=0;iColumn<numberCoefficients_;iColumn++) {
1297    which[column_[iColumn]]=1;
1298  }
1299  return numberCoefficients_;
1300}
1301// Number of coefficients
1302int 
1303ClpConstraintAmpl::numberCoefficients() const
1304{
1305  return numberCoefficients_;
1306}
1307// Say we have new primal solution - so may need to recompute
1308void 
1309ClpConstraintAmpl::newXValues() 
1310{
1311  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1312  info->conval_called_with_current_x_ = false;
1313  info->objval_called_with_current_x_ = false;
1314  info->jacval_called_with_current_x_ = false;
1315}
1316/* Load nonlinear part of problem from AMPL info
1317   Returns 0 if linear
1318   1 if quadratic objective
1319   2 if quadratic constraints
1320   3 if nonlinear objective
1321   4 if nonlinear constraints
1322   -1 on failure
1323*/
1324int 
1325ClpSimplex::loadNonLinear(void * amplInfo, int & numberConstraints, 
1326                          ClpConstraint ** & constraints)
1327{
1328  numberConstraints=0;
1329  constraints=NULL;
1330  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
1331  ASL_pfgh* asl = info->asl_;
1332  // For moment don't say quadratic
1333  int type=0;
1334  if (nlo+nlc) {
1335    // nonlinear
1336    if (!nlc) {
1337      type=3;
1338      delete objective_;
1339      objective_= new ClpAmplObjective(amplInfo);
1340    } else {
1341      type=4;
1342      numberConstraints = nlc;
1343      constraints = new ClpConstraint * [numberConstraints];
1344      if (nlo) {
1345        delete objective_;
1346        objective_= new ClpAmplObjective(amplInfo);
1347      }
1348      for (int i=0;i<numberConstraints;i++) {
1349        constraints[i] = new ClpConstraintAmpl(i,amplInfo);
1350      }
1351    }
1352  }
1353  return type;
1354}
1355#else
1356#include "ClpSimplex.hpp"
1357#include "ClpConstraint.hpp"
1358int 
1359ClpSimplex::loadNonLinear(void * amplInfo, int & numberConstraints, 
1360                          ClpConstraint ** & constraints)
1361{
1362  abort();
1363  return 0;
1364}
1365#endif
Note: See TracBrowser for help on using the repository browser.