source: branches/devel/Cbc/src/ClpAmplStuff.cpp @ 648

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

update branches/devel for threads

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