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