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