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