452 lines
17 KiB
C++
452 lines
17 KiB
C++
/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
|
|
/*************************************************************************
|
|
*
|
|
* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
|
|
*
|
|
* Copyright 2000, 2010 Oracle and/or its affiliates.
|
|
*
|
|
* OpenOffice.org - a multi-platform office productivity suite
|
|
*
|
|
* This file is part of OpenOffice.org.
|
|
*
|
|
* OpenOffice.org is free software: you can redistribute it and/or modify
|
|
* it under the terms of the GNU Lesser General Public License version 3
|
|
* only, as published by the Free Software Foundation.
|
|
*
|
|
* OpenOffice.org is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU Lesser General Public License version 3 for more details
|
|
* (a copy is included in the LICENSE file that accompanied this code).
|
|
*
|
|
* You should have received a copy of the GNU Lesser General Public License
|
|
* version 3 along with OpenOffice.org. If not, see
|
|
* <http://www.openoffice.org/license.html>
|
|
* for a copy of the LGPLv3 License.
|
|
*
|
|
* This file incorporates work covered by the following license notice:
|
|
*
|
|
* Licensed to the Apache Software Foundation (ASF) under one or more
|
|
* contributor license agreements. See the NOTICE file distributed
|
|
* with this work for additional information regarding copyright
|
|
* ownership. The ASF licenses this file to you under the Apache
|
|
* License, Version 2.0 (the "License"); you may not use this file
|
|
* except in compliance with the License. You may obtain a copy of
|
|
* the License at http://www.apache.org/licenses/LICENSE-2.0 .
|
|
*
|
|
************************************************************************/
|
|
|
|
#include <sal/config.h>
|
|
#include <sal/log.hxx>
|
|
|
|
#undef LANGUAGE_NONE
|
|
#if defined _WIN32
|
|
#define WINAPI __stdcall
|
|
#endif
|
|
#define LoadInverseLib FALSE
|
|
#define LoadLanguageLib FALSE
|
|
#ifdef SYSTEM_LPSOLVE
|
|
#include <lpsolve/lp_lib.h>
|
|
#else
|
|
#include <lp_lib.h>
|
|
#endif
|
|
#undef LANGUAGE_NONE
|
|
|
|
#include "SolverComponent.hxx"
|
|
#include <strings.hrc>
|
|
|
|
#include <com/sun/star/frame/XModel.hpp>
|
|
#include <com/sun/star/table/CellAddress.hpp>
|
|
#include <rtl/math.hxx>
|
|
#include <algorithm>
|
|
#include <memory>
|
|
#include <vector>
|
|
|
|
namespace com::sun::star::uno { class XComponentContext; }
|
|
|
|
using namespace com::sun::star;
|
|
|
|
namespace {
|
|
|
|
class LpsolveSolver : public SolverComponent
|
|
{
|
|
public:
|
|
LpsolveSolver() {}
|
|
|
|
private:
|
|
virtual void SAL_CALL solve() override;
|
|
virtual OUString SAL_CALL getImplementationName() override
|
|
{
|
|
return u"com.sun.star.comp.Calc.LpsolveSolver"_ustr;
|
|
}
|
|
virtual OUString SAL_CALL getComponentDescription() override
|
|
{
|
|
return SolverComponent::GetResourceString( RID_SOLVER_COMPONENT );
|
|
}
|
|
};
|
|
|
|
}
|
|
|
|
void SAL_CALL LpsolveSolver::solve()
|
|
{
|
|
uno::Reference<frame::XModel> xModel( mxDoc, uno::UNO_QUERY_THROW );
|
|
|
|
maStatus.clear();
|
|
mbSuccess = false;
|
|
|
|
if ( mnEpsilonLevel < EPS_TIGHT || mnEpsilonLevel > EPS_BAGGY )
|
|
{
|
|
maStatus = SolverComponent::GetResourceString( RID_ERROR_EPSILONLEVEL );
|
|
return;
|
|
}
|
|
|
|
xModel->lockControllers();
|
|
|
|
// collect variables in vector (?)
|
|
|
|
const auto & aVariableCells = maVariables;
|
|
size_t nVariables = aVariableCells.size();
|
|
size_t nVar = 0;
|
|
|
|
// Store all RHS values
|
|
sal_uInt32 nConstraints = maConstraints.size();
|
|
m_aConstrRHS.realloc(nConstraints);
|
|
|
|
// collect all dependent cells
|
|
|
|
ScSolverCellHashMap aCellsHash;
|
|
aCellsHash[maObjective].reserve( nVariables + 1 ); // objective function
|
|
|
|
for (const auto& rConstr : maConstraints)
|
|
{
|
|
table::CellAddress aCellAddr = rConstr.Left;
|
|
aCellsHash[aCellAddr].reserve( nVariables + 1 ); // constraints: left hand side
|
|
|
|
if ( rConstr.Right >>= aCellAddr )
|
|
aCellsHash[aCellAddr].reserve( nVariables + 1 ); // constraints: right hand side
|
|
}
|
|
|
|
// set all variables to zero
|
|
//! store old values?
|
|
//! use old values as initial values?
|
|
for ( const auto& rVarCell : aVariableCells )
|
|
{
|
|
SolverComponent::SetValue( mxDoc, rVarCell, 0.0 );
|
|
}
|
|
|
|
// read initial values from all dependent cells
|
|
for ( auto& rEntry : aCellsHash )
|
|
{
|
|
double fValue = SolverComponent::GetValue( mxDoc, rEntry.first );
|
|
rEntry.second.push_back( fValue ); // store as first element, as-is
|
|
}
|
|
|
|
// loop through variables
|
|
for ( const auto& rVarCell : aVariableCells )
|
|
{
|
|
SolverComponent::SetValue( mxDoc, rVarCell, 1.0 ); // set to 1 to examine influence
|
|
|
|
// read value change from all dependent cells
|
|
for ( auto& rEntry : aCellsHash )
|
|
{
|
|
double fChanged = SolverComponent::GetValue( mxDoc, rEntry.first );
|
|
double fInitial = rEntry.second.front();
|
|
rEntry.second.push_back( fChanged - fInitial );
|
|
}
|
|
|
|
SolverComponent::SetValue( mxDoc, rVarCell, 2.0 ); // minimal test for linearity
|
|
|
|
for ( const auto& rEntry : aCellsHash )
|
|
{
|
|
double fInitial = rEntry.second.front();
|
|
double fCoeff = rEntry.second.back(); // last appended: coefficient for this variable
|
|
double fTwo = SolverComponent::GetValue( mxDoc, rEntry.first );
|
|
|
|
bool bLinear = rtl::math::approxEqual( fTwo, fInitial + 2.0 * fCoeff ) ||
|
|
rtl::math::approxEqual( fInitial, fTwo - 2.0 * fCoeff );
|
|
// second comparison is needed in case fTwo is zero
|
|
if ( !bLinear )
|
|
maStatus = SolverComponent::GetResourceString( RID_ERROR_NONLINEAR );
|
|
}
|
|
|
|
SolverComponent::SetValue( mxDoc, rVarCell, 0.0 ); // set back to zero for examining next variable
|
|
}
|
|
|
|
xModel->unlockControllers();
|
|
|
|
if ( !maStatus.isEmpty() )
|
|
return;
|
|
|
|
|
|
// build lp_solve model
|
|
|
|
|
|
lprec* lp = make_lp( 0, nVariables );
|
|
if ( !lp )
|
|
return;
|
|
|
|
set_outputfile( lp, const_cast<char*>( "" ) ); // no output
|
|
|
|
// set objective function
|
|
|
|
const std::vector<double>& rObjCoeff = aCellsHash[maObjective];
|
|
std::unique_ptr<REAL[]> pObjVal(new REAL[nVariables+1]);
|
|
pObjVal[0] = 0.0; // ignored
|
|
for (nVar=0; nVar<nVariables; nVar++)
|
|
pObjVal[nVar+1] = rObjCoeff[nVar+1];
|
|
set_obj_fn( lp, pObjVal.get() );
|
|
pObjVal.reset();
|
|
set_rh( lp, 0, rObjCoeff[0] ); // constant term of objective
|
|
|
|
// add rows
|
|
|
|
set_add_rowmode(lp, TRUE);
|
|
|
|
sal_uInt32 nConstrCount(0);
|
|
double* pConstrRHS = m_aConstrRHS.getArray();
|
|
|
|
for (const auto& rConstr : maConstraints)
|
|
{
|
|
// integer constraints are set later
|
|
sheet::SolverConstraintOperator eOp = rConstr.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 = rConstr.Right;
|
|
if ( rRightAny >>= aRightAddr )
|
|
bRightCell = true; // cell specified as right-hand side
|
|
else
|
|
rRightAny >>= fDirectValue; // constant value
|
|
|
|
table::CellAddress aLeftAddr = rConstr.Left;
|
|
|
|
const std::vector<double>& rLeftCoeff = aCellsHash[aLeftAddr];
|
|
std::unique_ptr<REAL[]> pValues(new REAL[nVariables+1] );
|
|
pValues[0] = 0.0; // ignored?
|
|
for (nVar=0; nVar<nVariables; nVar++)
|
|
pValues[nVar+1] = 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+1] -= rRightCoeff[nVar+1];
|
|
|
|
fRightValue += rRightCoeff[0]; // constant term
|
|
}
|
|
else
|
|
fRightValue += fDirectValue;
|
|
|
|
// Remember the RHS value used for sensitivity analysis later
|
|
pConstrRHS[nConstrCount] = fRightValue;
|
|
|
|
int nConstrType = LE;
|
|
switch ( eOp )
|
|
{
|
|
case sheet::SolverConstraintOperator_LESS_EQUAL: nConstrType = LE; break;
|
|
case sheet::SolverConstraintOperator_GREATER_EQUAL: nConstrType = GE; break;
|
|
case sheet::SolverConstraintOperator_EQUAL: nConstrType = EQ; break;
|
|
default:
|
|
OSL_FAIL( "unexpected enum type" );
|
|
}
|
|
add_constraint( lp, pValues.get(), nConstrType, fRightValue );
|
|
nConstrCount++;
|
|
}
|
|
}
|
|
|
|
set_add_rowmode(lp, FALSE);
|
|
|
|
// apply settings to all variables
|
|
|
|
for (nVar=0; nVar<nVariables; nVar++)
|
|
{
|
|
if ( !mbNonNegative )
|
|
set_unbounded(lp, nVar+1); // allow negative (default is non-negative)
|
|
//! collect bounds from constraints?
|
|
if ( mbInteger )
|
|
set_int(lp, nVar+1, TRUE);
|
|
}
|
|
|
|
// apply single-var integer constraints
|
|
|
|
for (const auto& rConstr : maConstraints)
|
|
{
|
|
sheet::SolverConstraintOperator eOp = rConstr.Operator;
|
|
if ( eOp == sheet::SolverConstraintOperator_INTEGER ||
|
|
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 )
|
|
set_int(lp, nVar+1, TRUE);
|
|
else
|
|
set_binary(lp, nVar+1, TRUE);
|
|
}
|
|
}
|
|
}
|
|
|
|
if ( mbMaximize )
|
|
set_maxim(lp);
|
|
else
|
|
set_minim(lp);
|
|
|
|
if ( !mbLimitBBDepth )
|
|
set_bb_depthlimit( lp, 0 );
|
|
|
|
set_epslevel( lp, mnEpsilonLevel );
|
|
set_timeout( lp, mnTimeout );
|
|
|
|
// solve model
|
|
|
|
int nResult = ::solve( lp );
|
|
|
|
mbSuccess = ( nResult == OPTIMAL );
|
|
if ( mbSuccess )
|
|
{
|
|
// get solution
|
|
|
|
maSolution.realloc( nVariables );
|
|
|
|
REAL* pResultVar = nullptr;
|
|
get_ptr_variables( lp, &pResultVar );
|
|
std::copy_n(pResultVar, nVariables, maSolution.getArray());
|
|
|
|
mfResultValue = get_objective( lp );
|
|
|
|
// Initially set to false because getting the report might fail
|
|
m_aSensitivityReport.HasReport = false;
|
|
|
|
// Get sensitivity report if the user set SensitivityReport parameter to true
|
|
if (mbGenSensitivity)
|
|
{
|
|
// Get sensitivity data about the objective function
|
|
// LpSolve returns an interval for the coefficients of the objective function
|
|
// instead of returning an allowable increase/decrease (which is what we want to show
|
|
// in the sensitivity report; so we these from/till values are converted into increase
|
|
// and decrease values later)
|
|
REAL* pObjFrom = nullptr;
|
|
REAL* pObjTill = nullptr;
|
|
bool bHasObjReport = false;
|
|
bHasObjReport = get_ptr_sensitivity_obj(lp, &pObjFrom, &pObjTill);
|
|
|
|
// Get sensitivity data about constraints
|
|
// Similarly to the objective function, the sensitivity values returned for the
|
|
// constraints are in the form from/till and are later converted to increase and
|
|
// decrease values later
|
|
REAL* pConstrValue = nullptr;
|
|
REAL* pConstrDual = nullptr;
|
|
REAL* pConstrFrom = nullptr;
|
|
REAL* pConstrTill = nullptr;
|
|
bool bHasConstrReport = false;
|
|
bHasConstrReport = get_ptr_sensitivity_rhs(lp, &pConstrDual, &pConstrFrom, &pConstrTill);
|
|
|
|
// When successful, store sensitivity data in the solver component
|
|
if (bHasObjReport && bHasConstrReport)
|
|
{
|
|
m_aSensitivityReport.HasReport = true;
|
|
m_aObjDecrease.realloc(nVariables);
|
|
m_aObjIncrease.realloc(nVariables);
|
|
double* pObjDecrease = m_aObjDecrease.getArray();
|
|
double* pObjIncrease = m_aObjIncrease.getArray();
|
|
for (size_t i = 0; i < nVariables; i++)
|
|
{
|
|
// Allowed decrease. Note that the indices of rObjCoeff are offset by 1
|
|
// because of the objective function
|
|
if (static_cast<bool>(is_infinite(lp, pObjFrom[i])))
|
|
pObjDecrease[i] = get_infinite(lp);
|
|
else
|
|
pObjDecrease[i] = rObjCoeff[i + 1] - pObjFrom[i];
|
|
|
|
// Allowed increase
|
|
if (static_cast<bool>(is_infinite(lp, pObjTill[i])))
|
|
pObjIncrease[i] = get_infinite(lp);
|
|
else
|
|
pObjIncrease[i] = pObjTill[i] - rObjCoeff[i + 1];
|
|
}
|
|
|
|
// Save objective coefficients for the sensitivity report
|
|
double* pObjCoefficients(new double[nVariables]);
|
|
for (size_t i = 0; i < nVariables; i++)
|
|
pObjCoefficients[i] = rObjCoeff[i + 1];
|
|
m_aObjCoefficients.realloc(nVariables);
|
|
std::copy_n(pObjCoefficients, nVariables, m_aObjCoefficients.getArray());
|
|
|
|
// The reduced costs are in pConstrDual after the constraints
|
|
double* pObjRedCost(new double[nVariables]);
|
|
for (size_t i = 0; i < nVariables; i++)
|
|
pObjRedCost[i] = pConstrDual[nConstraints + i];
|
|
m_aObjRedCost.realloc(nVariables);
|
|
std::copy_n(pObjRedCost, nVariables, m_aObjRedCost.getArray());
|
|
|
|
// Final value of constraints
|
|
get_ptr_constraints(lp, &pConstrValue);
|
|
m_aConstrValue.realloc(nConstraints);
|
|
std::copy_n(pConstrValue, nConstraints, m_aConstrValue.getArray());
|
|
|
|
// The RHS contains information for each constraint
|
|
m_aConstrDual.realloc(nConstraints);
|
|
m_aConstrDecrease.realloc(nConstraints);
|
|
m_aConstrIncrease.realloc(nConstraints);
|
|
std::copy_n(pConstrDual, nConstraints, m_aConstrDual.getArray());
|
|
double* pConstrDecrease = m_aConstrDecrease.getArray();
|
|
double* pConstrIncrease = m_aConstrIncrease.getArray();
|
|
|
|
for (sal_uInt32 i = 0; i < nConstraints; i++)
|
|
{
|
|
// Allowed decrease
|
|
pConstrDecrease[i] = m_aConstrRHS[i] - pConstrFrom[i];
|
|
if (static_cast<bool>(is_infinite(lp, pConstrFrom[i]))
|
|
&& maConstraints[i].Operator == sheet::SolverConstraintOperator_LESS_EQUAL)
|
|
pConstrDecrease[i] = m_aConstrRHS[i] - m_aConstrValue[i];
|
|
|
|
// Allowed increase
|
|
pConstrIncrease[i] = pConstrTill[i] - m_aConstrRHS[i];
|
|
if (static_cast<bool>(is_infinite(lp, pConstrTill[i]))
|
|
&& maConstraints[i].Operator == sheet::SolverConstraintOperator_GREATER_EQUAL)
|
|
pConstrIncrease[i] = m_aConstrValue[i] - m_aConstrRHS[i];
|
|
}
|
|
|
|
// Set all values of the SensitivityReport object
|
|
m_aSensitivityReport.ObjCoefficients = m_aObjCoefficients;
|
|
m_aSensitivityReport.ObjReducedCosts = m_aObjRedCost;
|
|
m_aSensitivityReport.ObjAllowableDecreases = m_aObjDecrease;
|
|
m_aSensitivityReport.ObjAllowableIncreases = m_aObjIncrease;
|
|
m_aSensitivityReport.ConstrValues = m_aConstrValue;
|
|
m_aSensitivityReport.ConstrRHS = m_aConstrRHS;
|
|
m_aSensitivityReport.ConstrShadowPrices = m_aConstrDual;
|
|
m_aSensitivityReport.ConstrAllowableDecreases = m_aConstrDecrease;
|
|
m_aSensitivityReport.ConstrAllowableIncreases = m_aConstrIncrease;
|
|
}
|
|
}
|
|
}
|
|
else if ( nResult == INFEASIBLE )
|
|
maStatus = SolverComponent::GetResourceString( RID_ERROR_INFEASIBLE );
|
|
else if ( nResult == UNBOUNDED )
|
|
maStatus = SolverComponent::GetResourceString( RID_ERROR_UNBOUNDED );
|
|
else if ( nResult == TIMEOUT || nResult == SUBOPTIMAL )
|
|
maStatus = SolverComponent::GetResourceString( RID_ERROR_TIMEOUT );
|
|
// SUBOPTIMAL is assumed to be caused by a timeout, and reported as an error
|
|
|
|
delete_lp( lp );
|
|
}
|
|
|
|
extern "C" SAL_DLLPUBLIC_EXPORT css::uno::XInterface *
|
|
com_sun_star_comp_Calc_LpsolveSolver_get_implementation(
|
|
css::uno::XComponentContext *,
|
|
css::uno::Sequence<css::uno::Any> const &)
|
|
{
|
|
return cppu::acquire(new LpsolveSolver());
|
|
}
|
|
|
|
/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
|