source: branches/sandbox/Cbc/src/ClpAmplStuff.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: 43.4 KB
Line 
1// Copyright (C) 2007, International Business Machines
2// Corporation and others.  All Rights Reserved.
3/* $Id: ClpAmplStuff.cpp 1200 2009-07-25 08:44:13Z forrest $ */
4
5#include "ClpConfig.h"
6#include "CbcConfig.h"
7#ifdef COIN_HAS_ASL
8#include "CoinPragma.hpp"
9#include "CoinHelperFunctions.hpp"
10#include "CoinIndexedVector.hpp"
11#include "ClpFactorization.hpp"
12#include "ClpSimplex.hpp"
13#include "ClpAmplObjective.hpp"
14#include "ClpConstraintAmpl.hpp"
15#include "ClpMessage.hpp"
16#include "CoinUtilsConfig.h"
17#include "CoinHelperFunctions.hpp"
18#include "CoinWarmStartBasis.hpp"
19#include "OsiSolverInterface.hpp"
20#include "CbcSolver.hpp"
21#include "Cbc_ampl.h"
22#include "CoinTime.hpp"
23#include "CglStored.hpp"
24#include "CoinModel.hpp"
25#include "CbcLinked.hpp"
26class CbcAmpl  : public CbcUser {
27
28public:
29    ///@name usage methods
30    //@{
31    /// Solve (whatever that means)
32    virtual void solve(CbcSolver * model, const char * options);
33    /// Returns true if function knows about option
34    virtual bool canDo(const char * options) ;
35    /** Import - gets full command arguments
36        Returns -1 - no action
37                 0 - data read in without error
38             1 - errors
39    */
40    virtual int importData(CbcSolver * model, int & argc, char ** & argv);
41    /// Export 1 OsiClpSolver, 2 CbcModel - add 10 if infeasible from odd situation
42    virtual void exportSolution(CbcSolver * model, int mode, const char * message = NULL) ;
43    /// Export Data (i.e. at very end)
44    virtual void exportData(CbcSolver * model);
45    /// Get useful stuff
46    virtual void fillInformation(CbcSolver * model,
47                                 CbcSolverUsefulData & info);
48    //@}
49    ///@name Constructors and destructors etc
50    //@{
51    /// Default Constructor
52    CbcAmpl();
53
54    /** Copy constructor .
55     */
56    CbcAmpl(const CbcAmpl & rhs);
57
58    /// Assignment operator
59    CbcAmpl & operator=(const CbcAmpl& rhs);
60
61    /// Clone
62    virtual CbcUser * clone() const;
63
64    /// Destructor
65    virtual ~CbcAmpl ();
66    //@}
67private:
68    ///@name Private member data
69    //@{
70    /// AMPL info
71    ampl_info info_;
72    //@}
73};
74// Mechanicsburg stuff
75CbcAmpl::CbcAmpl()
76        : CbcUser()
77{
78    userName_ = "mech";
79    memset(&info_, 0, sizeof(info_));
80}
81CbcAmpl::~CbcAmpl()
82{
83}
84// Copy constructor
85CbcAmpl::CbcAmpl ( const CbcAmpl & rhs)
86        : CbcUser(rhs)
87{
88    info_ = rhs.info_;
89}
90// Assignment operator
91CbcAmpl &
92CbcAmpl::operator=(const CbcAmpl & rhs)
93{
94    if (this != &rhs) {
95        CbcUser::operator=(rhs);
96        info_ = rhs.info_;
97    }
98    return *this;
99}
100// Clone
101CbcUser *
102CbcAmpl::clone() const
103{
104    return new CbcAmpl(*this);
105}
106// Solve (whatever that means)
107void
108CbcAmpl::solve(CbcSolver * controlModel, const char * options)
109{
110    assert (controlModel->model());
111    //OsiClpSolverInterface * clpSolver = dynamic_cast< OsiClpSolverInterface*> (model->solver());
112    //ClpSimplex * lpSolver = clpSolver->getModelPtr();
113    if (!strcmp(options, "cbc_load")) {
114    } else if (!strcmp(options, "cbc_quit")) {
115    } else {
116        printf("unknown option for CbcAmpl is %s\n", options);
117        abort();
118    }
119}
120// Returns true if function knows about option
121bool
122CbcAmpl::canDo(const char * options)
123{
124    return (!strcmp(options, "cbc_load") || !strcmp(options, "cbc_quit"));
125}
126/* Import - gets full command arguments
127   Returns -1 - no action
128            0 - data read in without error
129            1 - errors
130*/
131int
132CbcAmpl::importData(CbcSolver * control, int &argc, char ** & argv)
133{
134    CbcModel * babModel = control->model();
135    assert (babModel);
136    CoinMessageHandler * generalMessageHandler = babModel->messageHandler();
137    OsiClpSolverInterface * solver = dynamic_cast< OsiClpSolverInterface*> (control->model()->solver());
138    assert (solver);
139    CoinMessages generalMessages = solver->getModelPtr()->messages();
140    char generalPrint[10000];
141    OsiSolverLink * si = NULL;
142    ClpSimplex * lpSolver = solver->getModelPtr();
143    if (argc > 2 && !strcmp(argv[2], "-AMPL")) {
144        // see if log in list
145        bool printing = false;
146        for (int i = 1; i < argc; i++) {
147            if (!strncmp(argv[i], "log", 3)) {
148                const char * equals = strchr(argv[i], '=');
149                if (equals && atoi(equals + 1) > 0) {
150                    printing = true;
151                    info_.logLevel = atoi(equals + 1);
152                    control->setIntValue(LOGLEVEL, info_.logLevel);
153                    // mark so won't be overWritten
154                    info_.numberRows = -1234567;
155                    break;
156                }
157            }
158        }
159        union {
160            void * voidModel;
161            CoinModel * model;
162        } coinModelStart;
163        coinModelStart.model = NULL;
164        int returnCode = readAmpl(&info_, argc, argv, & coinModelStart.voidModel);
165        CoinModel * coinModel = coinModelStart.model;
166        if (returnCode)
167            return returnCode;
168        control->setReadMode(3); // so will start with parameters
169        // see if log in list (including environment)
170        for (int i = 1; i < info_.numberArguments; i++) {
171            if (!strcmp(info_.arguments[i], "log")) {
172                if (i < info_.numberArguments - 1 && atoi(info_.arguments[i+1]) > 0)
173                    printing = true;
174                break;
175            }
176        }
177        control->setPrinting(printing);
178        if (printing)
179            printf("%d rows, %d columns and %d elements\n",
180                   info_.numberRows, info_.numberColumns, info_.numberElements);
181        if (!coinModel) {
182            solver->loadProblem(info_.numberColumns, info_.numberRows, info_.starts,
183                                info_.rows, info_.elements,
184                                info_.columnLower, info_.columnUpper, info_.objective,
185                                info_.rowLower, info_.rowUpper);
186            if (info_.numberSos) {
187                // SOS
188                solver->setSOSData(info_.numberSos, info_.sosType, info_.sosStart,
189                                   info_.sosIndices, info_.sosReference);
190            }
191        } else {
192            // save
193            control->setOriginalCoinModel(coinModel);
194            // load from coin model
195            OsiSolverLink solver1;
196            OsiSolverInterface * solver2 = solver1.clone();
197            babModel->assignSolver(solver2, false);
198            si = dynamic_cast<OsiSolverLink *>(babModel->solver()) ;
199            assert (si != NULL);
200            si->setDefaultMeshSize(0.001);
201            // need some relative granularity
202            si->setDefaultBound(100.0);
203            double dextra3 = control->doubleValue(DEXTRA3);
204            if (dextra3)
205                si->setDefaultMeshSize(dextra3);
206            si->setDefaultBound(100000.0);
207            si->setIntegerPriority(1000);
208            si->setBiLinearPriority(10000);
209            CoinModel * model2 = (CoinModel *) coinModel;
210            int logLevel = control->intValue(LOGLEVEL);
211            si->load(*model2, true, logLevel);
212            // redo
213            solver = dynamic_cast< OsiClpSolverInterface*> (control->model()->solver());
214            lpSolver = solver->getModelPtr();
215            solver->messageHandler()->setLogLevel(0) ;
216            control->setIntValue(TESTOSI, 0);
217            if (info_.cut) {
218                printf("Sorry - can't do cuts with LOS as ruins delicate row order\n");
219                abort();
220            }
221        }
222        if (info_.cut) {
223            int numberRows = info_.numberRows;
224            int * whichRow = new int [numberRows];
225            // Row copy
226            const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
227            const double * elementByRow = matrixByRow->getElements();
228            const int * column = matrixByRow->getIndices();
229            const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
230            const int * rowLength = matrixByRow->getVectorLengths();
231
232            const double * rowLower = solver->getRowLower();
233            const double * rowUpper = solver->getRowUpper();
234            int nDelete = 0;
235            CglStored storedAmpl;
236            for (int iRow = 0; iRow < numberRows; iRow++) {
237                if (info_.cut[iRow]) {
238                    whichRow[nDelete++] = iRow;
239                    int start = rowStart[iRow];
240                    storedAmpl.addCut(rowLower[iRow], rowUpper[iRow],
241                                      rowLength[iRow], column + start, elementByRow + start);
242                }
243            }
244            control->addCutGenerator(&storedAmpl);
245            solver->deleteRows(nDelete, whichRow);
246            // and special matrix
247            si->cleanMatrix()->deleteRows(nDelete, whichRow);
248            delete [] whichRow;
249        }
250        // If we had a solution use it
251        if (info_.primalSolution) {
252            solver->setColSolution(info_.primalSolution);
253        }
254        // status
255        if (info_.rowStatus) {
256            unsigned char * statusArray = lpSolver->statusArray();
257            memset(statusArray, 0, lpSolver->numberColumns() + lpSolver->numberRows());
258            int i;
259            for (i = 0; i < info_.numberColumns; i++)
260                statusArray[i] = (char)info_.columnStatus[i];
261            statusArray += info_.numberColumns;
262            for (i = 0; i < info_.numberRows; i++)
263                statusArray[i] = (char)info_.rowStatus[i];
264            CoinWarmStartBasis * basis = lpSolver->getBasis();
265            solver->setWarmStart(basis);
266            delete basis;
267        }
268        freeArrays1(&info_);
269        // modify objective if necessary
270        solver->setObjSense(info_.direction);
271        solver->setDblParam(OsiObjOffset, info_.offset);
272        if (info_.offset) {
273            sprintf(generalPrint, "Ampl objective offset is %g",
274                    info_.offset);
275            generalMessageHandler->message(CLP_GENERAL, generalMessages)
276            << generalPrint
277            << CoinMessageEol;
278        }
279        // Set integer variables (unless nonlinear when set)
280        if (!info_.nonLinear) {
281            for (int i = info_.numberColumns - info_.numberIntegers;
282                    i < info_.numberColumns; i++)
283                solver->setInteger(i);
284        }
285        // change argc etc
286        argc = info_.numberArguments;
287        argv = info_.arguments;
288        return 0;
289    } else {
290        return -1;
291    }
292    abort();
293    return -1;
294}
295// Export 1 OsiClpSolver, 2 CbcModel - add 10 if infeasible from odd situation
296void
297CbcAmpl::exportSolution(CbcSolver * model, int mode, const char * message)
298{
299    OsiClpSolverInterface * solver = model->originalSolver();
300    if (!solver) {
301        solver = dynamic_cast< OsiClpSolverInterface*> (model->model()->solver());
302        assert (solver);
303    }
304    ClpSimplex * lpSolver = solver->getModelPtr();
305    int numberColumns = lpSolver->numberColumns();
306    int numberRows = lpSolver->numberRows();
307    double totalTime = CoinCpuTime() - model->startTime();
308    if (mode == 1) {
309        double value = lpSolver->getObjValue() * lpSolver->getObjSense();
310        char buf[300];
311        int pos = 0;
312        int iStat = lpSolver->status();
313        if (iStat == 0) {
314            pos += sprintf(buf + pos, "optimal," );
315        } else if (iStat == 1) {
316            // infeasible
317            pos += sprintf(buf + pos, "infeasible,");
318        } else if (iStat == 2) {
319            // unbounded
320            pos += sprintf(buf + pos, "unbounded,");
321        } else if (iStat == 3) {
322            pos += sprintf(buf + pos, "stopped on iterations or time,");
323        } else if (iStat == 4) {
324            iStat = 7;
325            pos += sprintf(buf + pos, "stopped on difficulties,");
326        } else if (iStat == 5) {
327            iStat = 3;
328            pos += sprintf(buf + pos, "stopped on ctrl-c,");
329        } else {
330            pos += sprintf(buf + pos, "status unknown,");
331            iStat = 6;
332        }
333        info_.problemStatus = iStat;
334        info_.objValue = value;
335        pos += sprintf(buf + pos, " objective %.*g", ampl_obj_prec(),
336                       value);
337        sprintf(buf + pos, "\n%d iterations",
338                lpSolver->getIterationCount());
339        free(info_.primalSolution);
340        info_.primalSolution = (double *) malloc(numberColumns * sizeof(double));
341        CoinCopyN(lpSolver->primalColumnSolution(), numberColumns, info_.primalSolution);
342        free(info_.dualSolution);
343        info_.dualSolution = (double *) malloc(numberRows * sizeof(double));
344        CoinCopyN(lpSolver->dualRowSolution(), numberRows, info_.dualSolution);
345        CoinWarmStartBasis * basis = lpSolver->getBasis();
346        free(info_.rowStatus);
347        info_.rowStatus = (int *) malloc(numberRows * sizeof(int));
348        free(info_.columnStatus);
349        info_.columnStatus = (int *) malloc(numberColumns * sizeof(int));
350        // Put basis in
351        int i;
352        // free,basic,ub,lb are 0,1,2,3
353        for (i = 0; i < numberRows; i++) {
354            CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
355            info_.rowStatus[i] = status;
356        }
357        for (i = 0; i < numberColumns; i++) {
358            CoinWarmStartBasis::Status status = basis->getStructStatus(i);
359            info_.columnStatus[i] = status;
360        }
361        // put buffer into info_
362        strcpy(info_.buffer, buf);
363        delete basis;
364    } else if (mode == 2) {
365        CbcModel * babModel = model->model();
366        int iStat = babModel->status();
367        int iStat2 = babModel->secondaryStatus();
368        double value = babModel->getObjValue() * lpSolver->getObjSense();
369        char buf[300];
370        int pos = 0;
371        if (iStat == 0) {
372            if (babModel->getObjValue() < 1.0e40) {
373                pos += sprintf(buf + pos, "optimal," );
374            } else {
375                // infeasible
376                iStat = 1;
377                pos += sprintf(buf + pos, "infeasible,");
378            }
379        } else if (iStat == 1) {
380            if (iStat2 != 6)
381                iStat = 3;
382            else
383                iStat = 4;
384            std::string minor[] = {"", "", "gap", "nodes", "time", "", "solutions", "user ctrl-c"};
385            pos += sprintf(buf + pos, "stopped on %s,", minor[iStat2].c_str());
386        } else if (iStat == 2) {
387            iStat = 7;
388            pos += sprintf(buf + pos, "stopped on difficulties,");
389        } else if (iStat == 5) {
390            iStat = 3;
391            pos += sprintf(buf + pos, "stopped on ctrl-c,");
392        } else {
393            pos += sprintf(buf + pos, "status unknown,");
394            iStat = 6;
395        }
396        info_.problemStatus = iStat;
397        info_.objValue = value;
398        if (babModel->getObjValue() < 1.0e40) {
399            int precision = ampl_obj_prec();
400            if (precision > 0)
401                pos += sprintf(buf + pos, " objective %.*g", precision,
402                               value);
403            else
404                pos += sprintf(buf + pos, " objective %g", value);
405        }
406        sprintf(buf + pos, "\n%d nodes, %d iterations, %g seconds",
407                babModel->getNodeCount(),
408                babModel->getIterationCount(),
409                totalTime);
410        if (babModel->bestSolution()) {
411            free(info_.primalSolution);
412            info_.primalSolution = (double *) malloc(numberColumns * sizeof(double));
413            CoinCopyN(lpSolver->primalColumnSolution(), numberColumns, info_.primalSolution);
414            free(info_.dualSolution);
415            info_.dualSolution = (double *) malloc(numberRows * sizeof(double));
416            CoinCopyN(lpSolver->dualRowSolution(), numberRows, info_.dualSolution);
417        } else {
418            info_.primalSolution = NULL;
419            info_.dualSolution = NULL;
420        }
421        // put buffer into info
422        strcpy(info_.buffer, buf);
423    } else if (mode == 11 || mode == 12) {
424        // infeasible
425        info_.problemStatus = 1;
426        info_.objValue = 1.0e100;
427        sprintf(info_.buffer, "%s", message);
428        info_.primalSolution = NULL;
429        info_.dualSolution = NULL;
430    }
431}
432// Export Data (i.e. at very end)
433void
434CbcAmpl::exportData(CbcSolver * model)
435{
436    writeAmpl(&info_);
437    freeArrays2(&info_);
438    freeArgs(&info_);
439}
440// Get useful stuff
441void
442CbcAmpl::fillInformation(CbcSolver * model,
443                         CbcSolverUsefulData & info)
444{
445    memset(&info, 0, sizeof(info));
446    info.priorities_ = info_.priorities;
447    info.sosPriority_ = info_.sosPriority;
448    info.branchDirection_ = info_.branchDirection;
449    info.primalSolution_ = info_.primalSolution;
450    info.pseudoDown_ = info_.pseudoDown;
451    info.pseudoUp_ = info_.pseudoUp;
452}
453void addAmplToCbc(CbcSolver * control)
454{
455    CbcAmpl ampl;
456    control->addUserFunction(&ampl);
457}
458extern "C" {
459    //# include "getstub.h"
460# include "asl_pfgh.h"
461}
462//#############################################################################
463// Constructors / Destructor / Assignment
464//#############################################################################
465
466//-------------------------------------------------------------------
467// Default Constructor
468//-------------------------------------------------------------------
469ClpAmplObjective::ClpAmplObjective ()
470        : ClpObjective()
471{
472    type_ = 12;
473    objective_ = NULL;
474    amplObjective_ = NULL;
475    gradient_ = NULL;
476    offset_ = 0.0;
477}
478// stolen from IPopt with changes
479typedef struct {
480    double obj_sign_;
481    ASL_pfgh * asl_;
482    double * non_const_x_;
483    int * column_; // for jacobian
484    int * rowStart_;
485    double * gradient_;
486    double * constraintValues_;
487    int nz_h_full_; // number of nonzeros in hessian
488    int nerror_;
489    bool objval_called_with_current_x_;
490    bool conval_called_with_current_x_;
491    bool jacval_called_with_current_x_;
492} CbcAmplInfo;
493#if 0
494static bool get_nlp_info(void * amplInfo, int & n, int & m, int & nnz_jac_g,
495                         int & nnz_h_lag)
496{
497    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
498    ASL_pfgh* asl = info->asl_;
499
500    n = n_var; // # of variables
501    m = n_con; // # of constraints
502    nnz_jac_g = nzc; // # of non-zeros in the jacobian
503    nnz_h_lag = info->nz_h_full_; // # of non-zeros in the hessian
504
505    return true;
506}
507
508static bool get_bounds_info(void * amplInfo, int  n, double * x_l,
509                            double * x_u, int  m, double * g_l, double * g_u)
510{
511    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
512    ASL_pfgh* asl = info->asl_;
513    assert(n == n_var);
514    assert(m == n_con);
515    int i;
516    for (i = 0; i < n; i++) {
517        x_l[i] = LUv[2*i];
518        x_u[i] = LUv[2*i+1];
519    }
520
521    for (i = 0; i < m; i++) {
522        g_l[i] = LUrhs[2*i];
523        g_u[i] = LUrhs[2*i+1];
524    }
525    return true;
526}
527
528#endif
529bool get_constraints_linearity(void * amplInfo, int  n,
530                               int * const_types)
531{
532    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
533    ASL_pfgh* asl = info->asl_;
534    //check that n is good
535    assert(n == n_con);
536    // check that there are no network constraints
537    assert(nlnc == 0 && lnc == 0);
538    //the first nlc constraints are non linear the rest is linear
539    int i;
540    for (i = 0; i < nlc; i++) {
541        const_types[i] = 1;
542    }
543    // the rest is linear
544    for (i = nlc; i < n_con; i++)
545        const_types[i] = 0;
546    return true;
547}
548#if 0
549bool get_starting_point(int  n, bool init_x, double * x, bool init_z,
550                        double * z_L, double * z_U, int  m, bool init_lambda, double * lambda)
551{
552    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
553    ASL_pfgh* asl = info->asl_;
554    assert(n == n_var);
555    assert(m == n_con);
556    int i;
557
558    if (init_x) {
559        for (i = 0; i < n; i++) {
560            if (havex0[i]) {
561                x[i] = X0[i];
562            } else {
563                x[i] = 0.0;
564            }
565        }
566    }
567
568    if (init_z) {
569        for (i = 0; i < n; i++) {
570            z_L[i] = z_U[i] = 1.0;
571        }
572    }
573
574    if (init_lambda) {
575        for (i = 0; i < m; i++) {
576            if (havepi0[i]) {
577                lambda[i] = pi0[i];
578            } else {
579                lambda[i] = 0.0;
580            }
581        }
582    }
583
584    return true;
585}
586#endif
587static bool internal_objval(CbcAmplInfo * info , double & obj_val)
588{
589    ASL_pfgh* asl = info->asl_;
590    info->objval_called_with_current_x_ = false; // in case the call below fails
591
592    if (n_obj == 0) {
593        obj_val = 0;
594        info->objval_called_with_current_x_ = true;
595        return true;
596    }  else {
597        double  retval = objval(0, info->non_const_x_, (fint*)info->nerror_);
598        if (!info->nerror_) {
599            obj_val = info->obj_sign_ * retval;
600            info->objval_called_with_current_x_ = true;
601            return true;
602        } else {
603            abort();
604        }
605    }
606
607    return false;
608}
609static bool internal_conval(CbcAmplInfo * info , double * g)
610{
611    ASL_pfgh* asl = info->asl_;
612    info->conval_called_with_current_x_ = false; // in case the call below fails
613    assert (g);
614
615    conval(info->non_const_x_, g, (fint*)info->nerror_);
616
617    if (!info->nerror_) {
618        info->conval_called_with_current_x_ = true;
619        return true;
620    } else {
621        abort();
622    }
623    return false;
624}
625
626static bool apply_new_x(CbcAmplInfo * info  , bool new_x, int  n, const double * x)
627{
628    ASL_pfgh* asl = info->asl_;
629
630    if (new_x) {
631        // update the flags so these methods are called
632        // before evaluating the hessian
633        info->conval_called_with_current_x_ = false;
634        info->objval_called_with_current_x_ = false;
635        info->jacval_called_with_current_x_ = false;
636
637        //copy the data to the non_const_x_
638        if (!info->non_const_x_) {
639            info->non_const_x_ = new double [n];
640        }
641
642        for (int  i = 0; i < n; i++) {
643            info->non_const_x_[i] = x[i];
644        }
645
646        // tell ampl that we have a new x
647        xknowne(info->non_const_x_, (fint*)info->nerror_);
648        return info->nerror_ ? false : true;
649    }
650
651    return true;
652}
653static bool eval_f(void * amplInfo, int  n, const double * x, bool new_x, double & obj_value)
654{
655    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
656    if (!apply_new_x(info, new_x, n, x)) {
657        return false;
658    }
659
660    return internal_objval(info, obj_value);
661}
662
663static bool eval_grad_f(void * amplInfo, int  n, const double * x, bool new_x, double * grad_f)
664{
665    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
666    ASL_pfgh* asl = info->asl_;
667    if (!apply_new_x(info, new_x, n, x)) {
668        return false;
669    }
670    int i;
671
672    if (n_obj == 0) {
673        for (i = 0; i < n; i++) {
674            grad_f[i] = 0.;
675        }
676    } else {
677        objgrd(0, info->non_const_x_, grad_f, (fint*)info->nerror_);
678        if (info->nerror_) {
679            return false;
680        }
681
682        if (info->obj_sign_ == -1) {
683            for (i = 0; i < n; i++) {
684                grad_f[i] = -grad_f[i];
685            }
686        }
687    }
688    return true;
689}
690static bool eval_g(void * amplInfo, int  n, const double * x, bool new_x, double * g)
691{
692    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
693#ifndef NDEBUG
694    ASL_pfgh* asl = info->asl_;
695#endif
696    // warning: n_var is a macro that assumes we have a variable called asl
697    assert(n == n_var);
698
699    if (!apply_new_x(info, new_x, n, x)) {
700        return false;
701    }
702
703    return internal_conval(info, g);
704}
705
706static bool eval_jac_g(void * amplInfo, int  n, const double * x, bool new_x,
707                       double * values)
708{
709    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
710    ASL_pfgh* asl = info->asl_;
711    assert(n == n_var);
712
713    assert (values);
714    if (!apply_new_x(info, new_x, n, x)) {
715        return false;
716    }
717
718    jacval(info->non_const_x_, values, (fint*)info->nerror_);
719    if (!info->nerror_) {
720        return true;
721    } else {
722        abort();
723    }
724    return false;
725}
726#if 0
727
728static bool eval_h(void * amplInfo, int  n, const double * x, bool new_x,
729                   double  obj_factor, int  m, const double * lambda,
730                   bool new_lambda, int  nele_hess, int * iRow,
731                   int * jCol, double * values)
732{
733    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
734    ASL_pfgh* asl = info->asl_;
735    assert(n == n_var);
736    assert(m == n_con);
737    int i;
738
739    if (iRow && jCol && !values) {
740        // setup the structure
741        int k = 0;
742        for (int i = 0; i < n; i++) {
743            for (int j = sputinfo->hcolstarts[i]; j < sputinfo->hcolstarts[i+1]; j++) {
744                //iRow[k] = i + 1;
745                iRow[k] = i;
746                jCol[k] = sputinfo->hrownos[j] + 1;
747                k++;
748            }
749        }
750        assert(k == nele_hess);
751        return true;
752    } else if (!iRow & !jCol && values) {
753        if (!apply_new_x(info, new_x, n, x)) {
754            return false;
755        }
756        if (!info->objval_called_with_current_x_) {
757            double  dummy;
758            internal_objval(info, dummy);
759            internal_conval(info, m, NULL);
760        }
761        if (!info->conval_called_with_current_x_) {
762            internal_conval(info, m, NULL);
763        }
764        // copy lambda to non_const_lambda - note, we do not store a copy like
765        // we do with x since lambda is only used here and not in other calls
766        double * non_const_lambda = new double [m];
767        for (i = 0; i < m; i++) {
768            non_const_lambda[i] = lambda[i];
769        }
770
771        real ow = info->obj_sign_ * obj_factor;
772        sphes(values, -1, &ow, non_const_lambda);
773
774        delete [] non_const_lambda;
775        return true;
776    } else {
777        assert(false && "Invalid combination of iRow, jCol, and values pointers");
778    }
779
780    return false;
781}
782#endif
783//-------------------------------------------------------------------
784// Useful Constructor
785//-------------------------------------------------------------------
786ClpAmplObjective::ClpAmplObjective (void * amplInfo)
787        : ClpObjective()
788{
789    type_ = 12;
790    activated_ = 1;
791    gradient_ = NULL;
792    objective_ = NULL;
793    offset_ = 0.0;
794    amplObjective_ = amplInfo;
795}
796
797//-------------------------------------------------------------------
798// Copy constructor
799//-------------------------------------------------------------------
800ClpAmplObjective::ClpAmplObjective (const ClpAmplObjective & rhs)
801        : ClpObjective(rhs)
802{
803    amplObjective_ = rhs.amplObjective_;
804    offset_ = rhs.offset_;
805    type_ = rhs.type_;
806    if (!amplObjective_) {
807        objective_ = NULL;
808        gradient_ = NULL;
809    } else {
810        CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
811        ASL_pfgh* asl = info->asl_;
812
813        int numberColumns = n_var;;
814        if (rhs.objective_) {
815            objective_ = new double [numberColumns];
816            memcpy(objective_, rhs.objective_, numberColumns*sizeof(double));
817        } else {
818            objective_ = NULL;
819        }
820        if (rhs.gradient_) {
821            gradient_ = new double [numberColumns];
822            memcpy(gradient_, rhs.gradient_, numberColumns*sizeof(double));
823        } else {
824            gradient_ = NULL;
825        }
826    }
827}
828
829
830//-------------------------------------------------------------------
831// Destructor
832//-------------------------------------------------------------------
833ClpAmplObjective::~ClpAmplObjective ()
834{
835    delete [] objective_;
836    delete [] gradient_;
837}
838
839//----------------------------------------------------------------
840// Assignment operator
841//-------------------------------------------------------------------
842ClpAmplObjective &
843ClpAmplObjective::operator=(const ClpAmplObjective & rhs)
844{
845    if (this != &rhs) {
846        delete [] objective_;
847        delete [] gradient_;
848        amplObjective_ = rhs.amplObjective_;
849        offset_ = rhs.offset_;
850        type_ = rhs.type_;
851        if (!amplObjective_) {
852            objective_ = NULL;
853            gradient_ = NULL;
854        } else {
855            CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
856            ASL_pfgh* asl = info->asl_;
857
858            int numberColumns = n_var;;
859            if (rhs.objective_) {
860                objective_ = new double [numberColumns];
861                memcpy(objective_, rhs.objective_, numberColumns*sizeof(double));
862            } else {
863                objective_ = NULL;
864            }
865            if (rhs.gradient_) {
866                gradient_ = new double [numberColumns];
867                memcpy(gradient_, rhs.gradient_, numberColumns*sizeof(double));
868            } else {
869                gradient_ = NULL;
870            }
871        }
872    }
873    return *this;
874}
875
876// Returns gradient
877double *
878ClpAmplObjective::gradient(const ClpSimplex * model,
879                           const double * solution, double & offset, bool refresh,
880                           int includeLinear)
881{
882    if (model)
883        assert (model->optimizationDirection() == 1.0);
884    bool scaling = false;
885    if (model && (model->rowScale() ||
886                  model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0))
887        scaling = true;
888    const double * cost = NULL;
889    if (model)
890        cost = model->costRegion();
891    if (!cost) {
892        // not in solve
893        cost = objective_;
894        scaling = false;
895    }
896    assert (!scaling);
897    if (!amplObjective_ || !solution || !activated_) {
898        offset = offset_;
899        return objective_;
900    } else {
901        if (refresh || !gradient_) {
902            CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
903            ASL_pfgh* asl = info->asl_;
904            int numberColumns = n_var;;
905
906            if (!gradient_)
907                gradient_ = new double[numberColumns];
908            assert (solution);
909            eval_grad_f(amplObjective_, numberColumns, solution, true, gradient_);
910            // Is this best way?
911            double objValue = 0.0;
912            eval_f(amplObjective_, numberColumns, solution, false, objValue);
913            double objValue2 = 0.0;
914            for (int i = 0; i < numberColumns; i++)
915                objValue2 += gradient_[i] * solution[i];
916            offset_ = objValue2 - objValue; // or other way???
917            if (model && model->optimizationDirection() != 1.0) {
918                offset *= model->optimizationDirection();
919                for (int i = 0; i < numberColumns; i++)
920                    gradient_[i] *= -1.0;
921            }
922        }
923        offset = offset_;
924        return gradient_;
925    }
926}
927
928//-------------------------------------------------------------------
929// Clone
930//-------------------------------------------------------------------
931ClpObjective * ClpAmplObjective::clone() const
932{
933    return new ClpAmplObjective(*this);
934}
935// Resize objective
936void
937ClpAmplObjective::resize(int newNumberColumns)
938{
939    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
940    ASL_pfgh* asl = info->asl_;
941    int numberColumns = n_var;;
942    if (numberColumns != newNumberColumns) {
943        abort();
944    }
945
946}
947// Delete columns in  objective
948void
949ClpAmplObjective::deleteSome(int numberToDelete, const int * which)
950{
951    if (numberToDelete)
952        abort();
953}
954/* Returns reduced gradient.Returns an offset (to be added to current one).
955 */
956double
957ClpAmplObjective::reducedGradient(ClpSimplex * model, double * region,
958                                  bool useFeasibleCosts)
959{
960    int numberRows = model->numberRows();
961    int numberColumns = model->numberColumns();
962
963    //work space
964    CoinIndexedVector  * workSpace = model->rowArray(0);
965
966    CoinIndexedVector arrayVector;
967    arrayVector.reserve(numberRows + 1);
968
969    int iRow;
970#ifdef CLP_DEBUG
971    workSpace->checkClear();
972#endif
973    double * array = arrayVector.denseVector();
974    int * index = arrayVector.getIndices();
975    int number = 0;
976    const double * costNow = gradient(model, model->solutionRegion(), offset_,
977                                      true, useFeasibleCosts ? 2 : 1);
978    double * cost = model->costRegion();
979    const int * pivotVariable = model->pivotVariable();
980    for (iRow = 0; iRow < numberRows; iRow++) {
981        int iPivot = pivotVariable[iRow];
982        double value;
983        if (iPivot < numberColumns)
984            value = costNow[iPivot];
985        else if (!useFeasibleCosts)
986            value = cost[iPivot];
987        else
988            value = 0.0;
989        if (value) {
990            array[iRow] = value;
991            index[number++] = iRow;
992        }
993    }
994    arrayVector.setNumElements(number);
995
996    // Btran basic costs
997    model->factorization()->updateColumnTranspose(workSpace, &arrayVector);
998    double * work = workSpace->denseVector();
999    ClpFillN(work, numberRows, 0.0);
1000    // now look at dual solution
1001    double * rowReducedCost = region + numberColumns;
1002    double * dual = rowReducedCost;
1003    const double * rowCost = cost + numberColumns;
1004    for (iRow = 0; iRow < numberRows; iRow++) {
1005        dual[iRow] = array[iRow];
1006    }
1007    double * dj = region;
1008    ClpDisjointCopyN(costNow, numberColumns, dj);
1009
1010    model->transposeTimes(-1.0, dual, dj);
1011    for (iRow = 0; iRow < numberRows; iRow++) {
1012        // slack
1013        double value = dual[iRow];
1014        value += rowCost[iRow];
1015        rowReducedCost[iRow] = value;
1016    }
1017    return offset_;
1018}
1019/* Returns step length which gives minimum of objective for
1020   solution + theta * change vector up to maximum theta.
1021
1022   arrays are numberColumns+numberRows
1023*/
1024double
1025ClpAmplObjective::stepLength(ClpSimplex * model,
1026                             const double * solution,
1027                             const double * change,
1028                             double maximumTheta,
1029                             double & currentObj,
1030                             double & predictedObj,
1031                             double & thetaObj)
1032{
1033    // Assume convex
1034    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1035    ASL_pfgh* asl = info->asl_;
1036
1037    int numberColumns = n_var;;
1038    double * tempSolution = new double [numberColumns];
1039    double * tempGradient = new double [numberColumns];
1040    // current
1041    eval_f(amplObjective_, numberColumns, solution, true, currentObj);
1042    double objA = currentObj;
1043    double thetaA = 0.0;
1044    // at maximum
1045    int i;
1046    for (i = 0; i < numberColumns; i++)
1047        tempSolution[i] = solution[i] + maximumTheta * change[i];
1048    eval_f(amplObjective_, numberColumns, tempSolution, true, thetaObj);
1049    double objC = thetaObj;
1050    double thetaC = maximumTheta;
1051    double objB = 0.5 * (objA + objC);
1052    double thetaB = 0.5 * maximumTheta;
1053    double gradientNorm = 1.0e6;
1054    while (gradientNorm > 1.0e-6 && thetaC - thetaA > 1.0e-8) {
1055        for (i = 0; i < numberColumns; i++)
1056            tempSolution[i] = solution[i] + thetaB * change[i];
1057        eval_grad_f(amplObjective_, numberColumns, tempSolution, true, tempGradient);
1058        eval_f(amplObjective_, numberColumns, tempSolution, false, objB);
1059        double changeObj = 0.0;
1060        gradientNorm = 0.0;
1061        for (i = 0; i < numberColumns; i++) {
1062            changeObj += tempGradient[i] * change[i];
1063            gradientNorm += tempGradient[i] * tempGradient[i];
1064        }
1065        gradientNorm = fabs(changeObj) / sqrt(gradientNorm);
1066        // Should try and get quadratic convergence by interpolation
1067        if (changeObj < 0.0) {
1068            // increasing is good
1069            thetaA = thetaB;
1070        } else {
1071            // decreasing is good
1072            thetaC = thetaB;
1073        }
1074        thetaB = 0.5 * (thetaA + thetaC);
1075    }
1076    delete [] tempSolution;
1077    delete [] tempGradient;
1078    predictedObj = objB;
1079    return thetaB;
1080}
1081// Return objective value (without any ClpModel offset) (model may be NULL)
1082double
1083ClpAmplObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
1084{
1085    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1086    ASL_pfgh* asl = info->asl_;
1087
1088    int numberColumns = n_var;;
1089    // current
1090    double currentObj = 0.0;
1091    eval_f(amplObjective_, numberColumns, solution, true, currentObj);
1092    return currentObj;
1093}
1094// Scale objective
1095void
1096ClpAmplObjective::reallyScale(const double * columnScale)
1097{
1098    abort();
1099}
1100/* Given a zeroed array sets nonlinear columns to 1.
1101   Returns number of nonlinear columns
1102*/
1103int
1104ClpAmplObjective::markNonlinear(char * which)
1105{
1106    int iColumn;
1107    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1108    ASL_pfgh* asl = info->asl_;
1109    int nonLinear = CoinMax(nlvc, nlvo);
1110    for (iColumn = 0; iColumn < nonLinear; iColumn++) {
1111        which[iColumn] = 1;
1112    }
1113    int numberNonLinearColumns = 0;
1114    int numberColumns = n_var;;
1115    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1116        if (which[iColumn])
1117            numberNonLinearColumns++;
1118    }
1119    return numberNonLinearColumns;
1120}
1121// Say we have new primal solution - so may need to recompute
1122void
1123ClpAmplObjective::newXValues()
1124{
1125    CbcAmplInfo * info = (CbcAmplInfo *) amplObjective_;
1126    info->conval_called_with_current_x_ = false;
1127    info->objval_called_with_current_x_ = false;
1128    info->jacval_called_with_current_x_ = false;
1129}
1130//#############################################################################
1131// Constructors / Destructor / Assignment
1132//#############################################################################
1133//-------------------------------------------------------------------
1134// Default Constructor
1135//-------------------------------------------------------------------
1136ClpConstraintAmpl::ClpConstraintAmpl ()
1137        : ClpConstraint()
1138{
1139    type_ = 3;
1140    column_ = NULL;
1141    coefficient_ = NULL;
1142    numberCoefficients_ = 0;
1143    amplInfo_ = NULL;
1144}
1145
1146//-------------------------------------------------------------------
1147// Useful Constructor
1148//-------------------------------------------------------------------
1149ClpConstraintAmpl::ClpConstraintAmpl (int row, void * amplInfo)
1150        : ClpConstraint()
1151{
1152    type_ = 3;
1153    rowNumber_ = row;
1154    amplInfo_ = amplInfo;
1155    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1156#ifndef NDEBUG
1157    ASL_pfgh* asl = info->asl_;
1158#endif
1159    // warning: nlc is a macro that assumes we have a variable called asl
1160    assert (rowNumber_ < nlc);
1161    numberCoefficients_ = info->rowStart_[rowNumber_+1] - info->rowStart_[rowNumber_];
1162    column_ = CoinCopyOfArray(info->column_ + info->rowStart_[rowNumber_], numberCoefficients_);
1163    coefficient_ = new double [numberCoefficients_];;
1164}
1165
1166//-------------------------------------------------------------------
1167// Copy constructor
1168//-------------------------------------------------------------------
1169ClpConstraintAmpl::ClpConstraintAmpl (const ClpConstraintAmpl & rhs)
1170        : ClpConstraint(rhs)
1171{
1172    numberCoefficients_ = rhs.numberCoefficients_;
1173    column_ = CoinCopyOfArray(rhs.column_, numberCoefficients_);
1174    coefficient_ = CoinCopyOfArray(rhs.coefficient_, numberCoefficients_);
1175}
1176
1177
1178//-------------------------------------------------------------------
1179// Destructor
1180//-------------------------------------------------------------------
1181ClpConstraintAmpl::~ClpConstraintAmpl ()
1182{
1183    delete [] column_;
1184    delete [] coefficient_;
1185}
1186
1187//----------------------------------------------------------------
1188// Assignment operator
1189//-------------------------------------------------------------------
1190ClpConstraintAmpl &
1191ClpConstraintAmpl::operator=(const ClpConstraintAmpl & rhs)
1192{
1193    if (this != &rhs) {
1194        delete [] column_;
1195        delete [] coefficient_;
1196        numberCoefficients_ = rhs.numberCoefficients_;
1197        column_ = CoinCopyOfArray(rhs.column_, numberCoefficients_);
1198        coefficient_ = CoinCopyOfArray(rhs.coefficient_, numberCoefficients_);
1199    }
1200    return *this;
1201}
1202//-------------------------------------------------------------------
1203// Clone
1204//-------------------------------------------------------------------
1205ClpConstraint * ClpConstraintAmpl::clone() const
1206{
1207    return new ClpConstraintAmpl(*this);
1208}
1209
1210// Returns gradient
1211int
1212ClpConstraintAmpl::gradient(const ClpSimplex * model,
1213                            const double * solution,
1214                            double * gradient,
1215                            double & functionValue,
1216                            double & offset,
1217                            bool useScaling,
1218                            bool refresh) const
1219{
1220    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1221    ASL_pfgh* asl = info->asl_;
1222    int numberColumns = n_var;;
1223    // If not done then do all
1224    if (!info->jacval_called_with_current_x_) {
1225        bool getStuff = eval_g(amplInfo_, numberColumns, solution, true, info->constraintValues_);
1226        assert (getStuff);
1227        getStuff = eval_jac_g(amplInfo_, numberColumns, solution, false, info->gradient_);
1228        assert (getStuff);
1229        info->jacval_called_with_current_x_ = getStuff;
1230    }
1231    if (refresh || !lastGradient_) {
1232        functionValue_ = info->constraintValues_[rowNumber_];
1233        offset_ = functionValue_; // sign??
1234        if (!lastGradient_)
1235            lastGradient_ = new double[numberColumns];
1236        CoinZeroN(lastGradient_, numberColumns);
1237        assert (!(model && model->rowScale() && useScaling));
1238        int i;
1239        int start = info->rowStart_[rowNumber_];
1240        assert (numberCoefficients_ == info->rowStart_[rowNumber_+1] - start);
1241        for (i = 0; i < numberCoefficients_; i++) {
1242            int iColumn = column_[i];
1243            double valueS = solution[iColumn];
1244            double valueG = info->gradient_[start+i];
1245            lastGradient_[iColumn] = valueG;
1246            offset_ -= valueS * valueG;
1247        }
1248    }
1249    functionValue = functionValue_;
1250    offset = offset_;
1251    memcpy(gradient, lastGradient_, numberColumns*sizeof(double));
1252    return 0;
1253}
1254// Resize constraint
1255void
1256ClpConstraintAmpl::resize(int newNumberColumns)
1257{
1258    abort();
1259}
1260// Delete columns in  constraint
1261void
1262ClpConstraintAmpl::deleteSome(int numberToDelete, const int * which)
1263{
1264    if (numberToDelete) {
1265        abort();
1266    }
1267}
1268// Scale constraint
1269void
1270ClpConstraintAmpl::reallyScale(const double * columnScale)
1271{
1272    abort();
1273}
1274/* Given a zeroed array sets nonlinear columns to 1.
1275   Returns number of nonlinear columns
1276*/
1277int
1278ClpConstraintAmpl::markNonlinear(char * which) const
1279{
1280    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1281    ASL_pfgh* asl = info->asl_;
1282    int iColumn;
1283    int numberNon = 0;
1284    int nonLinear = CoinMax(nlvc, nlvo);
1285    for (iColumn = 0; iColumn < numberCoefficients_; iColumn++) {
1286        int jColumn = column_[iColumn];
1287        if (jColumn < nonLinear) {
1288            which[jColumn] = 1;
1289            numberNon++;
1290        }
1291    }
1292    return numberNon;
1293}
1294/* Given a zeroed array sets possible nonzero coefficients to 1.
1295   Returns number of nonzeros
1296*/
1297int
1298ClpConstraintAmpl::markNonzero(char * which) const
1299{
1300    int iColumn;
1301    for (iColumn = 0; iColumn < numberCoefficients_; iColumn++) {
1302        which[column_[iColumn]] = 1;
1303    }
1304    return numberCoefficients_;
1305}
1306// Number of coefficients
1307int
1308ClpConstraintAmpl::numberCoefficients() const
1309{
1310    return numberCoefficients_;
1311}
1312// Say we have new primal solution - so may need to recompute
1313void
1314ClpConstraintAmpl::newXValues()
1315{
1316    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo_;
1317    info->conval_called_with_current_x_ = false;
1318    info->objval_called_with_current_x_ = false;
1319    info->jacval_called_with_current_x_ = false;
1320}
1321/* Load nonlinear part of problem from AMPL info
1322   Returns 0 if linear
1323   1 if quadratic objective
1324   2 if quadratic constraints
1325   3 if nonlinear objective
1326   4 if nonlinear constraints
1327   -1 on failure
1328*/
1329int
1330ClpSimplex::loadNonLinear(void * amplInfo, int & numberConstraints,
1331                          ClpConstraint ** & constraints)
1332{
1333    numberConstraints = 0;
1334    constraints = NULL;
1335    CbcAmplInfo * info = (CbcAmplInfo *) amplInfo;
1336    ASL_pfgh* asl = info->asl_;
1337    // For moment don't say quadratic
1338    int type = 0;
1339    if (nlo + nlc) {
1340        // nonlinear
1341        if (!nlc) {
1342            type = 3;
1343            delete objective_;
1344            objective_ = new ClpAmplObjective(amplInfo);
1345        } else {
1346            type = 4;
1347            numberConstraints = nlc;
1348            constraints = new ClpConstraint * [numberConstraints];
1349            if (nlo) {
1350                delete objective_;
1351                objective_ = new ClpAmplObjective(amplInfo);
1352            }
1353            for (int i = 0; i < numberConstraints; i++) {
1354                constraints[i] = new ClpConstraintAmpl(i, amplInfo);
1355            }
1356        }
1357    }
1358    return type;
1359}
1360#else
1361#include "ClpSimplex.hpp"
1362#include "ClpConstraint.hpp"
1363int
1364ClpSimplex::loadNonLinear(void * , int & ,
1365                          ClpConstraint ** & )
1366{
1367    abort();
1368    return 0;
1369}
1370#endif
Note: See TracBrowser for help on using the repository browser.