source: trunk/Cbc/src/ClpAmplStuff.cpp @ 700

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

changes for postprocess loop and LOS

File size: 24.5 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 "CoinUtilsConfig.h"
15#include "CoinHelperFunctions.hpp"
16extern "C" {
17  //# include "getstub.h"
18# include "asl_pfgh.h"
19}
20//#############################################################################
21// Constructors / Destructor / Assignment
22//#############################################################################
23
24//-------------------------------------------------------------------
25// Default Constructor
26//-------------------------------------------------------------------
27ClpAmplObjective::ClpAmplObjective () 
28: ClpObjective()
29{
30  type_=12;
31  objective_=NULL;
32  amplObjective_=NULL;
33  gradient_ = NULL;
34  offset_ = 0.0;
35}
36// stolen from IPopt with changes
37typedef struct {
38  double obj_sign_;
39  ASL_pfgh * asl_;
40  double * non_const_x_;
41  int * column_; // for jacobian
42  int * rowStart_;
43  double * gradient_;
44  double * constraintValues_;
45  int nz_h_full_; // number of nonzeros in hessian
46  int nerror_;
47  bool objval_called_with_current_x_;
48  bool conval_called_with_current_x_;
49  bool jacval_called_with_current_x_;
50} CbcAmplInfo;
51#if 0
52static bool get_nlp_info(void * amplInfo,int & n, int & m, int & nnz_jac_g,
53                              int & nnz_h_lag)
54{
55  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
56  ASL_pfgh* asl = info->asl_;
57 
58  n = n_var; // # of variables
59  m = n_con; // # of constraints
60  nnz_jac_g = nzc; // # of non-zeros in the jacobian
61  nnz_h_lag = info->nz_h_full_; // # of non-zeros in the hessian
62 
63  return true;
64}
65
66static bool get_bounds_info(void * amplInfo,int  n, double * x_l,
67                            double * x_u, int  m, double * g_l, double * g_u)
68{
69  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
70  ASL_pfgh* asl = info->asl_;
71  assert(n == n_var);
72  assert(m == n_con);
73  int i;
74  for (i=0; i<n; i++) {
75    x_l[i] = LUv[2*i];
76    x_u[i] = LUv[2*i+1];
77  }
78 
79  for (i=0; i<m; i++) {
80    g_l[i] = LUrhs[2*i];
81    g_u[i] = LUrhs[2*i+1];
82  }
83  return true;
84}
85
86#endif
87bool get_constraints_linearity(void * amplInfo,int  n,
88      int * const_types)
89{
90  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
91  ASL_pfgh* asl = info->asl_;
92  //check that n is good
93  assert(n == n_con);
94  // check that there are no network constraints
95  assert(nlnc == 0 && lnc == 0);
96  //the first nlc constraints are non linear the rest is linear
97  int i;
98  for (i=0; i<nlc; i++) {
99    const_types[i]=1;
100  }
101  // the rest is linear
102  for (i=nlc; i<n_con; i++)
103    const_types[i]=0;
104  return true;
105}
106#if 0
107bool get_starting_point(int  n, bool init_x, double * x, bool init_z,
108                        double * z_L, double * z_U, int  m, bool init_lambda, double * lambda)
109{
110  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
111  ASL_pfgh* asl = info->asl_;
112  assert(n == n_var);
113  assert(m == n_con);
114  int i;
115 
116  if (init_x) {
117    for (i=0; i<n; i++) {
118      if (havex0[i]) {
119        x[i] = X0[i];
120      }
121      else {
122        x[i] = 0.0;
123      }
124    }
125  }
126 
127  if (init_z) {
128    for (i=0; i<n; i++) {
129      z_L[i] = z_U[i] = 1.0;
130    }
131  }
132 
133  if (init_lambda) {
134    for (i=0; i<m; i++) {
135      if (havepi0[i]) {
136        lambda[i] = pi0[i];
137      }
138      else {
139        lambda[i] = 0.0;
140      }
141    }
142  }
143 
144  return true;
145}
146#endif
147static bool internal_objval(CbcAmplInfo * info ,double & obj_val)
148{
149  ASL_pfgh* asl = info->asl_;
150  info->objval_called_with_current_x_ = false; // in case the call below fails
151
152  if (n_obj==0) {
153    obj_val = 0;
154    info->objval_called_with_current_x_ = true;
155    return true;
156  }  else {
157    double  retval = objval(0, info->non_const_x_, (fint*)info->nerror_);
158    if (!info->nerror_) {
159      obj_val = info->obj_sign_*retval;
160      info->objval_called_with_current_x_ = true;
161      return true;
162    } else {
163      abort();
164    }
165  }
166 
167  return false;
168}
169static bool internal_conval(CbcAmplInfo * info ,double * g)
170{
171  ASL_pfgh* asl = info->asl_;
172  info->conval_called_with_current_x_ = false; // in case the call below fails
173  assert (g);
174 
175  conval(info->non_const_x_, g, (fint*)info->nerror_);
176
177  if (!info->nerror_) {
178    info->conval_called_with_current_x_ = true;
179    return true;
180  } else {
181    abort();
182  }
183  return false;
184}
185
186static bool apply_new_x(CbcAmplInfo * info  ,bool new_x, int  n, const double * x)
187{
188  ASL_pfgh* asl = info->asl_;
189 
190  if (new_x) {
191    // update the flags so these methods are called
192    // before evaluating the hessian
193    info->conval_called_with_current_x_ = false;
194    info->objval_called_with_current_x_ = false;
195    info->jacval_called_with_current_x_ = false;
196
197    //copy the data to the non_const_x_
198    if (!info->non_const_x_) {
199      info->non_const_x_ = new double [n];
200    }
201
202    for (int  i=0; i<n; i++) {
203      info->non_const_x_[i] = x[i];
204    }
205   
206    // tell ampl that we have a new x
207    xknowne(info->non_const_x_, (fint*)info->nerror_);
208    return info->nerror_ ? false : true;
209  }
210 
211  return true;
212}
213static bool eval_f(void * amplInfo,int  n, const double * x, bool new_x, double & obj_value)
214{
215  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
216  if (!apply_new_x(info,new_x, n, x)) {
217    return false;
218  }
219 
220  return internal_objval(info,obj_value);
221}
222
223static bool eval_grad_f(void * amplInfo,int  n, const double * x, bool new_x, double * grad_f)
224{
225  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
226  ASL_pfgh* asl = info->asl_;
227  if (!apply_new_x(info,new_x, n, x)) {
228    return false;
229  }
230  int i;
231 
232  if (n_obj==0) {
233    for (i=0; i<n; i++) {
234      grad_f[i] = 0.;
235    }
236  }
237  else {
238    objgrd(0, info->non_const_x_, grad_f, (fint*)info->nerror_);
239    if (info->nerror_) {
240      return false;
241    }
242   
243    if (info->obj_sign_==-1) {
244      for (i=0; i<n; i++) {
245        grad_f[i] = -grad_f[i];
246      }
247    }
248  }
249  return true;
250}
251static bool eval_g(void * amplInfo,int  n, const double * x, bool new_x, double * g)
252{
253  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
254  ASL_pfgh* asl = info->asl_;
255  assert(n == n_var);
256 
257  if (!apply_new_x(info,new_x, n, x)) {
258    return false;
259  }
260 
261  return internal_conval(info, g);
262}
263
264static bool eval_jac_g(void * amplInfo,int  n, const double * x, bool new_x,
265                       double * values)
266{
267  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
268  ASL_pfgh* asl = info->asl_;
269  assert(n == n_var);
270 
271  assert (values);
272  if (!apply_new_x(info,new_x, n, x)) {
273    return false;
274  }
275   
276  jacval(info->non_const_x_, values, (fint*)info->nerror_);
277  if (!info->nerror_) {
278    return true;
279  } else {
280    abort();
281  }
282  return false;
283}
284#if 0
285
286static bool eval_h(void * amplInfo,int  n, const double * x, bool new_x,
287            double  obj_factor, int  m, const double * lambda,
288            bool new_lambda, int  nele_hess, int * iRow,
289            int * jCol, double * values)
290{
291  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
292  ASL_pfgh* asl = info->asl_;
293  assert(n == n_var);
294  assert(m == n_con);
295  int i;
296 
297  if (iRow && jCol && !values) {
298    // setup the structure
299    int k=0;
300    for (int i=0; i<n; i++) {
301      for (int j=sputinfo->hcolstarts[i]; j<sputinfo->hcolstarts[i+1]; j++) {
302        //iRow[k] = i + 1;
303        iRow[k] = i;
304        jCol[k] = sputinfo->hrownos[j]+1;
305        k++;
306      }
307    }
308    assert(k==nele_hess);
309    return true;
310  }
311  else if (!iRow & !jCol && values) {
312    if (!apply_new_x(info,new_x, n, x)) {
313      return false;
314    }
315    if (!info->objval_called_with_current_x_) {
316      double  dummy;
317      internal_objval(info,dummy);
318      internal_conval(info,m,NULL);
319    }
320    if (!info->conval_called_with_current_x_) {
321      internal_conval(info,m,NULL);
322    }
323    // copy lambda to non_const_lambda - note, we do not store a copy like
324    // we do with x since lambda is only used here and not in other calls
325    double * non_const_lambda = new double [m];
326    for (i=0; i<m; i++) {
327      non_const_lambda[i] = lambda[i];
328    }
329   
330    real ow=info->obj_sign_*obj_factor;
331    sphes(values, -1, &ow, non_const_lambda);
332   
333    delete [] non_const_lambda;
334    return true;
335  }
336  else {
337    assert(false && "Invalid combination of iRow, jCol, and values pointers");
338  }
339 
340  return false;
341}
342#endif
343//-------------------------------------------------------------------
344// Useful Constructor
345//-------------------------------------------------------------------
346ClpAmplObjective::ClpAmplObjective (void * amplInfo)
347  : ClpObjective()
348{
349  type_=12;
350  activated_=1;
351  gradient_ = NULL;
352  objective_ = NULL;
353  offset_ = 0.0;
354  amplObjective_ = amplInfo;
355}
356
357//-------------------------------------------------------------------
358// Copy constructor
359//-------------------------------------------------------------------
360ClpAmplObjective::ClpAmplObjective (const ClpAmplObjective & rhs) 
361: ClpObjective(rhs)
362{ 
363  amplObjective_ = rhs.amplObjective_;
364  offset_ = rhs.offset_;
365  type_=rhs.type_;
366  if (!amplObjective_) {
367    objective_=NULL;
368    gradient_ = NULL;
369  } else {
370    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
371    ASL_pfgh* asl = info->asl_;
372   
373    int numberColumns = n_var;;
374    if (rhs.objective_) {
375      objective_ = new double [numberColumns];
376      memcpy(objective_,rhs.objective_,numberColumns*sizeof(double));
377    } else {
378      objective_=NULL;
379    }
380    if (rhs.gradient_) {
381      gradient_ = new double [numberColumns];
382      memcpy(gradient_,rhs.gradient_,numberColumns*sizeof(double));
383    } else {
384      gradient_=NULL;
385    }
386  }
387}
388 
389
390//-------------------------------------------------------------------
391// Destructor
392//-------------------------------------------------------------------
393ClpAmplObjective::~ClpAmplObjective ()
394{
395  delete [] objective_;
396  delete [] gradient_;
397}
398
399//----------------------------------------------------------------
400// Assignment operator
401//-------------------------------------------------------------------
402ClpAmplObjective &
403ClpAmplObjective::operator=(const ClpAmplObjective& rhs)
404{
405  if (this != &rhs) {
406    delete [] objective_;
407    delete [] gradient_;
408    amplObjective_ = rhs.amplObjective_;
409    offset_ = rhs.offset_;
410    type_=rhs.type_;
411    if (!amplObjective_) {
412      objective_=NULL;
413      gradient_ = NULL;
414    } else {
415      CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
416      ASL_pfgh* asl = info->asl_;
417     
418      int numberColumns = n_var;;
419      if (rhs.objective_) {
420        objective_ = new double [numberColumns];
421        memcpy(objective_,rhs.objective_,numberColumns*sizeof(double));
422      } else {
423        objective_=NULL;
424      }
425      if (rhs.gradient_) {
426        gradient_ = new double [numberColumns];
427        memcpy(gradient_,rhs.gradient_,numberColumns*sizeof(double));
428      } else {
429        gradient_=NULL;
430      }
431    }
432  }
433  return *this;
434}
435
436// Returns gradient
437double * 
438ClpAmplObjective::gradient(const ClpSimplex * model,
439                                const double * solution, double & offset,bool refresh,
440                                int includeLinear)
441{
442  if (model)
443    assert (model->optimizationDirection()==1.0);
444  bool scaling=false;
445  if (model&&(model->rowScale()||
446              model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0))
447    scaling=true;
448  const double * cost = NULL;
449  if (model)
450    cost = model->costRegion();
451  if (!cost) {
452    // not in solve
453    cost = objective_;
454    scaling=false;
455  }
456  assert (!scaling);
457  if (!amplObjective_||!solution||!activated_) {
458    offset=offset_;
459    return objective_;
460  } else {
461    if (refresh||!gradient_) {
462      CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
463      ASL_pfgh* asl = info->asl_;
464      int numberColumns = n_var;;
465     
466      if (!gradient_) 
467        gradient_ = new double[numberColumns];
468      assert (solution);
469      eval_grad_f(amplObjective_,numberColumns,solution,true,gradient_);
470      // Is this best way?
471      double objValue=0.0;
472      eval_f(amplObjective_,numberColumns,solution,false,objValue);
473      double objValue2=0.0;
474      for (int i=0;i<numberColumns;i++)
475        objValue2 += gradient_[i]*solution[i];
476      offset_ = objValue2 - objValue; // or other way???
477      if (model&&model->optimizationDirection()!=1.0) {
478        offset *= model->optimizationDirection();
479        for (int i=0;i<numberColumns;i++)
480          gradient_[i] *= -1.0;
481      }
482    }
483    offset=offset_;
484    return gradient_;
485  }
486}
487 
488//-------------------------------------------------------------------
489// Clone
490//-------------------------------------------------------------------
491ClpObjective * ClpAmplObjective::clone() const
492{
493  return new ClpAmplObjective(*this);
494}
495// Resize objective
496void 
497ClpAmplObjective::resize(int newNumberColumns)
498{
499  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
500  ASL_pfgh* asl = info->asl_;
501  int numberColumns = n_var;;
502  if (numberColumns!=newNumberColumns) {
503    abort();
504  } 
505 
506}
507// Delete columns in  objective
508void 
509ClpAmplObjective::deleteSome(int numberToDelete, const int * which) 
510{
511  if (numberToDelete)
512    abort();
513}
514/* Returns reduced gradient.Returns an offset (to be added to current one).
515 */
516double 
517ClpAmplObjective::reducedGradient(ClpSimplex * model, double * region,
518                                       bool useFeasibleCosts)
519{
520  int numberRows = model->numberRows();
521  int numberColumns=model->numberColumns();
522 
523  //work space
524  CoinIndexedVector  * workSpace = model->rowArray(0);
525 
526  CoinIndexedVector arrayVector;
527  arrayVector.reserve(numberRows+1);
528 
529  int iRow;
530#ifdef CLP_DEBUG
531  workSpace->checkClear();
532#endif
533  double * array = arrayVector.denseVector();
534  int * index = arrayVector.getIndices();
535  int number=0;
536  const double * costNow = gradient(model,model->solutionRegion(),offset_,
537                                    true,useFeasibleCosts ? 2 : 1);
538  double * cost = model->costRegion();
539  const int * pivotVariable = model->pivotVariable();
540  for (iRow=0;iRow<numberRows;iRow++) {
541    int iPivot=pivotVariable[iRow];
542    double value;
543    if (iPivot<numberColumns)
544      value = costNow[iPivot];
545    else if (!useFeasibleCosts) 
546      value = cost[iPivot];
547    else 
548      value=0.0;
549    if (value) {
550      array[iRow]=value;
551      index[number++]=iRow;
552    }
553  }
554  arrayVector.setNumElements(number);
555
556  // Btran basic costs
557  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
558  double * work = workSpace->denseVector();
559  ClpFillN(work,numberRows,0.0);
560  // now look at dual solution
561  double * rowReducedCost = region+numberColumns;
562  double * dual = rowReducedCost;
563  const double * rowCost = cost+numberColumns;
564  for (iRow=0;iRow<numberRows;iRow++) {
565    dual[iRow]=array[iRow];
566  }
567  double * dj = region;
568  ClpDisjointCopyN(costNow,numberColumns,dj);
569 
570  model->transposeTimes(-1.0,dual,dj);
571  for (iRow=0;iRow<numberRows;iRow++) {
572    // slack
573    double value = dual[iRow];
574    value += rowCost[iRow];
575    rowReducedCost[iRow]=value;
576  }
577  return offset_;
578}
579/* Returns step length which gives minimum of objective for
580   solution + theta * change vector up to maximum theta.
581   
582   arrays are numberColumns+numberRows
583*/
584double 
585ClpAmplObjective::stepLength(ClpSimplex * model,
586                                  const double * solution,
587                                  const double * change,
588                                  double maximumTheta,
589                                  double & currentObj,
590                                  double & predictedObj,
591                                  double & thetaObj)
592{
593  // Assume convex
594  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
595  ASL_pfgh* asl = info->asl_;
596 
597  int numberColumns = n_var;;
598  double * tempSolution = new double [numberColumns];
599  double * tempGradient = new double [numberColumns];
600  // current
601  eval_f(amplObjective_,numberColumns,solution,true,currentObj);
602  double objA = currentObj;
603  double thetaA = 0.0;
604  // at maximum
605  int i;
606  for (i=0;i<numberColumns;i++)
607    tempSolution[i] = solution[i] + maximumTheta*change[i];
608  eval_f(amplObjective_,numberColumns,tempSolution,true,thetaObj);
609  double objC = thetaObj;
610  double thetaC = maximumTheta;
611  double objB=0.5*(objA+objC);
612  double thetaB=0.5*maximumTheta;
613  double gradientNorm=1.0e6;
614  while (gradientNorm>1.0e-6&&thetaC-thetaA>1.0e-8) {
615    for (i=0;i<numberColumns;i++)
616      tempSolution[i] = solution[i] + thetaB*change[i];
617    eval_grad_f(amplObjective_,numberColumns,tempSolution,true,tempGradient);
618    eval_f(amplObjective_,numberColumns,tempSolution,false,objB);
619    double changeObj=0.0;
620    gradientNorm=0.0;
621    for (i=0;i<numberColumns;i++) {
622      changeObj += tempGradient[i]*change[i];
623      gradientNorm += tempGradient[i]*tempGradient[i];
624    }
625    gradientNorm = fabs(changeObj)/sqrt(gradientNorm);
626    // Should try and get quadratic convergence by interpolation
627    if (changeObj<0.0) {
628      // increasing is good
629      thetaA = thetaB;
630    } else {
631      // decreasing is good
632      thetaC = thetaB;
633    }
634    thetaB = 0.5*(thetaA+thetaC);
635  }
636  delete [] tempSolution;
637  delete [] tempGradient;
638  predictedObj = objB;
639  return thetaB;
640}
641// Return objective value (without any ClpModel offset) (model may be NULL)
642double 
643ClpAmplObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
644{
645  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
646  ASL_pfgh* asl = info->asl_;
647 
648  int numberColumns = n_var;;
649  // current
650  double currentObj=0.0;
651  eval_f(amplObjective_,numberColumns,solution,true,currentObj);
652  return currentObj;
653}
654// Scale objective
655void 
656ClpAmplObjective::reallyScale(const double * columnScale) 
657{
658  abort();
659}
660/* Given a zeroed array sets nonlinear columns to 1.
661   Returns number of nonlinear columns
662*/
663int 
664ClpAmplObjective::markNonlinear(char * which)
665{
666  int iColumn;
667  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
668  ASL_pfgh* asl = info->asl_;
669  int nonLinear = CoinMax(nlvc,nlvo);
670  for (iColumn=0;iColumn<nonLinear;iColumn++) {
671    which[iColumn]=1;
672  }
673  int numberNonLinearColumns = 0;
674  int numberColumns = n_var;;
675  for (iColumn=0;iColumn<numberColumns;iColumn++) {
676    if(which[iColumn])
677      numberNonLinearColumns++;
678  }
679  return numberNonLinearColumns;
680}
681// Say we have new primal solution - so may need to recompute
682void 
683ClpAmplObjective::newXValues() 
684{
685  CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
686  info->conval_called_with_current_x_ = false;
687  info->objval_called_with_current_x_ = false;
688  info->jacval_called_with_current_x_ = false;
689}
690//#############################################################################
691// Constructors / Destructor / Assignment
692//#############################################################################
693//-------------------------------------------------------------------
694// Default Constructor
695//-------------------------------------------------------------------
696ClpConstraintAmpl::ClpConstraintAmpl () 
697: ClpConstraint()
698{
699  type_=3;
700  column_=NULL;
701  coefficient_ = NULL;
702  numberCoefficients_=0;
703  amplInfo_ = NULL;
704}
705
706//-------------------------------------------------------------------
707// Useful Constructor
708//-------------------------------------------------------------------
709ClpConstraintAmpl::ClpConstraintAmpl (int row, void * amplInfo)
710  : ClpConstraint()
711{
712  type_=3;
713  rowNumber_=row;
714  amplInfo_ = amplInfo;
715  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
716  ASL_pfgh* asl = info->asl_;
717  assert (rowNumber_<nlc);
718  numberCoefficients_ = info->rowStart_[rowNumber_+1]-info->rowStart_[rowNumber_];
719  column_ = CoinCopyOfArray(info->column_+info->rowStart_[rowNumber_],numberCoefficients_);
720  coefficient_ = new double [numberCoefficients_];;
721}
722
723//-------------------------------------------------------------------
724// Copy constructor
725//-------------------------------------------------------------------
726ClpConstraintAmpl::ClpConstraintAmpl (const ClpConstraintAmpl & rhs)
727: ClpConstraint(rhs)
728{ 
729  numberCoefficients_=rhs.numberCoefficients_;
730  column_ = CoinCopyOfArray(rhs.column_,numberCoefficients_);
731  coefficient_ = CoinCopyOfArray(rhs.coefficient_,numberCoefficients_);
732}
733 
734
735//-------------------------------------------------------------------
736// Destructor
737//-------------------------------------------------------------------
738ClpConstraintAmpl::~ClpConstraintAmpl ()
739{
740  delete [] column_;
741  delete [] coefficient_;
742}
743
744//----------------------------------------------------------------
745// Assignment operator
746//-------------------------------------------------------------------
747ClpConstraintAmpl &
748ClpConstraintAmpl::operator=(const ClpConstraintAmpl& rhs)
749{
750  if (this != &rhs) {
751    delete [] column_;
752    delete [] coefficient_;
753    numberCoefficients_=rhs.numberCoefficients_;
754    column_ = CoinCopyOfArray(rhs.column_,numberCoefficients_);
755    coefficient_ = CoinCopyOfArray(rhs.coefficient_,numberCoefficients_);
756  }
757  return *this;
758}
759//-------------------------------------------------------------------
760// Clone
761//-------------------------------------------------------------------
762ClpConstraint * ClpConstraintAmpl::clone() const
763{
764  return new ClpConstraintAmpl(*this);
765}
766
767// Returns gradient
768int 
769ClpConstraintAmpl::gradient(const ClpSimplex * model,
770                              const double * solution, 
771                              double * gradient,
772                              double & functionValue, 
773                              double & offset,
774                              bool useScaling,
775                              bool refresh) const
776{
777  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
778  ASL_pfgh* asl = info->asl_;
779  int numberColumns = n_var;;
780  // If not done then do all
781  if (!info->jacval_called_with_current_x_) {
782    bool getStuff = eval_g(amplInfo_,numberColumns,solution,true,info->constraintValues_);
783    assert (getStuff);
784    getStuff = eval_jac_g(amplInfo_,numberColumns,solution,false,info->gradient_);
785    assert (getStuff);
786    info->jacval_called_with_current_x_=getStuff;
787  }
788  if (refresh||!lastGradient_) {
789    functionValue_ = info->constraintValues_[rowNumber_]; 
790    offset_=functionValue_; // sign??
791    if (!lastGradient_)
792      lastGradient_ = new double[numberColumns];
793    CoinZeroN(lastGradient_,numberColumns);
794    bool scaling=(model&&model->rowScale()&&useScaling);
795    assert (!scaling);
796    int i;
797    int start = info->rowStart_[rowNumber_];
798    assert (numberCoefficients_==info->rowStart_[rowNumber_+1]-start);
799    for (i=0;i<numberCoefficients_;i++) {
800      int iColumn= column_[i];
801      double valueS = solution[iColumn];
802      double valueG = info->gradient_[start+i];
803      lastGradient_[iColumn]=valueG;
804      offset_ -= valueS*valueG;
805    }
806  }
807  functionValue = functionValue_;
808  offset = offset_;
809  memcpy(gradient,lastGradient_,numberColumns*sizeof(double));
810  return 0;
811}
812// Resize constraint
813void 
814ClpConstraintAmpl::resize(int newNumberColumns)
815{
816  abort();
817}
818// Delete columns in  constraint
819void 
820ClpConstraintAmpl::deleteSome(int numberToDelete, const int * which) 
821{
822  if (numberToDelete) {
823    abort();
824  }
825}
826// Scale constraint
827void 
828ClpConstraintAmpl::reallyScale(const double * columnScale) 
829{
830  abort();
831}
832/* Given a zeroed array sets nonlinear columns to 1.
833   Returns number of nonlinear columns
834*/
835int 
836ClpConstraintAmpl::markNonlinear(char * which) const
837{
838  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
839  ASL_pfgh* asl = info->asl_;
840  int iColumn;
841  int numberNon=0;
842  int nonLinear = CoinMax(nlvc,nlvo);
843  for (iColumn=0;iColumn<numberCoefficients_;iColumn++) {
844    int jColumn = column_[iColumn];
845    if (jColumn<nonLinear) {
846      which[jColumn]=1;
847      numberNon++;
848    }
849  }
850  return numberNon;
851}
852/* Given a zeroed array sets possible nonzero coefficients to 1.
853   Returns number of nonzeros
854*/
855int 
856ClpConstraintAmpl::markNonzero(char * which) const
857{
858  int iColumn;
859  for (iColumn=0;iColumn<numberCoefficients_;iColumn++) {
860    which[column_[iColumn]]=1;
861  }
862  return numberCoefficients_;
863}
864// Number of coefficients
865int 
866ClpConstraintAmpl::numberCoefficients() const
867{
868  return numberCoefficients_;
869}
870// Say we have new primal solution - so may need to recompute
871void 
872ClpConstraintAmpl::newXValues() 
873{
874  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
875  info->conval_called_with_current_x_ = false;
876  info->objval_called_with_current_x_ = false;
877  info->jacval_called_with_current_x_ = false;
878}
879/* Load nonlinear part of problem from AMPL info
880   Returns 0 if linear
881   1 if quadratic objective
882   2 if quadratic constraints
883   3 if nonlinear objective
884   4 if nonlinear constraints
885   -1 on failure
886*/
887int 
888ClpSimplex::loadNonLinear(void * amplInfo, int & numberConstraints, 
889                          ClpConstraint ** & constraints)
890{
891  numberConstraints=0;
892  constraints=NULL;
893  CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
894  ASL_pfgh* asl = info->asl_;
895  // For moment don't say quadratic
896  int type=0;
897  if (nlo+nlc) {
898    // nonlinear
899    if (!nlc) {
900      type=3;
901      delete objective_;
902      objective_= new ClpAmplObjective(amplInfo);
903    } else {
904      type=4;
905      numberConstraints = nlc;
906      constraints = new ClpConstraint * [numberConstraints];
907      if (nlo) {
908        delete objective_;
909        objective_= new ClpAmplObjective(amplInfo);
910      }
911      for (int i=0;i<numberConstraints;i++) {
912        constraints[i] = new ClpConstraintAmpl(i,amplInfo);
913      }
914    }
915  }
916  return type;
917}
918#else
919#include "ClpSimplex.hpp"
920#include "ClpConstraint.hpp"
921int 
922ClpSimplex::loadNonLinear(void * amplInfo, int & numberConstraints, 
923                          ClpConstraint ** & constraints)
924{
925  abort();
926  return 0;
927}
928#endif
Note: See TracBrowser for help on using the repository browser.