source: stable/BSP/Cbc/src/ClpAmplStuff.cpp @ 1185

Last change on this file since 1185 was 1185, checked in by stefan, 11 years ago

fix compiler warnings

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