Repository: Grid2op/lightsim2grid (language: python)

Issue:
ImportError when importing ContingencyAnalysisCPP and grid2op is not installed
## Environment
 - Grid2op version: `x.y.z` (**eg** 1.4.0)
 - lightsim version: `0.8.0` (**eg** 0.5.0)
 - System: `GutHub Actions for Windows.`
 - Additional system information

## Bug description

As seen below in the log from pytest, if we try importing ContingencyAnalysisCPP (from lightsim2grid.securityAnalysis import ContingencyAnalysisCPP), and grid2op is not installed, there is an ImportError.

Can you please add grid2op to the required packages for lightsim2grid, or change lightsim2grid imports so that importing ContingencyAnalysisCPP does not require grid2op (preferred)?

Here the error traceback:

==================================== ERRORS ====================================
_______ ERROR collecting pandapower/test/contingency/test_contingency.py _______
ImportError while importing test module '/home/runner/work/pandapower/pandapower/pandapower/test/contingency/test_contingency.py'.
Hint: make sure your test modules/packages have valid Python names.
Traceback:
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/importlib/__init__.py:127: in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
pandapower/test/contingency/test_contingency.py:15: in <module>
    import pandapower.contingency
pandapower/contingency/__init__.py:1: in <module>
    from pandapower.contingency.contingency import *
pandapower/contingency/contingency.py:14: in <module>
    from lightsim2grid.securityAnalysis import ContingencyAnalysisCPP
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/lightsim2grid/securityAnalysis.py:9: in <module>
    from lightsim2grid.contingencyAnalysis import ContingencyAnalysisCPP, ContingencyAnalysis
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/lightsim2grid/contingencyAnalysis.py:18: in <module>
    from lightsim2grid.lightSimBackend import LightSimBackend
/opt/hostedtoolcache/Python/3.8.18/x64/lib/python3.8/site-packages/lightsim2grid/lightSimBackend.py:16: in <module>
    from grid2op.Action import CompleteAction
E   ModuleNotFoundError: No module named 'grid2op'




Discussion / hints:

Fix patch:
diff --git a/.circleci/config.yml b/.circleci/config.yml
index cc6fe65..8f6291b 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -55,6 +55,31 @@ executors:
             - image: silkeh/clang:9  # no c++ 11, does not work
 
 jobs:
+    test_legacy_grid2op:
+        executor: gcc_13
+        resource_class: small
+        steps:
+            - checkout
+            - run: apt-get update && apt-get install python3-full python3-dev python3-pip python3-virtualenv git -y
+            - run: python3 -m virtualenv venv_test
+            - run:
+                  name: "Set up virtual environment"
+                  command: |
+                      source venv_test/bin/activate
+                      pip install --upgrade pip setuptools wheel gym "numpy<1.22" grid2op==1.6.4
+                      pip install -U pybind11
+                      git submodule init
+                      git submodule update
+                      make
+                      CC=gcc python setup.py build
+                      python -m pip install -U .
+                      pip freeze
+            - run:
+                name: legacy test
+                command: |
+                      source venv_test/bin/activate
+                      python -m unittest lightsim2grid/tests/test_compat_legacy_grid2op.py
+
     compile_gcc12:
         executor: gcc_12
         resource_class: small
diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml
index 8ce9807..dadcd64 100644
--- a/.github/workflows/main.yml
+++ b/.github/workflows/main.yml
@@ -50,7 +50,7 @@ jobs:
     steps:
 
       - name: Checkout sources
-        uses: actions/checkout@v3
+        uses: actions/checkout@v3 # v4 not working...
         with:
           submodules: true
 
@@ -76,16 +76,22 @@ jobs:
         run: python setup.py sdist
 
       - name: Install wheel
-        run: pip3 install wheelhouse/*.whl --user
+        run: |
+          pip3 install wheelhouse/*.whl --user
+          pip freeze
 
       # - name: Install GDB
       #   run: yum install -y gdb
 
-      - name: Check package can be imported
+      - name: Check package can be imported (bare install)
         run: |
           python3 -c "import lightsim2grid"
           python3 -c "from lightsim2grid import *"
           python3 -c "from lightsim2grid.newtonpf import newtonpf"
+          python3 -c "from lightsim2grid.timeSerie import TimeSeriesCPP"
+          python3 -c "from lightsim2grid.contingencyAnalysis import ContingencyAnalysisCPP"
+          python3 -c "from lightsim2grid.securityAnalysis import SecurityAnalysisCPP"
+          python3 -c "from lightsim2grid.gridmodel import init, GridModel"
       
       - name: Fix urllib3 (python 3.7)
         if: matrix.python.name == 'cp37'
@@ -98,22 +104,26 @@ jobs:
           python3 -m pip install grid2op
           python3 -m pip freeze
 
-      - name: Check LightSimBackend can be imported 1
+      - name: Check extra can be imported can be imported (with grid2op)
         run:
           python3 -v -c "from lightsim2grid import LightSimBackend"
-
-      - name: Check LightSimBackend can be imported 2
+          python3 -c "from lightsim2grid.timeSerie import TimeSerie"
+          python3 -c "from lightsim2grid.contingencyAnalysis import ContingencyAnalysis"
+          python3 -c "from lightsim2grid.physical_law_checker import PhysicalLawChecker"
+          python3 -c "from lightsim2grid.securityAnalysis import SecurityAnalysis"
+      
+      - name: Check LightSimBackend can be used to create env
         run:
           python3 -v -c "from lightsim2grid import LightSimBackend; import grid2op; env = grid2op.make('l2rpn_case14_sandbox', test=True, backend=LightSimBackend())"
 
       - name: Upload wheel
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v3  # v4 broken
         with:
           name: lightsim2grid-wheel-linux-${{ matrix.python.name }}
           path: wheelhouse/*.whl
 
       - name: Upload source archive
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v3  # v4 broken
         if: matrix.python.name == 'cp311'
         with:
           name: lightsim2grid-sources
@@ -157,12 +167,12 @@ jobs:
     steps:
 
       - name: Checkout sources
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           submodules: true
 
       - name: Setup Python
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python.version }}
 
@@ -210,7 +220,7 @@ jobs:
           python -c "from lightsim2grid import LightSimBackend; import grid2op; env = grid2op.make('l2rpn_case14_sandbox', test=True, backend=LightSimBackend())"
 
       - name: Upload wheel
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: lightsim2grid-wheel-${{ matrix.config.name }}-${{ matrix.python.name }}
           path: dist/*.whl
@@ -228,12 +238,12 @@ jobs:
           }
     steps:
       - name: Checkout sources
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           submodules: true
 
       - name: Setup Python
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python.version }}
 
@@ -272,7 +282,7 @@ jobs:
           python -c "from lightsim2grid import LightSimBackend; import grid2op; env = grid2op.make('l2rpn_case14_sandbox', test=True, backend=LightSimBackend())"
         
       - name: Upload wheel
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: lightsim2grid-wheel-darwin-${{ matrix.python.name }}
           path: dist/*.whl
@@ -311,12 +321,12 @@ jobs:
     steps:
 
       - name: Checkout sources
-        uses: actions/checkout@v3
+        uses: actions/checkout@v4
         with:
           submodules: true
 
       - name: Setup Python
-        uses: actions/setup-python@v4
+        uses: actions/setup-python@v5
         with:
           python-version: ${{ matrix.python.version }}
 
@@ -346,7 +356,7 @@ jobs:
             python -c "from lightsim2grid import LightSimBackend; import grid2op; env = grid2op.make('l2rpn_case14_sandbox', test=True, backend=LightSimBackend())"
 
       - name: Upload wheel
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: wheels-darwin-${{ matrix.python.name }}
           path: ./wheelhouse/*.whl
@@ -363,7 +373,7 @@ jobs:
           path: download
 
       - name: Upload wheels
-        uses: actions/upload-artifact@v3
+        uses: actions/upload-artifact@v4
         with:
           name: lightsim2grid-wheels
           path: |
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 63f0aab..ca87f38 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -18,7 +18,39 @@ Change Log
 - maybe have a look at suitesparse "sliplu" tools ?
 - easier building (get rid of the "make" part)
 
-[0.8.0] 2023-03-18
+[0.8.1] 2024-03-26
+--------------------
+- [FIXED] a bug with shunts when `nb_busbar_per_sub` >= 2
+- [FIXED] some bugs preventing backward compatibility
+- [FIXED] an issue in the computation of gen_q when intialized with pypowsybl
+  (some overflow cpp side leading to infinite number in gen_q)
+- [FIXED] a bug in the "containers" cpp side (wrong bus was assigned)
+  when elements was disconnected, which lead to wrong computations for 
+  time series or contingency analysis.
+- [FIXED] another bug in ContingencyAnalysis (cpp side) leading to wrong computation
+  when a powerline was disconnected
+- [FIXED] some broken imports when grid2op was not installed
+- [FIXED] missing "typing_extension" as required when installation
+- [ADDED] some information of compilation directly in the cpp module
+- [ADDED] some information of compilation available in the python `compilation_options`
+  module python side
+- [ADDED] some convenient methods for `ContingencyAnalysis` python side (most 
+  notably the possibility to initialize it from a `LightSimBackend` and to
+  change the topology of the grid)
+- [ADDED] a "reward" module in lightsim2grid with custom reward
+  based on lightsim2grid.
+- [ADDED] a class `N1ContingencyReward` that can leverage lightsim2grid to 
+  assess the number of safe / unsafe N-1.
+- [IMPROVED] time measurments in python and c++
+- [IMPROVED] now test lightsim2grid with oldest grid2op version
+- [IMPROVED] speed, by accelerating the reading back of the data (now read only once and then
+  pointers are re used)
+- [IMPROVED] c++ side avoid allocating memory (which allow to gain speed python side too)
+- [IMPROVED] type hinting in `LightSimBackend` for all 'public' methods (most 
+  notably the one used by grid2op)
+- [IMPROVED] now the benchmarks are more verbose (detailing some compilation options)
+
+[0.8.0] 2024-03-18
 --------------------
 - [BREAKING] now able to retrieve `dcSbus` with a dedicated method (and not with the old `get_Sbus`).
   If you previously used `gridmodel.get_Sbus()` to retrieve the Sbus used for DC powerflow, please use
diff --git a/README.md b/README.md
index cf1c078..c989f34 100644
--- a/README.md
+++ b/README.md
@@ -222,7 +222,7 @@ For example: `export PATH_NICSLU=/home/user/Documents/cktso`
 
 #### Enable 03 optimization
 By default, at least on ubuntu, only the "-O2" compiler flags is used. To use the O3 optimization flag, you need
-to specify the `__COMPLILE_O3` environment variable: `set __COMPLILE_O3=1` before the compilation (so before
+to specify the `__O3_OPTIM` environment variable: `set __O3_OPTIM=1` before the compilation (so before
 `python3 setup.py build` or `python -m pip install -e .`)
 
 This compilation argument will increase the compilation time, but will make the package faster.
diff --git a/benchmarks/benchmark_grid_size.py b/benchmarks/benchmark_grid_size.py
index a2e2668..550c633 100644
--- a/benchmarks/benchmark_grid_size.py
+++ b/benchmarks/benchmark_grid_size.py
@@ -46,7 +46,7 @@
               "case300.json",
               "case1354pegase.json",
               "case1888rte.json",
-              "GBnetwork.json",  # 2224 buses
+            #   "GBnetwork.json",  # 2224 buses
               "case2848rte.json",
               "case2869pegase.json",
               "case3120sp.json",
diff --git a/benchmarks/utils_benchmark.py b/benchmarks/utils_benchmark.py
index 1286737..f1fe163 100644
--- a/benchmarks/utils_benchmark.py
+++ b/benchmarks/utils_benchmark.py
@@ -11,6 +11,7 @@
 import numpy as np
 from tqdm import tqdm
 import argparse
+import datetime
 from grid2op.Environment import MultiMixEnvironment
 import pdb
 
@@ -127,6 +128,8 @@ def str2bool(v):
 
 
 def print_configuration():
+    print()
+    print(f"- date: {datetime.datetime.now():%Y-%m-%d %H:%M %z} {time.localtime().tm_zone}")
     try:
         import platform
         print(f"- system: {platform.system()} {platform.release()}")
@@ -156,6 +159,14 @@ def print_configuration():
     print(f"- numpy version: {np.__version__}")
     print(f"- pandas version: {pd.__version__}")
     print(f"- pandapower version: {pp.__version__}")
-    print(f"- lightsim2grid version: {lightsim2grid.__version__}")
     print(f"- grid2op version: {grid2op.__version__}")
+    print(f"- lightsim2grid version: {lightsim2grid.__version__}")
+    if hasattr(lightsim2grid, "compilation_options"):
+        print(f"- lightsim2grid extra information: ")
+        print()
+        print(f"\t- klu_solver_available: {lightsim2grid.compilation_options.klu_solver_available} ")
+        print(f"\t- nicslu_solver_available: {lightsim2grid.compilation_options.nicslu_solver_available} ")
+        print(f"\t- cktso_solver_available: {lightsim2grid.compilation_options.cktso_solver_available} ")
+        print(f"\t- compiled_march_native: {lightsim2grid.compilation_options.compiled_march_native} ")
+        print(f"\t- compiled_o3_optim: {lightsim2grid.compilation_options.compiled_o3_optim} ")
     print()
diff --git a/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.10.0.rst b/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.10.0.rst
new file mode 100644
index 0000000..9015dad
--- /dev/null
+++ b/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.10.0.rst
@@ -0,0 +1,138 @@
+Lightsim2grid 0.8.0 and grid2op 1.10.0
+======================================
+
+l2rpn_case14_sandbox
+---------------------
+
+Configuration:
+
+- date: 2024-03-18 14:43  CET
+- system: Linux 5.15.0-56-generic
+- OS: ubuntu 20.04
+- processor: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
+- python version: 3.10.13.final.0 (64 bit)
+- numpy version: 1.23.5
+- pandas version: 2.2.1
+- pandapower version: 2.13.1
+- lightsim2grid version: 0.8.0
+- grid2op version: 1.10.0
+
+====================  ======================  ===================================  ============================
+case14_sandbox          grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
+====================  ======================  ===================================  ============================
+PP                                      46.1                               18.5                          6.67
+GS                                     712                                  0.539                        0.398
+GS synch                               735                                  0.487                        0.345
+NR single (SLU)                        916                                  0.223                        0.0809
+NR (SLU)                               897                                  0.236                        0.0842
+NR single (KLU)                        971                                  0.163                        0.0217
+NR (KLU)                               971                                  0.162                        0.0222
+NR single (NICSLU )                   979                                  0.16                         0.0212
+NR (NICSLU )                          973                                  0.162                        0.0217
+NR single (CKTSO )                    975                                  0.16                         0.0198
+NR (CKTSO )                           964                                  0.163                        0.0212
+FDPF XB (SLU)                          964                                  0.172                        0.0318
+FDPF BX (SLU)                          954                                  0.183                        0.0445
+FDPF XB (KLU)                          970                                  0.166                        0.0263
+FDPF BX (KLU)                          964                                  0.175                        0.0371
+FDPF XB (NICSLU )                     971                                  0.164                        0.026
+FDPF BX (NICSLU )                     939                                  0.18                         0.0381
+FDPF XB (CKTSO )                      964                                  0.169                        0.027
+FDPF BX (CKTSO )                      958                                  0.177                        0.0381
+====================  ======================  ===================================  ============================
+
+
+Differences:
+
+============================  ==============  ==============  ================
+case14_sandbox (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
+============================  ==============  ==============  ================
+PP (ref)                            0               0                 0
+GS                                  0.000122        7.63e-06          7.63e-06
+GS synch                            0.000122        7.63e-06          7.63e-06
+NR single (SLU)                     0.000122        7.63e-06          7.63e-06
+NR (SLU)                            0.000122        7.63e-06          7.63e-06
+NR single (KLU)                     0.000122        7.63e-06          7.63e-06
+NR (KLU)                            0.000122        7.63e-06          7.63e-06
+NR single (NICSLU )                0.000122        7.63e-06          7.63e-06
+NR (NICSLU )                       0.000122        7.63e-06          7.63e-06
+NR single (CKTSO )                 0.000122        7.63e-06          7.63e-06
+NR (CKTSO )                        0.000122        7.63e-06          7.63e-06
+FDPF XB (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (NICSLU )                  0.000122        7.63e-06          7.63e-06
+FDPF BX (NICSLU )                  0.000122        7.63e-06          7.63e-06
+FDPF XB (CKTSO )                   0.000122        7.63e-06          7.63e-06
+FDPF BX (CKTSO )                   0.000122        7.63e-06          7.63e-06
+============================  ==============  ==============  ================
+
+
+l2rpn_neurips_2020_track2_small
+---------------------------------
+
+Configuration: 
+
+- date: 2024-03-18 14:50  CET
+- system: Linux 5.15.0-56-generic
+- OS: ubuntu 20.04
+- processor: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
+- python version: 3.10.13.final.0 (64 bit)
+- numpy version: 1.23.5
+- pandas version: 2.2.1
+- pandapower version: 2.13.1
+- lightsim2grid version: 0.8.0
+- grid2op version: 1.10.0
+
+=====================  ======================  ===================================  ============================
+neurips_2020_track2      grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
+=====================  ======================  ===================================  ============================
+PP                                      41.4                                20.8                           8.67
+GS                                       3.56                              280                           279
+GS synch                                36                                  26.7                          26.5
+NR single (SLU)                        522                                   0.952                         0.776
+NR (SLU)                               508                                   0.996                         0.813
+NR single (KLU)                        786                                   0.314                         0.143
+NR (KLU)                               787                                   0.315                         0.144
+NR single (NICSLU )                   786                                   0.308                         0.135
+NR (NICSLU )                          782                                   0.308                         0.136
+NR single (CKTSO )                    785                                   0.301                         0.128
+NR (CKTSO )                           782                                   0.306                         0.134
+FDPF XB (SLU)                          736                                   0.403                         0.233
+FDPF BX (SLU)                          722                                   0.424                         0.252
+FDPF XB (KLU)                          753                                   0.363                         0.191
+FDPF BX (KLU)                          745                                   0.381                         0.209
+FDPF XB (NICSLU )                     754                                   0.363                         0.191
+FDPF BX (NICSLU )                     741                                   0.38                          0.207
+FDPF XB (CKTSO )                      748                                   0.364                         0.193
+FDPF BX (CKTSO )                      736                                   0.382                         0.21
+=====================  ======================  ===================================  ============================
+
+
+Differences:
+=================================  ==============  ==============  ================
+neurips_2020_track2 (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
+=================================  ==============  ==============  ================
+PP (ref)                                  0              0                 0
+GS                                        6.1e-05        3.81e-06          1.53e-05
+GS synch                                  6.1e-05        3.81e-06          1.53e-05
+NR single (SLU)                           6.1e-05        0                 9.54e-07
+NR (SLU)                                  6.1e-05        0                 9.54e-07
+NR single (KLU)                           6.1e-05        0                 9.54e-07
+NR (KLU)                                  6.1e-05        0                 9.54e-07
+NR single (NICSLU )                      6.1e-05        0                 9.54e-07
+NR (NICSLU )                             6.1e-05        0                 9.54e-07
+NR single (CKTSO )                       6.1e-05        0                 9.54e-07
+NR (CKTSO )                              6.1e-05        0                 9.54e-07
+FDPF XB (SLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (SLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (KLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (KLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (NICSLU )                        6.1e-05        1.91e-06          1.53e-05
+FDPF BX (NICSLU )                        6.1e-05        1.91e-06          7.63e-06
+FDPF XB (CKTSO )                         6.1e-05        1.91e-06          1.53e-05
+FDPF BX (CKTSO )                         6.1e-05        1.91e-06          7.63e-06
+=================================  ==============  ==============  ================
+
+
diff --git a/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.9.8.rst b/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.9.8.rst
new file mode 100644
index 0000000..69ca3a9
--- /dev/null
+++ b/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.9.8.rst
@@ -0,0 +1,136 @@
+Lightsim2grid 0.8.0 and grid2op 1.9.8
+======================================
+
+l2rpn_case14_sandbox
+----------------------
+
+Configuration:
+
+- date: 2024-03-18 14:52  CET
+- system: Linux 5.15.0-56-generic
+- OS: ubuntu 20.04
+- processor: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
+- python version: 3.10.13.final.0 (64 bit)
+- numpy version: 1.23.5
+- pandas version: 2.2.1
+- pandapower version: 2.13.1
+- lightsim2grid version: 0.8.0
+- grid2op version: 1.9.8
+
+====================  ======================  ===================================  ============================
+case14_sandbox          grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
+====================  ======================  ===================================  ============================
+PP                                      37.9                               18.7                          6.43
+GS                                     728                                  0.535                        0.395
+GS synch                               758                                  0.483                        0.345
+NR single (SLU)                        952                                  0.219                        0.0798
+NR (SLU)                               929                                  0.231                        0.0843
+NR single (KLU)                        979                                  0.164                        0.0222
+NR (KLU)                              1010                                  0.159                        0.022
+NR single (NICSLU)                  1000                                  0.16                         0.0213
+NR (NICSLU)                         1010                                  0.159                        0.0218
+NR single (CKTSO)                    992                                  0.159                        0.0204
+NR (CKTSO)                          1000                                  0.159                        0.0212
+FDPF XB (SLU)                         1010                                  0.169                        0.0323
+FDPF BX (SLU)                          989                                  0.182                        0.0455
+FDPF XB (KLU)                         1010                                  0.164                        0.0263
+FDPF BX (KLU)                         1000                                  0.174                        0.0376
+FDPF XB (NICSLU)                     999                                  0.164                        0.0263
+FDPF BX (NICSLU)                     985                                  0.176                        0.0379
+FDPF XB (CKTSO)                      994                                  0.166                        0.0266
+FDPF BX (CKTSO)                      990                                  0.173                        0.0374
+====================  ======================  ===================================  ============================
+
+
+Differences:
+
+============================  ==============  ==============  ================
+case14_sandbox (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
+============================  ==============  ==============  ================
+PP (ref)                            0               0                 0
+GS                                  0.000122        7.63e-06          7.63e-06
+GS synch                            0.000122        7.63e-06          7.63e-06
+NR single (SLU)                     0.000122        7.63e-06          7.63e-06
+NR (SLU)                            0.000122        7.63e-06          7.63e-06
+NR single (KLU)                     0.000122        7.63e-06          7.63e-06
+NR (KLU)                            0.000122        7.63e-06          7.63e-06
+NR single (NICSLU)                0.000122        7.63e-06          7.63e-06
+NR (NICSLU)                       0.000122        7.63e-06          7.63e-06
+NR single (CKTSO)                 0.000122        7.63e-06          7.63e-06
+NR (CKTSO)                        0.000122        7.63e-06          7.63e-06
+FDPF XB (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (NICSLU)                  0.000122        7.63e-06          7.63e-06
+FDPF BX (NICSLU)                  0.000122        7.63e-06          7.63e-06
+FDPF XB (CKTSO)                   0.000122        7.63e-06          7.63e-06
+FDPF BX (CKTSO)                   0.000122        7.63e-06          7.63e-06
+============================  ==============  ==============  ================
+
+l2rpn_neurips_2020_track2_small
+--------------------------------
+
+Configuration: 
+
+- date: 2024-03-18 15:00  CET
+- system: Linux 5.15.0-56-generic
+- OS: ubuntu 20.04
+- processor: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
+- python version: 3.10.13.final.0 (64 bit)
+- numpy version: 1.23.5
+- pandas version: 2.2.1
+- pandapower version: 2.13.1
+- lightsim2grid version: 0.8.0
+- grid2op version: 1.9.8
+
+=====================  ======================  ===================================  ============================
+neurips_2020_track2      grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
+=====================  ======================  ===================================  ============================
+PP                                       16.3                               21.2                           8.24
+GS                                        3.5                              284                           284
+GS synch                                 35.7                               26.9                          26.7
+NR single (SLU)                         531                                  0.952                         0.776
+NR (SLU)                                513                                  0.995                         0.812
+NR single (KLU)                         798                                  0.317                         0.145
+NR (KLU)                                796                                  0.317                         0.145
+NR single (NICSLU)                    806                                  0.307                         0.136
+NR (NICSLU)                           812                                  0.305                         0.134
+NR single (CKTSO)                     815                                  0.299                         0.129
+NR (CKTSO)                            809                                  0.302                         0.131
+FDPF XB (SLU)                           736                                  0.415                         0.24
+FDPF BX (SLU)                           735                                  0.425                         0.254
+FDPF XB (KLU)                           769                                  0.366                         0.195
+FDPF BX (KLU)                           761                                  0.381                         0.211
+FDPF XB (NICSLU)                      774                                  0.362                         0.191
+FDPF BX (NICSLU)                      766                                  0.378                         0.209
+FDPF XB (CKTSO)                       736                                  0.38                          0.202
+FDPF BX (CKTSO)                       729                                  0.4                           0.222
+=====================  ======================  ===================================  ============================
+
+
+Differences:
+=================================  ==============  ==============  ================
+neurips_2020_track2 (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
+=================================  ==============  ==============  ================
+PP (ref)                                  0              0                 0
+GS                                        6.1e-05        3.81e-06          1.53e-05
+GS synch                                  6.1e-05        3.81e-06          1.53e-05
+NR single (SLU)                           6.1e-05        0                 9.54e-07
+NR (SLU)                                  6.1e-05        0                 9.54e-07
+NR single (KLU)                           6.1e-05        0                 9.54e-07
+NR (KLU)                                  6.1e-05        0                 9.54e-07
+NR single (NICSLU)                      6.1e-05        0                 9.54e-07
+NR (NICSLU)                             6.1e-05        0                 9.54e-07
+NR single (CKTSO)                       6.1e-05        0                 9.54e-07
+NR (CKTSO)                              6.1e-05        0                 9.54e-07
+FDPF XB (SLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (SLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (KLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (KLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (NICSLU)                        6.1e-05        1.91e-06          1.53e-05
+FDPF BX (NICSLU)                        6.1e-05        1.91e-06          7.63e-06
+FDPF XB (CKTSO)                         6.1e-05        1.91e-06          1.53e-05
+FDPF BX (CKTSO)                         6.1e-05        1.91e-06          7.63e-06
+=================================  ==============  ==============  ================
+
diff --git a/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.1_glop1.10.0.rst b/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.1_glop1.10.0.rst
new file mode 100644
index 0000000..2f17aa4
--- /dev/null
+++ b/docs/benchmark_solver/ubuntu_2004_dell/ls0.8.1_glop1.10.0.rst
@@ -0,0 +1,149 @@
+Lightsim2grid 0.8.0 and grid2op 1.10.0
+======================================
+
+l2rpn_case14_sandbox
+---------------------
+
+Configuration:
+
+- date: 2024-03-25 17:53  CET
+- system: Linux 5.15.0-56-generic
+- OS: ubuntu 20.04
+- processor: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
+- python version: 3.10.13.final.0 (64 bit)
+- numpy version: 1.23.5
+- pandas version: 2.2.1
+- pandapower version: 2.13.1
+- grid2op version: 1.10.1
+- lightsim2grid version: 0.8.1
+- lightsim2grid extra information: 
+
+	- klu_solver_available: True 
+	- nicslu_solver_available: True 
+	- cktso_solver_available: True 
+	- compiled_march_native: True 
+	- compiled_o3_optim: True 
+
+====================  ======================  ===================================  ============================
+case14_sandbox          grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
+====================  ======================  ===================================  ============================
+PP                                      46.3                               18.4                          6.57
+GS                                     757                                  0.474                        0.378
+GS synch                               769                                  0.445                        0.348
+NR single (SLU)                        960                                  0.184                        0.0831
+NR (SLU)                               952                                  0.189                        0.0819
+NR single (KLU)                       1030                                  0.12                         0.0221
+NR (KLU)                              1030                                  0.118                        0.0202
+NR single (NICSLU *)                  1020                                  0.121                        0.022
+NR (NICSLU *)                         1020                                  0.119                        0.02
+NR single (CKTSO *)                   1020                                  0.119                        0.0211
+NR (CKTSO *)                           989                                  0.121                        0.0192
+FDPF XB (SLU)                         1010                                  0.13                         0.032
+FDPF BX (SLU)                         1010                                  0.143                        0.0451
+FDPF XB (KLU)                         1020                                  0.124                        0.0263
+FDPF BX (KLU)                         1010                                  0.134                        0.0377
+FDPF XB (NICSLU *)                    1010                                  0.126                        0.0267
+FDPF BX (NICSLU *)                    1020                                  0.134                        0.0383
+FDPF XB (CKTSO *)                     1010                                  0.125                        0.0268
+FDPF BX (CKTSO *)                     1000                                  0.136                        0.0381
+====================  ======================  ===================================  ============================
+
+Differences:
+
+============================  ==============  ==============  ================
+case14_sandbox (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
+============================  ==============  ==============  ================
+PP (ref)                            0               0                 0
+GS                                  0.000122        7.63e-06          7.63e-06
+GS synch                            0.000122        7.63e-06          7.63e-06
+NR single (SLU)                     0.000122        7.63e-06          7.63e-06
+NR (SLU)                            0.000122        7.63e-06          7.63e-06
+NR single (KLU)                     0.000122        7.63e-06          7.63e-06
+NR (KLU)                            0.000122        7.63e-06          7.63e-06
+NR single (NICSLU *)                0.000122        7.63e-06          7.63e-06
+NR (NICSLU *)                       0.000122        7.63e-06          7.63e-06
+NR single (CKTSO *)                 0.000122        7.63e-06          7.63e-06
+NR (CKTSO *)                        0.000122        7.63e-06          7.63e-06
+FDPF XB (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (NICSLU *)                  0.000122        7.63e-06          7.63e-06
+FDPF BX (NICSLU *)                  0.000122        7.63e-06          7.63e-06
+FDPF XB (CKTSO *)                   0.000122        7.63e-06          7.63e-06
+FDPF BX (CKTSO *)                   0.000122        7.63e-06          7.63e-06
+============================  ==============  ==============  ================
+
+
+l2rpn_neurips_2020_track2_small
+---------------------------------
+
+Configuration: 
+
+- date: 2024-03-25 17:59  CET
+- system: Linux 5.15.0-56-generic
+- OS: ubuntu 20.04
+- processor: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
+- python version: 3.10.13.final.0 (64 bit)
+- numpy version: 1.23.5
+- pandas version: 2.2.1
+- pandapower version: 2.13.1
+- grid2op version: 1.10.1
+- lightsim2grid version: 0.8.1
+- lightsim2grid extra information: 
+
+	- klu_solver_available: True 
+	- nicslu_solver_available: True 
+	- cktso_solver_available: True 
+	- compiled_march_native: True 
+	- compiled_o3_optim: True 
+
+=====================  ======================  ===================================  ============================
+neurips_2020_track2      grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
+=====================  ======================  ===================================  ============================
+PP                                      41.5                                20.7                           8.6
+GS                                       3.74                              266                           266
+GS synch                                35.8                                26.9                          26.8
+NR single (SLU)                        536                                   0.897                         0.767
+NR (SLU)                               505                                   0.959                         0.818
+NR single (KLU)                        811                                   0.268                         0.144
+NR (KLU)                               820                                   0.256                         0.131
+NR single (NICSLU *)                   813                                   0.259                         0.134
+NR (NICSLU *)                          827                                   0.243                         0.118
+NR single (CKTSO *)                    814                                   0.257                         0.131
+NR (CKTSO *)                           829                                   0.24                          0.116
+FDPF XB (SLU)                          762                                   0.352                         0.232
+FDPF BX (SLU)                          749                                   0.373                         0.252
+FDPF XB (KLU)                          786                                   0.307                         0.188
+FDPF BX (KLU)                          776                                   0.327                         0.206
+FDPF XB (NICSLU *)                     786                                   0.308                         0.188
+FDPF BX (NICSLU *)                     771                                   0.324                         0.204
+FDPF XB (CKTSO *)                      784                                   0.309                         0.19
+FDPF BX (CKTSO *)                      773                                   0.329                         0.209
+=====================  ======================  ===================================  ============================
+
+Differences:
+
+=================================  ==============  ==============  ================
+neurips_2020_track2 (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
+=================================  ==============  ==============  ================
+PP (ref)                                  0              0                 0
+GS                                        6.1e-05        3.81e-06          1.53e-05
+GS synch                                  6.1e-05        3.81e-06          1.53e-05
+NR single (SLU)                           6.1e-05        0                 9.54e-07
+NR (SLU)                                  6.1e-05        0                 9.54e-07
+NR single (KLU)                           6.1e-05        0                 9.54e-07
+NR (KLU)                                  6.1e-05        0                 9.54e-07
+NR single (NICSLU *)                      6.1e-05        0                 9.54e-07
+NR (NICSLU *)                             6.1e-05        0                 9.54e-07
+NR single (CKTSO *)                       6.1e-05        0                 9.54e-07
+NR (CKTSO *)                              6.1e-05        0                 9.54e-07
+FDPF XB (SLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (SLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (KLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (KLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (NICSLU *)                        6.1e-05        1.91e-06          1.53e-05
+FDPF BX (NICSLU *)                        6.1e-05        1.91e-06          7.63e-06
+FDPF XB (CKTSO *)                         6.1e-05        1.91e-06          1.53e-05
+FDPF BX (CKTSO *)                         6.1e-05        1.91e-06          7.63e-06
+=================================  ==============  ==============  ================
diff --git a/docs/benchmarks.rst b/docs/benchmarks.rst
index 4ad4aa2..f4f028d 100644
--- a/docs/benchmarks.rst
+++ b/docs/benchmarks.rst
@@ -1,6 +1,6 @@
 
-Benchmarks
-============
+Benchmarks (solvers)
+======================
 
 In this paragraph we will expose some brief benchmarks about the use of lightsim2grid in the grid2op settings.
 The code to run these benchmarks are given with this package int the [benchmark](./benchmarks) folder.
@@ -16,15 +16,23 @@ compared with pandapower when using grid2op.
 
 All of them has been run on a computer with a the following characteristics:
 
-- system: Linux 5.11.0-40-generic
+- date: 2024-03-25 17:53  CET
+- system: Linux 5.15.0-56-generic
 - OS: ubuntu 20.04
-- processor: Intel(R) Core(TM) i7-4790K CPU @ 4.00GHz
-- python version: 3.8.10.final.0 (64 bit)
-- numpy version: 1.18.5
-- pandas version: 1.1.4
-- pandapower version: 2.7.0
-- lightsim2grid version: 0.6.0
-- grid2op version: 1.6.4
+- processor: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70GHz
+- python version: 3.10.13.final.0 (64 bit)
+- numpy version: 1.23.5
+- pandas version: 2.2.1
+- pandapower version: 2.13.1
+- grid2op version: 1.10.1
+- lightsim2grid version: 0.8.1
+- lightsim2grid extra information: 
+
+	- klu_solver_available: True 
+	- nicslu_solver_available: True 
+	- cktso_solver_available: True 
+	- compiled_march_native: True 
+	- compiled_o3_optim: True 
 
 
 To run the benchmark `cd` in the [benchmark](./benchmarks) folder and type:
@@ -130,19 +138,29 @@ stricly greater, for all benchmarks, than **solver powerflow time** (the closer
 
 First on an environment based on the IEEE case 14 grid:
 
-==================  ======================  ===================================  ============================
-case14_sandbox        grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
-==================  ======================  ===================================  ============================
-PP                                    70.5                               11                            4.27
-LS+GS                                881                                  0.447                        0.327
-LS+GS S                              877                                  0.446                        0.327
-LS+SLU (single)                     1110                                  0.191                        0.0655
-LS+SLU                              1120                                  0.195                        0.0683
-LS+KLU (single)                     1200                                  0.138                        0.0176
-LS+KLU                              1180                                  0.141                        0.0188
-LS+NICSLU (single)                  1200                                  0.139                        0.0179
-LS+NICSLU                           1200                                  0.139                        0.0184
-==================  ======================  ===================================  ============================
+====================  ======================  ===================================  ============================
+case14_sandbox          grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
+====================  ======================  ===================================  ============================
+PP                                      46.3                               18.4                          6.57
+GS                                     757                                  0.474                        0.378
+GS synch                               769                                  0.445                        0.348
+NR single (SLU)                        960                                  0.184                        0.0831
+NR (SLU)                               952                                  0.189                        0.0819
+NR single (KLU)                       1030                                  0.12                         0.0221
+NR (KLU)                              1030                                  0.118                        0.0202
+NR single (NICSLU *)                  1020                                  0.121                        0.022
+NR (NICSLU *)                         1020                                  0.119                        0.02
+NR single (CKTSO *)                   1020                                  0.119                        0.0211
+NR (CKTSO *)                           989                                  0.121                        0.0192
+FDPF XB (SLU)                         1010                                  0.13                         0.032
+FDPF BX (SLU)                         1010                                  0.143                        0.0451
+FDPF XB (KLU)                         1020                                  0.124                        0.0263
+FDPF BX (KLU)                         1010                                  0.134                        0.0377
+FDPF XB (NICSLU *)                    1010                                  0.126                        0.0267
+FDPF BX (NICSLU *)                    1020                                  0.134                        0.0383
+FDPF XB (CKTSO *)                     1010                                  0.125                        0.0268
+FDPF BX (CKTSO *)                     1000                                  0.136                        0.0381
+====================  ======================  ===================================  ============================
 
 From a grid2op perspective, lightsim2grid allows to compute up to ~1200 steps each second on the case 14 and
 "only" 70 for the default PandaPower Backend, leading to a speed up of **~17** in this case
@@ -155,15 +173,25 @@ Then on an environment based on the IEEE case 118:
 =====================  ======================  ===================================  ============================
 neurips_2020_track2      grid2op speed (it/s)    grid2op 'backend.runpf' time (ms)    solver powerflow time (ms)
 =====================  ======================  ===================================  ============================
-PP                                       39.6                               13.3                           5.58
-LS+GS                                     5.3                              188                           188
-LS+GS S                                  36.5                               26.6                          26.4
-LS+SLU (single)                         642                                  0.775                         0.607
-LS+SLU                                  588                                  0.932                         0.769
-LS+KLU (single)                         945                                  0.277                         0.116
-LS+KLU                                  918                                  0.306                         0.144
-LS+NICSLU (single)                      947                                  0.274                         0.11
-LS+NICSLU                               929                                  0.298                         0.134
+PP                                      41.5                                20.7                           8.6
+GS                                       3.74                              266                           266
+GS synch                                35.8                                26.9                          26.8
+NR single (SLU)                        536                                   0.897                         0.767
+NR (SLU)                               505                                   0.959                         0.818
+NR single (KLU)                        811                                   0.268                         0.144
+NR (KLU)                               820                                   0.256                         0.131
+NR single (NICSLU *)                   813                                   0.259                         0.134
+NR (NICSLU *)                          827                                   0.243                         0.118
+NR single (CKTSO *)                    814                                   0.257                         0.131
+NR (CKTSO *)                           829                                   0.24                          0.116
+FDPF XB (SLU)                          762                                   0.352                         0.232
+FDPF BX (SLU)                          749                                   0.373                         0.252
+FDPF XB (KLU)                          786                                   0.307                         0.188
+FDPF BX (KLU)                          776                                   0.327                         0.206
+FDPF XB (NICSLU *)                     786                                   0.308                         0.188
+FDPF BX (NICSLU *)                     771                                   0.324                         0.204
+FDPF XB (CKTSO *)                      784                                   0.309                         0.19
+FDPF BX (CKTSO *)                      773                                   0.329                         0.209
 =====================  ======================  ===================================  ============================
 
 For an environment based on the IEEE 118, the speed up in using lightsim + KLU (LS+KLU) is **~24** time faster than
@@ -215,14 +243,24 @@ Here are the results for the IEEE case 14 (max over 1000 powerflows):
 case14_sandbox (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
 ============================  ==============  ==============  ================
 PP (ref)                            0               0                 0
-LS+GS                               0.000122        7.63e-06          7.63e-06
-LS+GS S                             0.000122        7.63e-06          7.63e-06
-LS+SLU (single)                     0.000122        7.63e-06          7.63e-06
-LS+SLU                              0.000122        7.63e-06          7.63e-06
-LS+KLU (single)                     0.000122        7.63e-06          7.63e-06
-LS+KLU                              0.000122        7.63e-06          7.63e-06
-LS+NICSLU (single)                  0.000122        7.63e-06          7.63e-06
-LS+NICSLU                           0.000122        7.63e-06          7.63e-06
+GS                                  0.000122        7.63e-06          7.63e-06
+GS synch                            0.000122        7.63e-06          7.63e-06
+NR single (SLU)                     0.000122        7.63e-06          7.63e-06
+NR (SLU)                            0.000122        7.63e-06          7.63e-06
+NR single (KLU)                     0.000122        7.63e-06          7.63e-06
+NR (KLU)                            0.000122        7.63e-06          7.63e-06
+NR single (NICSLU *)                0.000122        7.63e-06          7.63e-06
+NR (NICSLU *)                       0.000122        7.63e-06          7.63e-06
+NR single (CKTSO *)                 0.000122        7.63e-06          7.63e-06
+NR (CKTSO *)                        0.000122        7.63e-06          7.63e-06
+FDPF XB (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (SLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF BX (KLU)                       0.000122        7.63e-06          7.63e-06
+FDPF XB (NICSLU *)                  0.000122        7.63e-06          7.63e-06
+FDPF BX (NICSLU *)                  0.000122        7.63e-06          7.63e-06
+FDPF XB (CKTSO *)                   0.000122        7.63e-06          7.63e-06
+FDPF BX (CKTSO *)                   0.000122        7.63e-06          7.63e-06
 ============================  ==============  ==============  ================
 
 .. note::
@@ -237,14 +275,24 @@ Here are the results for the IEEE case 118 (max over 1000 powerflows):
 neurips_2020_track2 (1000 iter)      Δ aor (amps)    Δ gen_p (MW)    Δ gen_q (MVAr)
 =================================  ==============  ==============  ================
 PP (ref)                                  0              0                 0
-LS+GS                                     6.1e-05        3.81e-06          1.53e-05
-LS+GS S                                   6.1e-05        3.81e-06          1.53e-05
-LS+SLU (single)                           6.1e-05        0                 9.54e-07
-LS+SLU                                    6.1e-05        0                 9.54e-07
-LS+KLU (single)                           6.1e-05        0                 9.54e-07
-LS+KLU                                    6.1e-05        0                 9.54e-07
-LS+NICSLU (single)                        6.1e-05        0                 9.54e-07
-LS+NICSLU                                 6.1e-05        0                 9.54e-07
+GS                                        6.1e-05        3.81e-06          1.53e-05
+GS synch                                  6.1e-05        3.81e-06          1.53e-05
+NR single (SLU)                           6.1e-05        0                 9.54e-07
+NR (SLU)                                  6.1e-05        0                 9.54e-07
+NR single (KLU)                           6.1e-05        0                 9.54e-07
+NR (KLU)                                  6.1e-05        0                 9.54e-07
+NR single (NICSLU *)                      6.1e-05        0                 9.54e-07
+NR (NICSLU *)                             6.1e-05        0                 9.54e-07
+NR single (CKTSO *)                       6.1e-05        0                 9.54e-07
+NR (CKTSO *)                              6.1e-05        0                 9.54e-07
+FDPF XB (SLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (SLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (KLU)                             6.1e-05        1.91e-06          1.53e-05
+FDPF BX (KLU)                             6.1e-05        1.91e-06          7.63e-06
+FDPF XB (NICSLU *)                        6.1e-05        1.91e-06          1.53e-05
+FDPF BX (NICSLU *)                        6.1e-05        1.91e-06          7.63e-06
+FDPF XB (CKTSO *)                         6.1e-05        1.91e-06          1.53e-05
+FDPF BX (CKTSO *)                         6.1e-05        1.91e-06          7.63e-06
 =================================  ==============  ==============  ================
 
 As we can see on all the tables above, the difference when using lightsim and pandapower is rather
@@ -255,8 +303,24 @@ When using Newton Raphson solvers, the difference in absolute values when using
 with using PandaPowerBackend is neglectible: less than 1e-06 in all cases (and 0.00 when comparing the
 flows on the powerline for both environments).
 
-Other benchmarks
-----------------------------
+Other benchmark
+----------------
+
+We have at our disposal different computers with different software / hardware.
+
+From time to time, we benchmark grid2op and lightsim2grid. 
+The results can be found in:
+
+.. toctree::
+  :maxdepth: 1
+  :caption: For a laptop with a i7 of 2015 wth a frequency of 2.70 GHz
+
+  benchmark_solver/ubuntu_2004_dell/ls0.8.1_glop1.10.0
+  benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.10.0
+  benchmark_solver/ubuntu_2004_dell/ls0.8.0_glop1.9.8
+
+Benchmarks of other lightsim2grid functions
+--------------------------------------------
 
 With lightsim2grid 0.5.5 some new feature has been introduced, which are the "security analysis" and the "comptuation 
 of time series". 
diff --git a/docs/conf.py b/docs/conf.py
index c50cbef..3e4db7b 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -22,7 +22,7 @@
 author = 'Benjamin DONNOT'
 
 # The full version, including alpha/beta/rc tags
-release = "0.8.0"
+release = "0.8.1"
 version = '0.8'
 
 # -- General configuration ---------------------------------------------------
diff --git a/docs/index.rst b/docs/index.rst
index f13c36d..66f4dcd 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -30,6 +30,7 @@ As from version 0.5.3:
    use_with_grid2op
    benchmarks
    use_solver
+   rewards
    physical_law_checker
 
 
diff --git a/docs/install_from_source.rst b/docs/install_from_source.rst
index 10e8447..c44a00a 100644
--- a/docs/install_from_source.rst
+++ b/docs/install_from_source.rst
@@ -195,5 +195,3 @@ Now you simply need to install the lightsim2grid package this way, like any pyth
     pip install -U .
 
 And you are done :-)
-
-
diff --git a/docs/quickstart.rst b/docs/quickstart.rst
index 46b4ba4..7ab2224 100644
--- a/docs/quickstart.rst
+++ b/docs/quickstart.rst
@@ -34,7 +34,7 @@ this way:
     from grid2op.Agent import RandomAgent
 
     # create an environment
-    env_name = "rte_case14_realistic"  # for example, other environments might be usable
+    env_name = "l2rpn_case14_sandbox"  # for example, other environments might be usable
     env = grid2op.make(env_name,
                        backend=LightSimBackend()  # this is the only change you have to make!
                        )
diff --git a/docs/rewards.rst b/docs/rewards.rst
new file mode 100644
index 0000000..db2593b
--- /dev/null
+++ b/docs/rewards.rst
@@ -0,0 +1,13 @@
+Custom Rewards (doc in progress)
+=======================================
+
+Detailed usage
+--------------------------
+
+.. automodule:: lightsim2grid.rewards
+    :members:
+    :autosummary:
+
+* :ref:`genindex`
+* :ref:`modindex`
+* :ref:`search`
\ No newline at end of file
diff --git a/docs/security_analysis.rst b/docs/security_analysis.rst
index 9e568c2..a652c62 100644
--- a/docs/security_analysis.rst
+++ b/docs/security_analysis.rst
@@ -1,4 +1,4 @@
-Security Analysis (doc in progress)
+Contingency Analysis (doc in progress)
 =======================================
 
 The documentation of this section is in progress. It is rather incomplete for the moment, and only expose the most
@@ -52,8 +52,8 @@ to / from other data sources.
 
 .. _sa_benchmarks:
 
-Benchmarks
------------------
+Benchmarks (Contingency Analysis)
+----------------------------------
 
 Here are some benchmarks made with:
 
diff --git a/docs/time_series.rst b/docs/time_series.rst
index 426e95f..94a97be 100644
--- a/docs/time_series.rst
+++ b/docs/time_series.rst
@@ -60,8 +60,8 @@ Importantly, this method is around **13x** faster than simulating "do nothing" (
 
 .. _ts_benchmarks:
 
-Benchmarks
------------------
+Benchmarks (Time Series)
+-------------------------
 
 Here are some benchmarks made with:
 
diff --git a/lightsim2grid/__init__.py b/lightsim2grid/__init__.py
index ec429e7..85ece59 100644
--- a/lightsim2grid/__init__.py
+++ b/lightsim2grid/__init__.py
@@ -1,13 +1,13 @@
-# Copyright (c) 2020, RTE (https://www.rte-france.com)
+# Copyright (c) 2020-2024, RTE (https://www.rte-france.com)
 # See AUTHORS.txt
 # This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 # If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
 # you can obtain one at http://mozilla.org/MPL/2.0/.
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
-__version__ = "0.8.0"
+__version__ = "0.8.1"
 
-__all__ = ["newtonpf", "SolverType", "ErrorType", "solver"]
+__all__ = ["newtonpf", "SolverType", "ErrorType", "solver", "compilation_options"]
 
 # import directly from c++ module
 from lightsim2grid.solver import SolverType
@@ -46,3 +46,11 @@
     # grid2op is not installed, the SecurtiyAnalysis module will not be available
     pass
     print(f"ContingencyAnalysis import error: {exc_}")
+    
+try:
+    from lightsim2grid.rewards import *
+    __all__.append("rewards")
+except ImportError as exc_:
+    # grid2op is not installed, the SecurtiyAnalysis module will not be available
+    pass
+    print(f"rewards import error: {exc_}")
diff --git a/lightsim2grid/compilation_options/__init__.py b/lightsim2grid/compilation_options/__init__.py
new file mode 100644
index 0000000..cdcb8df
--- /dev/null
+++ b/lightsim2grid/compilation_options/__init__.py
@@ -0,0 +1,37 @@
+# Copyright (c) 2024, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+__all__ = ["klu_solver_available",
+           "nicslu_solver_available",
+           "cktso_solver_available",
+           "compiled_march_native",
+           "compiled_o3_optim",
+           "version",
+           "lightsim2grid_lib"]
+
+from lightsim2grid_cpp import klu_solver_available
+from lightsim2grid_cpp import nicslu_solver_available
+from lightsim2grid_cpp import cktso_solver_available
+from lightsim2grid_cpp import compiled_march_native
+from lightsim2grid_cpp import compiled_o3_optim
+from lightsim2grid_cpp import version
+from lightsim2grid_cpp import __file__ as lightsim2grid_lib
+
+try:
+    from lightsim2grid_cpp import nicslu_lib
+    __all__.append("nicslu_lib")
+except ImportError :
+    # NICSLU linear solver is not available
+    pass
+
+try:
+    from lightsim2grid_cpp import cktso_lib
+    __all__.append("cktso_lib")
+except ImportError :
+    # CKTSO linear solver is not available
+    pass
diff --git a/lightsim2grid/contingencyAnalysis.py b/lightsim2grid/contingencyAnalysis.py
index f16a504..fcd0066 100644
--- a/lightsim2grid/contingencyAnalysis.py
+++ b/lightsim2grid/contingencyAnalysis.py
@@ -6,21 +6,25 @@
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
 
-__all__ = ["ContingencyAnalysisCPP", "ContingencyAnalysis",
-           # deprecated
-           "SecurityAnalysisCPP", "SecurityAnalysis",
-           ]
+__all__ = ["ContingencyAnalysisCPP"]
 
 import copy
 import numpy as np
 from collections.abc import Iterable
 
-from lightsim2grid.lightSimBackend import LightSimBackend
 from lightsim2grid.solver import SolverType
 from lightsim2grid_cpp import ContingencyAnalysisCPP
 
+try:
+    from lightsim2grid.lightSimBackend import LightSimBackend
+    __all__.append("ContingencyAnalysis")
+    GRID2OP_INSTALLED = True
+except ImportError as exc_:
+    # grid2op is not installed
+    GRID2OP_INSTALLED = False
 
-class ContingencyAnalysis(object):
+
+class __ContingencyAnalysis(object):
     """
     This class allows to perform a "security analysis" from a given grid state.
 
@@ -76,10 +80,26 @@ class ContingencyAnalysis(object):
     STR_TYPES = (str, np.str_)  # np.str deprecated in numpy 1.20 and earlier versions not supported anyway
         
     def __init__(self, grid2op_env):
-        if not isinstance(grid2op_env.backend, LightSimBackend):
-            raise RuntimeError("This class only works with LightSimBackend")
-        self.grid2op_env = grid2op_env.copy()
-        self.computer = ContingencyAnalysisCPP(self.grid2op_env.backend._grid)
+        if not GRID2OP_INSTALLED:
+            raise RuntimeError("Impossible to use the python wrapper `ContingencyAnalysis` "
+                               "when grid2op is not installed. Please fall back to the "
+                               "c++ version (available in python) with:\n"
+                               "\tfrom lightsim2grid.contingencyAnalysis import ContingencyAnalysisCPP\n"
+                               "and refer to the appropriate documentation.")
+        from grid2op.Environment import Environment
+        if isinstance(grid2op_env, Environment):    
+            if not isinstance(grid2op_env.backend, LightSimBackend):
+                raise RuntimeError("This class only works with LightSimBackend")
+            self._ls_backend = grid2op_env.backend.copy()
+        elif isinstance(grid2op_env, LightSimBackend):
+            if hasattr(grid2op_env, "_is_loaded") and not grid2op_env._is_loaded:
+                raise RuntimeError("Impossible to init a `ContingencyAnalysis` "
+                                   "with a backend that has not been initialized.")
+            self._ls_backend = grid2op_env.copy()
+        else:
+            raise RuntimeError("`ContingencyAnalysis` can only be created "
+                               "with a grid2op `Environment` or a `LightSimBackend`")
+        self.computer = ContingencyAnalysisCPP(self._ls_backend._grid)
         self._contingency_order = {}  # key: contingency (as tuple), value: order in which it is entered
         self._all_contingencies = []
         self.__computed = False
@@ -100,16 +120,29 @@ def all_contingencies(self, val):
         raise RuntimeError("Impossible to add new topologies like this. Please use `add_single_contingency` "
                            "or `add_multiple_contingencies`.")
 
-    def clear(self):
+    # TODO implement that !
+    def __update_grid(self, backend_act):
+        self.clear(with_contlist=False)
+        self._ls_backend.apply_action(backend_act)
+        # run the powerflow
+        self._ls_backend.runpf()
+        # update the computer
+        # self.computer = ContingencyAnalysisCPP(self._ls_backend._grid)
+        # self.computer.update_grid(...)  # not implemented
+        
+    def clear(self, with_contlist=True):
         """
         Clear the list of contingencies to simulate
         """
-        self.computer.clear()
-        self._contingency_order = {}
         self.__computed = False
         self._vs = None
         self._ampss = None
-        self._all_contingencies = []
+        if with_contlist:
+            self.computer.clear()
+            self._contingency_order = {}
+            self._all_contingencies = []
+        else:
+            self.computer.clear_results_only()
 
     def _single_cont_to_li_int(self, single_cont):
         li_disc = []
@@ -118,7 +151,7 @@ def _single_cont_to_li_int(self, single_cont):
 
         for stuff in single_cont:
             if isinstance(stuff, type(self).STR_TYPES):
-                stuff = np.where(self.grid2op_env.name_line == stuff)
+                stuff = np.where(type(self._ls_backend).name_line == stuff)
                 stuff = stuff[0]
                 if stuff.size == 0:
                     # name is not found
@@ -233,7 +266,7 @@ def add_all_n1_contingencies(self):
             for single_cont_id in range(env.n_line):
                 self.add_single_contingency(single_cont_id)
         """
-        for single_cont_id in range(self.grid2op_env.n_line):
+        for single_cont_id in range(type(self._ls_backend).n_line):
             self.add_single_contingency(single_cont_id)
 
     def get_flows(self, *args):
@@ -302,10 +335,10 @@ def compute_V(self):
             been entered. Please use `get_flows()` method for easier reading back of the results
 
         """
-        v_init = self.grid2op_env.backend.V
+        v_init = 1. * self._ls_backend.V
         self.computer.compute(v_init,
-                              self.grid2op_env.backend.max_it,
-                              self.grid2op_env.backend.tol)
+                              self._ls_backend.max_it,
+                              self._ls_backend.tol)
         self._vs = self.computer.get_voltages()
         self.__computed = True
         return self._vs
@@ -342,6 +375,10 @@ def compute_P(self):
 
     def close(self):
         """permanently close the object"""
-        self.grid2op_env.close()
+        self._ls_backend.close()
         self.clear()
         self.computer.close()
+        
+        
+if GRID2OP_INSTALLED:
+    ContingencyAnalysis = __ContingencyAnalysis
diff --git a/lightsim2grid/gridmodel/from_pypowsybl.py b/lightsim2grid/gridmodel/from_pypowsybl.py
index 6501ead..798742c 100644
--- a/lightsim2grid/gridmodel/from_pypowsybl.py
+++ b/lightsim2grid/gridmodel/from_pypowsybl.py
@@ -78,8 +78,8 @@ def init(net : pypo.network,
     # to handle encoding in 32 bits and overflow when "splitting" the Q values among 
     min_q = df_gen["min_q"].values.astype(np.float32)
     max_q = df_gen["max_q"].values.astype(np.float32)
-    min_q[~np.isfinite(min_q)] = np.finfo(np.float32).min * 0.5 + 1.
-    max_q[~np.isfinite(max_q)] = np.finfo(np.float32).max * 0.5 - 1.
+    min_q[~np.isfinite(min_q)] = np.finfo(np.float32).min * 1e-4 + 1.
+    max_q[~np.isfinite(max_q)] = np.finfo(np.float32).max * 1e-4 - 1.
     gen_bus, gen_disco = _aux_get_bus(bus_df, df_gen)
     model.init_generators_full(df_gen["target_p"].values,
                                df_gen["target_v"].values / voltage_levels.loc[df_gen["voltage_level_id"].values]["nominal_v"].values,
diff --git a/lightsim2grid/gridmodel/initGridModel.py b/lightsim2grid/gridmodel/initGridModel.py
index 78293cd..220f4e1 100644
--- a/lightsim2grid/gridmodel/initGridModel.py
+++ b/lightsim2grid/gridmodel/initGridModel.py
@@ -15,6 +15,9 @@
 import numpy as np
 from numbers import Number
 import warnings
+
+import pandapower
+import lightsim2grid
 from lightsim2grid_cpp import GridModel, PandaPowerConverter
 from lightsim2grid.gridmodel._aux_add_sgen import _aux_add_sgen
 from lightsim2grid.gridmodel._aux_add_load import _aux_add_load
@@ -28,7 +31,7 @@
 from lightsim2grid.gridmodel._aux_add_dc_line import _aux_add_dc_line
 
 
-def init(pp_net):
+def init(pp_net: "pandapower.auxiliary.pandapowerNet") -> GridModel:
     """
     Convert a pandapower network as input into a GridModel.
 
@@ -53,7 +56,7 @@ def init(pp_net):
 
     Parameters
     ----------
-    pp_net: :class:`pandapower.grid`
+    pp_net: :class:`pandapower.auxiliary.pandapowerNet`
         The initial pandapower network you want to convert
 
     Returns
diff --git a/lightsim2grid/lightSimBackend.py b/lightsim2grid/lightSimBackend.py
index a7c6319..7119b80 100644
--- a/lightsim2grid/lightSimBackend.py
+++ b/lightsim2grid/lightSimBackend.py
@@ -6,13 +6,21 @@
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
 
+import os
 import copy
-from typing import Optional, Union
+from typing import Tuple, Optional, Any, Dict, Union
+try:
+    from typing import Self
+except ImportError:
+    # python version is probably bellow 3.11
+    from typing_extensions import Self
+    
 import warnings
 import numpy as np
 import pandas as pd
 import time
 
+import grid2op
 from grid2op.Action import CompleteAction
 from grid2op.Backend import Backend
 from grid2op.Exceptions import BackendError
@@ -22,6 +30,12 @@
 except ImportError as exc_:
     from grid2op.Action._BackendAction import _BackendAction
 
+try:
+    from grid2op.Space import DEFAULT_N_BUSBAR_PER_SUB
+except ImportError:
+    # for backward compatibility with grid2op <= 1.9.8
+    DEFAULT_N_BUSBAR_PER_SUB = 2
+    
 try:
     from typing import Literal
 except ImportError:
@@ -94,6 +108,10 @@ def __init__(self,
                              loader_kwargs,
                              stop_if_load_disco,
                              stop_if_gen_disco)
+        
+        # backward compat: need to define it if not done by grid2op
+        if not hasattr(self, "_can_be_copied"):
+            self._can_be_copied = can_be_copied
 
         #: .. versionadded:: 0.8.0
         #: Which type of grid format can be read by your backend.
@@ -130,6 +148,9 @@ def __init__(self,
         self._timer_preproc = 0.
         self._timer_postproc = 0.
         self._timer_solver = 0.
+        self._timer_read_data_back = 0.
+        self._timer_fetch_data_cpp = 0.
+        self._timer_apply_act = 0.
 
         self.next_prod_p = None  # this vector is updated with the action that will modify the environment
         # it is done to keep track of the redispatching
@@ -175,6 +196,7 @@ def __init__(self,
         self.sh_p = None
         self.sh_q = None
         self.sh_v = None
+        self.sh_theta = None
         self.sh_bus = None
 
         # voltage angle
@@ -201,14 +223,17 @@ def __init__(self,
         # it takes only into account the time spend in the powerflow algorithm
         self.comp_time = 0.  
 
-        # computation time of the powerflow
-        # it takes into account everything in the gridmodel, including the mapping 
-        # to the solver, building of Ybus and Sbus AND the time to solve the powerflow
+        #: computation time of the powerflow
+        #: it takes into account everything in the gridmodel, including the mapping 
+        #: to the solver, building of Ybus and Sbus AND the time to solve the powerflow
         self.timer_gridmodel_xx_pf = 0.
         
         self._timer_postproc = 0.
         self._timer_preproc = 0.
         self._timer_solver = 0.
+        self._timer_read_data_back = 0.
+        self._timer_fetch_data_cpp = 0.
+        self._timer_apply_act = 0.
 
         # hack for the storage unit:
         # in grid2op, for simplicity, I suppose that if a storage is alone on a busbar, and
@@ -221,6 +246,17 @@ def __init__(self,
         # backend SHOULD not do these kind of stuff
         self._idx_hack_storage = []
         
+        # speed optimization
+        self._lineor_res = None
+        self._lineex_res = None
+        self._load_res = None
+        self._gen_res = None
+        self._shunt_res = None
+        self._trafo_hv_res = None
+        self._trafo_lv_res = None
+        self._storage_res = None
+        self._reset_res_pointers()
+        
     def _aux_init_super(self, 
                         detailed_infos_for_cascading_failures,
                         can_be_copied,
@@ -255,7 +291,11 @@ def _aux_init_super(self,
                           "you cannot set max_iter, tol nor solver_type arguments.")
             Backend.__init__(self,
                              detailed_infos_for_cascading_failures=detailed_infos_for_cascading_failures)
-                 
+            
+        if hasattr(type(self), "can_handle_more_than_2_busbar"):
+            # do not forget to propagate this if needed
+            self.can_handle_more_than_2_busbar()
+            
     def turnedoff_no_pv(self):
         self._turned_off_pv = False
         self._grid.turnedoff_no_pv()
@@ -263,25 +303,8 @@ def turnedoff_no_pv(self):
     def turnedoff_pv(self):
         self._turned_off_pv = True
         self._grid.turnedoff_pv()
-        
-    def _fill_theta(self):
-        # line_or_theta = np.empty(self.n_line)
-        self.line_or_theta[:self.__nb_powerline] = self._grid.get_lineor_theta()
-        self.line_or_theta[self.__nb_powerline:] = self._grid.get_trafohv_theta()
 
-        # line_ex_theta = np.empty(self.n_line)
-        self.line_ex_theta[:self.__nb_powerline] = self._grid.get_lineex_theta()
-        self.line_ex_theta[self.__nb_powerline:] = self._grid.get_trafolv_theta()
-
-        # line_or_theta = np.concatenate((self._grid.get_lineor_theta(), self._grid.get_trafohv_theta()))
-        # line_ex_theta = np.concatenate((self._grid.get_lineex_theta(), self._grid.get_trafolv_theta()))
-        self.load_theta[:] = self._grid.get_load_theta()
-        self.gen_theta[:] = self._grid.get_gen_theta()
-
-        if self.__has_storage:
-            self.storage_theta[:] = self._grid.get_storage_theta()
-
-    def get_theta(self):
+    def get_theta(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
         """
 
         Returns
@@ -329,7 +352,7 @@ def get_solver_types(self) -> Union[SolverType, SolverType]:
         """
         return self._grid.get_solver_type(), self._grid.get_dc_solver_type()
         
-    def set_solver_type(self, solver_type: SolverType):
+    def set_solver_type(self, solver_type: SolverType) -> None:
         """
         Change the type of solver you want to use.
 
@@ -371,7 +394,7 @@ def _check_suitable_solver_type(self, solver_type, check_in_avail_solver=True):
             raise BackendError(f"The solver type provided \"{solver_type}\" is not available on your system. Available"
                                f"solvers are {self.available_solvers}")
             
-    def set_solver_max_iter(self, max_iter):
+    def set_solver_max_iter(self, max_iter: int) -> None:
         """
         Set the maximum number of iteration the solver is allowed to perform.
 
@@ -404,7 +427,7 @@ def set_solver_max_iter(self, max_iter):
             raise BackendError("max_iter should be a strictly positive integer (integer >= 1)")
         self.max_it = max_iter
 
-    def set_tol(self, new_tol):
+    def set_tol(self, new_tol: float) -> None:
         """
         Set the tolerance of the powerflow. This means that the powerflow will stop when the Kirchhoff's Circuit Laws
         are met up to a tolerance of "new_tol".
@@ -454,7 +477,9 @@ def _assign_right_solver(self):
             else:
                 self._grid.change_solver(SolverType.SparseLU)
                 
-    def load_grid(self, path=None, filename=None):    
+    def load_grid(self,
+                  path : Union[os.PathLike, str],
+                  filename : Optional[Union[os.PathLike, str]]=None) -> None: 
         if hasattr(type(self), "can_handle_more_than_2_busbar"):
             # grid2op version >= 1.10.0 then we use this
             self.can_handle_more_than_2_busbar()
@@ -466,6 +491,7 @@ def load_grid(self, path=None, filename=None):
         else:
             raise BackendError(f"Impossible to initialize the backend with '{self._loader_method}'")
         self._grid.tell_solver_need_reset()
+        self._reset_res_pointers()  # force the re reading of the accessors at next powerflow
     
     def _should_not_have_to_do_this(self, path=None, filename=None):
         # included in grid2op now !
@@ -600,6 +626,7 @@ def _load_grid_pypowsybl(self, path=None, filename=None):
             self.name_line = np.concatenate((lor_sub.index, tor_sub.index))
             self.name_storage = np.array(batt_sub.index)
             self.name_shunt = np.array(sh_sub.index)
+            self.name_sub = np.array(buses_sub_id.index)
         
         if "reconnect_disco_gen" in loader_kwargs and loader_kwargs["reconnect_disco_gen"]:
             for el in self._grid.get_generators():
@@ -663,11 +690,45 @@ def _aux_setup_right_after_grid_init(self):
             self._grid._max_nb_bus_per_sub = self.n_busbar_per_sub
             
         self._grid.tell_solver_need_reset()
-            
+    
+    def init_from_loaded_pandapower(self, pp_net):
+        if hasattr(type(self), "can_handle_more_than_2_busbar"):
+            type(self.init_pp_backend).n_busbar_per_sub = self.n_busbar_per_sub
+        self.init_pp_backend = pp_net.copy()
+        self._aux_init_pandapower()
+        
+        # handles redispatching
+        if type(pp_net).redispatching_unit_commitment_availble:
+            self.redispatching_unit_commitment_availble = True
+            for attr_nm in ["gen_type", "gen_pmin", "gen_pmax",
+                            "gen_redispatchable", "gen_max_ramp_up",
+                            "gen_max_ramp_down", "gen_min_uptime",
+                            "gen_min_downtime", "gen_cost_per_MW",
+                            "gen_startup_cost", "gen_shutdown_cost",
+                            "gen_renewable"
+                            ]:
+                setattr(self, attr_nm, copy.deepcopy(getattr( type(pp_net), attr_nm)))
+        
+        # handles storages
+        for attr_nm in ["storage_type",
+                        "storage_Emax",
+                        "storage_Emin",
+                        "storage_max_p_prod" ,
+                        "storage_max_p_absorb",
+                        "storage_marginal_cost",
+                        "storage_loss",
+                        "storage_charging_efficiency",
+                        "storage_discharging_efficiency",
+                        ]:
+            setattr(self, attr_nm, copy.deepcopy(getattr( type(pp_net), attr_nm)))
+        
     def _load_grid_pandapower(self, path=None, filename=None):
         if hasattr(type(self), "can_handle_more_than_2_busbar"):
             type(self.init_pp_backend).n_busbar_per_sub = self.n_busbar_per_sub
         self.init_pp_backend.load_grid(path, filename)
+        self._aux_init_pandapower()
+    
+    def _aux_init_pandapower(self):
         self.can_output_theta = True  # i can compute the "theta" and output it to grid2op
 
         self._grid = init(self.init_pp_backend._grid)    
@@ -815,12 +876,13 @@ def _aux_finish_setup_after_reading(self):
 
         self.topo_vect = np.ones(cls.dim_topo, dtype=dt_int).reshape(-1)
         
-        if type(self).shunts_data_available:
+        if cls.shunts_data_available:
             self.shunt_topo_vect = np.ones(cls.n_shunt, dtype=dt_int)
              # shunts
             self.sh_p = np.full(cls.n_shunt, dtype=dt_float, fill_value=np.NaN).reshape(-1)
             self.sh_q = np.full(cls.n_shunt, dtype=dt_float, fill_value=np.NaN).reshape(-1)
             self.sh_v = np.full(cls.n_shunt, dtype=dt_float, fill_value=np.NaN).reshape(-1)
+            self.sh_theta = np.full(cls.n_shunt, dtype=dt_float, fill_value=np.NaN).reshape(-1)
             self.sh_bus = np.full(cls.n_shunt, dtype=dt_int, fill_value=-1).reshape(-1)
 
         self.p_or = np.full(cls.n_line, dtype=dt_float, fill_value=np.NaN).reshape(-1)
@@ -857,8 +919,9 @@ def _aux_finish_setup_after_reading(self):
         self.__me_at_init = self._grid.copy()
         self.__init_topo_vect = np.ones(cls.dim_topo, dtype=dt_int)
         self.__init_topo_vect[:] = self.topo_vect
+        self.sh_bus[:] = 1
         
-    def assert_grid_correct_after_powerflow(self):
+    def assert_grid_correct_after_powerflow(self) -> None:
         """
         This method is called by the environment. It ensure that the backend remains consistent even after a powerflow
         has be run with :func:`Backend.runpf` method.
@@ -926,14 +989,15 @@ def _count_object_per_bus(self):
             arr_ = self.shunt_to_subid[is_connected] + self.__nb_bus_before * (arr_[is_connected] - 1)
             self.nb_obj_per_bus[arr_] += 1
 
-    def close(self):
+    def close(self) -> None:
         self.init_pp_backend.close()
         self._grid = None
 
-    def apply_action(self, backendAction):
+    def apply_action(self, backendAction: Union["grid2op.Action._backendAction._BackendAction", None]) -> None:
         """
         Specific implementation of the method to apply an action modifying a powergrid in the pandapower format.
         """
+        tick = time.perf_counter()
         active_bus, *_, topo__, shunts__ = backendAction()
 
         # change the overall topology
@@ -977,18 +1041,41 @@ def apply_action(self, backendAction):
                         self._grid.change_bus_shunt(sh_id, type(self).local_bus_to_global_int(new_bus, self.shunt_to_subid[sh_id]))
                     else:
                         self._grid.change_bus_shunt(sh_id, self.shunt_to_subid[sh_id] + (new_bus == 2) * type(self).n_sub)
+            # remember the topology not to need to read it back from the grid
+            self.sh_bus[shunt_bus.changed] = shunt_bus.values[shunt_bus.changed]
             for sh_id, new_p in shunt_p:
                 self._grid.change_p_shunt(sh_id, new_p)
             for sh_id, new_q in shunt_q:
                 self._grid.change_q_shunt(sh_id, new_q)
         
         self._handle_dist_slack()
-    
+        self._timer_apply_act += time.perf_counter() - tick
+        
     def _handle_dist_slack(self):
         if self._dist_slack_non_renew:
             self._grid.update_slack_weights(type(self).gen_redispatchable)
 
-    def runpf(self, is_dc=False):
+    def _fetch_grid_data(self):
+        beg_test = time.perf_counter()
+        if self._lineor_res is None:
+            self._lineor_res = self._grid.get_lineor_res_full()
+        if self._lineex_res is None:
+            self._lineex_res = self._grid.get_lineex_res_full()
+        if self._load_res is None:
+            self._load_res = self._grid.get_loads_res_full()
+        if self._gen_res is None:
+            self._gen_res = self._grid.get_gen_res_full()
+        if self._trafo_hv_res is None:
+            self._trafo_hv_res = self._grid.get_trafohv_res_full()
+        if self._trafo_lv_res is None:
+            self._trafo_lv_res = self._grid.get_trafolv_res_full()
+        if self._storage_res is None:
+            self._storage_res = self._grid.get_storages_res_full()
+        if self._shunt_res is None:
+            self._shunt_res = self._grid.get_shunts_res_full()
+        self._timer_fetch_data_cpp += time.perf_counter() - beg_test
+            
+    def runpf(self, is_dc : bool=False) -> Tuple[bool, Union[Exception, None]]:
         my_exc_ = None
         res = False
         try:
@@ -1012,9 +1099,7 @@ def runpf(self, is_dc=False):
                     self._grid.deactivate_result_computation()
                     # if I init with dc values, it should depends on previous state
                     self.V[:] = self._grid.get_init_vm_pu()  # see issue 30
-                    # print(f"{self.V[:14] = }")
                     Vdc = self._grid.dc_pf(copy.deepcopy(self.V), self.max_it, self.tol)
-                    # print(f"{Vdc[:14] = }")
                     self._grid.reactivate_result_computation()
                     if Vdc.shape[0] == 0:
                         raise BackendError(f"Divergence of DC powerflow (non connected grid) at the initialization of AC powerflow. Detailed error: {self._grid.get_dc_solver().get_error()}")
@@ -1043,24 +1128,31 @@ def runpf(self, is_dc=False):
                 
                 self.timer_gridmodel_xx_pf += self._grid.timer_last_ac_pf
                 # timer_gridmodel_xx_pf takes all the time within the gridmodel "ac_pf"
-                   
+            
+            beg_readback = time.perf_counter()
             self.V[:] = V
+            self._fetch_grid_data()
+        
             (self.p_or[:self.__nb_powerline],
              self.q_or[:self.__nb_powerline],
              self.v_or[:self.__nb_powerline],
-             self.a_or[:self.__nb_powerline]) = self._grid.get_lineor_res()
+             self.a_or[:self.__nb_powerline],
+             self.line_or_theta[:self.__nb_powerline]) = self._lineor_res
             (self.p_or[self.__nb_powerline:],
              self.q_or[self.__nb_powerline:],
              self.v_or[self.__nb_powerline:],
-             self.a_or[self.__nb_powerline:]) = self._grid.get_trafohv_res()
+             self.a_or[self.__nb_powerline:],
+             self.line_or_theta[self.__nb_powerline:]) = self._trafo_hv_res
             (self.p_ex[:self.__nb_powerline],
              self.q_ex[:self.__nb_powerline],
              self.v_ex[:self.__nb_powerline],
-             self.a_ex[:self.__nb_powerline]) = self._grid.get_lineex_res()
+             self.a_ex[:self.__nb_powerline],
+             self.line_ex_theta[:self.__nb_powerline]) = self._lineex_res
             (self.p_ex[self.__nb_powerline:],
              self.q_ex[self.__nb_powerline:],
              self.v_ex[self.__nb_powerline:],
-             self.a_ex[self.__nb_powerline:]) = self._grid.get_trafolv_res()
+             self.a_ex[self.__nb_powerline:],
+             self.line_ex_theta[self.__nb_powerline:]) = self._trafo_lv_res
 
             self.a_or *= 1000.  # kA in lightsim, A expected in grid2op
             self.a_ex *= 1000.  # kA in lightsim, A expected in grid2op
@@ -1070,11 +1162,12 @@ def runpf(self, is_dc=False):
             self.a_ex[~np.isfinite(self.a_ex)] = 0.
             self.v_ex[~np.isfinite(self.v_ex)] = 0.
 
-            self.load_p[:], self.load_q[:], self.load_v[:] = self._grid.get_loads_res()
-            self.prod_p[:], self.prod_q[:], self.prod_v[:] = self._grid.get_gen_res()
+            self.load_p[:], self.load_q[:], self.load_v[:], self.load_theta[:] = self._load_res
+            self.prod_p[:], self.prod_q[:], self.prod_v[:], self.gen_theta[:] = self._gen_res
             if self.__has_storage:
-                self.storage_p[:], self.storage_q[:], self.storage_v[:] = self._grid.get_storages_res()
+                self.storage_p[:], self.storage_q[:], self.storage_v[:], self.storage_theta[:] = self._storage_res
             self.storage_v[self.storage_v == -1.] = 0.  # voltage is 0. for disconnected elements in grid2op
+            self._timer_read_data_back +=  time.perf_counter() - beg_readback
             
             self.next_prod_p[:] = self.prod_p
 
@@ -1092,8 +1185,6 @@ def runpf(self, is_dc=False):
 
             if type(self).shunts_data_available:
                 self._set_shunt_info()
-            
-            self._fill_theta()
 
             if (self.line_or_theta >= 1e6).any() or (self.line_ex_theta >= 1e6).any():
                 raise BackendError(f"Some theta are above 1e6 which should not be happening !")
@@ -1111,6 +1202,7 @@ def runpf(self, is_dc=False):
             if is_dc:
                 # set back the solver to its previous state
                 self._grid.change_solver(self.__current_solver_type)
+            self._grid.tell_solver_need_reset()
 
         # TODO grid2op compatibility ! (was a single returned element before storage were introduced)
         if self.__has_storage:
@@ -1143,6 +1235,7 @@ def _fill_nans(self):
             self.sh_p[:] = np.NaN
             self.sh_q[:] = np.NaN
             self.sh_v[:] = np.NaN
+            self.sh_theta[:] = np.NaN
             self.sh_bus[:] = -1
 
         if self.__has_storage:
@@ -1151,13 +1244,24 @@ def _fill_nans(self):
             self.storage_v[:] = np.NaN
             self.storage_theta[:] = np.NaN
         self.V[:] = self._grid.get_init_vm_pu()  # reset the V to its "original" value (see issue 30)
-
+        self._reset_res_pointers()
+    
+    def _reset_res_pointers(self):
+        self._lineor_res  = None
+        self._lineex_res = None
+        self._load_res = None
+        self._gen_res = None
+        self._trafo_hv_res = None
+        self._trafo_lv_res = None
+        self._storage_res = None
+        self._shunt_res = None
+        
     def __deepcopy__(self, memo):
         result = self.copy()
         memo[id(self)] = result
         return result
 
-    def copy(self):
+    def copy(self) -> Self:
         # i can perform a regular copy, everything has been initialized
         mygrid = self._grid
         __me_at_init = self.__me_at_init
@@ -1188,6 +1292,11 @@ def copy(self):
                             self._loader_kwargs,
                             self._stop_if_load_disco,
                             self._stop_if_gen_disco)
+        
+        # for backward compat (attribute was not necessarily present in early grid2op)
+        if not hasattr(res, "_can_be_copied"):
+            res._can_be_copied = self._can_be_copied
+            
         res.comp_time = self.comp_time
         res.timer_gridmodel_xx_pf = self.timer_gridmodel_xx_pf
 
@@ -1202,10 +1311,12 @@ def copy(self):
                            "_big_topo_to_obj", "dim_topo",
                            "_idx_hack_storage",
                            "_timer_preproc", "_timer_postproc", "_timer_solver",
+                           "_timer_read_data_back", "_timer_apply_act",
                            "supported_grid_format", 
                            "max_it", "tol", "_turned_off_pv", "_dist_slack_non_renew",
                            "_use_static_gen", "_loader_method", "_loader_kwargs",
-                           "_stop_if_load_disco", "_stop_if_gen_disco"
+                           "_stop_if_load_disco", "_stop_if_gen_disco",
+                           "_timer_fetch_data_cpp"
                            ]
         for attr_nm in li_regular_attr:
             if hasattr(self, attr_nm):
@@ -1223,7 +1334,7 @@ def copy(self):
                        "load_p", "load_q", "load_v",
                        "prod_p", "prod_q", "prod_v",
                        "storage_p", "storage_q", "storage_v",
-                       "sh_p", "sh_q", "sh_v", "sh_bus",
+                       "sh_p", "sh_q", "sh_v", "sh_bus", "sh_theta",
                        "line_or_theta", "line_ex_theta", "load_theta", "gen_theta", "storage_theta",                   
                        ]
         for attr_nm in li_attr_npy:
@@ -1265,63 +1376,70 @@ def copy(self):
         res._backend_action_class = self._backend_action_class  # this is const
         res.__init_topo_vect = self.__init_topo_vect
         res.available_solvers = self.available_solvers
-
+        res._reset_res_pointers()
+        res._fetch_grid_data()
+        
+        # assign back "self" attributes
         self._grid = mygrid
         self.init_pp_backend = inippbackend
         self.__me_at_init = __me_at_init
         return res
 
-    def get_line_status(self):
+    def get_line_status(self) -> np.ndarray:
         l_s = self._grid.get_lines_status()
         t_s = self._grid.get_trafo_status()
         return np.concatenate((l_s, t_s)).astype(dt_bool)
 
-    def get_line_flow(self):
+    def get_line_flow(self) -> np.ndarray:
         return self.a_or
+    
+    def get_topo_vect(self) -> np.ndarray:
+        return 1 * self.topo_vect
 
-    def _grid2op_bus_from_klu_bus(self, klu_bus):
-        res = 0
-        if klu_bus != 0:
-            # object is connected
-            res = 1 if klu_bus < self.__nb_bus_before else 2
-        return res
-
-    def _klu_bus_from_grid2op_bus(self, grid2op_bus, grid2op_bus_init):
-        return grid2op_bus_init[grid2op_bus - 1]
-
-    def get_topo_vect(self):
-        return self.topo_vect
-
-    def generators_info(self):
+    def generators_info(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
         return self.cst_1 * self.prod_p, self.cst_1 * self.prod_q, self.cst_1 * self.prod_v
 
-    def loads_info(self):
+    def loads_info(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
         return self.cst_1 * self.load_p, self.cst_1 * self.load_q, self.cst_1 * self.load_v
 
-    def lines_or_info(self):
+    def lines_or_info(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
         return self.cst_1 * self.p_or, self.cst_1 * self.q_or, self.cst_1 * self.v_or, self.cst_1 * self.a_or
 
-    def lines_ex_info(self):
+    def lines_ex_info(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
         return self.cst_1 * self.p_ex, self.cst_1 * self.q_ex, self.cst_1 * self.v_ex, self.cst_1 * self.a_ex
 
-    def storages_info(self):
+    def storages_info(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
         if not self.__has_storage:
             raise RuntimeError("Storage units are not supported with your grid2op version. Please upgrade to "
                                "grid2op >1.5")
         return self.cst_1 * self.storage_p, self.cst_1 * self.storage_q, self.cst_1 * self.storage_v
 
-    def shunt_info(self):
+    def shunt_info(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
         return self.cst_1 * self.sh_p, self.cst_1 * self.sh_q, self.cst_1 * self.sh_v, self.sh_bus
 
+    def _compute_shunt_bus_with_compat(self, shunt_bus):
+        cls = type(self)
+        if hasattr(cls, "global_bus_to_local"):
+            self.sh_bus[:] = cls.global_bus_to_local(shunt_bus, cls.shunt_to_subid)
+        else:
+            res = (1 * shunt_bus).astype(dt_int)  # make a copy
+            if hasattr(cls, "n_busbar_per_sub"):
+                n_busbar_per_sub = cls.n_busbar_per_sub
+            else:
+                # backward compat when this was not defined:
+                n_busbar_per_sub = DEFAULT_N_BUSBAR_PER_SUB
+            for i in range(n_busbar_per_sub):
+                res[(i * cls.n_sub <= shunt_bus) & (shunt_bus < (i+1) * cls.n_sub)] = i + 1
+            res[shunt_bus == -1] = -1
+            self.sh_bus[:] = res
+        
     def _set_shunt_info(self):
-        self.sh_p[:], self.sh_q[:], self.sh_v[:]  = self._grid.get_shunts_res()
-        shunt_bus = np.array([self._grid.get_bus_shunt(i) for i in range(type(self).n_shunt)], dtype=dt_int)
-        res_bus = np.ones(shunt_bus.shape[0], dtype=dt_int)  # by default all shunts are on bus one
-        res_bus[shunt_bus >= self.__nb_bus_before] = 2  # except the one that are connected to bus 2
-        res_bus[shunt_bus == -1] = -1  # or the one that are disconnected
-        self.sh_bus[:] = res_bus
+        tick = time.perf_counter()
+        self.sh_p[:], self.sh_q[:], self.sh_v[:], self.sh_theta[:]  = self._shunt_res
         self.sh_v[self.sh_v == -1.] = 0.  # in grid2op disco element have voltage of 0. and -1.
-
+        self._timer_read_data_back += time.perf_counter() - tick
+        # self._timer_fetch_data_cpp += time.perf_counter() - tick
+        
     def _disconnect_line(self, id_):
         self.topo_vect[self.line_ex_pos_topo_vect[id_]] = -1
         self.topo_vect[self.line_or_pos_topo_vect[id_]] = -1
@@ -1330,19 +1448,26 @@ def _disconnect_line(self, id_):
         else:
             self._grid.deactivate_trafo(id_ - self.__nb_powerline)
 
-    def get_current_solver_type(self):
+    def get_current_solver_type(self) -> SolverType:
         return self.__current_solver_type
 
-    def reset(self, grid_path, grid_filename=None):
+    def reset(self,
+              path : Union[os.PathLike, str],
+              grid_filename : Optional[Union[os.PathLike, str]]=None) -> None:
         self._fill_nans()
         self._grid = self.__me_at_init.copy()
         self._grid.unset_changes()
         self._grid.change_solver(self.__current_solver_type)
         self._handle_turnedoff_pv()
-        self.topo_vect[:] = self.__init_topo_vect
         self.comp_time = 0.
         self.timer_gridmodel_xx_pf = 0.
         self._timer_postproc = 0.
         self._timer_preproc = 0.
         self._timer_solver = 0.
+        self._timer_read_data_back = 0.
+        self._timer_fetch_data_cpp = 0.
+        self._timer_apply_act = 0.
         self._grid.tell_solver_need_reset()
+        self._reset_res_pointers()
+        self.sh_bus[:] = 1  # TODO self._compute_shunt_bus_with_compat(self._grid.get_all_shunt_buses())
+        self.topo_vect[:] = self.__init_topo_vect  # TODO
diff --git a/lightsim2grid/rewards/__init__.py b/lightsim2grid/rewards/__init__.py
new file mode 100644
index 0000000..be7dc2a
--- /dev/null
+++ b/lightsim2grid/rewards/__init__.py
@@ -0,0 +1,11 @@
+# Copyright (c) 2020-2024, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+__all__ = ["N1ContingencyReward"]
+
+from lightsim2grid.rewards.n1ContingencyReward import N1ContingencyReward
diff --git a/lightsim2grid/rewards/n1ContingencyReward.py b/lightsim2grid/rewards/n1ContingencyReward.py
new file mode 100644
index 0000000..e2db2a0
--- /dev/null
+++ b/lightsim2grid/rewards/n1ContingencyReward.py
@@ -0,0 +1,185 @@
+# Copyright (c) 2020-2024, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
+
+import time
+import numpy as np
+
+import grid2op
+from grid2op.Reward import BaseReward
+from grid2op.Action._backendAction import _BackendAction
+
+from lightsim2grid import LightSimBackend, ContingencyAnalysis
+from lightsim2grid.compilation_options import klu_solver_available
+from lightsim2grid.solver import SolverType
+
+
+class N1ContingencyReward(BaseReward):
+    """
+    This class implements a reward that leverage the :class:`lightsim2grid.ContingencyAnalysis`
+    to compute the number of unsafe contingency at any given time.
+
+    Examples
+    --------
+
+    This can be used as:
+
+    .. code-block:: python
+
+        import grid2op
+        from lightsim2grid.rewards import N1ContingencyReward
+        l_ids = [0, 1, 7]
+        env = grid2op.make("l2rpn_case14_sandbox",
+                           reward_class=N1ContingencyReward(l_ids=l_ids)
+                          )
+        obs = env.reset()
+        obs, reward, *_ = env.step(env.action_space())
+        print(f"reward: {reward:.3f}")
+        
+    """
+
+    def __init__(self,
+                 l_ids=None,
+                 threshold_margin=1.,
+                 dc=False,
+                 normalize=False,
+                 logger=None,
+                 tol=1e-8,
+                 nb_iter=10):
+        BaseReward.__init__(self, logger=logger)
+        self._backend : LightSimBackend = None
+        self._backend_action = None
+        self._l_ids = None
+        self._dc : bool = dc
+        self._normalize : bool = normalize
+        if l_ids is not None:
+            self._l_ids = [int(el) for el in l_ids]
+        self._threshold_margin :float = float(threshold_margin)
+        if klu_solver_available:
+            if self._dc:
+                self._solver_type = SolverType.KLUDC
+            else:
+                self._solver_type = SolverType.KLU
+        else:
+            if self._dc:
+                self._solver_type = SolverType.DC
+            else:
+                self._solver_type = SolverType.SparseLU
+        self._backend_ls = False
+        self._tol = tol
+        self._nb_iter = nb_iter
+        self._timer_call = 0.
+        self._timer_pre_proc = 0.
+        self._timer_compute = 0.
+        self._timer_post_proc = 0.
+            
+    def initialize(self, env: "grid2op.Environment.Environment"):
+        from grid2op.Environment import BaseEnv
+        from grid2op.Backend import PandaPowerBackend  # lazy import because grid2op -> pandapower-> lightsim2grid -> grid2op
+        if not isinstance(env, BaseEnv):
+            raise RuntimeError("You can only initialize this reward with a "
+                               "proper grid2op environment (`BaseEnv`)")
+             
+        if not isinstance(env.backend, (PandaPowerBackend, LightSimBackend)):
+            raise RuntimeError("Impossible to use the `N1ContingencyReward` with "
+                               "a environment with a backend that is not "
+                               "``PandaPowerBackend` nor `LightSimBackend`."
+                               )
+        if isinstance(env.backend, LightSimBackend):
+            self._backend : LightSimBackend = env.backend.copy()
+            self._backend_ls :bool  = True
+        elif isinstance(env.backend, PandaPowerBackend):
+            self._backend = LightSimBackend.init_grid(type(env.backend))()
+            self._backend.init_from_loaded_pandapower(env.backend)
+            self._backend.is_loaded = True
+        else:
+            raise NotImplementedError()
+        
+        self._backend.set_solver_type(self._solver_type)
+        conv, exc_ = self._backend.runpf()
+        if not conv:
+            raise RuntimeError(f"The reward N1ContingencyReward diverge with error {exc_}")
+        bk_act_cls = _BackendAction.init_grid(type(env.backend))
+        self._backend_action = bk_act_cls()
+        if self._l_ids is None:
+            self._l_ids = list(range(type(env).n_line))
+        
+        if len(self._l_ids) == 0:
+            raise RuntimeError("Impossible to use the N1ContingencyReward "
+                               "without any contingencies !")
+        self.reward_min = 0.
+        self.reward_max = len(self._l_ids) if not self._normalize else 1.
+        # self._contingecy_analyzer = ContingencyAnalysis(self._backend)
+        # self._contingecy_analyzer.add_multiple_contingencies(self._l_ids)
+
+    def __call__(self, action, env, has_error, is_done, is_illegal, is_ambiguous):
+        if is_done:
+            return self.reward_min
+        beg = time.perf_counter()
+        # retrieve the state of the grid
+        self._backend_action.reset()
+        act = env.backend.get_action_to_set()
+        th_lim_a = env.get_thermal_limit()
+        th_lim_a[th_lim_a <= 1] = 1  # assign 1 for the thermal limit
+        
+        # apply it to the backend
+        self._backend_action += act
+        self._backend.apply_action(self._backend_action)
+        conv, exc_ = self._backend.runpf()
+        if not conv:
+            self.logger.warn("Cannot set the backend of the `N1ContingencyReward` => divergence")
+            return self.reward_min
+        
+        # synch the contingency analyzer
+        # self._contingecy_analyzer.update_grid(self._backend_action)
+        contingecy_analyzer = ContingencyAnalysis(self._backend)
+        contingecy_analyzer.computer.change_solver(self._solver_type)
+        contingecy_analyzer.add_multiple_contingencies(*self._l_ids)
+        now_ = time.perf_counter()
+        self._timer_pre_proc += now_ - beg
+        tmp = contingecy_analyzer.get_flows()
+        self.logger.info(f"{contingecy_analyzer.computer.nb_solved()} converging contingencies")
+        now_2 = time.perf_counter()
+        self._timer_compute += now_2 - now_
+        if self._dc:
+            # In DC is study p, but take into account q in the limits
+            tmp_res = np.abs(tmp[0])  # this is Por
+            # now transform the limits in A in MW
+            por, qor, vor, aor = env.backend.lines_or_info()
+            p_sq = (1e-3*th_lim_a)**2 * 3. * vor**2 - qor**2
+            p_sq[p_sq <= 0.] = 0.
+            limits = np.sqrt(p_sq)
+        else:
+            tmp_res = tmp[1]
+            limits = th_lim_a
+        # print("Reward:")
+        # print(tmp_res)        
+        # print(self._threshold_margin * limits)
+        res = ((tmp_res > self._threshold_margin * limits) | (~np.isfinite(tmp_res))).any(axis=1)  # whether one powerline is above its limit, per cont
+        res |=  (np.abs(tmp_res) <= self._tol).all(axis=1)  # other type of divergence: all 0.
+        # print(res.nonzero())
+        res = res.sum()  # count total of n-1 unsafe 
+        res = len(self._l_ids) - res  # reward = things to maximise
+        if self._normalize:
+            res /= len(self._l_ids)
+        now_3 = time.perf_counter()
+        self._timer_post_proc += now_3 - now_2
+        self._timer_call += time.perf_counter() - beg
+        return res
+
+    def reset(self, env: "grid2op.Environment.BaseEnv") -> None:
+        self._timer_call = 0.
+        self._timer_pre_proc = 0.
+        self._timer_compute = 0.
+        self._timer_post_proc = 0.
+        return super().reset(env)
+    
+    def close(self):
+        if self._backend is not None:
+            self._backend.close()
+        del self._backend
+        self._backend = None
diff --git a/lightsim2grid/securityAnalysis.py b/lightsim2grid/securityAnalysis.py
index f9db2b1..2940aea 100644
--- a/lightsim2grid/securityAnalysis.py
+++ b/lightsim2grid/securityAnalysis.py
@@ -6,8 +6,26 @@
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
 
-from lightsim2grid.contingencyAnalysis import ContingencyAnalysisCPP, ContingencyAnalysis
-
+import warnings
 # Deprecated now, will be removed
+warnings.warn("You are using old names. Please upgrade to SecurityAnalysisCPP > ContingencyAnalysisCPP"
+              " and SecurityAnalysis > ContingencyAnalysis instead.",
+              category=DeprecationWarning)
+__all__ = ["SecurityAnalysisCPP"]
+
+
+from lightsim2grid.contingencyAnalysis import ContingencyAnalysisCPP
+
+try:
+    from lightsim2grid.contingencyAnalysis import ContingencyAnalysis
+    GRID2OP_INSTALLED = True
+except ImportError as exc_:
+    GRID2OP_INSTALLED = False
+
+
 SecurityAnalysisCPP = ContingencyAnalysisCPP
-SecurityAnalysis = ContingencyAnalysis
+
+
+if GRID2OP_INSTALLED:
+    SecurityAnalysis = ContingencyAnalysis
+    __all__.append("SecurityAnalysis")
diff --git a/lightsim2grid/timeSerie.py b/lightsim2grid/timeSerie.py
index 2942441..72b6445 100644
--- a/lightsim2grid/timeSerie.py
+++ b/lightsim2grid/timeSerie.py
@@ -6,16 +6,23 @@
 # SPDX-License-Identifier: MPL-2.0
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
 
-__all__ = ["TimeSerieCPP", "TimeSerie",
+__all__ = ["TimeSeriesCPP", 
            # deprecated
            "Computers"]
 
 import numpy as np
 import warnings
 
-from grid2op.Chronics import Multifolder, GridStateFromFile
 
-from lightsim2grid.lightSimBackend import LightSimBackend
+try:
+    from grid2op.Chronics import Multifolder, GridStateFromFile
+    from lightsim2grid.lightSimBackend import LightSimBackend
+    __all__.append("TimeSerie")
+    GRID2OP_INSTALLED = True
+except ImportError as exc_:
+    # grid2Op is not installed
+    GRID2OP_INSTALLED = False
+
 from lightsim2grid.solver import SolverType
 from lightsim2grid_cpp import TimeSeriesCPP
 
@@ -23,7 +30,7 @@
 Computers = TimeSeriesCPP
 
 
-class TimeSerie:
+class ___TimeSerie:
     """
     This helper class, that only works with grid2op when using a LightSimBackend allows to compute
     the flows (at the origin side of the powerline / transformers). It is roughly equivalent to the
@@ -77,6 +84,13 @@ class TimeSerie:
 
     """
     def __init__(self, grid2op_env):
+        if not GRID2OP_INSTALLED:
+            raise RuntimeError("Impossible to use the python wrapper `TimeSerie` "
+                               "when grid2op is not installed. Please fall back to the "
+                               "c++ version (available in python) with:\n"
+                               "\tfrom lightsim2grid.timeSerie import TimeSerieCPP\n"
+                               "and refer to the appropriate documentation.")
+            
         from grid2op.Environment import Environment  # otherwise i got issues...
         if not isinstance(grid2op_env.backend, LightSimBackend):
             raise RuntimeError("This class only works with LightSimBackend")
@@ -263,3 +277,8 @@ def _extract_inj(self):
         self.load_p = 1.0 * data_loader.load_p
         self.load_q = 1.0 * data_loader.load_q
         return self.prod_p, self.load_p, self.load_q
+
+
+if GRID2OP_INSTALLED:
+    TimeSerie = ___TimeSerie
+    
\ No newline at end of file
diff --git a/setup.py b/setup.py
index cb4269d..a8ca33b 100644
--- a/setup.py
+++ b/setup.py
@@ -14,7 +14,7 @@
 from pybind11.setup_helpers import Pybind11Extension, build_ext
 
 
-__version__ = "0.8.0"
+__version__ = "0.8.1"
 KLU_SOLVER_AVAILABLE = False
 
 # Try to link against SuiteSparse (if available)
@@ -136,6 +136,7 @@
 extra_compile_args += [f"-DVERSION_MAJOR={VERSION_MAJOR}",
                        f"-DVERSION_MEDIUM={VERSION_MEDIUM}",
                        f"-DVERSION_MINOR={VERSION_MINOR}"]
+extra_compile_args += [f"-DVERSION=\"{__version__}\""]
 src_files = ['src/main.cpp',
              "src/powerflow_algorithm/GaussSeidelAlgo.cpp",
              "src/powerflow_algorithm/GaussSeidelSynchAlgo.cpp",
@@ -210,6 +211,7 @@
         include_dirs.append(os.path.join(path_nicslu, "include"))
         src_files.append("src/linear_solvers/NICSLUSolver.cpp")
         extra_compile_args.append("-DNICSLU_SOLVER_AVAILABLE")
+        extra_compile_args += [f"-DNICSLU_PATH=\"{os.path.join(path_nicslu, libnicslu_path)}\""]
         print("INFO: Using NICSLU package")
 
 # Try to locate the CKTSO sparse linear solver
@@ -257,6 +259,7 @@
         include_dirs.append(os.path.join(path_cktso, "include"))
         src_files.append("src/linear_solvers/CKTSOSolver.cpp")
         extra_compile_args.append("-DCKTSO_SOLVER_AVAILABLE")
+        extra_compile_args += [f"-DCKTSO_PATH=\"{os.path.join(path_cktso, libcktso_path)}\""]
         print("INFO: Using CKTSO package")
         
         
@@ -269,6 +272,7 @@
 if "__COMPILE_MARCHNATIVE" in os.environ:
     if os.environ["__COMPILE_MARCHNATIVE"] == "1":
         extra_compile_args.append("-march=native")
+        extra_compile_args.append("-D__COMPILE_MARCHNATIVE")
         print("INFO: Using \"-march=native\" compiler flag")
 
 # $Env:_READ_THE_DOCS = "1" in powershell
@@ -299,6 +303,7 @@
             print("INFO: Using \"-O3\" compiler flag")
         elif IS_WINDOWS:
             extra_compile_args.append("/O2")
+        extra_compile_args.append("-D__O3_OPTIM")
 
 ext_modules = [
     Pybind11Extension(
@@ -316,7 +321,8 @@
         "pandapower" if sys.version_info < (3, 10) else "pandapower>=2.8",
         "pytest",  # for pandapower see https://github.com/e2nIEE/pandapower/issues/1988
     ]
-
+if sys.version_info < (3, 11):
+    req_pkgs.append("typing_extensions")
 if sys.version_info.major == 3 and sys.version_info.minor <= 7:
     # typing "Literal" not available on python 3.7
     req_pkgs.append("typing_extensions")
diff --git a/src/BaseConstants.cpp b/src/BaseConstants.cpp
index ba02112..fe1ba35 100644
--- a/src/BaseConstants.cpp
+++ b/src/BaseConstants.cpp
@@ -14,3 +14,4 @@ const real_type BaseConstants::my_one_ = 1.0;
 const real_type BaseConstants::my_two_ = 2.0;
 const real_type BaseConstants::my_half_ = 0.5;
 const real_type BaseConstants::my_zero_ = 0.;
+const real_type BaseConstants::my_180_pi_ = 180. / M_PI;
diff --git a/src/BaseConstants.h b/src/BaseConstants.h
index 7d6f268..329aabf 100644
--- a/src/BaseConstants.h
+++ b/src/BaseConstants.h
@@ -24,6 +24,7 @@ class BaseConstants
         static const real_type my_one_;
         static const real_type my_two_;
         static const real_type my_zero_;
+        static const real_type my_180_pi_;
 };
 
 enum class FDPFMethod {XB, BX};  // Different type of FDPF powerflow
diff --git a/src/ChooseSolver.h b/src/ChooseSolver.h
index 9c17e43..383f12b 100644
--- a/src/ChooseSolver.h
+++ b/src/ChooseSolver.h
@@ -286,6 +286,7 @@ class ChooseSolver
             Eigen::SparseMatrix<real_type> res = get_J();
             return res;
         }
+
         double get_computation_time() const
         {
             auto p_solver = get_prt_solver("get_computation_time", true);
@@ -293,6 +294,18 @@ class ChooseSolver
             return std::get<3>(res);
         }
 
+        std::tuple<double, double, double, double> get_timers() const
+        {
+            auto p_solver = get_prt_solver("get_timers", true);
+            return p_solver -> get_timers();
+        }
+
+        TimerJacType get_timers_jacobian() const
+        {
+            const BaseAlgo * p_solver = get_prt_solver("get_timers_jacobian", true);
+            return p_solver -> get_timers_jacobian();
+        }
+
         ErrorType get_error() const{
             auto p_solver = get_prt_solver("get_error", true);
             return p_solver -> get_error();
diff --git a/src/GridModel.cpp b/src/GridModel.cpp
index 67b4853..1ca5fe0 100644
--- a/src/GridModel.cpp
+++ b/src/GridModel.cpp
@@ -371,7 +371,7 @@ CplxVect GridModel::ac_pf(const CplxVect & Vinit,
     bool conv = false;
     CplxVect res = CplxVect();
 
-    reset_results();  // reset the results
+    // reset_results();  // clear the results  No need to do it, results are neceassirly set or reset in post process
 
     // pre process the data to define a proper jacobian matrix, the proper voltage vector etc.
     bool is_ac = true;
@@ -387,20 +387,17 @@ CplxVect GridModel::ac_pf(const CplxVect & Vinit,
                                     solver_control_);
 
     // start the solver
-    // std::cout << "\tbefore get_slack_weights (ac)" << std::endl;
     if(solver_control_.need_reset_solver() || 
        solver_control_.has_dimension_changed() ||
        solver_control_.has_slack_participate_changed() || 
        solver_control_.has_pv_changed() || 
        solver_control_.has_slack_weight_changed()){
-        // std::cout << "get_slack_weights" << std::endl;
         slack_weights_ = generators_.get_slack_weights(Ybus_ac_.rows(), id_me_to_ac_solver_); 
     }
     // std::cout << "\tbefore compute_pf" << std::endl;
     conv = _solver.compute_pf(Ybus_ac_, V, acSbus_, slack_bus_id_ac_solver_, slack_weights_, bus_pv_, bus_pq_, max_iter, tol / sn_mva_);
 
     // store results (in ac mode) 
-    // std::cout << "\tbefore process_results" << std::endl;
     process_results(conv, res, Vinit, true, id_me_to_ac_solver_);
 
     timer_last_ac_pf_ = timer.duration();
@@ -535,47 +532,33 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit,
     if (solver_control.need_reset_solver() || 
         solver_control.has_dimension_changed() ||
         solver_control.has_slack_participate_changed()){
-            // std::cout << "\t\tslack_bus_id_solver\n";
             slack_bus_id_me = generators_.get_slack_bus_id();
             // this is the slack bus ids with the gridmodel ordering, not the solver ordering.
             // conversion to solver ordering is done in init_slack_bus
-
-            // std::cout << "slack_bus_id_solver 1: ";
-            // for(auto el: slack_bus_id_solver) std::cout << el << ", ";
-            // std::cout << std::endl;
         }
     if (solver_control.need_reset_solver() ||
         solver_control.ybus_change_sparsity_pattern() || 
         solver_control.has_dimension_changed()){
-            // std::cout << "\t\tinit_Ybus;" << std::endl;
             init_Ybus(Ybus, id_me_to_solver, id_solver_to_me);
-            // std::cout << "init_Ybus;" << std::endl;
         }
     if (solver_control.need_reset_solver() ||
         solver_control.ybus_change_sparsity_pattern() || 
         solver_control.has_dimension_changed() || 
         solver_control.need_recompute_ybus()){
-            // std::cout << "\t\tfillYbus;" << std::endl;
             fillYbus(Ybus, is_ac, id_me_to_solver);
-            // std::cout << "fillYbus;" << std::endl;
         }
     if (solver_control.need_reset_solver() || 
         solver_control.has_dimension_changed()) {
             // init Sbus
-            // std::cout << "\t\tinit_Sbus;" << std::endl;
             Sbus = CplxVect::Constant(id_solver_to_me.size(), 0.);
-            // std::cout << "init_Sbus;" << std::endl;
         }
     if (solver_control.need_reset_solver() || 
         solver_control.has_dimension_changed() ||
         solver_control.has_slack_participate_changed() || 
         solver_control.has_pv_changed() || 
         solver_control.has_pq_changed()) {
-            // std::cout << "\t\tinit_slack_bus;" << std::endl;
             init_slack_bus(Sbus, id_me_to_solver, id_solver_to_me, slack_bus_id_me, slack_bus_id_solver);
-            // std::cout << "\t\tfillpv_pq;" << std::endl;
             fillpv_pq(id_me_to_solver, id_solver_to_me, slack_bus_id_solver, solver_control);
-            // std::cout << "fillpv_pq;" << std::endl;
         }
     
     if (is_ac && (solver_control.need_reset_solver() || 
@@ -583,14 +566,12 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit,
                   solver_control.need_recompute_sbus() ||  // TODO do we need it ?
                   solver_control.has_pq_changed())  // TODO do we need it ?
         ){
-        // std::cout << "\t\tq vector" << std::endl;
         int nb_bus_total = static_cast<int>(bus_vn_kv_.size());
         total_q_min_per_bus_ = RealVect::Constant(nb_bus_total, 0.);
         total_q_max_per_bus_ = RealVect::Constant(nb_bus_total, 0.);
         total_gen_per_bus_ = Eigen::VectorXi::Constant(nb_bus_total, 0);
         generators_.init_q_vector(nb_bus_total, total_gen_per_bus_, total_q_min_per_bus_, total_q_max_per_bus_);
         dc_lines_.init_q_vector(nb_bus_total, total_gen_per_bus_, total_q_min_per_bus_, total_q_max_per_bus_);
-        // std::cout << "total_gen_per_bus_;" << std::endl;
     }
 
     if (solver_control.need_reset_solver() || 
@@ -598,11 +579,8 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit,
         solver_control.has_slack_participate_changed() || 
         solver_control.has_pq_changed() ||
         solver_control.need_recompute_sbus()) {
-            // std::cout << "\t\tfillSbus_me;" << std::endl;
             fillSbus_me(Sbus, is_ac, id_me_to_solver);
-            // std::cout << "fillSbus_me;" << std::endl;
         }
-    // std::cout << "\t\tstatic_cast<int>(id_solver_to_me.size());" << std::endl;
     const int nb_bus_solver = static_cast<int>(id_solver_to_me.size());
     CplxVect V = CplxVect::Constant(nb_bus_solver, init_vm_pu_);
     for(int bus_solver_id = 0; bus_solver_id < nb_bus_solver; ++bus_solver_id){
@@ -612,11 +590,6 @@ CplxVect GridModel::pre_process_solver(const CplxVect & Vinit,
     }
     generators_.set_vm(V, id_me_to_solver);
     dc_lines_.set_vm(V, id_me_to_solver);
-    // std::cout << "pre_process_solver: V result: "<<std::endl;
-    // for(auto el: V) std::cout << el << ", ";
-    // std::cout << std::endl;
-    // std::cout << "nb_bus_solver " << nb_bus_solver << std::endl;
-    // std::cout << "\t\tend pre_process_solver" << std::endl;
     return V;
 }
 
@@ -901,7 +874,7 @@ CplxVect GridModel::dc_pf(const CplxVect & Vinit,
     bool conv = false;
     CplxVect res = CplxVect();
 
-    reset_results();  // reset the results
+    // reset_results();  // clear the results  No need to do it, results are neceassirly set or reset in post process
 
     // pre process the data to define a proper jacobian matrix, the proper voltage vector etc.
     bool is_ac = false;
diff --git a/src/GridModel.h b/src/GridModel.h
index 7a4b69c..2ffd4f3 100644
--- a/src/GridModel.h
+++ b/src/GridModel.h
@@ -494,6 +494,21 @@ class GridModel : public GenericContainer
         Eigen::Ref<const RealVect> get_dclineor_theta() const {return dc_lines_.get_theta_or();}
         Eigen::Ref<const RealVect> get_dclineex_theta() const {return dc_lines_.get_theta_ex();}
 
+        Eigen::Ref<const IntVect> get_all_shunt_buses() const {return shunts_.get_buses();}
+
+        // complete results (with theta)
+        tuple4d get_loads_res_full() const {return loads_.get_res_full();}
+        tuple4d get_shunts_res_full() const {return shunts_.get_res_full();}
+        tuple4d get_gen_res_full() const {return generators_.get_res_full();}
+        tuple5d get_lineor_res_full() const {return powerlines_.get_res_or_full();}
+        tuple5d get_lineex_res_full() const {return powerlines_.get_res_ex_full();}
+        tuple5d get_trafohv_res_full() const {return trafos_.get_res_hv_full();}
+        tuple5d get_trafolv_res_full() const {return trafos_.get_res_lv_full();}
+        tuple4d get_storages_res_full() const {return storages_.get_res_full();}
+        tuple4d get_sgens_res_full() const {return sgens_.get_res_full();}
+        tuple4d get_dclineor_res_full() const {return dc_lines_.get_res_or_full();}
+        tuple4d get_dclineex_res_full() const {return dc_lines_.get_res_ex_full();}
+
         // get some internal information, be cerafull the ID of the buses might not be the same
         // TODO convert it back to this ID, that will make copies, but who really cares ?
         Eigen::SparseMatrix<cplx_type> get_Ybus(){
diff --git a/src/Utils.h b/src/Utils.h
index eeaec05..4cb448e 100644
--- a/src/Utils.h
+++ b/src/Utils.h
@@ -22,6 +22,10 @@ typedef double real_type;  // type for real numbers: can be changed if installed
 typedef std::complex<real_type> cplx_type;  // type for complex number
 
 typedef Eigen::Matrix<real_type, Eigen::Dynamic, 1> EigenPythonNumType;  // Eigen::VectorXd
+typedef Eigen::Matrix<int, Eigen::Dynamic, 1> IntVect;
+typedef Eigen::Matrix<real_type, Eigen::Dynamic, 1> RealVect;
+typedef Eigen::Matrix<cplx_type, Eigen::Dynamic, 1> CplxVect;
+
 typedef std::tuple<Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType> > tuple3d;
@@ -29,9 +33,11 @@ typedef std::tuple<Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType>,
                    Eigen::Ref<const EigenPythonNumType> > tuple4d;
-typedef Eigen::Matrix<int, Eigen::Dynamic, 1> IntVect;
-typedef Eigen::Matrix<real_type, Eigen::Dynamic, 1> RealVect;
-typedef Eigen::Matrix<cplx_type, Eigen::Dynamic, 1> CplxVect;
+typedef std::tuple<Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType> > tuple5d;
 
 typedef Eigen::Matrix<real_type, Eigen::Dynamic, Eigen::Dynamic> RealMat;
 typedef Eigen::Matrix<cplx_type, Eigen::Dynamic, Eigen::Dynamic> CplxMat;
diff --git a/src/batch_algorithm/BaseBatchSolverSynch.h b/src/batch_algorithm/BaseBatchSolverSynch.h
index 8256c62..42c1ae5 100644
--- a/src/batch_algorithm/BaseBatchSolverSynch.h
+++ b/src/batch_algorithm/BaseBatchSolverSynch.h
@@ -61,6 +61,16 @@ class BaseBatchSolverSynch
         std::vector<SolverType> available_solvers() const {return _solver.available_solvers(); }
         SolverType get_solver_type() const {return _solver.get_type(); }
 
+        // TODO
+        // void change_gridmodel(const GridModel & new_grid_model){
+        //     if(_grid_model.total_bus() != new_grid_model.total_bus()){
+        //         throw std::runtime_error("BaseBatchSolverSynch: impossible to change the gridmodel. It appears it's not the same grid (not the same number of buses).");
+        //     }
+        //     CplxVect V = CplxVect::Constant(new_grid_model.total_bus(), 1.04);
+        //     _grid_model = new_grid_model.copy(); // not implemented !!!
+
+        // }
+
         // utlities informations
         double amps_computation_time() const {return _timer_compute_A;}
         double solver_time() const {return _timer_solver;}
@@ -71,9 +81,9 @@ class BaseBatchSolverSynch
             _active_power_flows = RealMat();
             _voltages = CplxMat();
             _nb_solved = 0;
-            _timer_compute_A = 0;
-            _timer_compute_P = 0;
-            _timer_solver = 0;
+            _timer_compute_A = 0.;
+            _timer_compute_P = 0.;
+            _timer_solver = 0.;
         }
 
         // results
diff --git a/src/batch_algorithm/ContingencyAnalysis.cpp b/src/batch_algorithm/ContingencyAnalysis.cpp
index 1118289..abedbd4 100644
--- a/src/batch_algorithm/ContingencyAnalysis.cpp
+++ b/src/batch_algorithm/ContingencyAnalysis.cpp
@@ -52,6 +52,7 @@ void ContingencyAnalysis::init_li_coeffs(bool ac_solver_used){
     const auto & id_me_to_solver = ac_solver_used ? _grid_model.id_me_to_ac_solver(): _grid_model.id_me_to_dc_solver();
     Eigen::Index bus_1_id, bus_2_id;
     cplx_type y_ff, y_ft, y_tf, y_tt;
+    bool status;
     for(const auto & this_cont_id: _li_defaults){
         std::vector<Coeff> this_cont_coeffs;
         this_cont_coeffs.reserve(this_cont_id.size() * 4);  // usually there are 4 coeffs per powerlines / trafos
@@ -61,6 +62,7 @@ void ContingencyAnalysis::init_li_coeffs(bool ac_solver_used){
                 // this is a powerline
                 bus_1_id = id_me_to_solver[powerlines.get_bus_from()[line_id]];
                 bus_2_id = id_me_to_solver[powerlines.get_bus_to()[line_id]];
+                status = powerlines.get_status()[line_id];
                 if(ac_solver_used){
                     y_ff = powerlines.yac_ff()[line_id];
                     y_ft = powerlines.yac_ft()[line_id];
@@ -75,6 +77,7 @@ void ContingencyAnalysis::init_li_coeffs(bool ac_solver_used){
             }else{
                 // this is a trafo
                 const auto trafo_id = line_id - n_line_;
+                status = trafos.get_status()[trafo_id];
                 bus_1_id = id_me_to_solver[trafos.get_bus_from()[trafo_id]];
                 bus_2_id = id_me_to_solver[trafos.get_bus_to()[trafo_id]];
                 if(ac_solver_used){
@@ -90,7 +93,7 @@ void ContingencyAnalysis::init_li_coeffs(bool ac_solver_used){
                 }
             }
 
-            if(bus_1_id != GenericContainer::_deactivated_bus_id && bus_2_id != GenericContainer::_deactivated_bus_id)
+            if(status && bus_1_id != GenericContainer::_deactivated_bus_id && bus_2_id != GenericContainer::_deactivated_bus_id)
             {
                 // element is connected
                 this_cont_coeffs.push_back({bus_1_id, bus_1_id, y_ff});
@@ -103,15 +106,15 @@ void ContingencyAnalysis::init_li_coeffs(bool ac_solver_used){
     }
 }
 
-
 bool ContingencyAnalysis::remove_from_Ybus(Eigen::SparseMatrix<cplx_type> & Ybus,
-                                        const std::vector<Coeff> & coeffs) const
+                                           const std::vector<Coeff> & coeffs) const
 {
     for(const auto & coeff_to_remove: coeffs){
         Ybus.coeffRef(coeff_to_remove.row_id, coeff_to_remove.col_id) -= coeff_to_remove.value;
     }
     return check_invertible(Ybus);
 }
+
 void ContingencyAnalysis::readd_to_Ybus(Eigen::SparseMatrix<cplx_type> & Ybus,
                                      const std::vector<Coeff> & coeffs) const
 {
@@ -180,18 +183,12 @@ void ContingencyAnalysis::compute(const CplxVect & Vinit, int max_iter, real_typ
     Eigen::Index cont_id = 0;
     bool conv;
     CplxVect V;
-    // int contingency = 0;
     for(const auto & coeffs_modif: _li_coeffs){
         auto timer_modif_Ybus = CustTimer();
         bool invertible = remove_from_Ybus(Ybus, coeffs_modif);
         _timer_modif_Ybus += timer_modif_Ybus.duration();
         conv = false;
 
-        // I have absolutely no idea why, but if i add this "if"
-        // which reduces the computation time, it somehow increase it by A LOT !
-        // 5.2ms without it vs 81.9ms with it (for the iee 118)
-        // So better make the computation, even if it's not used...
-
         if(invertible)
         {
             V = Vinit_solver; // Vinit is reused for each contingencies
diff --git a/src/batch_algorithm/ContingencyAnalysis.h b/src/batch_algorithm/ContingencyAnalysis.h
index 052722e..86dcfb9 100644
--- a/src/batch_algorithm/ContingencyAnalysis.h
+++ b/src/batch_algorithm/ContingencyAnalysis.h
@@ -74,6 +74,12 @@ class ContingencyAnalysis: public BaseBatchSolverSynch
             _timer_modif_Ybus = 0.;
             _timer_pre_proc = 0.;
         }
+        void clear_results_only(){
+            BaseBatchSolverSynch::clear();
+            _timer_total = 0.;
+            _timer_modif_Ybus = 0.;
+            _timer_pre_proc = 0.;
+        }
         
         bool remove_n1(int line_id){
             check_ok_el(line_id);
diff --git a/src/element_container/DCLineContainer.cpp b/src/element_container/DCLineContainer.cpp
index 01a683b..7992598 100644
--- a/src/element_container/DCLineContainer.cpp
+++ b/src/element_container/DCLineContainer.cpp
@@ -36,6 +36,7 @@ void DCLineContainer::set_state(DCLineContainer::StateRes & my_state){
     status_ = status;
     loss_percent_ = RealVect::Map(&loss_percent[0], loss_percent.size());
     loss_mw_ = RealVect::Map(&loss_mw[0], loss_percent.size());
+    reset_results();
 }
 
 void DCLineContainer::init(const Eigen::VectorXi & branch_from_id,
diff --git a/src/element_container/DCLineContainer.h b/src/element_container/DCLineContainer.h
index a64fb1f..400c1a5 100644
--- a/src/element_container/DCLineContainer.h
+++ b/src/element_container/DCLineContainer.h
@@ -294,6 +294,9 @@ class DCLineContainer : public GenericContainer
     
     tuple3d get_or_res() const {return from_gen_.get_res();}
     tuple3d get_ex_res() const {return to_gen_.get_res();}
+    tuple4d get_res_or_full() const {return from_gen_.get_res_full();}
+    tuple4d get_res_ex_full() const {return to_gen_.get_res_full();}
+    
     Eigen::Ref<const RealVect> get_theta_or() const {return from_gen_.get_theta();}
     Eigen::Ref<const RealVect> get_theta_ex() const {return to_gen_.get_theta();}
 
diff --git a/src/element_container/GeneratorContainer.cpp b/src/element_container/GeneratorContainer.cpp
index 7598301..f5049ec 100644
--- a/src/element_container/GeneratorContainer.cpp
+++ b/src/element_container/GeneratorContainer.cpp
@@ -54,6 +54,7 @@ void GeneratorContainer::init(const RealVect & generators_p,
     turnedoff_gen_pv_ = true;
     voltage_regulator_on_ = std::vector<bool>(generators_p.size(), true);
     q_mvar_ = RealVect::Zero(generators_p.size());
+    reset_results();
 }
 
 void GeneratorContainer::init_full(const RealVect & generators_p,
@@ -94,7 +95,6 @@ GeneratorContainer::StateRes GeneratorContainer::get_state() const
 
 void GeneratorContainer::set_state(GeneratorContainer::StateRes & my_state)
 {
-    reset_results();
     names_ = std::get<0>(my_state);
     turnedoff_gen_pv_ = std::get<1>(my_state);
 
@@ -122,6 +122,7 @@ void GeneratorContainer::set_state(GeneratorContainer::StateRes & my_state)
     status_ = status;
     gen_slackbus_ = slack_bus;
     gen_slack_weight_ = slack_weight;
+    reset_results();
 }
 
 RealVect GeneratorContainer::get_slack_weights(Eigen::Index nb_bus_solver, const std::vector<int> & id_grid_to_solver){
@@ -218,15 +219,20 @@ void GeneratorContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
     const int nb_gen = nb();
     v_kv_from_vpu(Va, Vm, status_, nb_gen, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_);
     v_deg_from_va(Va, Vm, status_, nb_gen, bus_id_, id_grid_to_solver, bus_vn_kv, res_theta_);
-    res_p_ = p_mw_;
+    for(int gen_id = 0; gen_id < nb_gen; ++gen_id){
+        if(!status_[gen_id]){
+            res_p_(gen_id) = 0.;
+            continue;
+        }
+        res_p_(gen_id) = p_mw_(gen_id);
+    }
 }
 
 void GeneratorContainer::reset_results(){
-    res_p_ = RealVect();  // in MW
-    res_q_ = RealVect();  // in MVar
-    res_v_ = RealVect();  // in kV
-    res_theta_ = RealVect();  // in deg
-    // bus_slack_weight_ = RealVect();
+    res_p_ = RealVect(nb());  // in MW
+    res_q_ = RealVect(nb());  // in MVar
+    res_v_ = RealVect(nb());  // in kV
+    res_theta_ = RealVect(nb());  // in deg
 }
 
 void GeneratorContainer::get_vm_for_dc(RealVect & Vm){
@@ -380,8 +386,6 @@ void GeneratorContainer::set_p_slack(const RealVect& node_mismatch,
         // TODO DEBUG MODE: check bus_slack_weight_[bus_id_solver] > 0
         const auto total_contrib_slack = bus_slack_weight_(bus_id_solver);
         const auto my_contrib_slack = gen_slack_weight_[gen_id];
-        // now take "my part"
-        // std::cout << "gen_id " << gen_id << " my_contrib_slack " << my_contrib_slack << ", " << total_contrib_slack << " node_mismatch " << node_mismatch(bus_id_solver) << std::endl;
         res_p_(gen_id) += node_mismatch(bus_id_solver) * my_contrib_slack / total_contrib_slack;
     }
 }
@@ -414,16 +418,31 @@ void GeneratorContainer::set_q(const RealVect & reactive_mismatch,
                                const RealVect & total_q_max_per_bus)
 {
     const int nb_gen = nb();
-    res_q_ = RealVect::Constant(nb_gen, 0.);
-    if(!ac) return;  // do not consider Q values in dc mode
+    if(!ac){
+        // do not consider Q values in dc mode
+        for(int gen_id = 0; gen_id < nb_gen; ++gen_id) res_q_(gen_id) = 0.;
+        return;
+    }
+    
     real_type eps_q = 1e-8;
     for(int gen_id = 0; gen_id < nb_gen; ++gen_id)
     {
+        if(!status_[gen_id]){
+            // set at 0 for disconnected generators
+            res_q_(gen_id) = 0.;
+            continue;  
+        }
         real_type real_q = 0.;
-        if(!status_[gen_id]) continue;  // set at 0 for disconnected generators
-
-        if (!voltage_regulator_on_[gen_id]) continue;  // gen is purposedly not pv
-        if ((!turnedoff_gen_pv_) && p_mw_(gen_id) == 0.) continue;  // in this case turned off generators are not pv
+        if (!voltage_regulator_on_[gen_id]){
+            // gen is purposedly not pv
+            res_q_(gen_id) = 0.;
+            continue;
+        } 
+        if ((!turnedoff_gen_pv_) && p_mw_(gen_id) == 0.) {
+            // in this case turned off generators are not pv
+            res_q_(gen_id) = 0.;
+            continue;
+        }  
 
         int bus_id = bus_id_(gen_id);
         const auto bus_solver = id_grid_to_solver[bus_id];
@@ -444,7 +463,6 @@ void GeneratorContainer::set_q(const RealVect & reactive_mismatch,
     }
 }
 
-
 void GeneratorContainer::update_slack_weights(Eigen::Ref<Eigen::Array<bool, Eigen::Dynamic, Eigen::RowMajor> > could_be_slack,
                                               SolverControl & solver_control)
 {
diff --git a/src/element_container/GeneratorContainer.h b/src/element_container/GeneratorContainer.h
index c8c86e6..8f5e0a0 100644
--- a/src/element_container/GeneratorContainer.h
+++ b/src/element_container/GeneratorContainer.h
@@ -60,7 +60,7 @@ class GeneratorContainer: public GenericContainer
             id(-1),
             name(""),
             connected(false),
-            bus_id(-1),
+            bus_id(_deactivated_bus_id),
             is_slack(false),
             slack_weight(-1.0),
             voltage_regulator_on(false),
@@ -82,7 +82,7 @@ class GeneratorContainer: public GenericContainer
                         name = r_data_gen.names_[my_id];
                     }
                     connected = r_data_gen.status_[my_id];
-                    bus_id = r_data_gen.bus_id_[my_id];
+                    if(connected) bus_id = r_data_gen.bus_id_[my_id];
                     is_slack = r_data_gen.gen_slackbus_[my_id];
                     slack_weight = r_data_gen.gen_slack_weight_[my_id];
 
@@ -324,6 +324,7 @@ class GeneratorContainer: public GenericContainer
     virtual void gen_p_per_bus(std::vector<real_type> & res) const;
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
+    tuple4d get_res_full() const {return tuple4d(res_p_, res_q_, res_v_, res_theta_);}
     Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
 
     const std::vector<bool>& get_status() const {return status_;}
diff --git a/src/element_container/GenericContainer.cpp b/src/element_container/GenericContainer.cpp
index cf6dc51..e7c7ccd 100644
--- a/src/element_container/GenericContainer.cpp
+++ b/src/element_container/GenericContainer.cpp
@@ -14,7 +14,7 @@
 const int GenericContainer::_deactivated_bus_id = -1;
 
 // TODO all functions bellow are generic ! Make a base class for that
-void GenericContainer::_get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v){
+void GenericContainer::_get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v) const {
     const real_type _1_sqrt_3 = 1.0 / std::sqrt(3.);
     RealVect p2q2 = p.array() * p.array() + q.array() * q.array();
     p2q2 = p2q2.array().cwiseSqrt();
@@ -96,7 +96,7 @@ void GenericContainer::_change_bus(int el_id, int new_bus_me_id, Eigen::VectorXi
     bus_me_id = new_bus_me_id;
 }
 
-int GenericContainer::_get_bus(int el_id, const std::vector<bool> & status_, const Eigen::VectorXi & bus_id_)
+int GenericContainer::_get_bus(int el_id, const std::vector<bool> & status_, const Eigen::VectorXi & bus_id_) const
 {
     int res;
     bool val = status_.at(el_id);  // also check if the el_id is out of bound
@@ -114,12 +114,14 @@ void GenericContainer::v_kv_from_vpu(const Eigen::Ref<const RealVect> & Va,
                                      const Eigen::VectorXi & bus_me_id,
                                      const std::vector<int> & id_grid_to_solver,
                                      const RealVect & bus_vn_kv,
-                                     RealVect & v)
+                                     RealVect & v) const
 {
-    v = RealVect::Constant(nb_element, -1.0);
     for(int el_id = 0; el_id < nb_element; ++el_id){
         // if the element is disconnected, i leave it like that
-        if(!status[el_id]) continue;
+        if(!status[el_id]) {
+            v(el_id) = -1;
+            continue;
+        }
         int el_bus_me_id = bus_me_id(el_id);
         int bus_solver_id = id_grid_to_solver[el_bus_me_id];
         if(bus_solver_id == _deactivated_bus_id){
@@ -135,7 +137,6 @@ void GenericContainer::v_kv_from_vpu(const Eigen::Ref<const RealVect> & Va,
     }
 }
 
-
 void GenericContainer::v_deg_from_va(const Eigen::Ref<const RealVect> & Va,
                                      const Eigen::Ref<const RealVect> & Vm,
                                      const std::vector<bool> & status,
@@ -143,11 +144,14 @@ void GenericContainer::v_deg_from_va(const Eigen::Ref<const RealVect> & Va,
                                      const Eigen::VectorXi & bus_me_id,
                                      const std::vector<int> & id_grid_to_solver,
                                      const RealVect & bus_vn_kv,
-                                     RealVect & theta){
-    theta = RealVect::Constant(nb_element, 0.0);
+                                     RealVect & theta) const
+{
     for(int el_id = 0; el_id < nb_element; ++el_id){
         // if the element is disconnected, i leave it like that
-        if(!status[el_id]) continue;
+        if(!status[el_id]) {
+            theta(el_id) = -1;
+            continue;
+        }
         int el_bus_me_id = bus_me_id(el_id);
         int bus_solver_id = id_grid_to_solver[el_bus_me_id];
         if(bus_solver_id == _deactivated_bus_id){
@@ -158,6 +162,6 @@ void GenericContainer::v_deg_from_va(const Eigen::Ref<const RealVect> & Va,
             exc_ << " is connected to a disconnected bus";
             throw std::runtime_error(exc_.str());
         }
-        theta(el_id) = Va(bus_solver_id) * 180. / my_pi;
+        theta(el_id) = Va(bus_solver_id) * my_180_pi_;
     }
 }
diff --git a/src/element_container/GenericContainer.h b/src/element_container/GenericContainer.h
index b1fa679..98050a0 100644
--- a/src/element_container/GenericContainer.h
+++ b/src/element_container/GenericContainer.h
@@ -83,7 +83,9 @@ class GenericContainer : public BaseConstants
                                      int nb_line,
                                      bool transpose) const {};
 
-        virtual void fillYbus(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver) {};
+        // no more used !
+        virtual void fillYbus(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver) const {};
+
         virtual void fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const {};
         virtual void fillpv(std::vector<int>& bus_pv,
                             std::vector<bool> & has_bus_been_added,
@@ -120,12 +122,12 @@ class GenericContainer : public BaseConstants
         void _reactivate(int el_id, std::vector<bool> & status);
         void _deactivate(int el_id, std::vector<bool> & status);
         void _change_bus(int el_id, int new_bus_me_id, Eigen::VectorXi & el_bus_ids, SolverControl & solver_control, int nb_bus);
-        int _get_bus(int el_id, const std::vector<bool> & status_, const Eigen::VectorXi & bus_id_);
+        int _get_bus(int el_id, const std::vector<bool> & status_, const Eigen::VectorXi & bus_id_) const;
 
         /**
         compute the amps from the p, the q and the v (v should NOT be pair unit)
         **/
-        void _get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v);
+        void _get_amps(RealVect & a, const RealVect & p, const RealVect & q, const RealVect & v) const;
 
         /**
         convert v from pu to v in kv (and assign it to the right element...)
@@ -137,7 +139,7 @@ class GenericContainer : public BaseConstants
                            const Eigen::VectorXi & bus_me_id,
                            const std::vector<int> & id_grid_to_solver,
                            const RealVect & bus_vn_kv,
-                           RealVect & v);
+                           RealVect & v) const;
 
 
         /**
@@ -150,13 +152,13 @@ class GenericContainer : public BaseConstants
                            const Eigen::VectorXi & bus_me_id,
                            const std::vector<int> & id_grid_to_solver,
                            const RealVect & bus_vn_kv,
-                           RealVect & v);
+                           RealVect & v) const;
 
         /**
         check the size of the elements
         **/
         template<class T, class intType>
-        void check_size(const T & container, intType size, const std::string & container_name)
+        void check_size(const T & container, intType size, const std::string & container_name) const
         {
             if(static_cast<intType>(container.size()) != size) throw std::runtime_error(container_name + " do not have the proper size");
         }
diff --git a/src/element_container/LineContainer.cpp b/src/element_container/LineContainer.cpp
index efeac7d..fdcf6bc 100644
--- a/src/element_container/LineContainer.cpp
+++ b/src/element_container/LineContainer.cpp
@@ -44,6 +44,7 @@ void LineContainer::init(const RealVect & branch_r,
     powerlines_x_ = branch_x;
     status_ = std::vector<bool>(branch_r.size(), true); // by default everything is connected
     _update_model_coeffs();
+    reset_results();
 }
 
 void LineContainer::init(const RealVect & branch_r,
@@ -82,6 +83,7 @@ void LineContainer::init(const RealVect & branch_r,
     powerlines_x_ = branch_x;
     status_ = std::vector<bool>(branch_r.size(), true); // by default everything is connected
     _update_model_coeffs();
+    reset_results();
 }
 
 LineContainer::StateRes LineContainer::get_state() const
@@ -98,7 +100,6 @@ LineContainer::StateRes LineContainer::get_state() const
 }
 void LineContainer::set_state(LineContainer::StateRes & my_state)
 {
-    reset_results();
     names_ = std::get<0>(my_state);
     std::vector<real_type> & branch_r = std::get<1>(my_state);
     std::vector<real_type> & branch_x = std::get<2>(my_state);
@@ -121,6 +122,7 @@ void LineContainer::set_state(LineContainer::StateRes & my_state)
     status_ = status;
 
     _update_model_coeffs();
+    reset_results();
 }
 
 void LineContainer::_update_model_coeffs()
@@ -355,14 +357,16 @@ void LineContainer::fillBf_for_PTDF(std::vector<Eigen::Triplet<real_type> > & Bf
 
 void LineContainer::reset_results()
 {
-    res_powerline_por_ = RealVect();  // in MW
-    res_powerline_qor_ = RealVect();  // in MVar
-    res_powerline_vor_ = RealVect();  // in kV
-    res_powerline_aor_ = RealVect();  // in kA
-    res_powerline_pex_ = RealVect();  // in MW
-    res_powerline_qex_ = RealVect();  // in MVar
-    res_powerline_vex_ = RealVect();  // in kV
-    res_powerline_aex_ = RealVect();  // in kA
+    res_powerline_por_ = RealVect(nb());  // in MW
+    res_powerline_qor_ = RealVect(nb());  // in MVar
+    res_powerline_vor_ = RealVect(nb());  // in kV
+    res_powerline_aor_ = RealVect(nb());  // in kA
+    res_powerline_pex_ = RealVect(nb());  // in MW
+    res_powerline_qex_ = RealVect(nb());  // in MVar
+    res_powerline_vex_ = RealVect(nb());  // in kV
+    res_powerline_aex_ = RealVect(nb());  // in kA
+    res_powerline_thetaor_ = RealVect(nb());
+    res_powerline_thetaex_ = RealVect(nb());
 }
 
 void LineContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
@@ -375,19 +379,21 @@ void LineContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
 {
     // it needs to be initialized at 0.
     Eigen::Index nb_element = nb();
-    res_powerline_por_ = RealVect::Constant(nb_element, my_zero_);  // in MW
-    res_powerline_qor_ = RealVect::Constant(nb_element, my_zero_);  // in MVar
-    res_powerline_vor_ = RealVect::Constant(nb_element, my_zero_);  // in kV
-    res_powerline_aor_ = RealVect::Constant(nb_element, my_zero_);  // in kA
-    res_powerline_pex_ = RealVect::Constant(nb_element, my_zero_);  // in MW
-    res_powerline_qex_ = RealVect::Constant(nb_element, my_zero_);  // in MVar
-    res_powerline_vex_ = RealVect::Constant(nb_element, my_zero_);  // in kV
-    res_powerline_aex_ = RealVect::Constant(nb_element, my_zero_);  // in kA
-    res_powerline_thetaor_ = RealVect::Constant(nb_element, my_zero_);  // in kV
-    res_powerline_thetaex_ = RealVect::Constant(nb_element, my_zero_);  // in kV
     for(Eigen::Index line_id = 0; line_id < nb_element; ++line_id){
         // don't do anything if the element is disconnected
-        if(!status_[line_id]) continue;
+        if(!status_[line_id]) {
+            res_powerline_por_(line_id) = 0.0;  // in MW
+            res_powerline_qor_(line_id) = 0.0;  // in MVar
+            res_powerline_vor_(line_id) = 0.0;  // in kV
+            res_powerline_aor_(line_id) = 0.0;  // in kA
+            res_powerline_pex_(line_id) = 0.0;  // in MW
+            res_powerline_qex_(line_id) = 0.0;  // in MVar
+            res_powerline_vex_(line_id) = 0.0;  // in kV
+            res_powerline_aex_(line_id) = 0.0;  // in kA
+            res_powerline_thetaor_(line_id) = 0.0;  // in kV
+            res_powerline_thetaex_(line_id) = 0.0;  // in kV
+            continue;
+        }
 
         // connectivity
         int bus_or_id_me = bus_or_id_(line_id);
@@ -420,8 +426,8 @@ void LineContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
         res_powerline_vex_(line_id) = v_ex * bus_vn_kv_ex;
 
         // retrieve the voltage angle in degree (instead of radian)
-        res_powerline_thetaor_(line_id) = Va(bus_or_solver_id) * 180. / my_pi;
-        res_powerline_thetaex_(line_id) = Va(bus_ex_solver_id) * 180. / my_pi;
+        res_powerline_thetaor_(line_id) = Va(bus_or_solver_id) * my_180_pi_;
+        res_powerline_thetaex_(line_id) = Va(bus_ex_solver_id) * my_180_pi_;
 
         // results of the powerflow
         cplx_type Eor = V(bus_or_solver_id);
@@ -447,6 +453,8 @@ void LineContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
             res_powerline_por_(line_id) = (std::real(ydc_ff_(line_id)) * Va(bus_or_solver_id) + std::real(ydc_ft_(line_id)) * Va(bus_ex_solver_id)) * sn_mva;
             res_powerline_pex_(line_id) = (std::real(ydc_tt_(line_id)) * Va(bus_ex_solver_id) + std::real(ydc_tf_(line_id)) * Va(bus_or_solver_id)) * sn_mva;   
 
+            res_powerline_qor_(line_id) = 0.;
+            res_powerline_qex_(line_id) = 0.;
             // for the voltage (by hypothesis vm = 1)
             // res_powerline_vor_(line_id) = bus_vn_kv_or;
             // res_powerline_vex_(line_id) = bus_vn_kv_ex;        
diff --git a/src/element_container/LineContainer.h b/src/element_container/LineContainer.h
index 7e6acf4..1fc1470 100644
--- a/src/element_container/LineContainer.h
+++ b/src/element_container/LineContainer.h
@@ -57,8 +57,8 @@ class LineContainer : public GenericContainer
                 id(my_id),
                 name(""),
                 connected(false),
-                bus_or_id(-1),
-                bus_ex_id(-1),
+                bus_or_id(_deactivated_bus_id),
+                bus_ex_id(_deactivated_bus_id),
                 r_pu(-1.0),
                 x_pu(-1.0),
                 h_pu(0., 0.),
@@ -83,8 +83,8 @@ class LineContainer : public GenericContainer
                             name = r_data_line.names_[my_id];
                         }
                         connected = r_data_line.status_[my_id];
-                        bus_or_id = r_data_line.bus_or_id_.coeff(my_id);
-                        bus_ex_id = r_data_line.bus_ex_id_.coeff(my_id);
+                        if(connected)  bus_or_id = r_data_line.bus_or_id_.coeff(my_id);
+                        if(connected)  bus_ex_id = r_data_line.bus_ex_id_.coeff(my_id);
                         r_pu = r_data_line.powerlines_r_.coeff(my_id);
                         x_pu = r_data_line.powerlines_x_.coeff(my_id);
                         h_or_pu = r_data_line.powerlines_h_or_.coeff(my_id);
@@ -247,7 +247,8 @@ class LineContainer : public GenericContainer
 
     tuple4d get_lineor_res() const {return tuple4d(res_powerline_por_, res_powerline_qor_, res_powerline_vor_, res_powerline_aor_);}
     tuple4d get_lineex_res() const {return tuple4d(res_powerline_pex_, res_powerline_qex_, res_powerline_vex_, res_powerline_aex_);}
-
+    tuple5d get_res_or_full() const {return tuple5d(res_powerline_por_, res_powerline_qor_, res_powerline_vor_, res_powerline_aor_, res_powerline_thetaor_);}
+    tuple5d get_res_ex_full() const {return tuple5d(res_powerline_pex_, res_powerline_qex_, res_powerline_vex_, res_powerline_aex_, res_powerline_thetaex_);}
     Eigen::Ref<const RealVect> get_theta_or() const {return res_powerline_thetaor_;}
     Eigen::Ref<const RealVect> get_theta_ex() const {return res_powerline_thetaex_;}
     const std::vector<bool>& get_status() const {return status_;}
diff --git a/src/element_container/LoadContainer.cpp b/src/element_container/LoadContainer.cpp
index d9d0a1b..e7c7483 100644
--- a/src/element_container/LoadContainer.cpp
+++ b/src/element_container/LoadContainer.cpp
@@ -8,6 +8,7 @@
 
 #include "LoadContainer.h"
 #include <sstream>
+#include <iostream>
 
 void LoadContainer::init(const RealVect & loads_p,
                          const RealVect & loads_q,
@@ -22,9 +23,9 @@ void LoadContainer::init(const RealVect & loads_p,
     q_mvar_ = loads_q;
     bus_id_ = loads_bus_id;
     status_ = std::vector<bool>(loads_p.size(), true);
+    reset_results();
 }
 
-
 LoadContainer::StateRes LoadContainer::get_state() const
 {
      std::vector<real_type> p_mw(p_mw_.begin(), p_mw_.end());
@@ -37,7 +38,6 @@ LoadContainer::StateRes LoadContainer::get_state() const
 
 void LoadContainer::set_state(LoadContainer::StateRes & my_state )
 {
-    reset_results();
     names_ = std::get<0>(my_state);
     std::vector<real_type> & p_mw = std::get<1>(my_state);
     std::vector<real_type> & q_mvar = std::get<2>(my_state);
@@ -50,9 +50,9 @@ void LoadContainer::set_state(LoadContainer::StateRes & my_state )
     q_mvar_ = RealVect::Map(&q_mvar[0], q_mvar.size());
     bus_id_ = Eigen::VectorXi::Map(&bus_id[0], bus_id.size());
     status_ = status;
+    reset_results();
 }
 
-
 void LoadContainer::fillSbus(CplxVect & Sbus,
                              const std::vector<int> & id_grid_to_solver,
                              bool ac) const
@@ -87,17 +87,23 @@ void LoadContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
                                     real_type sn_mva,
                                     bool ac)
 {
-    int nb_load = nb();
+    const int nb_load = nb();
     v_kv_from_vpu(Va, Vm, status_, nb_load, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_);
     v_deg_from_va(Va, Vm, status_, nb_load, bus_id_, id_grid_to_solver, bus_vn_kv, res_theta_);
     res_p_ = p_mw_;
-    res_q_ = q_mvar_;
+    if(ac) res_q_ = q_mvar_;
+    else{
+        // no q in DC mode
+        for(int load_id = 0; load_id < nb_load; ++load_id) res_q_(load_id) = 0.;
+    }
 }
 
 void LoadContainer::reset_results(){
-    res_p_ = RealVect();  // in MW
-    res_q_ =  RealVect();  // in MVar
-    res_v_ = RealVect();  // in kV
+    // std::cout << "Loads reset_results \n";
+    res_p_ = RealVect(nb());  // in MW
+    res_q_ =  RealVect(nb());  // in MVar
+    res_v_ = RealVect(nb());  // in kV
+    res_theta_ = RealVect(nb());  // in deg
 }
 
 void LoadContainer::change_p(int load_id, real_type new_p, SolverControl & solver_control)
diff --git a/src/element_container/LoadContainer.h b/src/element_container/LoadContainer.h
index 0483f3f..0dc86cd 100644
--- a/src/element_container/LoadContainer.h
+++ b/src/element_container/LoadContainer.h
@@ -60,7 +60,7 @@ class LoadContainer : public GenericContainer
                 id(-1),
                 name(""),
                 connected(false),
-                bus_id(-1),
+                bus_id(_deactivated_bus_id),
                 target_p_mw(0.),
                 target_q_mvar(0.),
                 has_res(false),
@@ -76,7 +76,7 @@ class LoadContainer : public GenericContainer
                             name = r_data_load.names_[my_id];
                         }
                         connected = r_data_load.status_[my_id];
-                        bus_id = r_data_load.bus_id_[my_id];
+                        if(connected) bus_id = r_data_load.bus_id_[my_id];
 
                         target_p_mw = r_data_load.p_mw_.coeff(my_id);
                         target_q_mvar = r_data_load.q_mvar_.coeff(my_id);
@@ -177,6 +177,8 @@ class LoadContainer : public GenericContainer
     void reset_results();
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
+    tuple4d get_res_full() const {return tuple4d(res_p_, res_q_, res_v_, res_theta_);}
+    
     Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
     const std::vector<bool>& get_status() const {return status_;}
     const Eigen::VectorXi & get_bus_id() const {return bus_id_;}
diff --git a/src/element_container/SGenContainer.cpp b/src/element_container/SGenContainer.cpp
index 6300727..8690b74 100644
--- a/src/element_container/SGenContainer.cpp
+++ b/src/element_container/SGenContainer.cpp
@@ -34,6 +34,7 @@ void SGenContainer::init(const RealVect & sgen_p,
     q_max_mvar_ = sgen_qmax;
     bus_id_ = sgen_bus_id;
     status_ = std::vector<bool>(sgen_p.size(), true);
+    reset_results();
 }
 
 SGenContainer::StateRes SGenContainer::get_state() const
@@ -51,9 +52,7 @@ SGenContainer::StateRes SGenContainer::get_state() const
 }
 
 void SGenContainer::set_state(SGenContainer::StateRes & my_state )
-{
-    reset_results();
-    
+{    
     names_ = std::get<0>(my_state);
     std::vector<real_type> & p_mw = std::get<1>(my_state);
     std::vector<real_type> & q_mvar = std::get<2>(my_state);
@@ -82,6 +81,7 @@ void SGenContainer::set_state(SGenContainer::StateRes & my_state )
     q_max_mvar_ = RealVect::Map(&q_max[0], size);
     bus_id_ = Eigen::VectorXi::Map(&bus_id[0], bus_id.size());
     status_ = status;
+    reset_results();
 }
 
 void SGenContainer::fillSbus(CplxVect & Sbus, const std::vector<int> & id_grid_to_solver, bool ac) const {
@@ -119,13 +119,17 @@ void SGenContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
     v_deg_from_va(Va, Vm, status_, nb_sgen, bus_id_, id_grid_to_solver, bus_vn_kv, res_theta_);
     res_p_ = p_mw_;
     if(ac) res_q_ = q_mvar_;
-    else res_q_ = RealVect::Zero(nb_sgen);
+    else{
+        // no q in DC mode
+        for(int sgen_id = 0; sgen_id < nb_sgen; ++sgen_id) res_q_(sgen_id) = 0.;
+    }
 }
 
 void SGenContainer::reset_results(){
-    res_p_ = RealVect();  // in MW
-    res_q_ =  RealVect();  // in MVar
-    res_v_ = RealVect();  // in kV
+    res_p_ = RealVect(nb());  // in MW
+    res_q_ =  RealVect(nb());  // in MVar
+    res_v_ = RealVect(nb());  // in kV
+    res_theta_ = RealVect(nb());  // in deg
 }
 
 void SGenContainer::change_p(int sgen_id, real_type new_p, SolverControl & solver_control)
diff --git a/src/element_container/SGenContainer.h b/src/element_container/SGenContainer.h
index b50a269..6a3f831 100644
--- a/src/element_container/SGenContainer.h
+++ b/src/element_container/SGenContainer.h
@@ -64,7 +64,7 @@ class SGenContainer: public GenericContainer
                 id(-1),
                 name(""),
                 connected(false),
-                bus_id(-1),
+                bus_id(_deactivated_bus_id),
                 min_q_mvar(0.),
                 max_q_mvar(0.),
                 min_p_mw(0.),
@@ -84,7 +84,7 @@ class SGenContainer: public GenericContainer
                             name = r_data_sgen.names_[my_id];
                         }
                         connected = r_data_sgen.status_[my_id];
-                        bus_id = r_data_sgen.bus_id_[my_id];
+                        if(connected)  bus_id = r_data_sgen.bus_id_[my_id];
 
                         min_q_mvar = r_data_sgen.q_min_mvar_(my_id);
                         max_q_mvar = r_data_sgen.q_max_mvar_(my_id);
@@ -198,6 +198,7 @@ class SGenContainer: public GenericContainer
     void reset_results();
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
+    tuple4d get_res_full() const {return tuple4d(res_p_, res_q_, res_v_, res_theta_);}
     Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
     const std::vector<bool>& get_status() const {return status_;}
     const Eigen::VectorXi & get_bus_id() const {return bus_id_;}
diff --git a/src/element_container/ShuntContainer.cpp b/src/element_container/ShuntContainer.cpp
index 4b5687b..4bd38a5 100644
--- a/src/element_container/ShuntContainer.cpp
+++ b/src/element_container/ShuntContainer.cpp
@@ -23,6 +23,7 @@ void ShuntContainer::init(const RealVect & shunt_p_mw,
     q_mvar_ = shunt_q_mvar;
     bus_id_ = shunt_bus_id;
     status_ = std::vector<bool>(p_mw_.size(), true); // by default everything is connected
+    reset_results();
 }
 
 ShuntContainer::StateRes ShuntContainer::get_state() const
@@ -37,7 +38,6 @@ ShuntContainer::StateRes ShuntContainer::get_state() const
 
 void ShuntContainer::set_state(ShuntContainer::StateRes & my_state )
 {
-    reset_results();
     names_ = std::get<0>(my_state);
     std::vector<real_type> & p_mw = std::get<1>(my_state);
     std::vector<real_type> & q_mvar = std::get<2>(my_state);
@@ -50,6 +50,7 @@ void ShuntContainer::set_state(ShuntContainer::StateRes & my_state )
     q_mvar_ = RealVect::Map(&q_mvar[0], q_mvar.size());
     bus_id_ = Eigen::VectorXi::Map(&bus_id[0], bus_id.size());
     status_ = status;
+    reset_results();
 }
 
 void ShuntContainer::fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
@@ -147,10 +148,14 @@ void ShuntContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
     const int nb_shunt = static_cast<int>(p_mw_.size());
     v_kv_from_vpu(Va, Vm, status_, nb_shunt, bus_id_, id_grid_to_solver, bus_vn_kv, res_v_);
     v_deg_from_va(Va, Vm, status_, nb_shunt, bus_id_, id_grid_to_solver, bus_vn_kv, res_theta_);
-    res_p_ = RealVect::Constant(nb_shunt, my_zero_);
-    res_q_ = RealVect::Constant(nb_shunt, my_zero_);
+    // res_p_ = RealVect::Constant(nb_shunt, my_zero_);
+    // res_q_ = RealVect::Constant(nb_shunt, my_zero_);
     for(int shunt_id = 0; shunt_id < nb_shunt; ++shunt_id){
-        if(!status_[shunt_id]) continue;
+        if(!status_[shunt_id]) {
+            res_p_(shunt_id) = my_zero_;
+            res_q_(shunt_id) = my_zero_;
+            continue;
+        }
         int bus_id_me = bus_id_(shunt_id);
         int bus_solver_id = id_grid_to_solver[bus_id_me];
         if(bus_solver_id == _deactivated_bus_id){
@@ -163,13 +168,15 @@ void ShuntContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
         cplx_type s = E * I;
         res_p_(shunt_id) = std::real(s) * sn_mva;
         if(ac) res_q_(shunt_id) = std::imag(s) * sn_mva;
+        else res_q_(shunt_id) = my_zero_;
     }
 }
 
 void ShuntContainer::reset_results(){
-    res_p_ = RealVect();  // in MW
-    res_q_ = RealVect();  // in MVar
-    res_v_ = RealVect();  // in kV
+    res_p_ = RealVect(nb());  // in MW
+    res_q_ = RealVect(nb());  // in MVar
+    res_v_ = RealVect(nb());  // in kV
+    res_theta_ = RealVect(nb());  // in deg
 }
 
 void ShuntContainer::change_p(int shunt_id, real_type new_p, SolverControl & solver_control)
diff --git a/src/element_container/ShuntContainer.h b/src/element_container/ShuntContainer.h
index 8458acf..59c4a8a 100644
--- a/src/element_container/ShuntContainer.h
+++ b/src/element_container/ShuntContainer.h
@@ -52,7 +52,7 @@ class ShuntContainer : public GenericContainer
                 id(-1),
                 name(""),
                 connected(false),
-                bus_id(-1),
+                bus_id(_deactivated_bus_id),
                 target_p_mw(0.),
                 target_q_mvar(0.),
                 has_res(false),
@@ -68,7 +68,7 @@ class ShuntContainer : public GenericContainer
                             name = r_data_shunt.names_[my_id];
                         }
                         connected = r_data_shunt.status_[my_id];
-                        bus_id = r_data_shunt.bus_id_[my_id];
+                        if(connected)  bus_id = r_data_shunt.bus_id_[my_id];
 
                         target_p_mw = r_data_shunt.p_mw_.coeff(my_id);
                         target_q_mvar = r_data_shunt.q_mvar_.coeff(my_id);
@@ -146,7 +146,9 @@ class ShuntContainer : public GenericContainer
     void change_bus(int shunt_id, int new_bus_id, SolverControl & solver_control, int nb_bus) {_change_bus(shunt_id, new_bus_id, bus_id_, solver_control, nb_bus);}
     void change_p(int shunt_id, real_type new_p, SolverControl & solver_control);
     void change_q(int shunt_id, real_type new_q, SolverControl & solver_control);
-    int get_bus(int shunt_id) {return _get_bus(shunt_id, status_, bus_id_);}
+    int get_bus(int shunt_id) const {return _get_bus(shunt_id, status_, bus_id_);}
+    Eigen::Ref<const IntVect> get_buses() const {return bus_id_;}
+
     virtual void reconnect_connected_buses(std::vector<bool> & bus_status) const;
     virtual void disconnect_if_not_in_main_component(std::vector<bool> & busbar_in_main_component);
     
@@ -180,6 +182,8 @@ class ShuntContainer : public GenericContainer
     void reset_results();
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
+    tuple4d get_res_full() const {return tuple4d(res_p_, res_q_, res_v_, res_theta_);}
+    
     Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
     const std::vector<bool>& get_status() const {return status_;}
 
diff --git a/src/element_container/TrafoContainer.cpp b/src/element_container/TrafoContainer.cpp
index f33e1cd..488c2e7 100644
--- a/src/element_container/TrafoContainer.cpp
+++ b/src/element_container/TrafoContainer.cpp
@@ -46,12 +46,13 @@ void TrafoContainer::init(const RealVect & trafo_r,
     x_ = trafo_x;
     h_ = trafo_b;
     ratio_ = ratio;
-    shift_ = trafo_shift_degree / 180. * my_pi;  // do not forget conversion degree / rad here !
+    shift_ = trafo_shift_degree / my_180_pi_;  // do not forget conversion degree / rad here !
     bus_hv_id_ = trafo_hv_id;
     bus_lv_id_ = trafo_lv_id;
     is_tap_hv_side_ = trafo_tap_hv;
     status_ = std::vector<bool>(trafo_r.size(), true);
     _update_model_coeffs();
+    reset_results();
 
 }
 
@@ -74,8 +75,6 @@ TrafoContainer::StateRes TrafoContainer::get_state() const
 
 void TrafoContainer::set_state(TrafoContainer::StateRes & my_state)
 {
-    reset_results();
-
     names_ = std::get<0>(my_state);
     std::vector<real_type> & branch_r = std::get<1>(my_state);
     std::vector<real_type> & branch_x = std::get<2>(my_state);
@@ -111,6 +110,7 @@ void TrafoContainer::set_state(TrafoContainer::StateRes & my_state)
     shift_  = RealVect::Map(&shift[0], size);
     is_tap_hv_side_ = is_tap_hv_side;
     _update_model_coeffs();
+    reset_results();
 }
 
 
@@ -271,19 +271,21 @@ void TrafoContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
 {
     // it needs to be initialized at 0.
     const int nb_element = nb();
-    res_p_hv_ = RealVect::Constant(nb_element, 0.0);  // in MW
-    res_q_hv_ = RealVect::Constant(nb_element, 0.0);  // in MVar
-    res_v_hv_ = RealVect::Constant(nb_element, 0.0);  // in kV
-    res_a_hv_ = RealVect::Constant(nb_element, 0.0);  // in kA
-    res_p_lv_ = RealVect::Constant(nb_element, 0.0);  // in MW
-    res_q_lv_ = RealVect::Constant(nb_element, 0.0);  // in MVar
-    res_v_lv_ = RealVect::Constant(nb_element, 0.0);  // in kV
-    res_a_lv_ = RealVect::Constant(nb_element, 0.0);  // in kA
-    res_theta_hv_ = RealVect::Constant(nb_element, 0.0);  // in degree
-    res_theta_lv_ = RealVect::Constant(nb_element, 0.0);  // in degree
     for(int trafo_id = 0; trafo_id < nb_element; ++trafo_id){
         // don't do anything if the element is disconnected
-        if(!status_[trafo_id]) continue;
+        if(!status_[trafo_id]) {
+            res_p_hv_(trafo_id) = 0.0;  // in MW
+            res_q_hv_(trafo_id) = 0.0;  // in MVar
+            res_v_hv_(trafo_id) = 0.0;  // in kV
+            res_a_hv_(trafo_id) = 0.0;  // in kA
+            res_p_lv_(trafo_id) = 0.0;  // in MW
+            res_q_lv_(trafo_id) = 0.0;  // in MVar
+            res_v_lv_(trafo_id) = 0.0;  // in kV
+            res_a_lv_(trafo_id) = 0.0;  // in kA
+            res_theta_hv_(trafo_id) = 0.0;  // in degree
+            res_theta_lv_(trafo_id) = 0.0;  // in degree            
+            continue;
+        }
 
         // connectivity
         int bus_hv_id_me = bus_hv_id_(trafo_id);
@@ -315,8 +317,8 @@ void TrafoContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
         res_v_hv_(trafo_id) = v_hv * bus_vn_kv_hv;
         res_v_lv_(trafo_id) = v_lv * bus_vn_kv_lv;
 
-        res_theta_hv_(trafo_id) = Va(bus_hv_solver_id) * 180. / my_pi;
-        res_theta_lv_(trafo_id) = Va(bus_lv_solver_id) * 180. / my_pi;
+        res_theta_hv_(trafo_id) = Va(bus_hv_solver_id) * my_180_pi_;
+        res_theta_lv_(trafo_id) = Va(bus_lv_solver_id) * my_180_pi_;
 
         if(ac){
             // results of the ac powerflow
@@ -341,7 +343,10 @@ void TrafoContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
             // result of the dc powerflow
             res_p_hv_(trafo_id) = (std::real(ydc_ff_(trafo_id)) * Va(bus_hv_solver_id) + std::real(ydc_ft_(trafo_id)) * Va(bus_lv_solver_id) - dc_x_tau_shift_(trafo_id) ) * sn_mva;
             res_p_lv_(trafo_id) = (std::real(ydc_tt_(trafo_id)) * Va(bus_lv_solver_id) + std::real(ydc_tf_(trafo_id)) * Va(bus_hv_solver_id) + dc_x_tau_shift_(trafo_id) ) * sn_mva; 
-
+           
+            res_q_hv_(trafo_id) = 0.;
+            res_q_lv_(trafo_id) = 0.;
+            
             // for voltages, because vm = 1. pu by hypothesis
             // res_v_hv_(trafo_id) = bus_vn_kv_hv;
             // res_v_lv_(trafo_id) = bus_vn_kv_lv;
@@ -353,14 +358,16 @@ void TrafoContainer::compute_results(const Eigen::Ref<const RealVect> & Va,
 }
 
 void TrafoContainer::reset_results(){
-    res_p_hv_ = RealVect();  // in MW
-    res_q_hv_ = RealVect();  // in MVar
-    res_v_hv_ = RealVect();  // in kV
-    res_a_hv_ = RealVect();  // in kA
-    res_p_lv_ = RealVect();  // in MW
-    res_q_lv_ = RealVect();  // in MVar
-    res_v_lv_ = RealVect();  // in kV
-    res_a_lv_ = RealVect();  // in kA
+    res_p_hv_ = RealVect(nb());  // in MW
+    res_q_hv_ = RealVect(nb());  // in MVar
+    res_v_hv_ = RealVect(nb());  // in kV
+    res_a_hv_ = RealVect(nb());  // in kA
+    res_p_lv_ = RealVect(nb());  // in MW
+    res_q_lv_ = RealVect(nb());  // in MVar
+    res_v_lv_ = RealVect(nb());  // in kV
+    res_a_lv_ = RealVect(nb());  // in kA
+    res_theta_hv_ = RealVect(nb());
+    res_theta_lv_ = RealVect(nb());
 }
 
 
diff --git a/src/element_container/TrafoContainer.h b/src/element_container/TrafoContainer.h
index fbb07f8..50338e1 100644
--- a/src/element_container/TrafoContainer.h
+++ b/src/element_container/TrafoContainer.h
@@ -64,8 +64,8 @@ class TrafoContainer : public GenericContainer
                 id(-1),
                 name(""),
                 connected(false),
-                bus_hv_id(-1),
-                bus_lv_id(-1),
+                bus_hv_id(_deactivated_bus_id),
+                bus_lv_id(_deactivated_bus_id),
                 r_pu(-1.0),
                 x_pu(-1.0),
                 h_pu(0., 0.),
@@ -91,8 +91,8 @@ class TrafoContainer : public GenericContainer
                             name = r_data_trafo.names_[my_id];
                         }
                         connected = r_data_trafo.status_[my_id];
-                        bus_hv_id = r_data_trafo.bus_hv_id_.coeff(my_id);
-                        bus_lv_id = r_data_trafo.bus_lv_id_.coeff(my_id);
+                        if(connected)  bus_hv_id = r_data_trafo.bus_hv_id_.coeff(my_id);
+                        if(connected)  bus_lv_id = r_data_trafo.bus_lv_id_.coeff(my_id);
                         r_pu = r_data_trafo.r_.coeff(my_id);
                         x_pu = r_data_trafo.x_.coeff(my_id);
                         h_pu = r_data_trafo.h_.coeff(my_id);
@@ -237,6 +237,8 @@ class TrafoContainer : public GenericContainer
 
     tuple4d get_res_hv() const {return tuple4d(res_p_hv_, res_q_hv_, res_v_hv_, res_a_hv_);}
     tuple4d get_res_lv() const {return tuple4d(res_p_lv_, res_q_lv_, res_v_lv_, res_a_lv_);}
+    tuple5d get_res_hv_full() const {return tuple5d(res_p_hv_, res_q_hv_, res_v_hv_, res_a_hv_, res_theta_hv_);}
+    tuple5d get_res_lv_full() const {return tuple5d(res_p_lv_, res_q_lv_, res_v_lv_, res_a_lv_, res_theta_lv_);}
     Eigen::Ref<const RealVect> get_theta_hv() const {return res_theta_hv_;}
     Eigen::Ref<const RealVect> get_theta_lv() const {return res_theta_lv_;}
     Eigen::Ref<const Eigen::VectorXi> get_bus_from() const {return bus_hv_id_;}
diff --git a/src/main.cpp b/src/main.cpp
index dd001e1..b8aac5b 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,4 +1,4 @@
-// Copyright (c) 2020, RTE (https://www.rte-france.com)
+// Copyright (c) 2020-2024, RTE (https://www.rte-france.com)
 // See AUTHORS.txt
 // This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
 // If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
@@ -19,11 +19,62 @@
 
 #include "help_fun_msg.h"
 
+#ifndef KLU_SOLVER_AVAILABLE
+#define this_KLU_SOLVER_AVAILABLE 0
+#else 
+#define this_KLU_SOLVER_AVAILABLE 1
+#endif
+#ifndef NICSLU_SOLVER_AVAILABLE
+#define this_NICSLU_SOLVER_AVAILABLE 0
+#else 
+#define this_NICSLU_SOLVER_AVAILABLE 1
+#endif
+#ifndef CKTSO_SOLVER_AVAILABLE
+#define this_CKTSO_SOLVER_AVAILABLE 0
+#else 
+#define this_CKTSO_SOLVER_AVAILABLE 1
+#endif
+#ifndef __COMPILE_MARCHNATIVE
+#define this__COMPILE_MARCHNATIVE 0
+#else 
+#define this__COMPILE_MARCHNATIVE 1
+#endif
+#ifndef __O3_OPTIM
+#define this__O3_OPTIM 0
+#else 
+#define this__O3_OPTIM 1
+#endif
+#ifndef VERSION
+#define this_VERSION "unknown"
+#else 
+#define this_VERSION VERSION
+#endif
+#ifdef NICSLU_PATH
+#define this_NICSLU_PATH NICSLU_PATH
+#endif
+#ifdef CKTSO_PATH
+#define this_CKTSO_PATH CKTSO_PATH
+#endif
+
 namespace py = pybind11;
 
 PYBIND11_MODULE(lightsim2grid_cpp, m)
 {
 
+    // constant and compilation information
+    m.attr("klu_solver_available") = py::bool_(this_KLU_SOLVER_AVAILABLE);
+    m.attr("nicslu_solver_available") = py::bool_(this_NICSLU_SOLVER_AVAILABLE);
+    m.attr("cktso_solver_available") = py::bool_(this_CKTSO_SOLVER_AVAILABLE);
+    m.attr("compiled_march_native") = py::bool_(this__COMPILE_MARCHNATIVE);
+    m.attr("compiled_o3_optim") = py::bool_(this__O3_OPTIM);
+    m.attr("version") = py::str(this_VERSION);
+    #ifdef NICSLU_PATH
+    m.attr("nicslu_lib") = py::str(this_NICSLU_PATH);
+    #endif
+    #ifdef CKTSO_PATH
+    m.attr("cktso_lib") = py::str(this_CKTSO_PATH);
+    #endif
+    
     // solver method for FDPF
     py::enum_<FDPFMethod>(m, "FDPFMethod", "This enum controls the type of method you can use for Fast Decoupled Powerflow (XB or BX)")
         .value("XB", FDPFMethod::XB, "denotes the XB method")
@@ -393,6 +444,8 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("get_nb_iter", &ChooseSolver::get_nb_iter, DocSolver::get_nb_iter.c_str()) 
         .def("converged", &ChooseSolver::converged, DocSolver::converged.c_str()) 
         .def("get_computation_time", &ChooseSolver::get_computation_time, DocSolver::get_computation_time.c_str())
+        .def("get_timers", &ChooseSolver::get_timers, "TODO")
+        .def("get_timers_jacobian", &ChooseSolver::get_timers_jacobian, "TODO")
         .def("get_fdpf_xb_lu", &ChooseSolver::get_fdpf_xb_lu, py::return_value_policy::reference, DocGridModel::_internal_do_not_use.c_str())  // TODO this for all solver !
         .def("get_fdpf_bx_lu", &ChooseSolver::get_fdpf_bx_lu, py::return_value_policy::reference, DocGridModel::_internal_do_not_use.c_str());
 
@@ -815,6 +868,19 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("get_trafohv_theta", &GridModel::get_trafohv_theta, DocGridModel::_internal_do_not_use.c_str())
         .def("get_trafolv_theta", &GridModel::get_trafolv_theta, DocGridModel::_internal_do_not_use.c_str())
 
+        .def("get_all_shunt_buses", &GridModel::get_all_shunt_buses, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_loads_res_full", &GridModel::get_loads_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_shunts_res_full", &GridModel::get_shunts_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_gen_res_full", &GridModel::get_gen_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_lineor_res_full", &GridModel::get_lineor_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_lineex_res_full", &GridModel::get_lineex_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_trafohv_res_full", &GridModel::get_trafohv_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_trafolv_res_full", &GridModel::get_trafolv_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_storages_res_full", &GridModel::get_storages_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_sgens_res_full", &GridModel::get_sgens_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_dclineor_res_full", &GridModel::get_dclineor_res_full, DocGridModel::_internal_do_not_use.c_str())
+        .def("get_dclineex_res_full", &GridModel::get_dclineex_res_full, DocGridModel::_internal_do_not_use.c_str())
+        
         // do something with the grid
         // .def("init_Ybus", &DataModel::init_Ybus) // temporary
         .def("deactivate_result_computation", &GridModel::deactivate_result_computation, DocGridModel::deactivate_result_computation.c_str())
@@ -913,6 +979,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         // remove some defaults (TODO)
         .def("reset", &ContingencyAnalysis::clear, DocSecurityAnalysis::clear.c_str())
         .def("clear", &ContingencyAnalysis::clear, DocSecurityAnalysis::clear.c_str())
+        .def("clear_results_only", &ContingencyAnalysis::clear_results_only, DocSecurityAnalysis::clear.c_str())
         .def("close", &ContingencyAnalysis::clear, DocComputers::clear.c_str())
         .def("remove_n1", &ContingencyAnalysis::remove_n1, DocSecurityAnalysis::remove_n1.c_str())
         .def("remove_nk", &ContingencyAnalysis::remove_nk, DocSecurityAnalysis::remove_nk.c_str())
diff --git a/src/powerflow_algorithm/BaseAlgo.h b/src/powerflow_algorithm/BaseAlgo.h
index 8556609..fc7ce46 100644
--- a/src/powerflow_algorithm/BaseAlgo.h
+++ b/src/powerflow_algorithm/BaseAlgo.h
@@ -30,6 +30,9 @@
 
 class GridModel;
 
+typedef std::tuple<double, double, double, double, 
+                   double, double, double, double, 
+                   double> TimerJacType;
 
 /**
 This class represents a algorithm to compute powerflow.
@@ -85,6 +88,22 @@ class BaseAlgo : public BaseConstants
               timer_Fx_, timer_solve_, timer_check_, timer_total_nr_);
             return res;
         }
+        
+        virtual TimerJacType  get_timers_jacobian() const
+        {
+            TimerJacType res = {
+                timer_Fx_,
+                timer_solve_,
+                -1.,  // not available for non NR solver, so I put -1
+                timer_check_,
+                -1.,  // not available for non NR solver, so I put -1
+                -1.,  // not available for non NR solver, so I put -1
+                -1.,  // not available for non NR solver, so I put -1
+                -1.,  // not available for non NR solver, so I put -1
+                timer_total_nr_
+            };
+            return res;
+        }
 
         virtual
         bool compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
diff --git a/src/powerflow_algorithm/BaseNRAlgo.h b/src/powerflow_algorithm/BaseNRAlgo.h
index d593092..2eda817 100644
--- a/src/powerflow_algorithm/BaseNRAlgo.h
+++ b/src/powerflow_algorithm/BaseNRAlgo.h
@@ -18,7 +18,14 @@ template<class LinearSolver>
 class BaseNRAlgo : public BaseAlgo
 {
     public:
-        BaseNRAlgo():BaseAlgo(true), need_factorize_(true), timer_initialize_(0.), timer_dSbus_(0.), timer_fillJ_(0.) {}
+        BaseNRAlgo():
+            BaseAlgo(true),
+            need_factorize_(true),
+            timer_initialize_(0.),
+            timer_dSbus_(0.),
+            timer_fillJ_(0.),
+            timer_Va_Vm_(0.),
+            timer_pre_proc_(0.){}
 
         virtual
         Eigen::Ref<const Eigen::SparseMatrix<real_type> > get_J() const {
@@ -32,11 +39,18 @@ class BaseNRAlgo : public BaseAlgo
         }
 
         virtual
-        std::tuple<double, double, double, double, double, double, double> get_timers_jacobian()
+        TimerJacType get_timers_jacobian() const
         {
             // TODO refacto that, and change the order
-            auto res = std::tuple<double, double, double, double, double, double, double>(
-              timer_Fx_, timer_solve_, timer_initialize_, timer_check_, timer_dSbus_, timer_fillJ_, timer_total_nr_);
+            auto res = TimerJacType(timer_Fx_,
+                                    timer_solve_,
+                                    timer_initialize_,
+                                    timer_check_,
+                                    timer_dSbus_,
+                                    timer_fillJ_,
+                                    timer_Va_Vm_,
+                                    timer_pre_proc_,
+                                    timer_total_nr_);
             return res;
         }
 
@@ -59,6 +73,8 @@ class BaseNRAlgo : public BaseAlgo
             BaseAlgo::reset_timer();
             timer_dSbus_ = 0.;
             timer_fillJ_ = 0.;
+            timer_Va_Vm_ = 0.;
+            timer_pre_proc_ = 0.;
             timer_initialize_ = 0.;
         }
         virtual
@@ -171,6 +187,8 @@ class BaseNRAlgo : public BaseAlgo
         double timer_initialize_;
         double timer_dSbus_;
         double timer_fillJ_;
+        double timer_Va_Vm_;
+        double timer_pre_proc_;
 
 
     Eigen::SparseMatrix<real_type>
diff --git a/src/powerflow_algorithm/BaseNRAlgo.tpp b/src/powerflow_algorithm/BaseNRAlgo.tpp
index 33ae20e..e889771 100644
--- a/src/powerflow_algorithm/BaseNRAlgo.tpp
+++ b/src/powerflow_algorithm/BaseNRAlgo.tpp
@@ -13,15 +13,15 @@
 
 template<class LinearSolver>
 bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
-                                            CplxVect & V,
-                                            const CplxVect & Sbus,
-                                            const Eigen::VectorXi & slack_ids,
-                                            const RealVect & slack_weights,
-                                            const Eigen::VectorXi & pv,
-                                            const Eigen::VectorXi & pq,
-                                            int max_iter,
-                                            real_type tol
-                                            )
+                                          CplxVect & V,
+                                          const CplxVect & Sbus,
+                                          const Eigen::VectorXi & slack_ids,
+                                          const RealVect & slack_weights,
+                                          const Eigen::VectorXi & pv,
+                                          const Eigen::VectorXi & pq,
+                                          int max_iter,
+                                          real_type tol
+                                          )
 {
     /**
     This method uses the newton raphson algorithm to compute voltage angles and magnitudes at each bus
@@ -54,6 +54,7 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
     reset_if_needed();
     err_ = ErrorType::NoError;  // reset the error if previous error happened
     auto timer = CustTimer();
+    auto timer_pre_proc = CustTimer();
 
     Eigen::VectorXi my_pv = retrieve_pv_with_slack(slack_ids, pv);  // retrieve_pv_with_slack (not all), add_slack_to_pv (all)
     real_type slack_absorbed = std::real(Sbus.sum());  // initial guess for slack_absorbed
@@ -77,6 +78,7 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
     V_ = V;
     Vm_ = V_.array().abs();  // update Vm and Va again in case
     Va_ = V_.array().arg();  // we "wrapped around" with a negative Vm
+    timer_pre_proc_ += timer_pre_proc.duration();
 
     // first check, if the problem is already solved, i stop there
     // compute a first time the mismatch to initialize the slack bus
@@ -132,20 +134,31 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
         }
         // const auto dx = -F;  // removed for speed optimization (-= used below)
 
+        auto timer_va_vm = CustTimer();
+        slack_absorbed -= F(0); // by convention in fill_jacobian_matrix the slack bus is the first component
         // update voltage (this should be done consistently with "_evaluate_Fx")
         if (n_pv > 0) Va_(my_pv) -= F.segment(1, n_pv);
         if (n_pq > 0){
             Va_(pq) -= F.segment(n_pv + 1, n_pq);
             Vm_(pq) -= F.segment(n_pv + n_pq + 1, n_pq);
         }
-        slack_absorbed -= F(0); // by convention in fill_jacobian_matrix the slack bus is the first component
+        
 
         // std::cout << "iter " << nr_iter_ << " dx(0): " << -F(0) << " dx(1): " << -F(1) << std::endl;
         // std::cout << "slack_absorbed " << slack_absorbed << std::endl;
         // TODO change here for not having to cast all the time ... maybe
         V_ = Vm_.array() * (Va_.array().cos().template cast<cplx_type>() + my_i * Va_.array().sin().template cast<cplx_type>() );
-        Vm_ = V_.array().abs();  // update Vm and Va again in case
-        Va_ = V_.array().arg();  // we wrapped around with a negative Vm TODO more efficient way maybe ?
+        // V_ = Vm_.array() *  (my_i * Va_.array().template cast<cplx_type>()).exp() ;
+        if(Vm_.minCoeff() < 0.)
+        {
+            // update Vm and Va again in case
+            // we wrapped around with a negative Vm 
+            // TODO more efficient way maybe ?
+            Vm_ = V_.array().abs();  
+            Va_ = V_.array().arg();  
+            // TODO do I need to change V here ???
+        }
+        timer_Va_Vm_ += timer_va_vm.duration();
 
         F = _evaluate_Fx(Ybus, V_, Sbus, slack_bus_id, slack_absorbed, slack_weights, my_pv, pq);
         bool tmp = F.allFinite();
@@ -162,7 +175,7 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
     timer_total_nr_ += timer.duration();
     #ifdef __COUT_TIMES
         std::cout << "Computation time: " << "\n\t timer_initialize_: " << timer_initialize_
-                  << "\n\t timer_dSbus_ (called in _fillJ_): " << timer_dSbus_
+                  << "\n\t timer_dSbus_: " << timer_dSbus_
                   << "\n\t timer_fillJ_: " << timer_fillJ_
                   << "\n\t timer_Fx_: " << timer_Fx_
                   << "\n\t timer_check_: " << timer_check_
@@ -170,6 +183,8 @@ bool BaseNRAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<cplx_type> &
                   << "\n\t timer_total_nr_: " << timer_total_nr_
                   << "\n\n";
     #endif // __COUT_TIMES
+    Vm_ = V_.array().abs();  
+    Va_ = V_.array().arg();  
     _solver_control.tell_none_changed();
     return res;
 }
@@ -345,10 +360,9 @@ void BaseNRAlgo<LinearSolver>::fill_jacobian_matrix(const Eigen::SparseMatrix<cp
     `slack` is the representation of the equation connecting together the slack buses (represented by slack_weights)
     the remaining pq components are all 0.
     **/
-
-    auto timer = CustTimer();
     _dSbus_dV(Ybus, V);
 
+    auto timer = CustTimer();
     const auto n_pvpq = pvpq.size();
     const auto n_pq = pq.size();
     const auto size_j = n_pvpq + n_pq + 1;   // +1 because i add the slack bus
diff --git a/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp b/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp
index 612f963..efd86b5 100644
--- a/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp
+++ b/src/powerflow_algorithm/BaseNRSingleSlackAlgo.tpp
@@ -48,6 +48,7 @@ bool BaseNRSingleSlackAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<c
     BaseNRAlgo<LinearSolver>::err_ = ErrorType::NoError;  // reset the error if previous error happened
     
     auto timer = CustTimer();
+    auto timer_pre_proc = CustTimer();
     // initialize once and for all the "inverse" of these vectors
     // Eigen::VectorXi my_pv = BaseNRAlgo<LinearSolver>::retrieve_pv_with_slack(slack_ids, pv);
     Eigen::VectorXi my_pv = pv;
@@ -66,6 +67,7 @@ bool BaseNRSingleSlackAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<c
     BaseNRAlgo<LinearSolver>::V_ = V;
     BaseNRAlgo<LinearSolver>::Vm_ = BaseNRAlgo<LinearSolver>::V_.array().abs();  // update Vm and Va again in case
     BaseNRAlgo<LinearSolver>::Va_ = BaseNRAlgo<LinearSolver>::V_.array().arg();  // we wrapped around with a negative Vm
+    BaseNRAlgo<LinearSolver>::timer_pre_proc_ += timer_pre_proc.duration();
 
     // first check, if the problem is already solved, i stop there
     RealVect F = BaseNRAlgo<LinearSolver>::_evaluate_Fx(Ybus, V, Sbus, my_pv, pq);
@@ -124,7 +126,7 @@ bool BaseNRSingleSlackAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<c
             break;
         }
         // auto dx = -F;
-
+        auto timer_va_vm = CustTimer();
         BaseNRAlgo<LinearSolver>::Vm_ = BaseNRAlgo<LinearSolver>::V_.array().abs();  // update Vm and Va again in case
         BaseNRAlgo<LinearSolver>::Va_ = BaseNRAlgo<LinearSolver>::V_.array().arg();  // we wrapped around with a negative Vm
 
@@ -139,6 +141,14 @@ bool BaseNRSingleSlackAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<c
         const RealVect & Vm = BaseNRAlgo<LinearSolver>::Vm_;  // I am forced to redefine the type for it to compile properly
         const RealVect & Va = BaseNRAlgo<LinearSolver>::Va_;
         BaseNRAlgo<LinearSolver>::V_ = Vm.array() * (Va.array().cos().cast<cplx_type>() + m_i * Va.array().sin().cast<cplx_type>() );
+        if(BaseNRAlgo<LinearSolver>::Vm_.minCoeff() < 0.)
+        {
+            // update Vm and Va again in case
+            // we wrapped around with a negative Vm TODO more efficient way maybe ?
+            BaseNRAlgo<LinearSolver>::Vm_ = BaseNRAlgo<LinearSolver>::V_.array().abs();  
+            BaseNRAlgo<LinearSolver>::Va_ = BaseNRAlgo<LinearSolver>::V_.array().arg();  
+        }
+        BaseNRAlgo<LinearSolver>::timer_Va_Vm_ += timer_va_vm.duration();
 
         F = BaseNRAlgo<LinearSolver>::_evaluate_Fx(Ybus, BaseNRAlgo<LinearSolver>::V_, Sbus, my_pv, pq);
         bool tmp = F.allFinite();
@@ -164,6 +174,10 @@ bool BaseNRSingleSlackAlgo<LinearSolver>::compute_pf(const Eigen::SparseMatrix<c
                   << "\n\t timer_total_nr_: " << BaseNRAlgo<LinearSolver>::timer_total_nr_
                   << "\n\n";
     #endif // __COUT_TIMES
+    // update Vm and Va again in case
+    // we wrapped around with a negative Vm TODO more efficient way maybe ?
+    BaseNRAlgo<LinearSolver>::Vm_ = BaseNRAlgo<LinearSolver>::V_.array().abs();  
+    BaseNRAlgo<LinearSolver>::Va_ = BaseNRAlgo<LinearSolver>::V_.array().arg();  
     BaseNRAlgo<LinearSolver>::_solver_control.tell_none_changed();
     return res;
 }
@@ -189,9 +203,9 @@ void BaseNRSingleSlackAlgo<LinearSolver>::fill_jacobian_matrix(const Eigen::Spar
     J22 = dS_dVm[array([pq]).T, pq].imag
     **/
 
-    auto timer = CustTimer();
     BaseNRAlgo<LinearSolver>::_dSbus_dV(Ybus, V);
 
+    auto timer = CustTimer();
     const int n_pvpq = static_cast<int>(pvpq.size());
     const int n_pq = static_cast<int>(pq.size());
     const int size_j = n_pvpq + n_pq;


Test patch:
diff --git a/lightsim2grid/tests/test_DCPF.py b/lightsim2grid/tests/test_DCPF.py
index 3e910d7..6faf256 100644
--- a/lightsim2grid/tests/test_DCPF.py
+++ b/lightsim2grid/tests/test_DCPF.py
@@ -196,8 +196,9 @@ def _aux_test(self, pn_net):
         load_p, load_q, load_v = backend.loads_info()
         max_mis = np.max(np.abs(load_p - load_p_pp))
         assert max_mis <= self.tol, f"Error: load_p do not match, maximum absolute error is {max_mis:.5f} MW"
-        max_mis = np.max(np.abs(load_q - load_q_pp))
-        assert max_mis <= self.tol, f"Error: load_q do not match, maximum absolute error is {max_mis:.5f} MVAr"
+        # PP does not set "load_q" to 0. in DC
+        # max_mis = np.max(np.abs(load_q - load_q_pp))
+        # assert max_mis <= self.tol, f"Error: load_q do not match, maximum absolute error is {max_mis:.5f} MVAr"
         max_mis = np.max(np.abs(load_v - load_v_pp))
         assert max_mis <= self.tol, f"Error: load_v do not match, maximum absolute error is {max_mis:.5f} kV"
 
diff --git a/lightsim2grid/tests/test_compat_legacy_grid2op.py b/lightsim2grid/tests/test_compat_legacy_grid2op.py
new file mode 100644
index 0000000..4688610
--- /dev/null
+++ b/lightsim2grid/tests/test_compat_legacy_grid2op.py
@@ -0,0 +1,188 @@
+# Copyright (c) 2019-2023, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of Grid2Op, Grid2Op a testbed platform to model sequential decision making in power systems.
+
+import warnings
+import unittest
+import numpy as np
+
+try:
+    import grid2op
+    from grid2op.Action import PlayableAction
+    from grid2op.Environment import Environment
+    from grid2op.Runner import Runner
+    GLOP_AVAIL = True
+except ImportError:
+    GLOP_AVAIL = False
+
+try:
+    from grid2op.gym_compat import GymEnv, BoxGymActSpace, BoxGymObsSpace, DiscreteActSpace, MultiDiscreteActSpace
+    GYM_AVAIL = True
+except ImportError as exc_:
+    print(exc_)
+    GYM_AVAIL = False
+    
+from lightsim2grid import LightSimBackend
+            
+            
+class TestEnvironmentBasic(unittest.TestCase):          
+    def setUp(self) -> None:
+        if not GLOP_AVAIL:
+            self.skipTest("grid2op not installed")
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            self.env = grid2op.make("educ_case14_storage",
+                                    test=True,
+                                    action_class=PlayableAction,
+                                    _add_to_name=type(self).__name__,
+                                    backend=LightSimBackend())
+        self.line_id = 3
+        th_lim = self.env.get_thermal_limit() * 2.  # avoid all problem in general
+        th_lim[self.line_id] /= 10.  # make sure to get trouble in line 3
+        self.env.set_thermal_limit(th_lim)
+        
+        TestEnvironmentBasic._init_env(self.env)
+        
+    @staticmethod  
+    def _init_env(env):
+        env.set_id(0)
+        env.seed(0)
+        env.reset()
+        
+    def tearDown(self) -> None:
+        self.env.close()
+        return super().tearDown()
+    
+    def test_right_type(self):
+        assert isinstance(self.env, Environment)
+        
+    def test_ok(self):
+        act = self.env.action_space()
+        for i in range(10):
+            obs_in, reward, done, info = self.env.step(act)
+            if i < 3:
+                assert obs_in.timestep_overflow[self.line_id] == i + 1, f"error for step {i}: {obs_in.timestep_overflow[self.line_id]}"
+            else:
+                assert not obs_in.line_status[self.line_id]
+    
+    def test_reset(self):
+        # timestep_overflow should be 0 initially even if the flow is too high
+        obs = self.env.reset()
+        assert obs.timestep_overflow[self.line_id] == 0
+        assert obs.rho[self.line_id] > 1.
+        
+
+class TestEnvironmentBasicCpy(TestEnvironmentBasic):
+    def setUp(self) -> None:
+        super().setUp()
+        init_int = self.env
+        self.env = self.env.copy()
+        init_int.close()
+        
+
+class TestBasicEnvironmentRunner(unittest.TestCase):    
+    def setUp(self) -> None:
+        TestEnvironmentBasic.setUp(self)
+        self.max_iter = 10
+
+    def tearDown(self) -> None:
+        self.env.close()
+        return super().tearDown()
+        
+    def test_runner_can_make(self):
+        runner = Runner(**self.env.get_params_for_runner())
+        env2 = runner.init_env()
+        assert isinstance(env2, Environment)
+    
+    def test_runner(self):
+        # create the runner
+        runner_in = Runner(**self.env.get_params_for_runner())
+        res_in, *_ = runner_in.run(nb_episode=1, max_iter=self.max_iter, env_seeds=[0], episode_id=[0], add_detailed_output=True)
+        res_in2, *_ = runner_in.run(nb_episode=1, max_iter=self.max_iter, env_seeds=[0], episode_id=[0])
+        # check correct results are obtained when agregated
+        assert res_in[3] == 10
+        assert res_in2[3] == 10
+        ref_val = 164.1740722
+        assert np.allclose(res_in[2], ref_val), f"{res_in[2]} vs {ref_val}"
+        assert np.allclose(res_in2[2], ref_val), f"{res_in2[2]} vs {ref_val}"
+        
+        # check detailed results
+        ep_data_in = res_in[-1]
+        for i in range(1, self.max_iter + 1):
+            # there is a bug in grid2op 1.6.4 for i = 0...
+            obs_in = ep_data_in.observations[i]
+            assert obs_in is not None, f"error for step {i}"
+            if i < 4:
+                assert obs_in.timestep_overflow[self.line_id] == i, f"error for step {i}: {obs_in.timestep_overflow[self.line_id]}"
+            else:
+                assert not obs_in.line_status[self.line_id], f"error for step {i}: line is not disconnected"
+    
+        
+class TestBasicEnvironmentGym(unittest.TestCase):
+    def setUp(self) -> None:
+        if not GYM_AVAIL:
+            self.skipTest("gym is not installed (gymnasium did not exist in 1.6.4)")
+        TestEnvironmentBasic.setUp(self)
+
+    def tearDown(self) -> None:
+        self.env.close()
+        return super().tearDown()
+    
+    def _aux_run_envs(self, act, env_gym):
+        for i in range(10):
+            obs_in, reward, done, info = env_gym.step(act)
+            if i < 3:
+                assert obs_in["timestep_overflow"][self.line_id] == i + 1, f"error for step {i}: {obs_in['timestep_overflow'][self.line_id]}"
+            else:
+                assert not obs_in["line_status"][self.line_id]
+    
+    def test_gym_with_step(self):
+        """test the step function also disconnects (or not) the lines"""
+        env_gym = GymEnv(self.env)
+        act = {}
+        self._aux_run_envs(act, env_gym)
+        env_gym.reset()
+        self._aux_run_envs(act, env_gym)
+            
+    def test_gym_normal(self):
+        """test I can create the gym env"""
+        env_gym = GymEnv(self.env)
+        env_gym.reset()
+    
+    def test_gym_box(self):
+        """test I can create the gym env with box ob space and act space"""
+        env_gym = GymEnv(self.env)
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            env_gym.action_space = BoxGymActSpace(self.env.action_space)
+            env_gym.observation_space = BoxGymObsSpace(self.env.observation_space)
+        env_gym.reset()
+    
+    def test_gym_discrete(self):
+        """test I can create the gym env with discrete act space"""
+        env_gym = GymEnv(self.env)
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            env_gym.action_space = DiscreteActSpace(self.env.action_space)
+        env_gym.reset()
+        act = 0
+        self._aux_run_envs(act, env_gym)
+        
+    def test_gym_multidiscrete(self):
+        """test I can create the gym env with multi discrete act space"""
+        env_gym = GymEnv(self.env)
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            env_gym.action_space = MultiDiscreteActSpace(self.env.action_space)
+        env_gym.reset()
+        act = env_gym.action_space.sample()
+        act[:] = 0
+        self._aux_run_envs(act, env_gym)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/lightsim2grid/tests/test_n1contingencyrewards.py b/lightsim2grid/tests/test_n1contingencyrewards.py
new file mode 100644
index 0000000..11b4479
--- /dev/null
+++ b/lightsim2grid/tests/test_n1contingencyrewards.py
@@ -0,0 +1,202 @@
+# Copyright (c) 2024, RTE (https://www.rte-france.com)
+# See AUTHORS.txt
+# This Source Code Form is subject to the terms of the Mozilla Public License, version 2.0.
+# If a copy of the Mozilla Public License, version 2.0 was not distributed with this file,
+# you can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+# This file is part of LightSim2grid, LightSim2grid a implements a c++ backend targeting the Grid2Op platform.
+import unittest
+import warnings
+import copy
+import numpy as np
+
+import grid2op
+from grid2op.Backend import PandaPowerBackend
+from grid2op.Action import CompleteAction
+from grid2op.Reward import EpisodeDurationReward
+
+from lightsim2grid import LightSimBackend
+from lightsim2grid.rewards import N1ContingencyReward
+
+
+TH_LIM_A_REF = np.array([
+        541.0,
+        450.0,
+        375.0,
+        636.0,
+        175.0,
+        285.0,
+        335.0,
+        657.0,
+        496.0,
+        827.0,
+        442.0,
+        641.0,
+        840.0,
+        156.0,
+        664.0,
+        235.0,
+        119.0,
+        179.0,
+        1986.0,
+        1572.0,
+    ])
+ 
+ 
+class TestN1ContingencyReward_Base(unittest.TestCase):
+    def get_env_nm(self):
+        return "educ_case14_storage"
+    
+    def init_backend(self):
+        return LightSimBackend()
+    
+    def is_dc(self):
+        return False
+    
+    def threshold_margin(self):
+        return 1.
+    
+    def l_ids(self):
+        return None
+    
+    def setUp(self) -> None:
+        reward = N1ContingencyReward(dc=self.is_dc(),
+                                     threshold_margin=self.threshold_margin(),
+                                     l_ids=self.l_ids())
+        
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            self.env = grid2op.make(self.get_env_nm(), 
+                                    test=True,
+                                    backend=self.init_backend(),
+                                    reward_class=reward,
+                                    action_class=CompleteAction,
+                                    _add_to_name=type(self).__name__)
+        params = self.env.parameters
+        params.MAX_LINE_STATUS_CHANGED = 999999
+        params.MAX_SUB_CHANGED = 999999
+        params.NB_TIMESTEP_COOLDOWN_LINE = 0
+        params.NB_TIMESTEP_COOLDOWN_SUB = 0
+        self.env.change_parameters(params)
+        if self.is_dc():
+            params = copy.deepcopy(params)
+            params.ENV_DC = True
+            params.FORECAST_DC = True
+        self.env.change_forecast_parameters(params)
+        assert (self.env.get_thermal_limit() == TH_LIM_A_REF).all()
+        self.my_ids = self.l_ids()
+        if self.my_ids is None:
+            self.my_ids = list(range(type(self.env).n_line))
+        self.env.observation_space.change_reward(EpisodeDurationReward)
+        return super().setUp()
+    
+    def tearDown(self) -> None:
+        self.env.close()
+        return super().tearDown()
+    
+    def _aux_test_reward(self, obs, reward):
+        unsafe_cont = 0
+        if self.is_dc():
+            th_lim = obs._thermal_limit
+            # A = sqrt(p**2 + q**2) / (sqrt3) * v)
+            # A**2 = (p**2 + q**2) / (3. v**2)
+            # p**2 = A**2 * 3. * v**2 - q**2
+            # and now the units...
+            # W**2 = A * 3. * V**2 - VAr**2
+            # MW**2 = 1e12 * A**2 * 3. * v**2 - 1e12 * VAr**2
+            # MW**2 = 1e6*A**2 * 3. * 1e6*(V)**2 *  - MVAr**2
+            # MW**2 = kA**2 * 3. * kV**2 - MVAr**2
+            p_square = 3. * (1e-3*th_lim)**2 * (obs.v_or)**2 - (obs.q_or)**2
+            p_square[p_square <= 0.] = 0.
+            th_lim_p = np.sqrt(p_square) * self.threshold_margin()
+            
+        # print("test:")
+        for l_id in self.my_ids:
+            sim_obs, sim_r, sim_d, sim_i = obs.simulate(self.env.action_space({"set_line_status": [(l_id, -1)]}),
+                                                        time_step=0)
+            if not self.is_dc():
+                if np.any(sim_obs.a_or > obs._thermal_limit * self.threshold_margin()) or sim_d:
+                    unsafe_cont += 1       
+            else:
+                if np.any(np.abs(sim_obs.p_or) > th_lim_p) or sim_d:
+                    unsafe_cont += 1    
+                
+        assert reward == (len(self.my_ids) - unsafe_cont), f"{reward} vs {(len(self.my_ids) - unsafe_cont)}"
+        
+    def test_do_nothing(self):
+        obs = self.env.reset()
+        obs, reward, done, info = self.env.step(self.env.action_space())
+        self._aux_test_reward(obs, reward)
+        obs, reward, done, info = self.env.step(self.env.action_space())
+        self._aux_test_reward(obs, reward)
+        
+    def test_disconnected_line(self):
+        obs = self.env.reset()
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_line_status": [(0, -1)]}))
+        assert not done
+        self._aux_test_reward(obs, reward)
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_line_status": [(0, +1)]}))
+        self._aux_test_reward(obs, reward)
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_line_status": [(4, -1)]}))
+        assert not done
+        self._aux_test_reward(obs, reward)
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_line_status": [(4, +1)]}))
+        self._aux_test_reward(obs, reward)
+        
+    def test_modif_topo(self):
+        obs = self.env.reset()
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_bus": {"substations_id": [(1, (1, 2, 1, 2, 1, 2))]}}))
+        assert not done
+        self._aux_test_reward(obs, reward)
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_bus": {"substations_id": [(3, (1, 2, 1, 2, 1, 2))]}}))
+        assert not done
+        self._aux_test_reward(obs, reward)
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_bus": {"substations_id": [(1, (1, 1, 1, 1, 1, 1))]}}))
+        assert not done
+        self._aux_test_reward(obs, reward)
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_bus": {"substations_id": [(3, (1, 1, 1, 1, 1, 1))]}}))
+        assert not done
+        self._aux_test_reward(obs, reward)
+    
+    def test_copy(self):
+        obs = self.env.reset()
+        env_cpy = self.env.copy()
+        obs_cpy = env_cpy.copy()
+        assert self.env._reward_helper.template_reward._backend is not env_cpy._reward_helper.template_reward._backend
+        obs, reward, done, info = self.env.step(self.env.action_space(
+            {"set_bus": {"substations_id": [(3, (1, 2, 1, 2, 1, 2))]}}))
+        obs_cpy, reward_cpy, done_cpy, info_cpy = env_cpy.step(self.env.action_space(
+            {"set_line_status": [(0, -1)]}))
+        self._aux_test_reward(obs_cpy, reward_cpy)
+        self._aux_test_reward(obs, reward)
+        assert reward != reward_cpy
+
+
+class TestN1ContingencyReward_DC(TestN1ContingencyReward_Base):
+    def is_dc(self):
+        return True
+    # def l_ids(self):
+    #     return [18]
+    
+    
+class TestN1ContingencyReward_LIDS(TestN1ContingencyReward_Base):
+    def l_ids(self):
+        return [0, 1, 2, 18, 16]
+    
+    
+class TestN1ContingencyReward_Margins(TestN1ContingencyReward_Base):
+    def threshold_margin(self):
+        return 0.9
+    
+    
+class TestN1ContingencyReward_PP(TestN1ContingencyReward_Base):
+    def init_backend(self):
+        return PandaPowerBackend(with_numba=False, lightsim2grid=False)
diff --git a/lightsim2grid/tests/test_n_busbar_per_sub.py b/lightsim2grid/tests/test_n_busbar_per_sub.py
index aa51bb0..e67c5a2 100644
--- a/lightsim2grid/tests/test_n_busbar_per_sub.py
+++ b/lightsim2grid/tests/test_n_busbar_per_sub.py
@@ -14,6 +14,7 @@
 
 from lightsim2grid import LightSimBackend
  
+ 
 class TestLightSimBackend_3busbars(unittest.TestCase):
     def get_nb_bus(self):
         return 3
@@ -147,7 +148,7 @@ def test_move_gen(self):
             self.env.backend.apply_action(bk_act)
             global_bus = sub_id + (new_bus -1) * cls.n_sub 
             if new_bus >= 1:
-                assert self.env.backend._grid.get_generators()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_loads()[el_id].bus_id} vs {global_bus}"
+                assert self.env.backend._grid.get_generators()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_generators()[el_id].bus_id} vs {global_bus}"
                 if line_or_id is not None:
                     assert self.env.backend._grid.get_lines()[line_or_id].bus_or_id == global_bus
                 else:
@@ -179,7 +180,7 @@ def test_move_storage(self):
             self.env.backend.apply_action(bk_act)
             global_bus = sub_id + (new_bus -1) * cls.n_sub 
             if new_bus >= 1:
-                assert self.env.backend._grid.get_storages()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_loads()[el_id].bus_id} vs {global_bus}"
+                assert self.env.backend._grid.get_storages()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_sotrages()[el_id].bus_id} vs {global_bus}"
                 if line_or_id is not None:
                     assert self.env.backend._grid.get_lines()[line_or_id].bus_or_id == global_bus
                 else:
@@ -243,9 +244,12 @@ def test_move_shunt(self):
             bk_act = self.env._backend_action_class()
             bk_act += act
             self.env.backend.apply_action(bk_act)
+            self.env.backend._set_shunt_info()
+            sh_p, sh_q, sh_v, sh_bus = self.env.backend.shunt_info()
+            assert sh_bus[el_id] == new_bus
             global_bus = sub_id + (new_bus -1) * cls.n_sub 
             if new_bus >= 1:
-                assert self.env.backend._grid.get_shunts()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_loads()[el_id].bus_id} vs {global_bus}"
+                assert self.env.backend._grid.get_shunts()[el_id].bus_id == global_bus, f"error for new_bus {new_bus}: {self.env.backend._grid.get_shunts()[el_id].bus_id} vs {global_bus}"
                 if line_or_id is not None:
                     assert self.env.backend._grid.get_lines()[line_or_id].bus_or_id == global_bus
                 else: