source: stable/2.4/Cbc/src/ClpAmplStuff.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

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