Refactor hydro remix : try another version of hydro remix in order to…

… simplify it

src/solver/simulation/hydro-remix-new.cpp

@@ -14,12 +14,14 @@ int find_min_index(const std::vector<double>& G_plus_H,
                    const std::vector<int>& tried_creux,
                    const std::vector<double>& P_max,
                    const std::vector<bool>& filter_hours_remix,
-                   double top)
+                   const double top)
 {
     double min_val = top;
     int min_idx = -1;
     for (int i = 0; i < G_plus_H.size(); ++i)
     {
+        // Aren't conditions (new_D[i] > 0) and (new_H[i] < P_max[i]) always respected ?
+        // In which case, we should remove them from the following if statement
         if (new_D[i] > 0 && new_H[i] < P_max[i] && tried_creux[i] == 0 && filter_hours_remix[i])
         {
             if (G_plus_H[i] < min_val)
@@ -138,6 +140,64 @@ static void checkInputCorrectness(const std::vector<double>& G,
     }
 }
 
+struct Coordinates
+{
+    unsigned int hour;
+    double value;
+
+    bool operator==(const Coordinates& c)
+    {
+        return (hour == c.hour) && (value == c.value);
+    }
+};
+
+Coordinates find_minimum_coordinates(const std::vector<double>& G_plus_H,
+                                     const std::vector<double>& new_D,
+                                     const std::vector<double>& new_H,
+                                     const std::vector<double>& P_max,
+                                     const std::vector<bool>& filter_hours_remix)
+{
+    double min_val = G_plus_H[0];
+    unsigned int min_idx = 0;
+    for (int i = 0; i < G_plus_H.size(); ++i)
+    {
+        // Aren't conditions (new_D[i] > 0) and (new_H[i] < P_max[i]) always respected ?
+        // In which case, we should remove them from the following if statement
+        if (new_D[i] > 0 && new_H[i] < P_max[i] && filter_hours_remix[i])
+        {
+            if (G_plus_H[i] < min_val)
+            {
+                min_val = G_plus_H[i];
+                min_idx = i;
+            }
+        }
+    }
+    return {min_idx, min_val};
+}
+
+Coordinates find_maximum_coordinates(const std::vector<double>& G_plus_H,
+                                     const std::vector<double>& new_H,
+                                     const std::vector<double>& P_min,
+                                     const std::vector<bool>& filter_hours_remix,
+                                     double ref_value,
+                                     double eps)
+{
+    double max_val = G_plus_H[0];
+    unsigned int max_idx = 0;
+    for (int i = 0; i < G_plus_H.size(); ++i)
+    {
+        if (new_H[i] > P_min[i] && G_plus_H[i] >= ref_value + eps && filter_hours_remix[i])
+        {
+            if (G_plus_H[i] > max_val)
+            {
+                max_val = G_plus_H[i];
+                max_idx = i;
+            }
+        }
+    }
+    return {max_idx, max_val};
+}
+
 RemixHydroOutput new_remix_hydro(const std::vector<double>& G,
                                  const std::vector<double>& H,
                                  const std::vector<double>& D,
@@ -180,8 +240,8 @@ RemixHydroOutput new_remix_hydro(const std::vector<double>& G,
 
     int loop = 1000;
     double eps = 1e-3;
-    double top = *std::ranges::max_element(G) + *std::ranges::max_element(H)
-                 + *std::ranges::max_element(D) + 1;
+    const double top = *std::ranges::max_element(G) + *std::ranges::max_element(H)
+                       + *std::ranges::max_element(D) + 1;
 
     std::vector<bool> filter_hours_remix(G.size(), false);
     for (unsigned int h = 0; h < filter_hours_remix.size(); h++)
@@ -195,6 +255,69 @@ RemixHydroOutput new_remix_hydro(const std::vector<double>& G,
     std::vector<double> G_plus_H(G.size());
     std::transform(G.begin(), G.end(), new_H.begin(), G_plus_H.begin(), std::plus<>());
 
+    // ======= NEW =================================================================
+
+    double delta = std::numeric_limits<double>::max();
+
+    while (loop-- > 0 && delta > 0.)
+    {
+        Coordinates min = find_minimum_coordinates(G_plus_H,
+                                                   new_D,
+                                                   new_H,
+                                                   P_max,
+                                                   filter_hours_remix);
+
+        Coordinates max = find_maximum_coordinates(G_plus_H,
+                                                   new_H,
+                                                   P_min,
+                                                   filter_hours_remix,
+                                                   min.value,
+                                                   eps);
+
+        if (min == max)
+        {
+            break;
+        }
+
+        std::vector<double> intermediate_level(levels.begin() + std::min(min.hour, max.hour),
+                                               levels.begin() + std::max(min.hour, max.hour));
+        double max_pic, max_creux;
+        if (min.hour < max.hour)
+        {
+            max_pic = capa;
+            max_creux = *std::ranges::min_element(intermediate_level);
+        }
+        else
+        {
+            max_pic = capa - *std::ranges::max_element(intermediate_level);
+            max_creux = capa;
+        }
+
+        max_pic = std::min(new_H[max.hour] - P_min[max.hour], max_pic);
+        max_creux = std::min({P_max[min.hour] - new_H[min.hour], new_D[min.hour], max_creux});
+
+        double dif_pic_creux = std::max(G_plus_H[max.hour] - G_plus_H[min.hour], 0.);
+
+        delta = std::max(std::min({max_pic, max_creux, dif_pic_creux / 2.}), 0.);
+
+        if (delta > 0)
+        {
+            new_H[max.hour] -= delta;
+            new_H[min.hour] += delta;
+            new_D[max.hour] = H[max.hour] + D[max.hour] - new_H[max.hour];
+            new_D[min.hour] = H[min.hour] + D[min.hour] - new_H[min.hour];
+        }
+
+        std::transform(G.begin(), G.end(), new_H.begin(), G_plus_H.begin(), std::plus<>());
+        levels[0] = initial_level + inflows[0] - overflow[0] + pump[0] - new_H[0];
+        for (size_t i = 1; i < levels.size(); ++i)
+        {
+            levels[i] = levels[i - 1] + inflows[i] - overflow[i] + pump[i] - new_H[i];
+        }
+    }
+
+    // ========== OLD =====================================================================
+    /*
     while (loop-- > 0)
     {
         std::vector<int> tried_creux(G.size(), 0);
@@ -209,6 +332,11 @@ RemixHydroOutput new_remix_hydro(const std::vector<double>& G,
                                            P_max,
                                            filter_hours_remix,
                                            top);
+            // top : do need this to start min search ?
+            //       if we don't need it, remove it completely
+            // ==> removing top seems to lead to infinite loops in unit tests
+
+            // I tried to remove this if statement along with top, but it caused an infinite loop
             if (idx_creux == -1)
             {
                 break;
@@ -286,6 +414,8 @@ RemixHydroOutput new_remix_hydro(const std::vector<double>& G,
             levels[i] = levels[i - 1] + inflows[i] - overflow[i] + pump[i] - new_H[i];
         }
     }
+    */
+    // =================================================================================
     return {new_H, new_D, levels};
 }