diff --git a/src/calculate_reward_and_bellman_values.py b/src/calculate_reward_and_bellman_values.py
index a5aa2fe..3b6b32f 100644
--- a/src/calculate_reward_and_bellman_values.py
+++ b/src/calculate_reward_and_bellman_values.py
@@ -160,9 +160,6 @@ def get_disc(
             n_pts_below + n_pts_in + n_pts_above
         )  # Make sure total adds up
         if method == "lines":
-            # in_curve_pts = np.linspace(lbs, ubs, xNsteps-2) # Disc-2 * R
-            # all_pts = np.insert(in_curve_pts, 0, empty, axis=0) # (Disc-1) * R
-            # all_pts = np.insert(all_pts, all_pts.shape[0], full, axis=0) # Disc * R
             above_curve_pts = [
                 np.linspace(ubs[r], full[r], n_pts_above[r], endpoint=True)
                 for r in range(n_reservoirs)
@@ -183,7 +180,6 @@ def get_disc(
                     for r in range(n_reservoirs)
                 ]
             ).T
-            # all_pts = np.concatenate((below_curve_pts, in_curve_pts, above_curve_pts), axis=0) # Disc * R
             diffs_to_ref = all_pts[:, None] - reference_pt[None, :]  # Disc * R
             diffs_to_ref = (
                 diffs_to_ref[:, :, None] * np.eye(n_reservoirs)[None, :, :]
diff --git a/src/display.py b/src/display.py
index ee60fa0..8612dcf 100644
--- a/src/display.py
+++ b/src/display.py
@@ -31,8 +31,6 @@ def draw_usage_values(
                 closest_level
             ]
 
-    # z = np.maximum(np.zeros(reinterpolated_usage_values[area].T.shape), np.minimum(ub*np.ones(reinterpolated_usage_values[area].T.shape), reinterpolated_usage_values[area].T))
-    # z = np.maximum(np.zeros(reinterpolated_usage_values[area].T.shape) - ub, np.minimum(ub*np.ones(reinterpolated_usage_values[area].T.shape), reinterpolated_usage_values[area].T))
     usage_values_plot = go.Figure(
         data=[
             go.Heatmap(
diff --git a/src/estimation.py b/src/estimation.py
index 06d2cc3..c5bbed5 100644
--- a/src/estimation.py
+++ b/src/estimation.py
@@ -324,9 +324,6 @@ def __init__(
             costs:np.ndarray: Cost for every input,
             duals:np.ndarray: Duals for every input first dimension should be the same as inputs,
         """
-        # self.true_controls=controls
-        # self.true_costs=costs
-        # self.true_duals=duals
         self.estimators = np.array(
             [
                 [
@@ -335,7 +332,6 @@ def __init__(
                         costs=costs[week, scenario],
                         duals=duals[week, scenario],
                         correlations=correlations,
-                        # interp_mode=interp_mode,
                     )
                     for scenario in range(param.len_scenario)
                 ]
@@ -409,7 +405,6 @@ def update(
                     inputs=inputs,
                     costs=costs,
                     duals=duals,
-                    #    interp_mode=interp_mode,
                 )
 
     def enrich_estimator(self, n_splits: int = 3) -> None:
@@ -441,8 +436,6 @@ def cleanup_approximations(
         for week_estimators in self.estimators:
             for estimator in week_estimators:
                 estimator.remove_incoherence()
-                # controls=true_controls[week, scenario],
-                # real_costs=true_costs[week, scenario])
 
     def remove_redundants(
         self,
diff --git a/src/functions_iterative.py b/src/functions_iterative.py
index 32969d3..36ea514 100644
--- a/src/functions_iterative.py
+++ b/src/functions_iterative.py
@@ -274,7 +274,7 @@ def itr_control(
     ) = init_iterative_calculation(param, reservoir_management, output_path, X, solver)
     i = 0
 
-    while (gap >= tol_gap and gap >= 0) and i < N:  # and (i<3):
+    while (gap >= tol_gap and gap >= 0) and i < N:
         debut = time()
 
         initial_x, controls = compute_x_multi_scenario(
diff --git a/src/hyperplane_interpolation.py b/src/hyperplane_interpolation.py
index 2201deb..2745cb9 100644
--- a/src/hyperplane_interpolation.py
+++ b/src/hyperplane_interpolation.py
@@ -41,14 +41,11 @@ def BezierAv(
     if last_axify:
         p[:-1] = (p0[:-1] + p1[:-1]) / 2
         p[-1] = p0[-1] * diff_rate + p1[-1] * (1 - diff_rate)
-        # v_sum = np.sum(v[:-1])
-        # v[-1] /= v_sum
         v = v0 * (t) + v1 * (1 - t)
         v[-1] = v0[-1] * (t * (1 - partage) + partage * diff_rate) + v1[-1] * (
             (1 - t) * (1 - partage) + partage * (1 - diff_rate)
         )
     elif norm(v) > 1e-12:
-        # v = v / sum(v[:-1])
         v = v / norm(v) * (n0 + n1) / 2
     return np.array([p, v])
 
diff --git a/src/multi_stock_bellman_value_calculation.py b/src/multi_stock_bellman_value_calculation.py
index f99b5ee..f42eb98 100644
--- a/src/multi_stock_bellman_value_calculation.py
+++ b/src/multi_stock_bellman_value_calculation.py
@@ -11,7 +11,6 @@
 from estimation import Estimator, LinearCostEstimator, LinearInterpolator
 from read_antares_data import TimeScenarioParameter, TimeScenarioIndex
 
-# from julia_sddp import python_to_julia_data, manage_reservoirs_py
 from optimization import (
     AntaresProblem,
     WeeklyBellmanProblem,
@@ -23,8 +22,6 @@
 import numpy.typing as npt
 import numpy as np
 
-# from scipy.interpolate import interp1d
-# from functions_iterative import compute_upper_bound
 from display import ConvergenceProgressBar, draw_usage_values, draw_uvs_sddp
 from scipy.stats import random_correlation
 import pickle as pkl
@@ -383,7 +380,6 @@ def get_bellman_values_from_costs(
                     f"We were using the following future costs estimation: {[f'Cost(lvl) >= {cost} +  (lvl_0 - {input[0]})*{duals[0]} + (lvl_1 - {input[1]})*{duals[1]}' for input, cost, duals in zip(future_costs_approx.inputs, future_costs_approx.costs, future_costs_approx.duals)]}"
                 )
                 raise ValueError
-            # assert (np.min(np.abs(duals_wl))>.1 or np.random.rand()> .01)
 
             # Writing down results
             controls_w[lvl_id] = controls_wls
@@ -1025,7 +1021,6 @@ def select_controls_to_explore(
         * (2 * n_inputs + 1),
         size=(n_weeks, n_scenarios, n_reservoirs),
     )
-    # pseudo_opt_controls += np.random.normal(scale=np.exp(-relative_proximity) * max_var[: ,None ,:]/10, size=(n_weeks, n_scenarios, n_reservoirs))
     controls_to_explore[:, :, 0, :] = np.maximum(
         np.minimum(pseudo_opt_controls, max_gen[:, None, :]), max_pump[:, None, :]
     )
@@ -1174,11 +1169,8 @@ def compute_usage_values_from_costs(
             for j, levels in enumerate(levels_to_test[week, :, i]):
                 # Remove previous constraints / vars
                 problem.reset_solver()
-                # problem.solver.SetSolverSpecificParametersAsString("presolve off")
-
                 try:
                     # Rewrite problem
-                    # print(levels, future_costs_approx_l[week].costs)
                     problem.write_problem(
                         week=week,
                         level_init=levels,
@@ -1229,13 +1221,10 @@ def get_opt_gap(
         ]
     )
     # Computing the optimality gap
-    # pre_update_tot_costs = np.sum(pre_update_costs, axis=0)
-    # real_tot_costs = np.sum(costs, axis=0)
     costs = np.mean(costs, axis=-1)
     opt_gap = max(
         1e-10, min(1, np.max(np.mean((costs - pre_update_costs), axis=1) / max_gap))
     )
-    # print(f"Gap {np.max(np.mean((costs - pre_update_costs),axis=1)/max_gap)}")
     return opt_gap
 
 
@@ -1298,32 +1287,6 @@ def cutting_plane_method(
     max_gap = np.mean(np.max(costs, axis=2) - np.min(costs, axis=2), axis=1)
     rng = np.random.default_rng(1234115)
 
-    # hlpr = Caller(
-    #     param=param,
-    #     multi_stock_management=multi_stock_management,
-    #     name_solver=name_solver,
-    #     list_models=list_models,
-    #     starting_pt=starting_pt,
-    #     inflows=inflows,
-    #     trajectory=trajectory,
-    #     max_gap=max_gap,
-    #     costs_approx=costs_approx,
-    #     future_costs_approx=future_costs_approx,
-    #     nSteps_bellman=nSteps_bellman,
-    #     method=method,
-    #     correlations=correlations,
-    #     maxiter=maxiter,
-    #     precision=precision,
-    #     interp_mode=interp_mode,
-    #     opt_gap=opt_gap,
-    #     divisor=divisor,
-    #     rounding=rounding,
-    #     output_path=output_path,
-    #     rng=rng,
-    #     verbose=False,
-    #     prefix="cut_plan")
-    # prefix="cut_plan"
-
     # To display convergence:
     if verbose:
         pbar = ConvergenceProgressBar(
@@ -1335,8 +1298,6 @@ def cutting_plane_method(
         if verbose:
             pbar.describe("Dynamic Programming")
 
-        # levels, bellman_costs, _, _, future_costs_approx_l = hlpr(get_bellman_values_from_costs,\
-        #            ("levels", "bellman_costs", "bellman_duals", "bellman_controls", "future_costs_approx_l", "usage_values"))
         (
             levels,
             bellman_costs,
@@ -1358,11 +1319,8 @@ def cutting_plane_method(
             verbose=verbose,
         )
         future_costs_approx = future_costs_approx_l[0]
-        # hlpr.update(future_costs_approx = future_costs_approx)
 
         # Evaluate optimal
-        # trajectory, _, _ = hlpr(solve_for_optimal_trajectory,\
-        #                                     ("trajectory", "pseudo_opt_controls", "pseudo_opt_costs"))
         trajectory, pseudo_opt_controls, pseudo_opt_costs = (
             solve_for_optimal_trajectory(
                 param=param,
@@ -1378,7 +1336,6 @@ def cutting_plane_method(
         )
         # Beware, some trajectories seem to be overstep the bounds with values such as -2e-12
 
-        # controls_list = hlpr(select_controls_to_explore, ("controls_list"))
         controls_list = select_controls_to_explore(
             param=param,
             multi_stock_management=multi_stock_management,
@@ -1390,8 +1347,6 @@ def cutting_plane_method(
             pbar.describe("Simulation")
 
         # Evaluating this "optimal" trajectory
-        # hlpr.update(prefix=f"cut_plan_iter_{iter}_")
-        # controls_list, costs, slopes = hlpr(Lget_costs, ("controls_list", "costs", "slopes"))
         controls, costs, slopes = Lget_costs(
             param=param,
             multi_stock_management=multi_stock_management,
@@ -1404,7 +1359,6 @@ def cutting_plane_method(
             prefix=f"cut_plan_iter_{iter}_",
         )
 
-        # opt_gap = hlpr(get_opt_gap, "opt_gap")
         opt_gap = get_opt_gap(
             param=param,
             costs=costs,
@@ -1422,7 +1376,6 @@ def cutting_plane_method(
         # If we want to look at the usage values evolution
         if verbose:
             pbar.describe("Drawing")
-            # hlpr.update(optimal_trajectory=trajectory, n_weeks=param.len_week)
             usage_values, levels_uv = compute_usage_values_from_costs(
                 param=param,
                 multi_stock_management=multi_stock_management,
@@ -1435,7 +1388,6 @@ def cutting_plane_method(
                 divisor=divisor,
                 rounding=rounding,
             )
-            # hlpr(compute_usage_values_from_costs, ("usage_values", "levels_uv"))
             draw_usage_values(
                 usage_values=usage_values,
                 levels_uv=levels_uv,
@@ -1445,7 +1397,6 @@ def cutting_plane_method(
                 trajectory=trajectory,
                 ub=500,
             )
-            # hlpr(draw_usage_values)
             pbar.update(precision=opt_gap)
 
     if verbose:
@@ -1512,40 +1463,14 @@ def iter_bell_vals(
     Returns:
         tuple[np.ndarray, Estimator, list[LinearInterpolator], np.ndarray]: Bellman costs, costs approximation, future costs approximation, levels
     """
-    # Initialize the caller
-    # hlpr = Caller(
-    #     param=param,
-    #     multi_stock_management=multi_stock_management,
-    #     output_path=output_path,
-    #     n_controls_init=n_controls_init,
-    #     starting_pt=starting_pt,
-    #     nSteps_bellman=nSteps_bellman,
-    #     method=method,
-    #     name_solver=name_solver,
-    #     precision=precision,
-    #     maxiter=maxiter,
-    #     correlations=correlations,
-    #     interp_mode=interp_mode,
-    #     divisor=divisor,
-    #     rounding=rounding,
-    #     verbose=verbose,
-    #     already_init=already_init,
-    #     keep_intermed_res=keep_intermed_res)
-
     # Choose first controls to test
-    # controls = hlpr(initialize_controls, 'controls_list')
     controls_list = initialize_controls(
         param=param,
         multi_stock_management=multi_stock_management,
         n_controls_init=n_controls_init,
     )
 
-    # #Initialize the Antares problems
-    # hlpr(initialize_antares_problems, 'list_models')
-
-    # #Get hyperplanes resulting from initial controls
-    # _, duals = hlpr(get_all_costs, ('costs', 'slopes'))
-    # controls, costs, duals = hlpr(Lget_costs, ('controls_list', 'costs', 'slopes'))
+    # Get hyperplanes resulting from initial controls
     controls, costs, duals = Lget_costs(
         param=param,
         multi_stock_management=multi_stock_management,
@@ -1556,16 +1481,12 @@ def iter_bell_vals(
         load_from_protos=True,
         verbose=verbose,
     )
-    # hlpr.update(duals=duals, controls=controls)
-    # hlpr.update(list_models=[])
 
-    # costs_approx = hlpr(LinearCostEstimator, "costs_approx")
     costs_approx = LinearCostEstimator(
         param=param, controls=controls, costs=costs, duals=duals
     )
 
     # Initialize our approximation on future costs
-    # hlpr(initialize_future_costs, "future_costs_approx")
     future_costs_approx = initialize_future_costs(
         param=param,
         starting_pt=starting_pt,
@@ -1574,7 +1495,6 @@ def iter_bell_vals(
     )
 
     # Correlations matrix
-    # hlpr(get_correlation_matrix, "correlations")
     correlations = get_correlation_matrix(
         param=param,
         multi_stock_management=multi_stock_management,
@@ -1582,8 +1502,6 @@ def iter_bell_vals(
     )
 
     # Iterative part
-    # hlpr(cutting_plane_method, ("bellman_costs", "levels", "future_costs_approx_l",
-    #                                   "costs_approx", "optimal_trajectory"))
     bellman_costs, levels, future_costs_approx_l, costs_approx, optimal_trajectory = (
         cutting_plane_method(
             param=param,
@@ -1607,8 +1525,6 @@ def iter_bell_vals(
     )
 
     # Deducing usage values
-    # hlpr.update(nSteps_bellman=101)
-    # hlpr(compute_usage_values_from_costs, ('usage_values', 'levels_imposed'))
     usage_values, levels_imposed = compute_usage_values_from_costs(
         param=param,
         multi_stock_management=multi_stock_management,
@@ -1622,10 +1538,6 @@ def iter_bell_vals(
         rounding=rounding,
     )
 
-    # returns = ("bellman_costs", "costs_approx", "future_costs_approx_l", "levels", "optimal_trajectory",\
-    #       "usage_values")
-    # bellman_costs, costs_approx, future_costs_approx_l,\
-    #       levels, optimal_trajectory, usage_values = hlpr.get(returns)
     return (
         bellman_costs,
         costs_approx,
@@ -1801,24 +1713,6 @@ def iter_bell_vals_v2(
         n_controls_init=n_controls_init,
     )
 
-    # #Initialize the Antares problems
-    # list_models = initialize_antares_problems(
-    #     param=param,
-    #     multi_stock_management=multi_stock_management,
-    #     output_path=output_path,
-    #     name_solver=name_solver,
-    #     direct_bellman_calc=False,
-    #     verbose=verbose,
-    # )
-
-    # Get hyperplanes resulting from initial
-    # costs, duals = get_all_costs(
-    #     param=param,
-    #     list_models=list_models,
-    #     multi_stock_management=multi_stock_management,
-    #     controls_list=controls_list,
-    #     verbose=verbose,
-    # )
     controls_list, costs, duals = Lget_costs(
         param=param,
         multi_stock_management=multi_stock_management,
diff --git a/src/optimization.py b/src/optimization.py
index f331b79..62e0a04 100644
--- a/src/optimization.py
+++ b/src/optimization.py
@@ -1384,34 +1384,8 @@ def get_control_cost(
         # 0 <= control_cost
         control_cost = self.solver.NumVar(0, inf, name=f"control_cost")
 
-        # - max_pump * eff - max_generating <= control_diff <= max_pump * eff + max_generating
-        # control_diffs = [
-        #     [
-        #        [ self.solver.NumVar(
-        #             -mng.reservoir.max_pumping[week] * mng.reservoir.efficiency - mng.reservoir.max_generating[week],
-        #             mng.reservoir.max_pumping[week] * mng.reservoir.efficiency + mng.reservoir.max_generating[week],
-        #             name=f"control_diff_{area}_{s}_{c}"
-        #         )
-        #         for c, _ in enumerate(self.week_costs_estimation[week,s].inputs)]
-        #     for s in range(self.n_scenarios)]
-        # for area, mng in self.managements.items()]
-
         # ========= Constraints ==========
 
-        # Control diffs definition
-        # control_diff_constraints = [
-        #     [
-        #         [
-        #             self.solver.Add(
-        #                 control_diffs[r][s][c] == controls[r,s] - input[r],
-        #                 name=f"control_diff_{r}_{s}_{c}"
-        #             )
-        #         for c, input in enumerate(self.week_costs_estimation[week,s].inputs)]
-        #     for s in range(self.n_scenarios)]
-        # for r, _ in enumerate(self.managements)]
-
-        # Cuts on control costs
-        # control_cost_per_scenario >= cost + <X - X_ref | duals >
         control_cost_constraints = [
             [
                 self.solver.Add(
@@ -1719,7 +1693,6 @@ def solve(
         # Removing unwanted parts of the cost
         cost -= self.future_cost.solution_value() * remove_future_costs
         cost -= self.penalty_cost.solution_value() * remove_penalties
-        # penalty = self.total_penalty.solution_value()
         penalty_duals = np.array(
             [cstr.dual_value() for cstr in self.initial_level_csts]
         )
diff --git a/src/read_antares_data.py b/src/read_antares_data.py
index cc574eb..b24d88c 100644
--- a/src/read_antares_data.py
+++ b/src/read_antares_data.py
@@ -96,9 +96,9 @@ def read_rule_curves(self, dir_study: str) -> None:
             )[:, [0, 2]]
             * self.capacity
         )
-        # assert (
-        #     rule_curves[0, 0] == rule_curves[0, 1]
-        # ), "Initial level is not correctly defined by bottom and upper rule curves"
+        assert (
+            rule_curves[0, 0] == rule_curves[0, 1]
+        ), "Initial level is not correctly defined by bottom and upper rule curves"
         self.initial_level = rule_curves[0, 0]
         bottom_rule_curve = rule_curves[6:365:7, 0]
         upper_rule_curve = rule_curves[6:365:7, 1]
@@ -127,7 +127,7 @@ def read_efficiency(self, hydro_ini_file: ConfigParser) -> None:
 
 
 def generate_mps_file(study_path: str, antares_path: str) -> str:
-    # change_hydro_management_to_heuristic(dir_study=study_path)
+    change_hydro_management_to_heuristic(dir_study=study_path)
 
     name_solver = antares_path.split("/")[-1]
     assert "solver" in name_solver
@@ -162,9 +162,5 @@ def change_hydro_management_to_heuristic(dir_study: str) -> None:
             if "use heuristic" in hydro_ini.keys():
                 hydro_ini["use heuristic"][area] = "true"
 
-    # Solve the weird mps occurence of limited hydro power:
-    # for area in hydro_ini["intra-daily-modulation"].keys():
-    #     hydro_ini["intra-daily-modulation"][area] = "999999.000000"
-
     with open(dir_study + "/input/hydro/hydro.ini", "w") as configfile:  # save
         hydro_ini.write(configfile)
diff --git a/src/sddp.jl b/src/sddp.jl
index 0db5407..0cc476b 100644
--- a/src/sddp.jl
+++ b/src/sddp.jl
@@ -74,9 +74,7 @@ function generate_model(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.
         stages = 3*n_weeks,
         sense = :Min,
         lower_bound = 0.0,
-        optimizer = Clp.Optimizer,
-        # optimizer = HiGHS.Optimizer,
-        # cut_oracle = SDDP.LevelOneCutOracle()
+        optimizer = Clp.Optimizer
     ) do subproblem, stage
 
         # Declaring the state variable 
@@ -121,10 +119,6 @@ function generate_model(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.
         end)
         
         # Define scenarios for inflows
-        # inflows = [
-        #     [reservoirs[r]["inflows"][stage][scenario] for r in 1:n_reservoirs]
-        #     for scenario in 1:n_scenarios
-        # ]
         Ω = [
             (
                 inflows = [reservoirs[r].inflows[modulo_stage, scenario] for r in 1:n_reservoirs],
@@ -178,10 +172,6 @@ function manage_reservoirs(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Ma
     #Simulating
     simulation_results = get_trajectory(n_weeks, n_scenarios, reservoirs, model, norms)
     return simulation_results, model
-    # #Getting the usage values
-    # n_disc = 101
-    # VU, costs = get_usage_values(n_weeks, n_scenarios, reservoirs, model, norms, n_disc)
-    # return VU, costs, simulation_results
 end
 
 function get_trajectory(n_weeks::Int, n_scenarios::Int, reservoirs::Vector{Main.Jl_SDDP.Reservoir}, model, norms::Normalizer)
@@ -238,8 +228,3 @@ end
 
 export manage_reservoirs, get_usage_values, stability_report, reinit_cuts
 end
-# Call the function
-# results = manage_reservoirs(n_weeks, n_scenarios, reservoirs, costs_approx)
-
-# Print the results
-# println(results)