source: branches/sandbox/Cbc/src/ClpAmplObjective.cpp @ 1286

Last change on this file since 1286 was 1286, checked in by EdwinStraver, 10 years ago

Changed formatting using AStyle -A4 -p

File size: 30.8 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.