source: branches/heur/Cbc/src/ClpAmplStuff.cpp @ 885

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

for deterministic parallel

File size: 38.0 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    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  ASL_pfgh* asl = info->asl_;
694  assert(n == n_var);
695 
696  if (!apply_new_x(info,new_x, n, x)) {
697    return false;
698  }
699 
700  return internal_conval(info, g);
701}
702
703static bool eval_jac_g(void * amplInfo,int  n, const double * x, bool new_x,
704                       double * values)
705{
706  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
707  ASL_pfgh* asl = info->asl_;
708  assert(n == n_var);
709 
710  assert (values);
711  if (!apply_new_x(info,new_x, n, x)) {
712    return false;
713  }
714   
715  jacval(info->non_const_x_, values, (fint*)info->nerror_);
716  if (!info->nerror_) {
717    return true;
718  } else {
719    abort();
720  }
721  return false;
722}
723#if 0
724
725static bool eval_h(void * amplInfo,int  n, const double * x, bool new_x,
726            double  obj_factor, int  m, const double * lambda,
727            bool new_lambda, int  nele_hess, int * iRow,
728            int * jCol, double * values)
729{
730  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
731  ASL_pfgh* asl = info->asl_;
732  assert(n == n_var);
733  assert(m == n_con);
734  int i;
735 
736  if (iRow && jCol && !values) {
737    // setup the structure
738    int k=0;
739    for (int i=0; i<n; i++) {
740      for (int j=sputinfo->hcolstarts[i]; j<sputinfo->hcolstarts[i+1]; j++) {
741        //iRow[k] = i + 1;
742        iRow[k] = i;
743        jCol[k] = sputinfo->hrownos[j]+1;
744        k++;
745      }
746    }
747    assert(k==nele_hess);
748    return true;
749  }
750  else if (!iRow & !jCol && values) {
751    if (!apply_new_x(info,new_x, n, x)) {
752      return false;
753    }
754    if (!info->objval_called_with_current_x_) {
755      double  dummy;
756      internal_objval(info,dummy);
757      internal_conval(info,m,NULL);
758    }
759    if (!info->conval_called_with_current_x_) {
760      internal_conval(info,m,NULL);
761    }
762    // copy lambda to non_const_lambda - note, we do not store a copy like
763    // we do with x since lambda is only used here and not in other calls
764    double * non_const_lambda = new double [m];
765    for (i=0; i<m; i++) {
766      non_const_lambda[i] = lambda[i];
767    }
768   
769    real ow=info->obj_sign_*obj_factor;
770    sphes(values, -1, &ow, non_const_lambda);
771   
772    delete [] non_const_lambda;
773    return true;
774  }
775  else {
776    assert(false && "Invalid combination of iRow, jCol, and values pointers");
777  }
778 
779  return false;
780}
781#endif
782//-------------------------------------------------------------------
783// Useful Constructor
784//-------------------------------------------------------------------
785ClpAmplObjective::ClpAmplObjective (void * amplInfo)
786  : ClpObjective()
787{
788  type_=12;
789  activated_=1;
790  gradient_ = NULL;
791  objective_ = NULL;
792  offset_ = 0.0;
793  amplObjective_ = amplInfo;
794}
795
796//-------------------------------------------------------------------
797// Copy constructor
798//-------------------------------------------------------------------
799ClpAmplObjective::ClpAmplObjective (const ClpAmplObjective & rhs) 
800: ClpObjective(rhs)
801{ 
802  amplObjective_ = rhs.amplObjective_;
803  offset_ = rhs.offset_;
804  type_=rhs.type_;
805  if (!amplObjective_) {
806    objective_=NULL;
807    gradient_ = NULL;
808  } else {
809    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
810    ASL_pfgh* asl = info->asl_;
811   
812    int numberColumns = n_var;;
813    if (rhs.objective_) {
814      objective_ = new double [numberColumns];
815      memcpy(objective_,rhs.objective_,numberColumns*sizeof(double));
816    } else {
817      objective_=NULL;
818    }
819    if (rhs.gradient_) {
820      gradient_ = new double [numberColumns];
821      memcpy(gradient_,rhs.gradient_,numberColumns*sizeof(double));
822    } else {
823      gradient_=NULL;
824    }
825  }
826}
827 
828
829//-------------------------------------------------------------------
830// Destructor
831//-------------------------------------------------------------------
832ClpAmplObjective::~ClpAmplObjective ()
833{
834  delete [] objective_;
835  delete [] gradient_;
836}
837
838//----------------------------------------------------------------
839// Assignment operator
840//-------------------------------------------------------------------
841ClpAmplObjective &
842ClpAmplObjective::operator=(const ClpAmplObjective& rhs)
843{
844  if (this != &rhs) {
845    delete [] objective_;
846    delete [] gradient_;
847    amplObjective_ = rhs.amplObjective_;
848    offset_ = rhs.offset_;
849    type_=rhs.type_;
850    if (!amplObjective_) {
851      objective_=NULL;
852      gradient_ = NULL;
853    } else {
854      CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
855      ASL_pfgh* asl = info->asl_;
856     
857      int numberColumns = n_var;;
858      if (rhs.objective_) {
859        objective_ = new double [numberColumns];
860        memcpy(objective_,rhs.objective_,numberColumns*sizeof(double));
861      } else {
862        objective_=NULL;
863      }
864      if (rhs.gradient_) {
865        gradient_ = new double [numberColumns];
866        memcpy(gradient_,rhs.gradient_,numberColumns*sizeof(double));
867      } else {
868        gradient_=NULL;
869      }
870    }
871  }
872  return *this;
873}
874
875// Returns gradient
876double * 
877ClpAmplObjective::gradient(const ClpSimplex * model,
878                                const double * solution, double & offset,bool refresh,
879                                int includeLinear)
880{
881  if (model)
882    assert (model->optimizationDirection()==1.0);
883  bool scaling=false;
884  if (model&&(model->rowScale()||
885              model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0))
886    scaling=true;
887  const double * cost = NULL;
888  if (model)
889    cost = model->costRegion();
890  if (!cost) {
891    // not in solve
892    cost = objective_;
893    scaling=false;
894  }
895  assert (!scaling);
896  if (!amplObjective_||!solution||!activated_) {
897    offset=offset_;
898    return objective_;
899  } else {
900    if (refresh||!gradient_) {
901      CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
902      ASL_pfgh* asl = info->asl_;
903      int numberColumns = n_var;;
904     
905      if (!gradient_) 
906        gradient_ = new double[numberColumns];
907      assert (solution);
908      eval_grad_f(amplObjective_,numberColumns,solution,true,gradient_);
909      // Is this best way?
910      double objValue=0.0;
911      eval_f(amplObjective_,numberColumns,solution,false,objValue);
912      double objValue2=0.0;
913      for (int i=0;i<numberColumns;i++)
914        objValue2 += gradient_[i]*solution[i];
915      offset_ = objValue2 - objValue; // or other way???
916      if (model&&model->optimizationDirection()!=1.0) {
917        offset *= model->optimizationDirection();
918        for (int i=0;i<numberColumns;i++)
919          gradient_[i] *= -1.0;
920      }
921    }
922    offset=offset_;
923    return gradient_;
924  }
925}
926 
927//-------------------------------------------------------------------
928// Clone
929//-------------------------------------------------------------------
930ClpObjective * ClpAmplObjective::clone() const
931{
932  return new ClpAmplObjective(*this);
933}
934// Resize objective
935void 
936ClpAmplObjective::resize(int newNumberColumns)
937{
938  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
939  ASL_pfgh* asl = info->asl_;
940  int numberColumns = n_var;;
941  if (numberColumns!=newNumberColumns) {
942    abort();
943  } 
944 
945}
946// Delete columns in  objective
947void 
948ClpAmplObjective::deleteSome(int numberToDelete, const int * which) 
949{
950  if (numberToDelete)
951    abort();
952}
953/* Returns reduced gradient.Returns an offset (to be added to current one).
954 */
955double 
956ClpAmplObjective::reducedGradient(ClpSimplex * model, double * region,
957                                       bool useFeasibleCosts)
958{
959  int numberRows = model->numberRows();
960  int numberColumns=model->numberColumns();
961 
962  //work space
963  CoinIndexedVector  * workSpace = model->rowArray(0);
964 
965  CoinIndexedVector arrayVector;
966  arrayVector.reserve(numberRows+1);
967 
968  int iRow;
969#ifdef CLP_DEBUG
970  workSpace->checkClear();
971#endif
972  double * array = arrayVector.denseVector();
973  int * index = arrayVector.getIndices();
974  int number=0;
975  const double * costNow = gradient(model,model->solutionRegion(),offset_,
976                                    true,useFeasibleCosts ? 2 : 1);
977  double * cost = model->costRegion();
978  const int * pivotVariable = model->pivotVariable();
979  for (iRow=0;iRow<numberRows;iRow++) {
980    int iPivot=pivotVariable[iRow];
981    double value;
982    if (iPivot<numberColumns)
983      value = costNow[iPivot];
984    else if (!useFeasibleCosts) 
985      value = cost[iPivot];
986    else 
987      value=0.0;
988    if (value) {
989      array[iRow]=value;
990      index[number++]=iRow;
991    }
992  }
993  arrayVector.setNumElements(number);
994
995  // Btran basic costs
996  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
997  double * work = workSpace->denseVector();
998  ClpFillN(work,numberRows,0.0);
999  // now look at dual solution
1000  double * rowReducedCost = region+numberColumns;
1001  double * dual = rowReducedCost;
1002  const double * rowCost = cost+numberColumns;
1003  for (iRow=0;iRow<numberRows;iRow++) {
1004    dual[iRow]=array[iRow];
1005  }
1006  double * dj = region;
1007  ClpDisjointCopyN(costNow,numberColumns,dj);
1008 
1009  model->transposeTimes(-1.0,dual,dj);
1010  for (iRow=0;iRow<numberRows;iRow++) {
1011    // slack
1012    double value = dual[iRow];
1013    value += rowCost[iRow];
1014    rowReducedCost[iRow]=value;
1015  }
1016  return offset_;
1017}
1018/* Returns step length which gives minimum of objective for
1019   solution + theta * change vector up to maximum theta.
1020   
1021   arrays are numberColumns+numberRows
1022*/
1023double 
1024ClpAmplObjective::stepLength(ClpSimplex * model,
1025                                  const double * solution,
1026                                  const double * change,
1027                                  double maximumTheta,
1028                                  double & currentObj,
1029                                  double & predictedObj,
1030                                  double & thetaObj)
1031{
1032  // Assume convex
1033  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1034  ASL_pfgh* asl = info->asl_;
1035 
1036  int numberColumns = n_var;;
1037  double * tempSolution = new double [numberColumns];
1038  double * tempGradient = new double [numberColumns];
1039  // current
1040  eval_f(amplObjective_,numberColumns,solution,true,currentObj);
1041  double objA = currentObj;
1042  double thetaA = 0.0;
1043  // at maximum
1044  int i;
1045  for (i=0;i<numberColumns;i++)
1046    tempSolution[i] = solution[i] + maximumTheta*change[i];
1047  eval_f(amplObjective_,numberColumns,tempSolution,true,thetaObj);
1048  double objC = thetaObj;
1049  double thetaC = maximumTheta;
1050  double objB=0.5*(objA+objC);
1051  double thetaB=0.5*maximumTheta;
1052  double gradientNorm=1.0e6;
1053  while (gradientNorm>1.0e-6&&thetaC-thetaA>1.0e-8) {
1054    for (i=0;i<numberColumns;i++)
1055      tempSolution[i] = solution[i] + thetaB*change[i];
1056    eval_grad_f(amplObjective_,numberColumns,tempSolution,true,tempGradient);
1057    eval_f(amplObjective_,numberColumns,tempSolution,false,objB);
1058    double changeObj=0.0;
1059    gradientNorm=0.0;
1060    for (i=0;i<numberColumns;i++) {
1061      changeObj += tempGradient[i]*change[i];
1062      gradientNorm += tempGradient[i]*tempGradient[i];
1063    }
1064    gradientNorm = fabs(changeObj)/sqrt(gradientNorm);
1065    // Should try and get quadratic convergence by interpolation
1066    if (changeObj<0.0) {
1067      // increasing is good
1068      thetaA = thetaB;
1069    } else {
1070      // decreasing is good
1071      thetaC = thetaB;
1072    }
1073    thetaB = 0.5*(thetaA+thetaC);
1074  }
1075  delete [] tempSolution;
1076  delete [] tempGradient;
1077  predictedObj = objB;
1078  return thetaB;
1079}
1080// Return objective value (without any ClpModel offset) (model may be NULL)
1081double 
1082ClpAmplObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
1083{
1084  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1085  ASL_pfgh* asl = info->asl_;
1086 
1087  int numberColumns = n_var;;
1088  // current
1089  double currentObj=0.0;
1090  eval_f(amplObjective_,numberColumns,solution,true,currentObj);
1091  return currentObj;
1092}
1093// Scale objective
1094void 
1095ClpAmplObjective::reallyScale(const double * columnScale) 
1096{
1097  abort();
1098}
1099/* Given a zeroed array sets nonlinear columns to 1.
1100   Returns number of nonlinear columns
1101*/
1102int 
1103ClpAmplObjective::markNonlinear(char * which)
1104{
1105  int iColumn;
1106  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1107  ASL_pfgh* asl = info->asl_;
1108  int nonLinear = CoinMax(nlvc,nlvo);
1109  for (iColumn=0;iColumn<nonLinear;iColumn++) {
1110    which[iColumn]=1;
1111  }
1112  int numberNonLinearColumns = 0;
1113  int numberColumns = n_var;;
1114  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1115    if(which[iColumn])
1116      numberNonLinearColumns++;
1117  }
1118  return numberNonLinearColumns;
1119}
1120// Say we have new primal solution - so may need to recompute
1121void 
1122ClpAmplObjective::newXValues() 
1123{
1124  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1125  info->conval_called_with_current_x_ = false;
1126  info->objval_called_with_current_x_ = false;
1127  info->jacval_called_with_current_x_ = false;
1128}
1129//#############################################################################
1130// Constructors / Destructor / Assignment
1131//#############################################################################
1132//-------------------------------------------------------------------
1133// Default Constructor
1134//-------------------------------------------------------------------
1135ClpConstraintAmpl::ClpConstraintAmpl () 
1136: ClpConstraint()
1137{
1138  type_=3;
1139  column_=NULL;
1140  coefficient_ = NULL;
1141  numberCoefficients_=0;
1142  amplInfo_ = NULL;
1143}
1144
1145//-------------------------------------------------------------------
1146// Useful Constructor
1147//-------------------------------------------------------------------
1148ClpConstraintAmpl::ClpConstraintAmpl (int row, void * amplInfo)
1149  : ClpConstraint()
1150{
1151  type_=3;
1152  rowNumber_=row;
1153  amplInfo_ = amplInfo;
1154  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1155  ASL_pfgh* asl = info->asl_;
1156  assert (rowNumber_<nlc);
1157  numberCoefficients_ = info->rowStart_[rowNumber_+1]-info->rowStart_[rowNumber_];
1158  column_ = CoinCopyOfArray(info->column_+info->rowStart_[rowNumber_],numberCoefficients_);
1159  coefficient_ = new double [numberCoefficients_];;
1160}
1161
1162//-------------------------------------------------------------------
1163// Copy constructor
1164//-------------------------------------------------------------------
1165ClpConstraintAmpl::ClpConstraintAmpl (const ClpConstraintAmpl & rhs)
1166: ClpConstraint(rhs)
1167{ 
1168  numberCoefficients_=rhs.numberCoefficients_;
1169  column_ = CoinCopyOfArray(rhs.column_,numberCoefficients_);
1170  coefficient_ = CoinCopyOfArray(rhs.coefficient_,numberCoefficients_);
1171}
1172 
1173
1174//-------------------------------------------------------------------
1175// Destructor
1176//-------------------------------------------------------------------
1177ClpConstraintAmpl::~ClpConstraintAmpl ()
1178{
1179  delete [] column_;
1180  delete [] coefficient_;
1181}
1182
1183//----------------------------------------------------------------
1184// Assignment operator
1185//-------------------------------------------------------------------
1186ClpConstraintAmpl &
1187ClpConstraintAmpl::operator=(const ClpConstraintAmpl& rhs)
1188{
1189  if (this != &rhs) {
1190    delete [] column_;
1191    delete [] coefficient_;
1192    numberCoefficients_=rhs.numberCoefficients_;
1193    column_ = CoinCopyOfArray(rhs.column_,numberCoefficients_);
1194    coefficient_ = CoinCopyOfArray(rhs.coefficient_,numberCoefficients_);
1195  }
1196  return *this;
1197}
1198//-------------------------------------------------------------------
1199// Clone
1200//-------------------------------------------------------------------
1201ClpConstraint * ClpConstraintAmpl::clone() const
1202{
1203  return new ClpConstraintAmpl(*this);
1204}
1205
1206// Returns gradient
1207int 
1208ClpConstraintAmpl::gradient(const ClpSimplex * model,
1209                              const double * solution, 
1210                              double * gradient,
1211                              double & functionValue, 
1212                              double & offset,
1213                              bool useScaling,
1214                              bool refresh) const
1215{
1216  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1217  ASL_pfgh* asl = info->asl_;
1218  int numberColumns = n_var;;
1219  // If not done then do all
1220  if (!info->jacval_called_with_current_x_) {
1221    bool getStuff = eval_g(amplInfo_,numberColumns,solution,true,info->constraintValues_);
1222    assert (getStuff);
1223    getStuff = eval_jac_g(amplInfo_,numberColumns,solution,false,info->gradient_);
1224    assert (getStuff);
1225    info->jacval_called_with_current_x_=getStuff;
1226  }
1227  if (refresh||!lastGradient_) {
1228    functionValue_ = info->constraintValues_[rowNumber_]; 
1229    offset_=functionValue_; // sign??
1230    if (!lastGradient_)
1231      lastGradient_ = new double[numberColumns];
1232    CoinZeroN(lastGradient_,numberColumns);
1233    bool scaling=(model&&model->rowScale()&&useScaling);
1234    assert (!scaling);
1235    int i;
1236    int start = info->rowStart_[rowNumber_];
1237    assert (numberCoefficients_==info->rowStart_[rowNumber_+1]-start);
1238    for (i=0;i<numberCoefficients_;i++) {
1239      int iColumn= column_[i];
1240      double valueS = solution[iColumn];
1241      double valueG = info->gradient_[start+i];
1242      lastGradient_[iColumn]=valueG;
1243      offset_ -= valueS*valueG;
1244    }
1245  }
1246  functionValue = functionValue_;
1247  offset = offset_;
1248  memcpy(gradient,lastGradient_,numberColumns*sizeof(double));
1249  return 0;
1250}
1251// Resize constraint
1252void 
1253ClpConstraintAmpl::resize(int newNumberColumns)
1254{
1255  abort();
1256}
1257// Delete columns in  constraint
1258void 
1259ClpConstraintAmpl::deleteSome(int numberToDelete, const int * which) 
1260{
1261  if (numberToDelete) {
1262    abort();
1263  }
1264}
1265// Scale constraint
1266void 
1267ClpConstraintAmpl::reallyScale(const double * columnScale) 
1268{
1269  abort();
1270}
1271/* Given a zeroed array sets nonlinear columns to 1.
1272   Returns number of nonlinear columns
1273*/
1274int 
1275ClpConstraintAmpl::markNonlinear(char * which) const
1276{
1277  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1278  ASL_pfgh* asl = info->asl_;
1279  int iColumn;
1280  int numberNon=0;
1281  int nonLinear = CoinMax(nlvc,nlvo);
1282  for (iColumn=0;iColumn<numberCoefficients_;iColumn++) {
1283    int jColumn = column_[iColumn];
1284    if (jColumn<nonLinear) {
1285      which[jColumn]=1;
1286      numberNon++;
1287    }
1288  }
1289  return numberNon;
1290}
1291/* Given a zeroed array sets possible nonzero coefficients to 1.
1292   Returns number of nonzeros
1293*/
1294int 
1295ClpConstraintAmpl::markNonzero(char * which) const
1296{
1297  int iColumn;
1298  for (iColumn=0;iColumn<numberCoefficients_;iColumn++) {
1299    which[column_[iColumn]]=1;
1300  }
1301  return numberCoefficients_;
1302}
1303// Number of coefficients
1304int 
1305ClpConstraintAmpl::numberCoefficients() const
1306{
1307  return numberCoefficients_;
1308}
1309// Say we have new primal solution - so may need to recompute
1310void 
1311ClpConstraintAmpl::newXValues() 
1312{
1313  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1314  info->conval_called_with_current_x_ = false;
1315  info->objval_called_with_current_x_ = false;
1316  info->jacval_called_with_current_x_ = false;
1317}
1318/* Load nonlinear part of problem from AMPL info
1319   Returns 0 if linear
1320   1 if quadratic objective
1321   2 if quadratic constraints
1322   3 if nonlinear objective
1323   4 if nonlinear constraints
1324   -1 on failure
1325*/
1326int 
1327ClpSimplex::loadNonLinear(void * amplInfo, int & numberConstraints, 
1328                          ClpConstraint ** & constraints)
1329{
1330  numberConstraints=0;
1331  constraints=NULL;
1332  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
1333  ASL_pfgh* asl = info->asl_;
1334  // For moment don't say quadratic
1335  int type=0;
1336  if (nlo+nlc) {
1337    // nonlinear
1338    if (!nlc) {
1339      type=3;
1340      delete objective_;
1341      objective_= new ClpAmplObjective(amplInfo);
1342    } else {
1343      type=4;
1344      numberConstraints = nlc;
1345      constraints = new ClpConstraint * [numberConstraints];
1346      if (nlo) {
1347        delete objective_;
1348        objective_= new ClpAmplObjective(amplInfo);
1349      }
1350      for (int i=0;i<numberConstraints;i++) {
1351        constraints[i] = new ClpConstraintAmpl(i,amplInfo);
1352      }
1353    }
1354  }
1355  return type;
1356}
1357#else
1358#include "ClpSimplex.hpp"
1359#include "ClpConstraint.hpp"
1360int 
1361ClpSimplex::loadNonLinear(void * amplInfo, int & numberConstraints, 
1362                          ClpConstraint ** & constraints)
1363{
1364  abort();
1365  return 0;
1366}
1367#endif
Note: See TracBrowser for help on using the repository browser.