source: stable/2.3/Cbc/src/ClpAmplStuff.cpp @ 1203

Last change on this file since 1203 was 1203, checked in by tkr, 10 years ago

Merging r1179:1197 from stable/BSP/ stable/2.3

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