tdf#157519 Implement Sensitivity Report in LpSolve solver

This patch implements sensitivity analysis when using the LpSolve solver engine.

It also adds the infrastructure needed for future implementations in other solver engines via the css::sheet::SensitivityReport struct.

Change-Id: I74c2ed9c6201a0b2ffc29ef612d2b778d11a3bef
Reviewed-on: https://gerrit.libreoffice.org/c/core/+/173642
Tested-by: Jenkins
Reviewed-by: Mike Kaganski <mike.kaganski@collabora.com>

3 files changed, 150 insertions, 6 deletions

diff --git a/sccomp/source/solver/LpsolveSolver.cxx b/sccomp/source/solver/LpsolveSolver.cxx
index e80bd87e1a41..70a9239c5083 100644
--- a/sccomp/source/solver/LpsolveSolver.cxx
+++ b/sccomp/source/solver/LpsolveSolver.cxx
@@ -37,6 +37,7 @@
  ************************************************************************/
 
 #include <sal/config.h>
+#include <sal/log.hxx>
 
 #undef LANGUAGE_NONE
 #if defined _WIN32
@@ -107,6 +108,10 @@ void SAL_CALL LpsolveSolver::solve()
     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;
@@ -197,6 +202,9 @@ void SAL_CALL LpsolveSolver::solve()
 
     set_add_rowmode(lp, TRUE);
 
+    sal_uInt32 nConstrCount(0);
+    double* pConstrRHS = m_aConstrRHS.getArray();
+
     for (const auto& rConstr : maConstraints)
     {
         // integer constraints are set later
@@ -237,6 +245,9 @@ void SAL_CALL LpsolveSolver::solve()
             else
                 fRightValue += fDirectValue;
 
+            // Remember the RHS value used for sensitivity analysis later
+            pConstrRHS[nConstrCount] = fRightValue;
+
             int nConstrType = LE;
             switch ( eOp )
             {
@@ -247,6 +258,7 @@ void SAL_CALL LpsolveSolver::solve()
                     OSL_FAIL( "unexpected enum type" );
             }
             add_constraint( lp, pValues.get(), nConstrType, fRightValue );
+            nConstrCount++;
         }
     }
 
@@ -311,6 +323,112 @@ void SAL_CALL LpsolveSolver::solve()
         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 successfull, 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 );
diff --git a/sccomp/source/solver/SolverComponent.cxx b/sccomp/source/solver/SolverComponent.cxx
index b7038090e56c..a06b65f4c46a 100644
--- a/sccomp/source/solver/SolverComponent.cxx
+++ b/sccomp/source/solver/SolverComponent.cxx
@@ -37,6 +37,8 @@ constexpr OUStringLiteral STR_INTEGER = u"Integer";
 constexpr OUStringLiteral STR_TIMEOUT = u"Timeout";
 constexpr OUStringLiteral STR_EPSILONLEVEL = u"EpsilonLevel";
 constexpr OUStringLiteral STR_LIMITBBDEPTH = u"LimitBBDepth";
+constexpr OUStringLiteral STR_GEN_SENSITIVITY = u"GenSensitivityReport";
+constexpr OUStringLiteral STR_SENSITIVITY_REPORT = u"SensitivityReport";
 
 
 //  Resources from tools are used for translated strings
@@ -64,7 +66,9 @@ namespace
         PROP_INTEGER,
         PROP_TIMEOUT,
         PROP_EPSILONLEVEL,
-        PROP_LIMITBBDEPTH
+        PROP_LIMITBBDEPTH,
+        PROP_GEN_SENSITIVITY,
+        PROP_SENSITIVITY_REPORT
     };
 }
 
@@ -95,15 +99,20 @@ SolverComponent::SolverComponent() :
     mnTimeout( 100 ),
     mnEpsilonLevel( 0 ),
     mbLimitBBDepth( true ),
+    mbGenSensitivity(false),
     mbSuccess( false ),
     mfResultValue( 0.0 )
 {
     // for XPropertySet implementation:
-    registerProperty( STR_NONNEGATIVE,  PROP_NONNEGATIVE,  0, &mbNonNegative,  cppu::UnoType<decltype(mbNonNegative)>::get()  );
-    registerProperty( STR_INTEGER,      PROP_INTEGER,      0, &mbInteger,      cppu::UnoType<decltype(mbInteger)>::get()      );
-    registerProperty( STR_TIMEOUT,      PROP_TIMEOUT,      0, &mnTimeout,      cppu::UnoType<decltype(mnTimeout)>::get()      );
-    registerProperty( STR_EPSILONLEVEL, PROP_EPSILONLEVEL, 0, &mnEpsilonLevel, cppu::UnoType<decltype(mnEpsilonLevel)>::get() );
-    registerProperty( STR_LIMITBBDEPTH, PROP_LIMITBBDEPTH, 0, &mbLimitBBDepth, cppu::UnoType<decltype(mbLimitBBDepth)>::get() );
+    registerProperty(STR_NONNEGATIVE,  PROP_NONNEGATIVE,  0, &mbNonNegative,    cppu::UnoType<decltype(mbNonNegative)>::get());
+    registerProperty(STR_INTEGER,      PROP_INTEGER,      0, &mbInteger,        cppu::UnoType<decltype(mbInteger)>::get());
+    registerProperty(STR_TIMEOUT,      PROP_TIMEOUT,      0, &mnTimeout,        cppu::UnoType<decltype(mnTimeout)>::get());
+    registerProperty(STR_EPSILONLEVEL, PROP_EPSILONLEVEL, 0, &mnEpsilonLevel,   cppu::UnoType<decltype(mnEpsilonLevel)>::get());
+    registerProperty(STR_LIMITBBDEPTH, PROP_LIMITBBDEPTH, 0, &mbLimitBBDepth,   cppu::UnoType<decltype(mbLimitBBDepth)>::get());
+    registerProperty(STR_GEN_SENSITIVITY, PROP_GEN_SENSITIVITY, 0, &mbGenSensitivity, cppu::UnoType<decltype(mbGenSensitivity)>::get());
+
+    // Sensitivity report
+    registerProperty(STR_SENSITIVITY_REPORT, PROP_SENSITIVITY_REPORT, 0, &m_aSensitivityReport, cppu::UnoType<decltype(m_aSensitivityReport)>::get());
 }
 
 SolverComponent::~SolverComponent()
@@ -158,6 +167,9 @@ OUString SAL_CALL SolverComponent::getPropertyDescription( const OUString& rProp
         case PROP_LIMITBBDEPTH:
             pResId = RID_PROPERTY_LIMITBBDEPTH;
             break;
+        case PROP_GEN_SENSITIVITY:
+            pResId = RID_PROPERTY_SENSITIVITY;
+            break;
         default:
             {
                 // unknown - leave empty
diff --git a/sccomp/source/solver/SolverComponent.hxx b/sccomp/source/solver/SolverComponent.hxx
index 7b5ff1dd49f0..543eeedea282 100644
--- a/sccomp/source/solver/SolverComponent.hxx
+++ b/sccomp/source/solver/SolverComponent.hxx
@@ -21,6 +21,7 @@
 
 #include <com/sun/star/sheet/XSolver.hpp>
 #include <com/sun/star/sheet/XSolverDescription.hpp>
+#include <com/sun/star/sheet/SensitivityReport.hpp>
 #include <com/sun/star/table/CellAddress.hpp>
 #include <com/sun/star/lang/XServiceInfo.hpp>
 #include <cppuhelper/implbase.hxx>
@@ -76,12 +77,25 @@ protected:
     sal_Int32                                               mnTimeout;
     sal_Int32                                               mnEpsilonLevel;
     bool                                                    mbLimitBBDepth;
+    bool                                                    mbGenSensitivity;
     // results
     bool                                                    mbSuccess;
     double                                                  mfResultValue;
     css::uno::Sequence< double >                            maSolution;
     OUString                                                maStatus;
 
+    // Sensitivity report
+    css::uno::Sequence<double> m_aObjCoefficients;
+    css::uno::Sequence<double> m_aObjDecrease;
+    css::uno::Sequence<double> m_aObjIncrease;
+    css::uno::Sequence<double> m_aObjRedCost;
+    css::uno::Sequence<double> m_aConstrValue;
+    css::uno::Sequence<double> m_aConstrRHS;
+    css::uno::Sequence<double> m_aConstrDual;
+    css::uno::Sequence<double> m_aConstrIncrease;
+    css::uno::Sequence<double> m_aConstrDecrease;
+    css::sheet::SensitivityReport m_aSensitivityReport;
+
     static OUString GetResourceString(TranslateId aId);
     static css::uno::Reference<css::table::XCell> GetCell(
             const css::uno::Reference<css::sheet::XSpreadsheetDocument>& xDoc,