Ignore:
Timestamp:
Sep 21, 2007 4:43:51 PM (12 years ago)
Author:
forrest
Message:

trying to move ampl stuff away

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/Cbc/src/ClpAmplStuff.cpp

    r700 r789  
    1212#include "ClpAmplObjective.hpp"
    1313#include "ClpConstraintAmpl.hpp"
     14#include "ClpMessage.hpp"
    1415#include "CoinUtilsConfig.h"
    1516#include "CoinHelperFunctions.hpp"
     17#include "CoinWarmStartBasis.hpp"
     18#include "OsiSolverInterface.hpp"
     19#include "CbcSolver.hpp"
     20#include "Cbc_ampl.h"
     21#include "CoinTime.hpp"
     22#include "CglStored.hpp"
     23#include "CoinModel.hpp"
     24#include "CbcLinked.hpp"
     25class CbcAmpl  : public CbcUser {
     26 
     27public:
     28  ///@name usage methods
     29  //@{
     30  /// Solve (whatever that means)
     31  virtual void solve(CbcSolver * model, const char * options);
     32  /// Returns true if function knows about option
     33  virtual bool canDo(const char * options) ;
     34  /** Import - gets full command arguments
     35      Returns -1 - no action
     36               0 - data read in without error
     37               1 - errors
     38  */
     39  virtual int importData(CbcSolver * model, int & argc, char * argv[]);
     40  /// Export 1 OsiClpSolver, 2 CbcModel - add 10 if infeasible from odd situation
     41  virtual void exportSolution(CbcSolver * model, int mode,const char * message=NULL) ;
     42  /// Export Data (i.e. at very end)
     43  virtual void exportData(CbcSolver * model);
     44  /// Get useful stuff
     45  virtual void fillInformation(CbcSolver * model,
     46                               CbcSolverUsefulData & info);
     47  //@}
     48  ///@name Constructors and destructors etc
     49  //@{
     50  /// Default Constructor
     51  CbcAmpl();
     52 
     53  /** Copy constructor .
     54   */ 
     55  CbcAmpl(const CbcAmpl & rhs);
     56 
     57  /// Assignment operator
     58  CbcAmpl & operator=(const CbcAmpl& rhs);
     59
     60  /// Clone
     61  virtual CbcUser * clone() const;
     62 
     63  /// Destructor
     64  virtual ~CbcAmpl ();
     65  //@}
     66private:
     67  ///@name Private member data
     68  //@{
     69  /// AMPL info
     70  ampl_info info_;
     71  //@}
     72};
     73// Mechanicsburg stuff
     74CbcAmpl::CbcAmpl()
     75  : CbcUser()
     76{
     77  userName_ = "mech";
     78  memset(&info_,0,sizeof(info_));
     79}
     80CbcAmpl::~CbcAmpl()
     81{
     82}
     83// Copy constructor
     84CbcAmpl::CbcAmpl ( const CbcAmpl & rhs)
     85  : CbcUser(rhs)
     86{
     87  info_ = rhs.info_;
     88}
     89// Assignment operator
     90CbcAmpl &
     91CbcAmpl::operator=(const CbcAmpl & rhs)
     92{
     93  if (this != &rhs) {
     94    CbcUser::operator=(rhs);
     95    info_ = rhs.info_;
     96  }
     97  return *this;
     98}
     99// Clone
     100CbcUser *
     101CbcAmpl::clone() const
     102{
     103  return new CbcAmpl(*this);
     104}
     105// Solve (whatever that means)
     106void
     107CbcAmpl::solve(CbcSolver * controlModel, const char * options)
     108{
     109  CbcModel * model = controlModel->model();
     110  assert (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    CoinModel * coinModel=NULL;
     160    int returnCode = readAmpl(&info_,argc, argv,(void **) (& coinModel));
     161    if (returnCode)
     162      return returnCode;
     163    control->setReadMode(3); // so will start with parameters
     164    // see if log in list (including environment)
     165    for (int i=1;i<info_.numberArguments;i++) {
     166      if (!strcmp(info_.arguments[i],"log")) {
     167        if (i<info_.numberArguments-1&&atoi(info_.arguments[i+1])>0)
     168          printing=true;
     169        break;
     170      }
     171    }
     172    control->setPrinting(printing);
     173    if (printing)
     174      printf("%d rows, %d columns and %d elements\n",
     175             info_.numberRows,info_.numberColumns,info_.numberElements);
     176    if (!coinModel) {
     177      solver->loadProblem(info_.numberColumns,info_.numberRows,info_.starts,
     178                          info_.rows,info_.elements,
     179                          info_.columnLower,info_.columnUpper,info_.objective,
     180                          info_.rowLower,info_.rowUpper);
     181      if (info_.numberSos) {
     182        // SOS
     183        solver->setSOSData(info_.numberSos,info_.sosType,info_.sosStart,
     184                           info_.sosIndices,info_.sosReference);
     185      }
     186    } else {
     187      // save
     188      control->setOriginalCoinModel(coinModel);
     189      // load from coin model
     190      OsiSolverLink solver1;
     191      OsiSolverInterface * solver2 = solver1.clone();
     192      babModel->assignSolver(solver2,false);
     193      si = dynamic_cast<OsiSolverLink *>(babModel->solver()) ;
     194      assert (si != NULL);
     195      si->setDefaultMeshSize(0.001);
     196      // need some relative granularity
     197      si->setDefaultBound(100.0);
     198      double dextra3 = control->doubleValue(DEXTRA3);
     199      if (dextra3)
     200        si->setDefaultMeshSize(dextra3);
     201      si->setDefaultBound(100000.0);
     202      si->setIntegerPriority(1000);
     203      si->setBiLinearPriority(10000);
     204      CoinModel * model2 = (CoinModel *) coinModel;
     205      int logLevel = control->intValue(LOGLEVEL);
     206      si->load(*model2,true,logLevel);
     207      // redo
     208      solver = dynamic_cast< OsiClpSolverInterface*> (control->model()->solver());
     209      lpSolver = solver->getModelPtr();
     210      solver->messageHandler()->setLogLevel(0) ;
     211      control->setIntValue(TESTOSI,0);
     212      if (info_.cut) {
     213        printf("Sorry - can't do cuts with LOS as ruins delicate row order\n");
     214        abort();
     215      }
     216    }
     217    if (info_.cut) {
     218      int numberRows = info_.numberRows;
     219      int * whichRow = new int [numberRows];
     220      // Row copy
     221      const CoinPackedMatrix * matrixByRow = solver->getMatrixByRow();
     222      const double * elementByRow = matrixByRow->getElements();
     223      const int * column = matrixByRow->getIndices();
     224      const CoinBigIndex * rowStart = matrixByRow->getVectorStarts();
     225      const int * rowLength = matrixByRow->getVectorLengths();
     226     
     227      const double * rowLower = solver->getRowLower();
     228      const double * rowUpper = solver->getRowUpper();
     229      int nDelete=0;
     230      CglStored storedAmpl;
     231      for (int iRow=0;iRow<numberRows;iRow++) {
     232        if (info_.cut[iRow]) {
     233          whichRow[nDelete++]=iRow;
     234          int start = rowStart[iRow];
     235          storedAmpl.addCut(rowLower[iRow],rowUpper[iRow],
     236                            rowLength[iRow],column+start,elementByRow+start);
     237        }
     238      }
     239      control->addCutGenerator(&storedAmpl);
     240      solver->deleteRows(nDelete,whichRow);
     241      // and special matrix
     242      si->cleanMatrix()->deleteRows(nDelete,whichRow);
     243      delete [] whichRow;
     244    }
     245    // If we had a solution use it
     246    if (info_.primalSolution) {
     247      solver->setColSolution(info_.primalSolution);
     248    }
     249    // status
     250    if (info_.rowStatus) {
     251      unsigned char * statusArray = lpSolver->statusArray();
     252      memset(statusArray,0,lpSolver->numberColumns()+lpSolver->numberRows());
     253      int i;
     254      for (i=0;i<info_.numberColumns;i++)
     255        statusArray[i]=(char)info_.columnStatus[i];
     256      statusArray+=info_.numberColumns;
     257      for (i=0;i<info_.numberRows;i++)
     258        statusArray[i]=(char)info_.rowStatus[i];
     259      CoinWarmStartBasis * basis = lpSolver->getBasis();
     260      solver->setWarmStart(basis);
     261      delete basis;
     262    }
     263    freeArrays1(&info_);
     264    // modify objective if necessary
     265    solver->setObjSense(info_.direction);
     266    solver->setDblParam(OsiObjOffset,info_.offset);
     267    if (info_.offset) {
     268      sprintf(generalPrint,"Ampl objective offset is %g",
     269              info_.offset);
     270      generalMessageHandler->message(CLP_GENERAL,generalMessages)
     271        << generalPrint
     272        <<CoinMessageEol;
     273    }
     274    // Set integer variables (unless nonlinear when set)
     275    if (!info_.nonLinear) {
     276      for (int i=info_.numberColumns-info_.numberIntegers;
     277           i<info_.numberColumns;i++)
     278        solver->setInteger(i);
     279    }
     280    // change argc etc
     281    argc = info_.numberArguments;
     282    argv = info_.arguments;
     283    return 0;
     284  } else {
     285    return -1;
     286  }
     287  abort();
     288  return -1;
     289}
     290// Export 1 OsiClpSolver, 2 CbcModel - add 10 if infeasible from odd situation
     291void
     292CbcAmpl::exportSolution(CbcSolver * model, int mode,const char * message)
     293{
     294  OsiClpSolverInterface * solver = model->originalSolver();
     295  if (!solver) {
     296    solver = dynamic_cast< OsiClpSolverInterface*> (model->model()->solver());
     297    assert (solver);
     298  }
     299  ClpSimplex * lpSolver = solver->getModelPtr();
     300  int numberColumns=lpSolver->numberColumns();
     301  int numberRows = lpSolver->numberRows();
     302  double totalTime = CoinCpuTime()-model->startTime();
     303  if (mode==1) {
     304    double value = lpSolver->getObjValue()*lpSolver->getObjSense();
     305    char buf[300];
     306    int pos=0;
     307    int iStat = lpSolver->status();
     308    if (iStat==0) {
     309      pos += sprintf(buf+pos,"optimal," );
     310    } else if (iStat==1) {
     311      // infeasible
     312      pos += sprintf(buf+pos,"infeasible,");
     313    } else if (iStat==2) {
     314      // unbounded
     315      pos += sprintf(buf+pos,"unbounded,");
     316    } else if (iStat==3) {
     317      pos += sprintf(buf+pos,"stopped on iterations or time,");
     318    } else if (iStat==4) {
     319      iStat = 7;
     320      pos += sprintf(buf+pos,"stopped on difficulties,");
     321    } else if (iStat==5) {
     322      iStat = 3;
     323      pos += sprintf(buf+pos,"stopped on ctrl-c,");
     324    } else {
     325      pos += sprintf(buf+pos,"status unknown,");
     326      iStat=6;
     327    }
     328    info_.problemStatus=iStat;
     329    info_.objValue = value;
     330    pos += sprintf(buf+pos," objective %.*g",ampl_obj_prec(),
     331                   value);
     332    sprintf(buf+pos,"\n%d iterations",
     333            lpSolver->getIterationCount());
     334    free(info_.primalSolution);
     335    info_.primalSolution = (double *) malloc(numberColumns*sizeof(double));
     336    CoinCopyN(lpSolver->primalColumnSolution(),numberColumns,info_.primalSolution);
     337    free(info_.dualSolution);
     338    info_.dualSolution = (double *) malloc(numberRows*sizeof(double));
     339    CoinCopyN(lpSolver->dualRowSolution(),numberRows,info_.dualSolution);
     340    CoinWarmStartBasis * basis = lpSolver->getBasis();
     341    free(info_.rowStatus);
     342    info_.rowStatus = (int *) malloc(numberRows*sizeof(int));
     343    free(info_.columnStatus);
     344    info_.columnStatus = (int *) malloc(numberColumns*sizeof(int));
     345    // Put basis in
     346    int i;
     347    // free,basic,ub,lb are 0,1,2,3
     348    for (i=0;i<numberRows;i++) {
     349      CoinWarmStartBasis::Status status = basis->getArtifStatus(i);
     350      info_.rowStatus[i]=status;
     351    }
     352    for (i=0;i<numberColumns;i++) {
     353      CoinWarmStartBasis::Status status = basis->getStructStatus(i);
     354      info_.columnStatus[i]=status;
     355    }
     356    // put buffer into info_
     357    strcpy(info_.buffer,buf);
     358    delete basis;
     359  } else if (mode==2) {
     360    CbcModel * babModel = model->model();
     361    int iStat = babModel->status();
     362    int iStat2 = babModel->secondaryStatus();
     363    double value = babModel->getObjValue()*lpSolver->getObjSense();
     364    char buf[300];
     365    int pos=0;
     366    if (iStat==0) {
     367      if (babModel->getObjValue()<1.0e40) {
     368        pos += sprintf(buf+pos,"optimal," );
     369      } else {
     370        // infeasible
     371        iStat=1;
     372        pos += sprintf(buf+pos,"infeasible,");
     373      }
     374    } else if (iStat==1) {
     375      if (iStat2!=6)
     376        iStat=3;
     377      else
     378        iStat=4;
     379      std::string minor[]={"","","gap","nodes","time","","solutions","user ctrl-c"};
     380      pos += sprintf(buf+pos,"stopped on %s,",minor[iStat2].c_str());
     381    } else if (iStat==2) {
     382      iStat = 7;
     383      pos += sprintf(buf+pos,"stopped on difficulties,");
     384    } else if (iStat==5) {
     385      iStat = 3;
     386      pos += sprintf(buf+pos,"stopped on ctrl-c,");
     387    } else {
     388      pos += sprintf(buf+pos,"status unknown,");
     389      iStat=6;
     390    }
     391    info_.problemStatus=iStat;
     392    info_.objValue = value;
     393    if (babModel->getObjValue()<1.0e40) {
     394      int precision = ampl_obj_prec();
     395      if (precision>0)
     396        pos += sprintf(buf+pos," objective %.*g",precision,
     397                       value);
     398      else
     399        pos += sprintf(buf+pos," objective %g",value);
     400    }
     401    sprintf(buf+pos,"\n%d nodes, %d iterations, %g seconds",
     402            babModel->getNodeCount(),
     403            babModel->getIterationCount(),
     404            totalTime);
     405    if (babModel->bestSolution()) {
     406      free(info_.primalSolution);
     407      info_.primalSolution = (double *) malloc(numberColumns*sizeof(double));
     408      CoinCopyN(lpSolver->primalColumnSolution(),numberColumns,info_.primalSolution);
     409      free(info_.dualSolution);
     410      info_.dualSolution = (double *) malloc(numberRows*sizeof(double));
     411      CoinCopyN(lpSolver->dualRowSolution(),numberRows,info_.dualSolution);
     412    } else {
     413      info_.primalSolution=NULL;
     414      info_.dualSolution=NULL;
     415    }
     416    // put buffer into info
     417    strcpy(info_.buffer,buf);
     418  } else if (mode==11||mode==12) {
     419    // infeasible
     420    info_.problemStatus=1;
     421    info_.objValue = 1.0e100;
     422    sprintf(info_.buffer,"%s",message);
     423    info_.primalSolution=NULL;
     424    info_.dualSolution=NULL;
     425  }
     426}
     427// Export Data (i.e. at very end)
     428void
     429CbcAmpl::exportData(CbcSolver * model)
     430{
     431  writeAmpl(&info_);
     432  freeArrays2(&info_);
     433  freeArgs(&info_);
     434}
     435// Get useful stuff
     436void
     437CbcAmpl::fillInformation(CbcSolver * model,
     438                         CbcSolverUsefulData & info)
     439{
     440  memset(&info,0,sizeof(info));
     441  info.priorities_ = info_.priorities;
     442  info.sosPriority_ = info_.sosPriority;
     443  info.branchDirection_ = info_.branchDirection;
     444  info.primalSolution_ = info_.primalSolution;
     445  info.pseudoDown_ = info_.pseudoDown;
     446  info.pseudoUp_ = info_.pseudoUp;
     447}
     448void addAmplToCbc(CbcSolver * control)
     449{
     450  CbcAmpl ampl;
     451  control->addUserFunction(&ampl);
     452}
    16453extern "C" {
    17454  //# include "getstub.h"
Note: See TracChangeset for help on using the changeset viewer.