sccomp/source/solver/CoinMPSolver.cxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  *
   9  * This file incorporates work covered by the following license notice:
  10  *
  11  *   Licensed to the Apache Software Foundation (ASF) under one or more
  12  *   contributor license agreements. See the NOTICE file distributed
  13  *   with this work for additional information regarding copyright
  14  *   ownership. The ASF licenses this file to you under the Apache
  15  *   License, Version 2.0 (the "License"); you may not use this file
  16  *   except in compliance with the License. You may obtain a copy of
  17  *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
  18  */
  19
  20 #include <CoinMP.h>
  21 #include <CoinError.hpp>
  22
  23 #include "SolverComponent.hxx"
  24 #include "solver.hrc"
  25
  26 #include <com/sun/star/frame/XModel.hpp>
  27 #include <com/sun/star/table/CellAddress.hpp>
  28 #include <com/sun/star/uno/XComponentContext.hpp>
  29
  30 #include <rtl/math.hxx>
  31 #include <stdexcept>
  32 #include <vector>
  33
  34 using namespace com::sun::star;
  35
  36 class CoinMPSolver : public SolverComponent
  37 {
  38 public:
  39     CoinMPSolver() {}
  40     virtual ~CoinMPSolver() {}
  41
  42 private:
  43     virtual void SAL_CALL solve() throw(css::uno::RuntimeException, std::exception) SAL_OVERRIDE;
  44     virtual OUString SAL_CALL getImplementationName()
  45         throw(css::uno::RuntimeException, std::exception) SAL_OVERRIDE
  46     {
  47         return OUString("com.sun.star.comp.Calc.CoinMPSolver");
  48     }
  49     virtual OUString SAL_CALL getComponentDescription()
  50         throw (uno::RuntimeException, std::exception) SAL_OVERRIDE
  51     {
  52         return SolverComponent::GetResourceString( RID_COINMP_SOLVER_COMPONENT );
  53     }
  54 };
  55
  56 void SAL_CALL CoinMPSolver::solve() throw(uno::RuntimeException, std::exception)
  57 {
  58     uno::Reference<frame::XModel> xModel( mxDoc, uno::UNO_QUERY );
  59     if ( !xModel.is() )
  60         throw uno::RuntimeException();
  61
  62     maStatus.clear();
  63     mbSuccess = false;
  64
  65     xModel->lockControllers();
  66
  67     // collect variables in vector (?)
  68
  69     std::vector<table::CellAddress> aVariableCells;
  70     for (sal_Int32 nPos=0; nPos<maVariables.getLength(); nPos++)
  71         aVariableCells.push_back( maVariables[nPos] );
  72     size_t nVariables = aVariableCells.size();
  73     size_t nVar = 0;
  74
  75     // collect all dependent cells
  76
  77     ScSolverCellHashMap aCellsHash;
  78     aCellsHash[maObjective].reserve( nVariables + 1 );                  // objective function
  79
  80     for (sal_Int32 nConstrPos = 0; nConstrPos < maConstraints.getLength(); ++nConstrPos)
  81     {
  82         table::CellAddress aCellAddr = maConstraints[nConstrPos].Left;
  83         aCellsHash[aCellAddr].reserve( nVariables + 1 );                // constraints: left hand side
  84
  85         if ( maConstraints[nConstrPos].Right >>= aCellAddr )
  86             aCellsHash[aCellAddr].reserve( nVariables + 1 );            // constraints: right hand side
  87     }
  88
  89     // set all variables to zero
  90     //! store old values?
  91     //! use old values as initial values?
  92     std::vector<table::CellAddress>::const_iterator aVarIter;
  93     for ( aVarIter = aVariableCells.begin(); aVarIter != aVariableCells.end(); ++aVarIter )
  94     {
  95         SolverComponent::SetValue( mxDoc, *aVarIter, 0.0 );
  96     }
  97
  98     // read initial values from all dependent cells
  99     ScSolverCellHashMap::iterator aCellsIter;
 100     for ( aCellsIter = aCellsHash.begin(); aCellsIter != aCellsHash.end(); ++aCellsIter )
 101     {
 102         double fValue = SolverComponent::GetValue( mxDoc, aCellsIter->first );
 103         aCellsIter->second.push_back( fValue );                         // store as first element, as-is
 104     }
 105
 106     // loop through variables
 107     for ( aVarIter = aVariableCells.begin(); aVarIter != aVariableCells.end(); ++aVarIter )
 108     {
 109         SolverComponent::SetValue( mxDoc, *aVarIter, 1.0 );      // set to 1 to examine influence
 110
 111         // read value change from all dependent cells
 112         for ( aCellsIter = aCellsHash.begin(); aCellsIter != aCellsHash.end(); ++aCellsIter )
 113         {
 114             double fChanged = SolverComponent::GetValue( mxDoc, aCellsIter->first );
 115             double fInitial = aCellsIter->second.front();
 116             aCellsIter->second.push_back( fChanged - fInitial );
 117         }
 118
 119         SolverComponent::SetValue( mxDoc, *aVarIter, 2.0 );      // minimal test for linearity
 120
 121         for ( aCellsIter = aCellsHash.begin(); aCellsIter != aCellsHash.end(); ++aCellsIter )
 122         {
 123             double fInitial = aCellsIter->second.front();
 124             double fCoeff   = aCellsIter->second.back();       // last appended: coefficient for this variable
 125             double fTwo     = SolverComponent::GetValue( mxDoc, aCellsIter->first );
 126
 127             bool bLinear = rtl::math::approxEqual( fTwo, fInitial + 2.0 * fCoeff ) ||
 128                            rtl::math::approxEqual( fInitial, fTwo - 2.0 * fCoeff );
 129             // second comparison is needed in case fTwo is zero
 130             if ( !bLinear )
 131                 maStatus = SolverComponent::GetResourceString( RID_ERROR_NONLINEAR );
 132         }
 133
 134         SolverComponent::SetValue( mxDoc, *aVarIter, 0.0 );      // set back to zero for examining next variable
 135     }
 136
 137     xModel->unlockControllers();
 138
 139     if ( !maStatus.isEmpty() )
 140         return;
 141
 142     //
 143     // build parameter arrays for CoinMP
 144     //
 145
 146     // set objective function
 147
 148     const std::vector<double>& rObjCoeff = aCellsHash[maObjective];
 149     double* pObjectCoeffs = new double[nVariables];
 150     for (nVar=0; nVar<nVariables; nVar++)
 151         pObjectCoeffs[nVar] = rObjCoeff[nVar+1];
 152     double nObjectConst = rObjCoeff[0];             // constant term of objective
 153
 154     // add rows
 155
 156     size_t nRows = maConstraints.getLength();
 157     size_t nCompSize = nVariables * nRows;
 158     double* pCompMatrix = new double[nCompSize];    // first collect all coefficients, row-wise
 159     for (size_t i=0; i<nCompSize; i++)
 160         pCompMatrix[i] = 0.0;
 161
 162     double* pRHS = new double[nRows];
 163     char* pRowType = new char[nRows];
 164     for (size_t i=0; i<nRows; i++)
 165     {
 166         pRHS[i] = 0.0;
 167         pRowType[i] = 'N';
 168     }
 169
 170     for (sal_Int32 nConstrPos = 0; nConstrPos < maConstraints.getLength(); ++nConstrPos)
 171     {
 172         // integer constraints are set later
 173         sheet::SolverConstraintOperator eOp = maConstraints[nConstrPos].Operator;
 174         if ( eOp == sheet::SolverConstraintOperator_LESS_EQUAL ||
 175              eOp == sheet::SolverConstraintOperator_GREATER_EQUAL ||
 176              eOp == sheet::SolverConstraintOperator_EQUAL )
 177         {
 178             double fDirectValue = 0.0;
 179             bool bRightCell = false;
 180             table::CellAddress aRightAddr;
 181             const uno::Any& rRightAny = maConstraints[nConstrPos].Right;
 182             if ( rRightAny >>= aRightAddr )
 183                 bRightCell = true;                  // cell specified as right-hand side
 184             else
 185                 rRightAny >>= fDirectValue;         // constant value
 186
 187             table::CellAddress aLeftAddr = maConstraints[nConstrPos].Left;
 188
 189             const std::vector<double>& rLeftCoeff = aCellsHash[aLeftAddr];
 190             double* pValues = &pCompMatrix[nConstrPos * nVariables];
 191             for (nVar=0; nVar<nVariables; nVar++)
 192                 pValues[nVar] = rLeftCoeff[nVar+1];
 193
 194             // if left hand cell has a constant term, put into rhs value
 195             double fRightValue = -rLeftCoeff[0];
 196
 197             if ( bRightCell )
 198             {
 199                 const std::vector<double>& rRightCoeff = aCellsHash[aRightAddr];
 200                 // modify pValues with rhs coefficients
 201                 for (nVar=0; nVar<nVariables; nVar++)
 202                     pValues[nVar] -= rRightCoeff[nVar+1];
 203
 204                 fRightValue += rRightCoeff[0];      // constant term
 205             }
 206             else
 207                 fRightValue += fDirectValue;
 208
 209             switch ( eOp )
 210             {
 211                 case sheet::SolverConstraintOperator_LESS_EQUAL:    pRowType[nConstrPos] = 'L'; break;
 212                 case sheet::SolverConstraintOperator_GREATER_EQUAL: pRowType[nConstrPos] = 'G'; break;
 213                 case sheet::SolverConstraintOperator_EQUAL:         pRowType[nConstrPos] = 'E'; break;
 214                 default:
 215                     OSL_ENSURE( false, "unexpected enum type" );
 216             }
 217             pRHS[nConstrPos] = fRightValue;
 218         }
 219     }
 220
 221     // Find non-zero coefficients, column-wise
 222
 223     int* pMatrixBegin = new int[nVariables+1];
 224     int* pMatrixCount = new int[nVariables];
 225     double* pMatrix = new double[nCompSize];    // not always completely used
 226     int* pMatrixIndex = new int[nCompSize];
 227     int nMatrixPos = 0;
 228     for (nVar=0; nVar<nVariables; nVar++)
 229     {
 230         int nBegin = nMatrixPos;
 231         for (size_t nRow=0; nRow<nRows; nRow++)
 232         {
 233             double fCoeff = pCompMatrix[ nRow * nVariables + nVar ];    // row-wise
 234             if ( fCoeff != 0.0 )
 235             {
 236                 pMatrix[nMatrixPos] = fCoeff;
 237                 pMatrixIndex[nMatrixPos] = nRow;
 238                 ++nMatrixPos;
 239             }
 240         }
 241         pMatrixBegin[nVar] = nBegin;
 242         pMatrixCount[nVar] = nMatrixPos - nBegin;
 243     }
 244     pMatrixBegin[nVariables] = nMatrixPos;
 245     delete[] pCompMatrix;
 246     pCompMatrix = NULL;
 247
 248     // apply settings to all variables
 249
 250     double* pLowerBounds = new double[nVariables];
 251     double* pUpperBounds = new double[nVariables];
 252     for (nVar=0; nVar<nVariables; nVar++)
 253     {
 254         pLowerBounds[nVar] = mbNonNegative ? 0.0 : -DBL_MAX;
 255         pUpperBounds[nVar] = DBL_MAX;
 256
 257         // bounds could possibly be further restricted from single-cell constraints
 258     }
 259
 260     char* pColType = new char[nVariables];
 261     for (nVar=0; nVar<nVariables; nVar++)
 262         pColType[nVar] = mbInteger ? 'I' : 'C';
 263
 264     // apply single-var integer constraints
 265
 266     for (sal_Int32 nConstrPos = 0; nConstrPos < maConstraints.getLength(); ++nConstrPos)
 267     {
 268         sheet::SolverConstraintOperator eOp = maConstraints[nConstrPos].Operator;
 269         if ( eOp == sheet::SolverConstraintOperator_INTEGER ||
 270              eOp == sheet::SolverConstraintOperator_BINARY )
 271         {
 272             table::CellAddress aLeftAddr = maConstraints[nConstrPos].Left;
 273             // find variable index for cell
 274             for (nVar=0; nVar<nVariables; nVar++)
 275                 if ( AddressEqual( aVariableCells[nVar], aLeftAddr ) )
 276                 {
 277                     if ( eOp == sheet::SolverConstraintOperator_INTEGER )
 278                         pColType[nVar] = 'I';
 279                     else
 280                     {
 281                         pColType[nVar] = 'B';
 282                         pLowerBounds[nVar] = 0.0;
 283                         pUpperBounds[nVar] = 1.0;
 284                     }
 285                 }
 286         }
 287     }
 288
 289     int nObjectSense = mbMaximize ? SOLV_OBJSENS_MAX : SOLV_OBJSENS_MIN;
 290
 291     HPROB hProb = CoinCreateProblem("");
 292     int nResult = CoinLoadProblem( hProb, nVariables, nRows, nMatrixPos, 0,
 293                     nObjectSense, nObjectConst, pObjectCoeffs,
 294                     pLowerBounds, pUpperBounds, pRowType, pRHS, NULL,
 295                     pMatrixBegin, pMatrixCount, pMatrixIndex, pMatrix,
 296                     NULL, NULL, NULL );
 297     nResult = CoinLoadInteger( hProb, pColType );
 298
 299     delete[] pColType;
 300     delete[] pMatrixIndex;
 301     delete[] pMatrix;
 302     delete[] pMatrixCount;
 303     delete[] pMatrixBegin;
 304     delete[] pUpperBounds;
 305     delete[] pLowerBounds;
 306     delete[] pRowType;
 307     delete[] pRHS;
 308     delete[] pObjectCoeffs;
 309
 310     CoinSetRealOption( hProb, COIN_REAL_MAXSECONDS, mnTimeout );
 311     CoinSetRealOption( hProb, COIN_REAL_MIPMAXSEC, mnTimeout );
 312
 313     // TODO: handle (or remove) settings: epsilon, B&B depth
 314
 315     // solve model
 316
 317     nResult = CoinCheckProblem( hProb );
 318
 319     try
 320     {
 321         nResult = CoinOptimizeProblem( hProb, 0 );
 322     }
 323     catch (const CoinError& e)
 324     {
 325         throw std::runtime_error(e.message());
 326     }
 327
 328     mbSuccess = ( nResult == SOLV_CALL_SUCCESS );
 329     if ( mbSuccess )
 330     {
 331         // get solution
 332
 333         maSolution.realloc( nVariables );
 334         CoinGetSolutionValues( hProb, maSolution.getArray(), NULL, NULL, NULL );
 335         mfResultValue = CoinGetObjectValue( hProb );
 336     }
 337     else
 338     {
 339         int nSolutionStatus = CoinGetSolutionStatus( hProb );
 340         if ( nSolutionStatus == 1 )
 341             maStatus = SolverComponent::GetResourceString( RID_ERROR_INFEASIBLE );
 342         else if ( nSolutionStatus == 2 )
 343             maStatus = SolverComponent::GetResourceString( RID_ERROR_UNBOUNDED );
 344         // TODO: detect timeout condition and report as RID_ERROR_TIMEOUT
 345         // (currently reported as infeasible)
 346     }
 347
 348     CoinUnloadProblem( hProb );
 349 }
 350
 351 extern "C" SAL_DLLPUBLIC_EXPORT css::uno::XInterface * SAL_CALL
 352 com_sun_star_comp_Calc_CoinMPSolver_get_implementation(
 353     css::uno::XComponentContext *,
 354     css::uno::Sequence<css::uno::Any> const &)
 355 {
 356     return cppu::acquire(new CoinMPSolver());
 357 }
 358
 359 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */