source: stable/2.4/Cbc/src/ClpAmplObjective.cpp @ 1271

Last change on this file since 1271 was 1271, checked in by forrest, 10 years ago

Creating new stable branch 2.4 from trunk (rev 1270)

File size: 23.1 KB
Line 
1/* $Id: ClpAmplObjective.cpp 1173 2009-06-04 09:44:10Z forrest $ */
2// Copyright (C) 2007, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinPragma.hpp"
6#include "CoinHelperFunctions.hpp"
7#include "CoinIndexedVector.hpp"
8#include "ClpFactorization.hpp"
9#include "ClpSimplex.hpp"
10#include "ClpAmplObjective.hpp"
11//#############################################################################
12// Constructors / Destructor / Assignment
13//#############################################################################
14
15//-------------------------------------------------------------------
16// Default Constructor
17//-------------------------------------------------------------------
18ClpAmplObjective::ClpAmplObjective () 
19: ClpObjective()
20{
21  type_=12;
22  objective_=NULL;
23  amplObjective_=NULL;
24  gradient_ = NULL;
25}
26
27//-------------------------------------------------------------------
28// Useful Constructor
29//-------------------------------------------------------------------
30ClpAmplObjective::ClpAmplObjective (void * amplInfo)
31  : ClpObjective()
32{
33  type_=12;
34  loadAmplObjective(amplInfo);
35  numberColumns_ = numberColumns;
36  if (numberExtendedColumns>=0)
37    numberExtendedColumns_= CoinMax(numberColumns_,numberExtendedColumns);
38  else
39    numberExtendedColumns_= numberColumns_;
40  if (objective) {
41    objective_ = new double [numberExtendedColumns_];
42    memcpy(objective_,objective,numberColumns_*sizeof(double));
43    memset(objective_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
44  } else {
45    objective_ = new double [numberExtendedColumns_];
46    memset(objective_,0,numberExtendedColumns_*sizeof(double));
47  }
48  if (start) 
49    amplObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
50                                             start[numberColumns],element,column,start,NULL);
51  else
52  amplObjective_=NULL;
53  gradient_ = NULL;
54  activated_=1;
55  fullMatrix_=false;
56}
57
58//-------------------------------------------------------------------
59// Copy constructor
60//-------------------------------------------------------------------
61ClpAmplObjective::ClpAmplObjective (const ClpAmplObjective & rhs) 
62: ClpObjective(rhs)
63{ 
64  numberColumns_=rhs.numberColumns_;
65  numberExtendedColumns_=rhs.numberExtendedColumns_;
66  fullMatrix_=rhs.fullMatrix_;
67  if (rhs.objective_) {
68    objective_ = new double [numberExtendedColumns_];
69    memcpy(objective_,rhs.objective_,numberExtendedColumns_*sizeof(double));
70  } else {
71    objective_=NULL;
72  }
73  if (rhs.gradient_) {
74    gradient_ = new double [numberExtendedColumns_];
75    memcpy(gradient_,rhs.gradient_,numberExtendedColumns_*sizeof(double));
76  } else {
77    gradient_=NULL;
78  }
79  if (rhs.amplObjective_) {
80    // see what type of matrix wanted
81    if (type==0) {
82      // just copy
83      amplObjective_ = new CoinPackedMatrix(*rhs.amplObjective_);
84    } else if (type==1) {
85      // expand to full symmetric
86      fullMatrix_=true;
87      const int * columnAmpl1 = rhs.amplObjective_->getIndices();
88      const CoinBigIndex * columnAmplStart1 = rhs.amplObjective_->getVectorStarts();
89      const int * columnAmplLength1 = rhs.amplObjective_->getVectorLengths();
90      const double * amplElement1 = rhs.amplObjective_->getElements();
91      CoinBigIndex * columnAmplStart2 = new CoinBigIndex [numberExtendedColumns_+1];
92      int * columnAmplLength2 = new int [numberExtendedColumns_];
93      int iColumn;
94      int numberColumns = rhs.amplObjective_->getNumCols();
95      int numberBelow=0;
96      int numberAbove=0;
97      int numberDiagonal=0;
98      CoinZeroN(columnAmplLength2,numberExtendedColumns_);
99      for (iColumn=0;iColumn<numberColumns;iColumn++) {
100        for (CoinBigIndex j=columnAmplStart1[iColumn];
101             j<columnAmplStart1[iColumn]+columnAmplLength1[iColumn];j++) {
102          int jColumn = columnAmpl1[j];
103          if (jColumn>iColumn) {
104            numberBelow++;
105            columnAmplLength2[jColumn]++;
106            columnAmplLength2[iColumn]++;
107          } else if (jColumn==iColumn) {
108            numberDiagonal++;
109            columnAmplLength2[iColumn]++;
110          } else {
111            numberAbove++;
112          }
113        }
114      }
115      if (numberAbove>0) {
116        if (numberAbove==numberBelow) {
117          // already done
118          amplObjective_ = new CoinPackedMatrix(*rhs.amplObjective_);
119          delete [] columnAmplStart2;
120          delete [] columnAmplLength2;
121        } else {
122          printf("number above = %d, number below = %d, error\n",
123                 numberAbove,numberBelow);
124          abort();
125        }
126      } else {
127        int numberElements=numberDiagonal+2*numberBelow;
128        int * columnAmpl2 = new int [numberElements];
129        double * amplElement2 = new double [numberElements];
130        columnAmplStart2[0]=0;
131        numberElements=0;
132        for (iColumn=0;iColumn<numberColumns;iColumn++) {
133          int n=columnAmplLength2[iColumn];
134          columnAmplLength2[iColumn]=0;
135          numberElements += n;
136          columnAmplStart2[iColumn+1]=numberElements;
137        }
138        for (iColumn=0;iColumn<numberColumns;iColumn++) {
139          for (CoinBigIndex j=columnAmplStart1[iColumn];
140               j<columnAmplStart1[iColumn]+columnAmplLength1[iColumn];j++) {
141            int jColumn = columnAmpl1[j];
142            if (jColumn>iColumn) {
143              // put in two places
144              CoinBigIndex put=columnAmplLength2[jColumn]+columnAmplStart2[jColumn];
145              columnAmplLength2[jColumn]++;
146              amplElement2[put]=amplElement1[j];
147              columnAmpl2[put]=iColumn;
148              put=columnAmplLength2[iColumn]+columnAmplStart2[iColumn];
149              columnAmplLength2[iColumn]++;
150              amplElement2[put]=amplElement1[j];
151              columnAmpl2[put]=jColumn;
152            } else if (jColumn==iColumn) {
153              CoinBigIndex put=columnAmplLength2[iColumn]+columnAmplStart2[iColumn];
154              columnAmplLength2[iColumn]++;
155              amplElement2[put]=amplElement1[j];
156              columnAmpl2[put]=iColumn;
157            } else {
158              abort();
159            }
160          }
161        }
162        // Now create
163        amplObjective_ = 
164          new CoinPackedMatrix (true,
165                                rhs.numberExtendedColumns_,
166                                rhs.numberExtendedColumns_,
167                                numberElements,
168                                amplElement2,
169                                columnAmpl2,
170                                columnAmplStart2,
171                                columnAmplLength2,0.0,0.0);
172        delete [] columnAmplStart2;
173        delete [] columnAmplLength2;
174        delete [] columnAmpl2;
175        delete [] amplElement2;
176      }
177    } else {
178      fullMatrix_=false;
179      abort(); // code when needed
180    }
181           
182  } else {
183    amplObjective_=NULL;
184  }
185}
186 
187
188//-------------------------------------------------------------------
189// Destructor
190//-------------------------------------------------------------------
191ClpAmplObjective::~ClpAmplObjective ()
192{
193  delete [] objective_;
194  delete [] gradient_;
195}
196
197//----------------------------------------------------------------
198// Assignment operator
199//-------------------------------------------------------------------
200ClpAmplObjective &
201ClpAmplObjective::operator=(const ClpAmplObjective& rhs)
202{
203  if (this != &rhs) {
204    amplObjective_ = NULL;
205    ClpObjective::operator=(rhs);
206    numberColumns_=rhs.numberColumns_;
207    numberExtendedColumns_=rhs.numberExtendedColumns_;
208    if (rhs.objective_) {
209      objective_ = new double [numberExtendedColumns_];
210      memcpy(objective_,rhs.objective_,numberExtendedColumns_*sizeof(double));
211    } else {
212      objective_=NULL;
213    }
214    if (rhs.gradient_) {
215      gradient_ = new double [numberExtendedColumns_];
216      memcpy(gradient_,rhs.gradient_,numberExtendedColumns_*sizeof(double));
217    } else {
218      gradient_=NULL;
219    }
220    if (rhs.amplObjective_) {
221      amplObjective_ = new CoinPackedMatrix(*rhs.amplObjective_);
222    } else {
223      amplObjective_=NULL;
224    }
225  }
226  return *this;
227}
228
229// Returns gradient
230double * 
231ClpAmplObjective::gradient(const ClpSimplex * model,
232                                const double * solution, double & offset,bool refresh,
233                                int includeLinear)
234{
235  offset=0.0;
236  bool scaling=false;
237  if (model&&(model->rowScale()||
238              model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0))
239    scaling=true;
240  const double * cost = NULL;
241  if (model)
242    cost = model->costRegion();
243  if (!cost) {
244    // not in solve
245    cost = objective_;
246    scaling=false;
247  }
248  if (!scaling) {
249    if (!amplObjective_||!solution||!activated_) {
250      return objective_;
251    } else {
252      if (refresh||!gradient_) {
253        if (!gradient_) 
254          gradient_ = new double[numberExtendedColumns_];
255        const int * columnAmpl = amplObjective_->getIndices();
256        const CoinBigIndex * columnAmplStart = amplObjective_->getVectorStarts();
257        const int * columnAmplLength = amplObjective_->getVectorLengths();
258        const double * amplElement = amplObjective_->getElements();
259        offset=0.0;
260        // use current linear cost region
261        if (includeLinear==1)
262          memcpy(gradient_,cost,numberExtendedColumns_*sizeof(double));
263        else if (includeLinear==2)
264          memcpy(gradient_,objective_,numberExtendedColumns_*sizeof(double));
265        else
266          memset(gradient_,0,numberExtendedColumns_*sizeof(double));
267        if (activated_) {
268          if (!fullMatrix_) {
269            int iColumn;
270            for (iColumn=0;iColumn<numberColumns_;iColumn++) {
271              double valueI = solution[iColumn];
272              CoinBigIndex j;
273              for (j=columnAmplStart[iColumn];
274                   j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
275                int jColumn = columnAmpl[j];
276                double valueJ = solution[jColumn];
277                double elementValue = amplElement[j];
278                if (iColumn!=jColumn) {
279                  offset += valueI*valueJ*elementValue;
280                  //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
281                  //printf("%d %d %g %g %g -> %g\n",
282                  //       iColumn,jColumn,valueI,valueJ,elementValue,
283                  //       valueI*valueJ*elementValue);
284                  double gradientI = valueJ*elementValue;
285                  double gradientJ = valueI*elementValue;
286                  gradient_[iColumn] += gradientI;
287                  gradient_[jColumn] += gradientJ;
288                } else {
289                  offset += 0.5*valueI*valueI*elementValue;
290                  //if (fabs(valueI*valueI*elementValue)>1.0e-12)
291                  //printf("XX %d %g %g -> %g\n",
292                  //       iColumn,valueI,elementValue,
293                  //       0.5*valueI*valueI*elementValue);
294                  double gradientI = valueI*elementValue;
295                  gradient_[iColumn] += gradientI;
296                }
297              }
298            }
299          } else {
300            // full matrix
301            int iColumn;
302            offset *= 2.0;
303            for (iColumn=0;iColumn<numberColumns_;iColumn++) {
304              CoinBigIndex j;
305              double value=0.0;
306              double current = gradient_[iColumn];
307              for (j=columnAmplStart[iColumn];
308                   j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
309                int jColumn = columnAmpl[j];
310                double valueJ = solution[jColumn]*amplElement[j];
311                value += valueJ;
312              }
313              offset += value*solution[iColumn];
314              gradient_[iColumn] = current+value;
315            }
316            offset *= 0.5;
317          }
318        }
319      }
320      if (model)
321        offset *= model->optimizationDirection()*model->objectiveScale();
322      return gradient_;
323    }
324  } else {
325    // do scaling
326    assert (solution);
327    // for now only if half
328    assert (!fullMatrix_);
329    if (refresh||!gradient_) {
330      if (!gradient_) 
331        gradient_ = new double[numberExtendedColumns_];
332      double direction = model->optimizationDirection()*model->objectiveScale();
333      // direction is actually scale out not scale in
334      //if (direction)
335      //direction = 1.0/direction;
336      const int * columnAmpl = amplObjective_->getIndices();
337      const CoinBigIndex * columnAmplStart = amplObjective_->getVectorStarts();
338      const int * columnAmplLength = amplObjective_->getVectorLengths();
339      const double * amplElement = amplObjective_->getElements();
340      int iColumn;
341      const double * columnScale = model->columnScale();
342      // use current linear cost region (already scaled)
343      if (includeLinear==1) {
344        memcpy(gradient_,model->costRegion(),numberExtendedColumns_*sizeof(double));
345      } else if (includeLinear==2) {
346        memset(gradient_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
347        if (!columnScale) {
348          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
349            gradient_[iColumn]= objective_[iColumn]*direction;
350          }
351        } else {
352          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
353            gradient_[iColumn]= objective_[iColumn]*direction*columnScale[iColumn];
354          }
355        }
356      } else {
357        memset(gradient_,0,numberExtendedColumns_*sizeof(double));
358      }
359      if (!columnScale) {
360        if (activated_) {
361          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
362            double valueI = solution[iColumn];
363            CoinBigIndex j;
364            for (j=columnAmplStart[iColumn];
365                 j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
366              int jColumn = columnAmpl[j];
367              double valueJ = solution[jColumn];
368              double elementValue = amplElement[j];
369              elementValue *= direction;
370              if (iColumn!=jColumn) {
371                offset += valueI*valueJ*elementValue;
372                double gradientI = valueJ*elementValue;
373                double gradientJ = valueI*elementValue;
374                gradient_[iColumn] += gradientI;
375                gradient_[jColumn] += gradientJ;
376              } else {
377                offset += 0.5*valueI*valueI*elementValue;
378                double gradientI = valueI*elementValue;
379                gradient_[iColumn] += gradientI;
380              }
381            }
382          }
383        }
384      } else {
385        // scaling
386        if (activated_) {
387          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
388            double valueI = solution[iColumn];
389            double scaleI = columnScale[iColumn]*direction;
390            CoinBigIndex j;
391            for (j=columnAmplStart[iColumn];
392                 j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
393              int jColumn = columnAmpl[j];
394              double valueJ = solution[jColumn];
395              double elementValue = amplElement[j];
396              double scaleJ = columnScale[jColumn];
397              elementValue *= scaleI*scaleJ;
398              if (iColumn!=jColumn) {
399                offset += valueI*valueJ*elementValue;
400                double gradientI = valueJ*elementValue;
401                double gradientJ = valueI*elementValue;
402                gradient_[iColumn] += gradientI;
403                gradient_[jColumn] += gradientJ;
404              } else {
405                offset += 0.5*valueI*valueI*elementValue;
406                double gradientI = valueI*elementValue;
407                gradient_[iColumn] += gradientI;
408              }
409            }
410          }
411        }
412      }
413    }
414    if (model)
415      offset *= model->optimizationDirection();
416    return gradient_;
417  }
418}
419 
420//-------------------------------------------------------------------
421// Clone
422//-------------------------------------------------------------------
423ClpObjective * ClpAmplObjective::clone() const
424{
425  return new ClpAmplObjective(*this);
426}
427// Resize objective
428void 
429ClpAmplObjective::resize(int newNumberColumns)
430{
431  if (numberColumns_!=newNumberColumns) {
432    abort();
433  } 
434 
435}
436// Delete columns in  objective
437void 
438ClpAmplObjective::deleteSome(int numberToDelete, const int * which) 
439{
440  if (numberToDelete)
441    abort();
442}
443
444// Load up ampl objective
445void 
446ClpAmplObjective::loadAmplObjective(void * amplInfo)
447{
448  fullMatrix_=false;
449  delete amplObjective_;
450  amplObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
451                                             start[numberColumns],element,column,start,NULL);
452  numberColumns_=numberColumns;
453  if (numberExtended>numberExtendedColumns_) {
454    if (objective_) {
455      // make correct size
456      double * newArray = new double[numberExtended];
457      memcpy(newArray,objective_,numberColumns_*sizeof(double));
458      delete [] objective_;
459      objective_ = newArray;
460      memset(objective_+numberColumns_,0,(numberExtended-numberColumns_)*sizeof(double));
461    }
462    if (gradient_) {
463      // make correct size
464      double * newArray = new double[numberExtended];
465      memcpy(newArray,gradient_,numberColumns_*sizeof(double));
466      delete [] gradient_;
467      gradient_ = newArray;
468      memset(gradient_+numberColumns_,0,(numberExtended-numberColumns_)*sizeof(double));
469    }
470    numberExtendedColumns_ = numberExtended;
471  } else {
472    numberExtendedColumns_ = numberColumns_;
473  }
474}
475/* Returns reduced gradient.Returns an offset (to be added to current one).
476 */
477double 
478ClpAmplObjective::reducedGradient(ClpSimplex * model, double * region,
479                                       bool useFeasibleCosts)
480{
481  int numberRows = model->numberRows();
482  int numberColumns=model->numberColumns();
483 
484  //work space
485  CoinIndexedVector  * workSpace = model->rowArray(0);
486 
487  CoinIndexedVector arrayVector;
488  arrayVector.reserve(numberRows+1);
489 
490  int iRow;
491#ifdef CLP_DEBUG
492  workSpace->checkClear();
493#endif
494  double * array = arrayVector.denseVector();
495  int * index = arrayVector.getIndices();
496  int number=0;
497  const double * costNow = gradient(model,model->solutionRegion(),offset_,
498                                    true,useFeasibleCosts ? 2 : 1);
499  double * cost = model->costRegion();
500  const int * pivotVariable = model->pivotVariable();
501  for (iRow=0;iRow<numberRows;iRow++) {
502    int iPivot=pivotVariable[iRow];
503    double value;
504    if (iPivot<numberColumns)
505      value = costNow[iPivot];
506    else if (!useFeasibleCosts) 
507      value = cost[iPivot];
508    else 
509      value=0.0;
510    if (value) {
511      array[iRow]=value;
512      index[number++]=iRow;
513    }
514  }
515  arrayVector.setNumElements(number);
516
517  // Btran basic costs
518  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
519  double * work = workSpace->denseVector();
520  ClpFillN(work,numberRows,0.0);
521  // now look at dual solution
522  double * rowReducedCost = region+numberColumns;
523  double * dual = rowReducedCost;
524  const double * rowCost = cost+numberColumns;
525  for (iRow=0;iRow<numberRows;iRow++) {
526    dual[iRow]=array[iRow];
527  }
528  double * dj = region;
529  ClpDisjointCopyN(costNow,numberColumns,dj);
530 
531  model->transposeTimes(-1.0,dual,dj);
532  for (iRow=0;iRow<numberRows;iRow++) {
533    // slack
534    double value = dual[iRow];
535    value += rowCost[iRow];
536    rowReducedCost[iRow]=value;
537  }
538  return offset_;
539}
540/* Returns step length which gives minimum of objective for
541   solution + theta * change vector up to maximum theta.
542   
543   arrays are numberColumns+numberRows
544*/
545double 
546ClpAmplObjective::stepLength(ClpSimplex * model,
547                                  const double * solution,
548                                  const double * change,
549                                  double maximumTheta,
550                                  double & currentObj,
551                                  double & predictedObj,
552                                  double & thetaObj)
553{
554  const double * cost = model->costRegion();
555  bool inSolve=true;
556  if (!cost) {
557    // not in solve
558    cost = objective_;
559    inSolve=false;
560  }
561  double delta=0.0;
562  double linearCost =0.0;
563  int numberRows = model->numberRows();
564  int numberColumns = model->numberColumns();
565  int numberTotal = numberColumns;
566  if (inSolve)
567    numberTotal += numberRows;
568  currentObj=0.0;
569  thetaObj=0.0;
570  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
571    delta += cost[iColumn]*change[iColumn];
572    linearCost += cost[iColumn]*solution[iColumn];
573  }
574  if (!activated_||!amplObjective_) {
575    currentObj=linearCost;
576    thetaObj =currentObj + delta*maximumTheta;
577    if (delta<0.0) {
578      return maximumTheta;
579    } else {
580      printf("odd linear direction %g\n",delta);
581      return 0.0;
582    }
583  }
584  assert (model);
585  bool scaling=false;
586  if ((model->rowScale()||
587       model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0)&&inSolve)
588    scaling=true;
589  const int * columnAmpl = amplObjective_->getIndices();
590  const CoinBigIndex * columnAmplStart = amplObjective_->getVectorStarts();
591  const int * columnAmplLength = amplObjective_->getVectorLengths();
592  const double * amplElement = amplObjective_->getElements();
593  double a=0.0;
594  double b=delta;
595  double c=0.0;
596  if (!scaling) {
597    if (!fullMatrix_) {
598      int iColumn;
599      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
600        double valueI = solution[iColumn];
601        double changeI = change[iColumn];
602        CoinBigIndex j;
603        for (j=columnAmplStart[iColumn];
604             j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
605          int jColumn = columnAmpl[j];
606          double valueJ = solution[jColumn];
607          double changeJ = change[jColumn];
608          double elementValue = amplElement[j];
609          if (iColumn!=jColumn) {
610            a += changeI*changeJ*elementValue;
611            b += (changeI*valueJ+changeJ*valueI)*elementValue;
612            c += valueI*valueJ*elementValue;
613          } else {
614            a += 0.5*changeI*changeI*elementValue;
615            b += changeI*valueI*elementValue;
616            c += 0.5*valueI*valueI*elementValue;
617          }
618        }
619      }
620    } else {
621      // full matrix stored
622      int iColumn;
623      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
624        double valueI = solution[iColumn];
625        double changeI = change[iColumn];
626        CoinBigIndex j;
627        for (j=columnAmplStart[iColumn];
628             j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
629          int jColumn = columnAmpl[j];
630          double valueJ = solution[jColumn];
631          double changeJ = change[jColumn];
632          double elementValue = amplElement[j];
633          valueJ *= elementValue;
634          a += changeI*changeJ*elementValue;
635          b += changeI*valueJ;
636          c += valueI*valueJ;
637        }
638      }
639      a *= 0.5;
640      c *= 0.5;
641    }
642  } else {
643    // scaling
644    // for now only if half
645    assert (!fullMatrix_);
646    const double * columnScale = model->columnScale();
647    double direction = model->optimizationDirection()*model->objectiveScale();
648    // direction is actually scale out not scale in
649    if (direction)
650      direction = 1.0/direction;
651    if (!columnScale) {
652      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
653        double valueI = solution[iColumn];
654        double changeI = change[iColumn];
655        CoinBigIndex j;
656        for (j=columnAmplStart[iColumn];
657             j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
658          int jColumn = columnAmpl[j];
659          double valueJ = solution[jColumn];
660          double changeJ = change[jColumn];
661          double elementValue = amplElement[j];
662          elementValue *= direction;
663          if (iColumn!=jColumn) {
664            a += changeI*changeJ*elementValue;
665            b += (changeI*valueJ+changeJ*valueI)*elementValue;
666            c += valueI*valueJ*elementValue;
667          } else {
668            a += 0.5*changeI*changeI*elementValue;
669            b += changeI*valueI*elementValue;
670            c += 0.5*valueI*valueI*elementValue;
671          }
672        }
673      }
674    } else {
675      // scaling
676      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
677        double valueI = solution[iColumn];
678        double changeI = change[iColumn];
679        double scaleI = columnScale[iColumn]*direction;
680        CoinBigIndex j;
681        for (j=columnAmplStart[iColumn];
682             j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
683          int jColumn = columnAmpl[j];
684          double valueJ = solution[jColumn];
685          double changeJ = change[jColumn];
686          double elementValue = amplElement[j];
687          elementValue *= scaleI*columnScale[jColumn];
688          if (iColumn!=jColumn) {
689            a += changeI*changeJ*elementValue;
690            b += (changeI*valueJ+changeJ*valueI)*elementValue;
691            c += valueI*valueJ*elementValue;
692          } else {
693            a += 0.5*changeI*changeI*elementValue;
694            b += changeI*valueI*elementValue;
695            c += 0.5*valueI*valueI*elementValue;
696          }
697        }
698      }
699    }
700  }
701  double theta;
702  //printf("Current cost %g\n",c+linearCost);
703  currentObj = c+linearCost;
704  thetaObj = currentObj + a*maximumTheta*maximumTheta+b*maximumTheta;
705  // minimize a*x*x + b*x + c
706  if (a<=0.0) {
707    theta = maximumTheta;
708  } else {
709    theta = -0.5*b/a;
710  }
711  predictedObj = currentObj + a*theta*theta+b*theta;
712  if (b>0.0) {
713    if (model->messageHandler()->logLevel()&32)
714      printf("a %g b %g c %g => %g\n",a,b,c,theta); 
715    b=0.0;
716  }
717  return CoinMin(theta,maximumTheta);
718}
719// Scale objective
720void 
721ClpAmplObjective::reallyScale(const double * columnScale) 
722{
723  abort();
724}
725/* Given a zeroed array sets nonlinear columns to 1.
726   Returns number of nonlinear columns
727*/
728int 
729ClpAmplObjective::markNonlinear(char * which)
730{
731  int iColumn;
732  const int * columnAmpl = amplObjective_->getIndices();
733  const CoinBigIndex * columnAmplStart = amplObjective_->getVectorStarts();
734  const int * columnAmplLength = amplObjective_->getVectorLengths();
735  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
736    CoinBigIndex j;
737    for (j=columnAmplStart[iColumn];
738         j<columnAmplStart[iColumn]+columnAmplLength[iColumn];j++) {
739      int jColumn = columnAmpl[j];
740      which[jColumn]=1;
741      which[iColumn]=1;
742    }
743  }
744  int numberNonLinearColumns = 0;
745  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
746    if(which[iColumn])
747      numberNonLinearColumns++;
748  }
749  return numberNonLinearColumns;
750}
Note: See TracBrowser for help on using the repository browser.