Remove computation for number of non-zero terms in the constraint mat…

…rix [ANT-2258] (#2496)

The goal of this PR is to avoid having to compute the number of non-zero
terms in the constraint matrix.

Some C functions require a `T*` array. Previously, we used
`std::vector<T>::data` for `CoefficientsDeLaMatriceDesContraintes` and
`IndicesColonnes`.

Instead, we use a 2-step mechanism through templated class `VectorMap`
1. Fill a `std::vector<std::pair<int, T>>` (vector of pairs)
2. Extract the vector of pairs into an `std::vector`, use it's `data`.
This `std::vector` is kept alive as long as it's `data` is used,
preventing invalid memory reads.

src/libs/antares/InfoCollection/StudyInfoCollector.cpp

@@ -167,9 +167,6 @@ void SimulationInfoCollector::toFileContent(FileContent& file_content)
 {
     file_content.addItemToSection("optimization problem", "variables", opt_info_.nbVariables);
     file_content.addItemToSection("optimization problem", "constraints", opt_info_.nbConstraints);
-    file_content.addItemToSection("optimization problem",
-                                  "non-zero coefficients",
-                                  opt_info_.nbNonZeroCoeffs);
 }
 
 } // namespace Benchmarking

src/libs/antares/InfoCollection/include/antares/infoCollection/StudyInfoCollector.h

@@ -62,7 +62,6 @@ struct OptimizationInfo
 {
     unsigned int nbVariables = 0;
     unsigned int nbConstraints = 0;
-    unsigned int nbNonZeroCoeffs = 0;
 };
 
 class SimulationInfoCollector

src/solver/optimisation/CMakeLists.txt

@@ -5,6 +5,7 @@ set(RTESOLVER_OPT
         opt_optimisation_lineaire.cpp
         opt_chainage_intercos.cpp
         include/antares/solver/optimisation/opt_fonctions.h
+        include/antares/solver/optimisation/SparseVector.hxx
         opt_pilotage_optimisation_lineaire.cpp
         opt_pilotage_optimisation_quadratique.cpp
         include/antares/solver/optimisation/opt_structure_probleme_a_resoudre.h

src/solver/optimisation/HebdoProblemToLpsTranslator.cpp

@@ -43,6 +43,13 @@ void copy(const T& in, U& out)
 {
     std::ranges::copy(in, std::back_inserter(out));
 }
+
+template<class T, class U>
+void copy(const SparseVector<T>& in, U& out)
+{
+    copy(in.extract(), out);
+}
+
 } // namespace
 
 WeeklyDataFromAntares HebdoProblemToLpsTranslator::translate(

src/solver/optimisation/adequacy_patch_csr/adq_patch_curtailment_sharing.cpp

@@ -141,13 +141,9 @@ void HourlyCSRProblem::calculateCsrParameters()
 void HourlyCSRProblem::allocateProblem()
 {
     using namespace Antares::Data::AdequacyPatch;
-    int nbConst;
-
     problemeAResoudre_.NombreDeVariables = countVariables(problemeHebdo_);
-    nbConst = problemeAResoudre_.NombreDeContraintes = countConstraints(problemeHebdo_);
-    int nbTerms = 3 * nbConst; // This is a rough estimate, reallocations may happen later if it's
-                               // too low
-    OPT_AllocateFromNumberOfVariableConstraints(&problemeAResoudre_, nbTerms);
+    problemeAResoudre_.NombreDeContraintes = countConstraints(problemeHebdo_);
+    OPT_AllocateFromNumberOfVariableConstraints(&problemeAResoudre_);
 }
 
 void HourlyCSRProblem::buildProblemVariables()

src/solver/optimisation/constraints/ConsistenceNumberOfDispatchableUnits.cpp

@@ -60,7 +60,6 @@ void ConsistenceNumberOfDispatchableUnits::add(int pays, int index, int pdt)
     }
     else
     {
-        builder.data.NbTermesContraintesPourLesCoutsDeDemarrage += 4;
         builder.data.nombreDeContraintes++;
     }
 }

src/solver/optimisation/constraints/ConstraintBuilder.cpp

@@ -188,11 +188,6 @@ void ConstraintBuilder::OPT_ChargerLaContrainteDansLaMatriceDesContraintes()
           = data.Pi[i];
         data.IndicesColonnes[data.nombreDeTermesDansLaMatriceDeContrainte] = data.Colonne[i];
         data.nombreDeTermesDansLaMatriceDeContrainte++;
-        if (data.nombreDeTermesDansLaMatriceDeContrainte
-            == data.NombreDeTermesAllouesDansLaMatriceDesContraintes)
-        {
-            OPT_AugmenterLaTailleDeLaMatriceDesContraintes();
-        }
     }
     data.NombreDeTermesDesLignes[data.nombreDeContraintes] = nombreDeTermes_;
 
@@ -201,20 +196,3 @@ void ConstraintBuilder::OPT_ChargerLaContrainteDansLaMatriceDesContraintes()
 
     return;
 }
-
-void ConstraintBuilder::OPT_AugmenterLaTailleDeLaMatriceDesContraintes()
-{
-    int NbTermes = data.NombreDeTermesAllouesDansLaMatriceDesContraintes;
-    NbTermes += data.IncrementDAllocationMatriceDesContraintes;
-
-    logs.info();
-    logs.info() << " Expected Number of Non-zero terms in Problem Matrix : increased to : "
-                << NbTermes;
-    logs.info();
-
-    data.CoefficientsDeLaMatriceDesContraintes.resize(NbTermes);
-
-    data.IndicesColonnes.resize(NbTermes);
-
-    data.NombreDeTermesAllouesDansLaMatriceDesContraintes = NbTermes;
-}

src/solver/optimisation/constraints/MinDownTime.cpp

@@ -66,8 +66,6 @@ void MinDownTime::add(int pays, int index, int pdt)
     }
     else
     {
-        builder.data.NbTermesContraintesPourLesCoutsDeDemarrage
-          += 1 + DureeMinimaleDArretDUnGroupeDuPalierThermique;
         builder.data.nombreDeContraintes++;
     }
 }

src/solver/optimisation/constraints/NbDispUnitsMinBoundSinceMinUpTime.cpp

@@ -70,8 +70,6 @@ void NbDispUnitsMinBoundSinceMinUpTime::add(int pays, int index, int pdt)
     }
     else
     {
-        builder.data.NbTermesContraintesPourLesCoutsDeDemarrage
-          += 1 + 2 * DureeMinimaleDeMarcheDUnGroupeDuPalierThermique;
         builder.data.nombreDeContraintes++;
     }
 }

src/solver/optimisation/constraints/NbUnitsOutageLessThanNbUnitsStop.cpp

@@ -50,7 +50,6 @@ void NbUnitsOutageLessThanNbUnitsStop::add(int pays, int index, int pdt)
     }
     else
     {
-        builder.data.NbTermesContraintesPourLesCoutsDeDemarrage += 2;
         builder.data.nombreDeContraintes++;
     }
 }

src/solver/optimisation/constraints/PMaxDispatchableGeneration.cpp

@@ -48,7 +48,6 @@ void PMaxDispatchableGeneration::add(int pays, int index, int pdt)
     }
     else
     {
-        builder.data.NbTermesContraintesPourLesCoutsDeDemarrage += 2;
         builder.data.nombreDeContraintes++;
     }
 }

src/solver/optimisation/constraints/PMinDispatchableGeneration.cpp

@@ -49,7 +49,6 @@ void PMinDispatchableGeneration::add(int pays, int index, int pdt)
     }
     else
     {
-        builder.data.NbTermesContraintesPourLesCoutsDeDemarrage += 2;
         builder.data.nombreDeContraintes++;
     }
 }

src/solver/optimisation/constraints/constraint_builder_utils.cpp

@@ -32,7 +32,6 @@ ConstraintBuilderData NewGetConstraintBuilderFromProblemHebdoAndProblemAResoudre
             ProblemeAResoudre.IndicesDebutDeLigne,
             ProblemeAResoudre.CoefficientsDeLaMatriceDesContraintes,
             ProblemeAResoudre.IndicesColonnes,
-            ProblemeAResoudre.NombreDeTermesAllouesDansLaMatriceDesContraintes,
             ProblemeAResoudre.NombreDeTermesDesLignes,
             ProblemeAResoudre.Sens,
             ProblemeAResoudre.IncrementDAllocationMatriceDesContraintes,
@@ -44,6 +43,5 @@ ConstraintBuilderData NewGetConstraintBuilderFromProblemHebdoAndProblemAResoudre
             problemeHebdo->NamedProblems,
             problemeHebdo->NomsDesPays,
             problemeHebdo->weekInTheYear,
-            problemeHebdo->NombreDePasDeTemps,
-            problemeHebdo->NbTermesContraintesPourLesCoutsDeDemarrage};
+            problemeHebdo->NombreDePasDeTemps};
 }

src/solver/optimisation/include/antares/solver/optimisation/SparseVector.hxx

@@ -0,0 +1,44 @@
+#include <vector>
+
+template<class T>
+class SparseVector
+{
+public:
+    std::size_t size() const
+    {
+        return isExtracted ? v.size() : m.size();
+    }
+
+    T& operator[](std::size_t idx)
+    {
+        isExtracted = false;
+        m.push_back({idx, 0});
+        return m.back().second;
+    }
+
+    T* data() const
+    {
+        extract();
+        return v.data();
+    }
+
+    std::vector<T>& extract() const
+    {
+        if (!isExtracted)
+        {
+            isExtracted = true;
+            v.assign(m.size(), 0.);
+            for (auto [index, coeff]: m)
+            {
+                v[index] = coeff;
+            }
+            m.clear();
+        }
+        return v;
+    }
+
+private:
+    mutable bool isExtracted = false;
+    mutable std::vector<std::pair<int, T>> m;
+    mutable std::vector<T> v;
+};

src/solver/optimisation/include/antares/solver/optimisation/constraints/ConstraintBuilder.h

@@ -39,9 +39,8 @@ class ConstraintBuilderData
     int& nombreDeContraintes;
     int& nombreDeTermesDansLaMatriceDeContrainte;
     std::vector<int>& IndicesDebutDeLigne;
-    std::vector<double>& CoefficientsDeLaMatriceDesContraintes;
-    std::vector<int>& IndicesColonnes;
-    int& NombreDeTermesAllouesDansLaMatriceDesContraintes; // TODO Check if ref is needed
+    SparseVector<double>& CoefficientsDeLaMatriceDesContraintes;
+    SparseVector<int>& IndicesColonnes;
     std::vector<int>& NombreDeTermesDesLignes;
     std::string& Sens;
     int& IncrementDAllocationMatriceDesContraintes;
@@ -54,7 +53,6 @@ class ConstraintBuilderData
     const std::vector<const char*>& NomsDesPays;
     const uint32_t& weekInTheYear;
     const uint32_t& NombreDePasDeTemps;
-    uint32_t& NbTermesContraintesPourLesCoutsDeDemarrage;
 };
 
 /*! \verbatim
@@ -221,8 +219,6 @@ class ConstraintBuilder
 private:
     void OPT_ChargerLaContrainteDansLaMatriceDesContraintes();
 
-    void OPT_AugmenterLaTailleDeLaMatriceDesContraintes();
-
     unsigned int hourInWeek_ = 0;
 
     char operator_ = '=';

src/solver/optimisation/include/antares/solver/optimisation/opt_fonctions.h

@@ -95,8 +95,7 @@ double OPT_CalculerAireMaxPminJour(int, int, int, int, std::vector<int>&, std::v
 
 void OPT_ChainagesDesIntercoPartantDUnNoeud(PROBLEME_HEBDO*);
 
-void OPT_AllocateFromNumberOfVariableConstraints(PROBLEME_ANTARES_A_RESOUDRE* ProblemeAResoudre,
-                                                 int);
+void OPT_AllocateFromNumberOfVariableConstraints(PROBLEME_ANTARES_A_RESOUDRE* ProblemeAResoudre);
 void OPT_AllocDuProblemeAOptimiser(PROBLEME_HEBDO*);
 int OPT_DecompteDesVariablesEtDesContraintesDuProblemeAOptimiser(PROBLEME_HEBDO*);
 

src/solver/optimisation/include/antares/solver/optimisation/opt_structure_probleme_a_resoudre.h

@@ -27,6 +27,7 @@
 
 #include <antares/solver/utils/basis_status.h>
 
+#include "SparseVector.hxx"
 #include "opt_constants.h"
 
 /*--------------------------------------------------------------------------------------*/
@@ -45,9 +46,8 @@ struct PROBLEME_ANTARES_A_RESOUDRE
     std::string Sens;
     std::vector<int> IndicesDebutDeLigne;
     std::vector<int> NombreDeTermesDesLignes;
-    std::vector<double> CoefficientsDeLaMatriceDesContraintes;
-    std::vector<int> IndicesColonnes;
-    int NombreDeTermesAllouesDansLaMatriceDesContraintes;
+    SparseVector<double> CoefficientsDeLaMatriceDesContraintes;
+    SparseVector<int> IndicesColonnes;
     int IncrementDAllocationMatriceDesContraintes;
     int NombreDeTermesDansLaMatriceDesContraintes;
     /* Donnees variables de la matrice des contraintes */

src/solver/optimisation/opt_alloc_probleme_a_optimiser.cpp

@@ -20,13 +20,11 @@
 */
 
 #include <antares/logs/logs.h>
+#include "antares/solver/optimisation/opt_fonctions.h"
 #include "antares/solver/optimisation/opt_structure_probleme_a_resoudre.h"
 #include "antares/solver/simulation/sim_structure_probleme_economique.h"
 
-int OPT_DecompteDesVariablesEtDesContraintesDuProblemeAOptimiser(PROBLEME_HEBDO*);
-
-void OPT_AllocateFromNumberOfVariableConstraints(PROBLEME_ANTARES_A_RESOUDRE* ProblemeAResoudre,
-                                                 int NbTermes)
+void OPT_AllocateFromNumberOfVariableConstraints(PROBLEME_ANTARES_A_RESOUDRE* ProblemeAResoudre)
 {
     const size_t nbVariables = ProblemeAResoudre->NombreDeVariables;
     const size_t nbConstraints = ProblemeAResoudre->NombreDeContraintes;
@@ -35,12 +33,6 @@ void OPT_AllocateFromNumberOfVariableConstraints(PROBLEME_ANTARES_A_RESOUDRE* Pr
     ProblemeAResoudre->IndicesDebutDeLigne.assign(nbConstraints, 0);
     ProblemeAResoudre->NombreDeTermesDesLignes.assign(nbConstraints, 0);
 
-    ProblemeAResoudre->CoefficientsDeLaMatriceDesContraintes.assign(NbTermes, 0.);
-    ProblemeAResoudre->IndicesColonnes.assign(NbTermes, 0);
-
-    ProblemeAResoudre->NombreDeTermesAllouesDansLaMatriceDesContraintes = NbTermes;
-    ProblemeAResoudre->IncrementDAllocationMatriceDesContraintes = (int)(0.1 * NbTermes);
-
     ProblemeAResoudre->CoutQuadratique.assign(nbVariables, 0.);
     ProblemeAResoudre->CoutLineaire.assign(nbVariables, 0.);
     ProblemeAResoudre->TypeDeVariable.assign(nbVariables, 0);
@@ -70,62 +62,20 @@ void OPT_AllocateFromNumberOfVariableConstraints(PROBLEME_ANTARES_A_RESOUDRE* Pr
     ProblemeAResoudre->VariablesEntieres.resize(nbVariables);
 }
 
-static void optimisationAllocateProblem(PROBLEME_HEBDO* problemeHebdo, const int mxPaliers)
+static void optimisationAllocateProblem(PROBLEME_HEBDO* problemeHebdo)
 {
     const auto& ProblemeAResoudre = problemeHebdo->ProblemeAResoudre;
 
     int NombreDePasDeTempsPourUneOptimisation = problemeHebdo
                                                   ->NombreDePasDeTempsPourUneOptimisation;
 
-    int Sparsity = mxPaliers * problemeHebdo->NombreDePays;
-    Sparsity += problemeHebdo->NombreDInterconnexions;
-    if (Sparsity > 100)
-    {
-        Sparsity = 100;
-    }
-
-    int NbTermes = 0;
-    NbTermes += ProblemeAResoudre->NombreDeContraintes;
-
-    int Adder = mxPaliers;
-    Adder += 4;
-    Adder *= problemeHebdo->NombreDePays;
-    Adder += 2 * problemeHebdo->NombreDInterconnexions;
-    Adder *= NombreDePasDeTempsPourUneOptimisation;
-
-    NbTermes += Adder;
-
-    NbTermes += Adder;
-
-    Adder = 3 * problemeHebdo->NombreDInterconnexions * NombreDePasDeTempsPourUneOptimisation;
-    NbTermes += Adder;
-
-    Adder = Sparsity * problemeHebdo->NombreDeContraintesCouplantes;
-    Adder *= (NombreDePasDeTempsPourUneOptimisation);
-    Adder += Sparsity * (7 + 7) * problemeHebdo->NombreDeContraintesCouplantes;
-
-    NbTermes += Adder;
-
-    NbTermes += 3 * problemeHebdo->NombreDePays * NombreDePasDeTempsPourUneOptimisation;
-    NbTermes += problemeHebdo->NombreDePays * NombreDePasDeTempsPourUneOptimisation * 4;
-    NbTermes += problemeHebdo->NombreDePays * NombreDePasDeTempsPourUneOptimisation * 5;
-
-    NbTermes += problemeHebdo->NombreDePays * NombreDePasDeTempsPourUneOptimisation
-                * 2; /*inequality constraint on final hydros level*/
-    NbTermes += 1;   /* constraint includes hydro generation, pumping and final level */
-    NbTermes += 101; /* constraint expressing final level as a sum of stock layers */
-
-    NbTermes += problemeHebdo->NbTermesContraintesPourLesCoutsDeDemarrage;
-
     logs.info();
     logs.info()
       << " Starting Memory Allocation for a Weekly Optimization problem in Canonical form ";
     logs.info() << " ( Problem Size :" << ProblemeAResoudre->NombreDeVariables << " variables "
                 << ProblemeAResoudre->NombreDeContraintes << " Constraints) ";
-    logs.info() << " Expected Number of Non-zero terms in Problem Matrix : " << NbTermes;
-    logs.info();
 
-    OPT_AllocateFromNumberOfVariableConstraints(problemeHebdo->ProblemeAResoudre.get(), NbTermes);
+    OPT_AllocateFromNumberOfVariableConstraints(problemeHebdo->ProblemeAResoudre.get());
 
     int NbIntervalles = problemeHebdo->NombreDePasDeTemps / NombreDePasDeTempsPourUneOptimisation;
 
@@ -140,7 +90,7 @@ void OPT_AllocDuProblemeAOptimiser(PROBLEME_HEBDO* problemeHebdo)
 {
     problemeHebdo->ProblemeAResoudre = std::make_unique<PROBLEME_ANTARES_A_RESOUDRE>();
 
-    int mxPaliers = OPT_DecompteDesVariablesEtDesContraintesDuProblemeAOptimiser(problemeHebdo);
+    OPT_DecompteDesVariablesEtDesContraintesDuProblemeAOptimiser(problemeHebdo);
 
-    optimisationAllocateProblem(problemeHebdo, mxPaliers);
+    optimisationAllocateProblem(problemeHebdo);
 }

src/solver/simulation/adequacy.cpp

@@ -46,7 +46,6 @@ Benchmarking::OptimizationInfo Adequacy::getOptimizationInfo() const
 
     optInfo.nbVariables = Pb->NombreDeVariables;
     optInfo.nbConstraints = Pb->NombreDeContraintes;
-    optInfo.nbNonZeroCoeffs = Pb->NombreDeTermesAllouesDansLaMatriceDesContraintes;
     return optInfo;
 }
 

src/solver/simulation/economy.cpp

@@ -51,7 +51,6 @@ Benchmarking::OptimizationInfo Economy::getOptimizationInfo() const
 
     optInfo.nbVariables = Pb->NombreDeVariables;
     optInfo.nbConstraints = Pb->NombreDeContraintes;
-    optInfo.nbNonZeroCoeffs = Pb->NombreDeTermesAllouesDansLaMatriceDesContraintes;
     return optInfo;
 }
 

src/solver/simulation/include/antares/solver/simulation/sim_structure_probleme_economique.h

@@ -502,7 +502,6 @@ struct PROBLEME_HEBDO
     std::vector<::ShortTermStorage::AREA_INPUT> ShortTermStorage;
 
     /* Optimization problem */
-    uint32_t NbTermesContraintesPourLesCoutsDeDemarrage = 0;
     std::vector<bool> DefaillanceNegativeUtiliserPMinThermique;
     std::vector<bool> DefaillanceNegativeUtiliserHydro;
     std::vector<bool> DefaillanceNegativeUtiliserConsoAbattue;