ODFstudio/odfstudio
2
0
Fork 0
Code Issues Packages Activity

solvers: extract the CoinMP problem build and solve into a helper

Browse Source 
CoinMPSolver did the whole job in one function: sample the model by
perturbing each variable, build the CoinMP matrix from those
coefficients, solve, and read the result back. Move the matrix
build and solve into a standalone helper, coinmpSolveLinearModel,
that  takes the sampled coefficients plus the variable bounds and
column  types and returns the solution and status.

The helper treats the first entry of each cell's coefficient vector
as the value of an affine model at the origin and the rest as the
per variable slopes, which is what the existing sampling already
produces. That keeps the same convention while letting a future
caller supply coefficients sampled around a different point and its
own bounds, which is what a sequential linear programming loop
needs.

No behaviour change. CoinMPSolver still samples and sets up bounds
and column types exactly as before, then calls the helper. The
CoinMP solver unit test passes unchanged.

Signed-off-by: Tomaž Vajngerl <tomaz.vajngerl@collabora.co.uk>
Change-Id: I6d373a48caa9d67a932a599068368bae978f1f22
Reviewed-on: https://gerrit.collaboraoffice.com/c/online/+/4648
Tested-by: Tomaž Vajngerl <tomaz.vajngerl@collabora.com>
Reviewed-by: Miklos Vajna <vmiklos@collabora.com>
(cherry picked from commit 197f790a93880f56dd53db43e6bbb61c22ccfdc2)
Reviewed-on: https://gerrit.libreoffice.org/c/core/+/208337
Reviewed-by: Xisco Fauli <xiscofauli@libreoffice.org>
Tested-by: Xisco Fauli <xiscofauli@libreoffice.org>
This commit is contained in:
Tomaž Vajngerl
2026-07-21 11:39:34 +02:00
committed by Xisco Fauli
parent b45d468a3d
commit 2c87d392dc
4 changed files with 275 additions and 187 deletions
Download Patch File Download Diff File Expand all files Collapse all files
sc/Library_solver.mk
+1
View File
@@ -52,6 +52,7 @@ $(eval $(call gb_Library_add_exception_objects,solver,\
sc/source/solver/SwarmSolver \
sc/source/solver/SolverComponent \
$(if $(ENABLE_COINMP), sc/source/solver/CoinMPSolver) \
$(if $(ENABLE_COINMP), sc/source/solver/CoinMPModel) \
))
# vim: set noet sw=4 ts=4:
sc/source/solver/CoinMPModel.cxx
+207
View File
@@ -0,0 +1,207 @@
/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
* This file is part of the LibreOffice project.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
*/
#include <CoinMP.h>
#include <CoinError.hpp>
#include "CoinMPModel.hxx"
#include <memory>
#include <stdexcept>
using namespace com::sun::star;
CoinMpSolveResult coinmpSolveLinearModel(const uno::Sequence<table::CellAddress>& rVariables,
const uno::Sequence<sheet::SolverConstraint>& rConstraints,
const table::CellAddress& rObjective,
const ScSolverCellHashMap& rCells,
const std::vector<double>& rLowerBounds,
const std::vector<double>& rUpperBounds,
const std::vector<char>& rColumnType, bool bMaximize,
double fTimeoutSeconds)
{
CoinMpSolveResult aResult;
const size_t nVariables = rVariables.size();
size_t nVar = 0;
// objective: slope of each variable, plus the constant term at the origin
const std::vector<double>& rObjectCoefficients = rCells.at(rObjective);
std::unique_ptr<double[]> pObjectCoefficients(new double[nVariables]);
for (nVar = 0; nVar < nVariables; nVar++)
pObjectCoefficients[nVar] = rObjectCoefficients[nVar + 1];
double fObjectConstant = rObjectCoefficients[0];
// rows: collect all coefficients row-wise first
size_t nRows = rConstraints.getLength();
size_t nDenseSize = nVariables * nRows;
std::unique_ptr<double[]> pDenseMatrix(new double[nDenseSize]);
for (size_t i = 0; i < nDenseSize; i++)
pDenseMatrix[i] = 0.0;
std::unique_ptr<double[]> pRHS(new double[nRows]);
std::unique_ptr<char[]> pRowType(new char[nRows]);
for (size_t i = 0; i < nRows; i++)
{
pRHS[i] = 0.0;
pRowType[i] = 'N';
}
for (sal_Int32 nConstraintPosition = 0; nConstraintPosition < rConstraints.getLength();
++nConstraintPosition)
{
// integer and binary constraints are reflected in the column types
sheet::SolverConstraintOperator eOperator = rConstraints[nConstraintPosition].Operator;
if (eOperator == sheet::SolverConstraintOperator_LESS_EQUAL
|| eOperator == sheet::SolverConstraintOperator_GREATER_EQUAL
|| eOperator == sheet::SolverConstraintOperator_EQUAL)
{
double fDirectValue = 0.0;
bool bRightCell = false;
table::CellAddress aRightAddress;
const uno::Any& rRightAny = rConstraints[nConstraintPosition].Right;
if (rRightAny >>= aRightAddress)
bRightCell = true;
else
rRightAny >>= fDirectValue;
table::CellAddress aLeftAddress = rConstraints[nConstraintPosition].Left;
const std::vector<double>& rLeftCoefficients = rCells.at(aLeftAddress);
double* pValues = &pDenseMatrix[nConstraintPosition * nVariables];
for (nVar = 0; nVar < nVariables; nVar++)
pValues[nVar] = rLeftCoefficients[nVar + 1];
// if the left hand cell has a constant term, move it to the right
double fRightValue = -rLeftCoefficients[0];
if (bRightCell)
{
const std::vector<double>& rRightCoefficients = rCells.at(aRightAddress);
for (nVar = 0; nVar < nVariables; nVar++)
pValues[nVar] -= rRightCoefficients[nVar + 1];
fRightValue += rRightCoefficients[0];
}
else
fRightValue += fDirectValue;
switch (eOperator)
{
case sheet::SolverConstraintOperator_LESS_EQUAL:
pRowType[nConstraintPosition] = 'L';
break;
case sheet::SolverConstraintOperator_GREATER_EQUAL:
pRowType[nConstraintPosition] = 'G';
break;
case sheet::SolverConstraintOperator_EQUAL:
pRowType[nConstraintPosition] = 'E';
break;
default:
break;
}
pRHS[nConstraintPosition] = fRightValue;
}
}
// turn the row-wise matrix into the column-wise form CoinMP wants
std::unique_ptr<int[]> pMatrixBegin(new int[nVariables + 1]);
std::unique_ptr<int[]> pMatrixCount(new int[nVariables]);
std::unique_ptr<double[]> pMatrix(new double[nDenseSize]);
std::unique_ptr<int[]> pMatrixIndex(new int[nDenseSize]);
int nMatrixPosition = 0;
for (nVar = 0; nVar < nVariables; nVar++)
{
int nBegin = nMatrixPosition;
for (size_t nRow = 0; nRow < nRows; nRow++)
{
double fCoefficient = pDenseMatrix[nRow * nVariables + nVar];
if (fCoefficient != 0.0)
{
pMatrix[nMatrixPosition] = fCoefficient;
pMatrixIndex[nMatrixPosition] = nRow;
++nMatrixPosition;
}
}
pMatrixBegin[nVar] = nBegin;
pMatrixCount[nVar] = nMatrixPosition - nBegin;
}
pMatrixBegin[nVariables] = nMatrixPosition;
pDenseMatrix.reset();
std::unique_ptr<double[]> pLowerBounds(new double[nVariables]);
std::unique_ptr<double[]> pUpperBounds(new double[nVariables]);
std::unique_ptr<char[]> pColumnType(new char[nVariables]);
for (nVar = 0; nVar < nVariables; nVar++)
{
pLowerBounds[nVar] = rLowerBounds[nVar];
pUpperBounds[nVar] = rUpperBounds[nVar];
pColumnType[nVar] = rColumnType[nVar];
}
int nObjectSense = bMaximize ? SOLV_OBJSENS_MAX : SOLV_OBJSENS_MIN;
HPROB hProblem = CoinCreateProblem("");
int nResult = CoinLoadProblem(hProblem, nVariables, nRows, nMatrixPosition, 0, nObjectSense,
fObjectConstant, pObjectCoefficients.get(), pLowerBounds.get(),
pUpperBounds.get(), pRowType.get(), pRHS.get(), nullptr,
pMatrixBegin.get(), pMatrixCount.get(), pMatrixIndex.get(),
pMatrix.get(), nullptr, nullptr, nullptr);
if (nResult == SOLV_CALL_SUCCESS)
nResult = CoinLoadInteger(hProblem, pColumnType.get());
pColumnType.reset();
pMatrixIndex.reset();
pMatrix.reset();
pMatrixCount.reset();
pMatrixBegin.reset();
pUpperBounds.reset();
pLowerBounds.reset();
pRowType.reset();
pRHS.reset();
pObjectCoefficients.reset();
CoinSetRealOption(hProblem, COIN_REAL_MAXSECONDS, fTimeoutSeconds);
CoinSetRealOption(hProblem, COIN_REAL_MIPMAXSEC, fTimeoutSeconds);
if (nResult == SOLV_CALL_SUCCESS)
nResult = CoinCheckProblem(hProblem);
if (nResult == SOLV_CALL_SUCCESS)
{
try
{
nResult = CoinOptimizeProblem(hProblem, 0);
}
catch (const CoinError& rError)
{
CoinUnloadProblem(hProblem);
throw std::runtime_error(rError.message());
}
}
aResult.bSuccess = (nResult == SOLV_CALL_SUCCESS);
if (aResult.bSuccess)
{
aResult.aSolution.resize(nVariables);
CoinGetSolutionValues(hProblem, aResult.aSolution.data(), nullptr, nullptr, nullptr);
aResult.fObjective = CoinGetObjectValue(hProblem);
}
else
{
aResult.nStatus = CoinGetSolutionStatus(hProblem);
}
CoinUnloadProblem(hProblem);
return aResult;
}
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
sc/source/solver/CoinMPModel.hxx
+46
View File
@@ -0,0 +1,46 @@
/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
/*
* This file is part of the LibreOffice project.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
*/
#pragma once
#include <com/sun/star/sheet/SolverConstraint.hpp>
#include <com/sun/star/table/CellAddress.hpp>
#include <com/sun/star/uno/Sequence.hxx>
#include "SolverComponent.hxx"
#include <vector>
// Result of solving one linear problem with CoinMP.
struct CoinMpSolveResult
{
bool bSuccess = false;
// CoinGetSolutionStatus when the solve did not succeed: 1 infeasible,
// 2 unbounded, other values for the remaining cases.
int nStatus = 0;
std::vector<double> aSolution;
double fObjective = 0.0;
};
// Build a linear problem from already sampled coefficients and solve it with
// CoinMP. For the objective and for each constraint cell, rCells holds a vector
// whose first entry is the value of an affine model of that cell at the origin
// (all variables zero) and whose remaining entries are the slope of each
// variable. The caller supplies the variable bounds and column types (which
// encode integer and binary variables), so the same routine serves both a
// one-shot linear solve and a linearized step of an outer nonlinear loop.
CoinMpSolveResult coinmpSolveLinearModel(
const css::uno::Sequence<css::table::CellAddress>& rVariables,
const css::uno::Sequence<css::sheet::SolverConstraint>& rConstraints,
const css::table::CellAddress& rObjective, const ScSolverCellHashMap& rCells,
const std::vector<double>& rLowerBounds, const std::vector<double>& rUpperBounds,
const std::vector<char>& rColumnType, bool bMaximize, double fTimeoutSeconds);
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
sc/source/solver/CoinMPSolver.cxx
+21 -187
View File
@@ -17,9 +17,7 @@
* the License at http://www.apache.org/licenses/LICENSE-2.0 .
*/
#include <CoinMP.h>
#include <CoinError.hpp>
#include "CoinMPModel.hxx"
#include "SolverComponent.hxx"
#include "strings.hrc"
@@ -27,7 +25,7 @@
#include <com/sun/star/table/CellAddress.hpp>
#include <rtl/math.hxx>
#include <stdexcept>
#include <algorithm>
#include <vector>
#include <float.h>
@@ -136,128 +134,12 @@ void SAL_CALL CoinMPSolver::solve()
if ( !maStatus.isEmpty() )
return;
//
// build parameter arrays for CoinMP
//
// set objective function
const std::vector<double>& rObjCoeff = aCellsHash[maObjective];
std::unique_ptr<double[]> pObjectCoeffs(new double[nVariables]);
for (nVar=0; nVar<nVariables; nVar++)
pObjectCoeffs[nVar] = rObjCoeff[nVar+1];
double nObjectConst = rObjCoeff[0]; // constant term of objective
// add rows
size_t nRows = maConstraints.getLength();
size_t nCompSize = nVariables * nRows;
std::unique_ptr<double[]> pCompMatrix(new double[nCompSize]); // first collect all coefficients, row-wise
for (size_t i=0; i<nCompSize; i++)
pCompMatrix[i] = 0.0;
std::unique_ptr<double[]> pRHS(new double[nRows]);
std::unique_ptr<char[]> pRowType(new char[nRows]);
for (size_t i=0; i<nRows; i++)
{
pRHS[i] = 0.0;
pRowType[i] = 'N';
}
for (sal_Int32 nConstrPos = 0; nConstrPos < maConstraints.getLength(); ++nConstrPos)
{
// integer constraints are set later
sheet::SolverConstraintOperator eOp = maConstraints[nConstrPos].Operator;
if ( eOp == sheet::SolverConstraintOperator_LESS_EQUAL ||
eOp == sheet::SolverConstraintOperator_GREATER_EQUAL ||
eOp == sheet::SolverConstraintOperator_EQUAL )
{
double fDirectValue = 0.0;
bool bRightCell = false;
table::CellAddress aRightAddr;
const uno::Any& rRightAny = maConstraints[nConstrPos].Right;
if ( rRightAny >>= aRightAddr )
bRightCell = true; // cell specified as right-hand side
else
rRightAny >>= fDirectValue; // constant value
table::CellAddress aLeftAddr = maConstraints[nConstrPos].Left;
const std::vector<double>& rLeftCoeff = aCellsHash[aLeftAddr];
double* pValues = &pCompMatrix[nConstrPos * nVariables];
for (nVar=0; nVar<nVariables; nVar++)
pValues[nVar] = rLeftCoeff[nVar+1];
// if left hand cell has a constant term, put into rhs value
double fRightValue = -rLeftCoeff[0];
if ( bRightCell )
{
const std::vector<double>& rRightCoeff = aCellsHash[aRightAddr];
// modify pValues with rhs coefficients
for (nVar=0; nVar<nVariables; nVar++)
pValues[nVar] -= rRightCoeff[nVar+1];
fRightValue += rRightCoeff[0]; // constant term
}
else
fRightValue += fDirectValue;
switch ( eOp )
{
case sheet::SolverConstraintOperator_LESS_EQUAL: pRowType[nConstrPos] = 'L'; break;
case sheet::SolverConstraintOperator_GREATER_EQUAL: pRowType[nConstrPos] = 'G'; break;
case sheet::SolverConstraintOperator_EQUAL: pRowType[nConstrPos] = 'E'; break;
default:
OSL_ENSURE( false, "unexpected enum type" );
}
pRHS[nConstrPos] = fRightValue;
}
}
// Find non-zero coefficients, column-wise
std::unique_ptr<int[]> pMatrixBegin(new int[nVariables+1]);
std::unique_ptr<int[]> pMatrixCount(new int[nVariables]);
std::unique_ptr<double[]> pMatrix(new double[nCompSize]); // not always completely used
std::unique_ptr<int[]> pMatrixIndex(new int[nCompSize]);
int nMatrixPos = 0;
for (nVar=0; nVar<nVariables; nVar++)
{
int nBegin = nMatrixPos;
for (size_t nRow=0; nRow<nRows; nRow++)
{
double fCoeff = pCompMatrix[ nRow * nVariables + nVar ]; // row-wise
if ( fCoeff != 0.0 )
{
pMatrix[nMatrixPos] = fCoeff;
pMatrixIndex[nMatrixPos] = nRow;
++nMatrixPos;
}
}
pMatrixBegin[nVar] = nBegin;
pMatrixCount[nVar] = nMatrixPos - nBegin;
}
pMatrixBegin[nVariables] = nMatrixPos;
pCompMatrix.reset();
// apply settings to all variables
std::unique_ptr<double[]> pLowerBounds(new double[nVariables]);
std::unique_ptr<double[]> pUpperBounds(new double[nVariables]);
for (nVar=0; nVar<nVariables; nVar++)
{
pLowerBounds[nVar] = mbNonNegative ? 0.0 : -DBL_MAX;
pUpperBounds[nVar] = DBL_MAX;
// bounds could possibly be further restricted from single-cell constraints
}
std::unique_ptr<char[]> pColType(new char[nVariables]);
for (nVar=0; nVar<nVariables; nVar++)
pColType[nVar] = mbInteger ? 'I' : 'C';
// apply single-var integer constraints
// variable bounds and column types. The default range is unbounded unless
// non-negativity was requested. Integer and binary single-variable
// constraints refine these.
std::vector<double> aLowerBounds(nVariables, mbNonNegative ? 0.0 : -DBL_MAX);
std::vector<double> aUpperBounds(nVariables, DBL_MAX);
std::vector<char> aColumnType(nVariables, mbInteger ? 'I' : 'C');
for (const auto& rConstr : maConstraints)
{
@@ -266,92 +148,44 @@ void SAL_CALL CoinMPSolver::solve()
eOp == sheet::SolverConstraintOperator_BINARY )
{
table::CellAddress aLeftAddr = rConstr.Left;
// find variable index for cell
for (nVar=0; nVar<nVariables; nVar++)
if ( AddressEqual( aVariableCells[nVar], aLeftAddr ) )
{
if ( eOp == sheet::SolverConstraintOperator_INTEGER )
pColType[nVar] = 'I';
aColumnType[nVar] = 'I';
else
{
pColType[nVar] = 'B';
pLowerBounds[nVar] = 0.0;
pUpperBounds[nVar] = 1.0;
aColumnType[nVar] = 'B';
aLowerBounds[nVar] = 0.0;
aUpperBounds[nVar] = 1.0;
}
}
}
}
int nObjectSense = mbMaximize ? SOLV_OBJSENS_MAX : SOLV_OBJSENS_MIN;
HPROB hProb = CoinCreateProblem("");
int nResult = CoinLoadProblem( hProb, nVariables, nRows, nMatrixPos, 0,
nObjectSense, nObjectConst, pObjectCoeffs.get(),
pLowerBounds.get(), pUpperBounds.get(), pRowType.get(), pRHS.get(), nullptr,
pMatrixBegin.get(), pMatrixCount.get(), pMatrixIndex.get(), pMatrix.get(),
nullptr, nullptr, nullptr );
if (nResult == SOLV_CALL_SUCCESS)
{
nResult = CoinLoadInteger( hProb, pColType.get() );
}
pColType.reset();
pMatrixIndex.reset();
pMatrix.reset();
pMatrixCount.reset();
pMatrixBegin.reset();
pUpperBounds.reset();
pLowerBounds.reset();
pRowType.reset();
pRHS.reset();
pObjectCoeffs.reset();
CoinSetRealOption( hProb, COIN_REAL_MAXSECONDS, mnTimeout );
CoinSetRealOption( hProb, COIN_REAL_MIPMAXSEC, mnTimeout );
// TODO: handle (or remove) settings: epsilon, B&B depth
// solve model
CoinMpSolveResult aSolveResult
= coinmpSolveLinearModel(maVariables, maConstraints, maObjective, aCellsHash, aLowerBounds,
aUpperBounds, aColumnType, mbMaximize, mnTimeout);
if (nResult == SOLV_CALL_SUCCESS)
{
nResult = CoinCheckProblem( hProb );
}
if (nResult == SOLV_CALL_SUCCESS)
{
try
{
nResult = CoinOptimizeProblem( hProb, 0 );
}
catch (const CoinError& e)
{
CoinUnloadProblem(hProb);
throw std::runtime_error(e.message());
}
}
mbSuccess = ( nResult == SOLV_CALL_SUCCESS );
mbSuccess = aSolveResult.bSuccess;
if ( mbSuccess )
{
// get solution
maSolution.realloc( nVariables );
CoinGetSolutionValues( hProb, maSolution.getArray(), nullptr, nullptr, nullptr );
mfResultValue = CoinGetObjectValue( hProb );
std::copy(aSolveResult.aSolution.begin(), aSolveResult.aSolution.end(),
maSolution.getArray());
mfResultValue = aSolveResult.fObjective;
}
else
{
int nSolutionStatus = CoinGetSolutionStatus( hProb );
if ( nSolutionStatus == 1 )
if (aSolveResult.nStatus == 1)
maStatus = SolverComponent::GetResourceString( RID_ERROR_INFEASIBLE );
else if ( nSolutionStatus == 2 )
else if (aSolveResult.nStatus == 2)
maStatus = SolverComponent::GetResourceString( RID_ERROR_UNBOUNDED );
// TODO: detect timeout condition and report as RID_ERROR_TIMEOUT
// (currently reported as infeasible)
}
CoinUnloadProblem( hProb );
}
extern "C" SAL_DLLPUBLIC_EXPORT css::uno::XInterface *