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

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

need variable asl for ASL macro in order to have assert working

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.