Repository: Grid2op/lightsim2grid (language: python)

Issue:
The init voltage values could depends on the previous state even when "init_dc" is True
This was because the "init_dc" part did not reset "self.V" in the step function.

The major impact is that the backend was "stuck" in this state until it could converged again. And it might not be possible to make it converged this way. 
It entails major difference with pandapower, especially using the "obs.simulate" function.

This bug has been found while running this code


Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 3eba549..bcb6837 100644
--- a/.gitignore
+++ b/.gitignore
@@ -210,3 +210,6 @@ build_cmake/SuiteSparse/KLU/Lib
 build_cmake/SuiteSparse/KLU/Source
 build_cmake/SuiteSparse/SuiteSparse_config/SuiteSparse_config.c
 build_cmake/SuiteSparse/SuiteSparse_config/SuiteSparse_config.h
+lightsim2grid/tests/pn.case6470rte.json
+**/nicslu.lic
+benchmarks/debug_speed_improvments.py
\ No newline at end of file
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9cc1702..cfdfed2 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -9,16 +9,26 @@ Change Log
 - easier building (get rid of the "make" part)
 - code helm powerflow method
 - interface with gridpack (to enforce q limits for example)
-- have a solver with "NICSLU" (https://github.com/chenxm1986/nicslu) if available on the platform
 - maybe have a look at suitesparse "sliplu" tools ?
 
+[0.5.4] 2021-08-20
+------------------
+- [FIXED] a bug for static generator (wrong signed convention were used in some part of the c++ code). This has
+  no impact at all for provided grid2op environments.
+- [FIXED] An issue where the backend could get "stuck" in a wrong state because of the way the Vinit was computed (see
+  `Issue 30 <https://github.com/BDonnot/lightsim2grid/issues/30>`_)
+- [ADDED] experimental support for the `NICSLU` linear solver (requires a proper license and library, see
+  https://github.com/chenxm1986/nicslu for more information. Support does not include multi threaded at the moment.
+- [IMPROVED] minor performance improvements for the solvers based on Newton Raphson (faster filling of the Jacobian
+  matrix after the first iteration)
+
 [0.5.3] 2021-08-11
 -------------------
 - [FIXED] minor issues in the benchmark (some time measurments were wrong)
 - [ADDED] lightsim2grid package now can be distributed on pypi
 - [ADDED] compilation of SuiteSparse using cmake
 - [ADDED] compatibility with the KLU linear solver on windows based systems.
-- [IMPROVED] trying to make the package available on pypi
+- [IMPROVED] the package should now be available on pypi
 
 [0.5.2] 2021-07-26
 -------------------
@@ -72,7 +82,7 @@ Change Log
   using the SuiteSparse KLU linear solver but rather the Eigen SparseLU one)
 - [ADDED] start of the documentation.
 
-[0.2.4] - 2020-08-xx
+[0.2.4] - 2020-08-20
 --------------------
 - [FIXED] issue for copying environment
 
diff --git a/README.md b/README.md
index c8fe12d..ecc8068 100644
--- a/README.md
+++ b/README.md
@@ -64,6 +64,7 @@ get some benefit (in terms of performances) to install it from your on your mach
 You need to:
 - clone this repository and get the code of Eigen (mandatory for compilation) and SparseSuite (optional)
 - (optional) compile a piece of SparseSuite
+- (optional) [experimental] retrieve and get a proper license for the NICSLU linear solver (see https://github.com/chenxm1986/nicslu)
 - install the package
 
 ### Important note
@@ -141,6 +142,28 @@ The main steps (for windows, somme commands needs to be adapted on linux / macos
 
 For more information, feel free to read the dedicated [README](build_cmake/README.md).
 
+### (optional) Include NICSLU linear solver (experimental)
+Another linear solver that can be used with lighsim2grid is the "NICSLU" linear solver that might, in some cases, be
+even faster than the KLU linear solver. This can lead to more speed up if using lighsim2grid.
+
+To use it, you need to:
+
+1) retrieve the sources (only available as a freeware) from https://github.com/chenxm1986/nicslu and save
+   it on your machine. Say you clone this github repository in `NICSLU_GIT` 
+   (*eg* NICSLU_GIT="/home/user/Documents/nicslu/"). Also note that you need to check that your usage
+   is compliant with their license !
+2) define the "PATH_NICSLU" environment variable **before** compiling lightsim2grid, on linux you can do
+   `export PATH_NICSLU=NICSLU_GIT/nicsluDATE` 
+   (for example `export PATH_NICSLU=/home/user/Documents/nicslu/nicslu202103` if you cloned the repository 
+   as the example of `step 1)` and use the version of nicslu compiled by the author on March 2021 [version 
+   distributed at time of writing the readme] )
+
+And this is it. Lightsim will be able to use this linear solver.
+
+Be carefull though, you require a license file in order to use it. As of now, the best way is to copy paste the license
+file at the same location that the one you execute python from (*ie* you need to copy paste it each time). We will 
+try to find another solution.
+
 ### 2. Installation of the python package
 Now you simply need to install the lightsim2grid package this way, like any python package:
 
diff --git a/SuiteSparse b/SuiteSparse
index 1869379..538273c 160000
--- a/SuiteSparse
+++ b/SuiteSparse
@@ -1,1 +1,1 @@
-Subproject commit 1869379f464f0f8dac471edb4e6d010b2b0e639d
+Subproject commit 538273cfd53720a10e34a3d80d3779b607e1ac26
diff --git a/benchmarks/benchmark_solvers.py b/benchmarks/benchmark_solvers.py
index 5441a8e..6fa24d9 100644
--- a/benchmarks/benchmark_solvers.py
+++ b/benchmarks/benchmark_solvers.py
@@ -23,7 +23,7 @@
 from grid2op.Parameters import Parameters
 import lightsim2grid
 from lightsim2grid.LightSimBackend import LightSimBackend
-from utils_benchmark import print_res, run_env, str2bool, get_env_name_displayed
+from utils_benchmark import print_res, run_env, str2bool, get_env_name_displayed, print_configuration
 TABULATE_AVAIL = False
 try:
     from tabulate import tabulate
@@ -34,6 +34,8 @@
 MAX_TS = 1000
 ENV_NAME = "rte_case14_realistic"
 
+NICSLU_LICENSE_AVAIL = os.path.exists("./nicslu.lic") and os.path.isfile("./nicslu.lic")
+
 
 def main(max_ts, env_name_input, test=True,
          no_gs=False, no_gs_synch=False):
@@ -75,6 +77,17 @@ def main(max_ts, env_name_input, test=True,
         klu_comp_time = env_lightsim.backend.comp_time
         klu_time_pf = env_lightsim._time_powerflow
 
+    if lightsim2grid.SolverType.NICSLU in solver_types and NICSLU_LICENSE_AVAIL:
+        env_lightsim.backend.set_solver_type(lightsim2grid.SolverType.NICSLU)
+        env_lightsim.backend.set_solver_max_iter(10)
+        nb_ts_nicslu, time_nicslu, aor_nicslu, gen_p_nicslu, gen_q_nicslu = run_env(env_lightsim,
+                                                                                    max_ts,
+                                                                                    agent, chron_id=0,
+                                                                                    with_type_solver=wst,
+                                                                                    env_seed=0)
+        nicslu_comp_time = env_lightsim.backend.comp_time
+        nicslu_time_pf = env_lightsim._time_powerflow
+
     if lightsim2grid.SolverType.SparseLU in solver_types:
         env_lightsim.backend.set_solver_type(lightsim2grid.SolverType.SparseLU)
         env_lightsim.backend.set_solver_max_iter(10)
@@ -100,8 +113,11 @@ def main(max_ts, env_name_input, test=True,
         gsa_time_pf = env_lightsim._time_powerflow
 
     # NOW PRINT THE RESULTS
+    print("Configuration:")
+    print_configuration()
+
     env_name = get_env_name_displayed(env_name_input)
-    hds = [f"{env_name}", f"grid2op speed (it/s)", f"grid2op powerflow time (ms)", f"solver powerflow time (ms)"]
+    hds = [f"{env_name}", f"grid2op speed (it/s)", f"grid2op 'backend.runpf' time (ms)", f"solver powerflow time (ms)"]
     tab = [["PP", f"{nb_ts_pp/time_pp:.2e}",
             f"{1000.*pp_time_pf/nb_ts_pp:.2e}",
             f"{1000.*pp_comp_time/nb_ts_pp:.2e}"]]
@@ -121,6 +137,10 @@ def main(max_ts, env_name_input, test=True,
         tab.append(["LS+KLU", f"{nb_ts_klu/time_klu:.2e}",
                     f"{1000.*klu_time_pf/nb_ts_klu:.2e}",
                     f"{1000.*klu_comp_time/nb_ts_klu:.2e}"])
+    if lightsim2grid.SolverType.NICSLU in solver_types:
+        tab.append(["LS+NICSLU", f"{nb_ts_nicslu/time_nicslu:.2e}",
+                    f"{1000.*nicslu_time_pf/nb_ts_nicslu:.2e}",
+                    f"{1000.*nicslu_comp_time/nb_ts_nicslu:.2e}"])
 
     if TABULATE_AVAIL:
         res_use_with_grid2op_1 = tabulate(tab, headers=hds,  tablefmt="rst")
@@ -137,7 +157,7 @@ def main(max_ts, env_name_input, test=True,
     print()
 
     hds = [f"{env_name} ({nb_ts_pp} iter)", f"Δ aor (amps)", f"Δ gen_p (MW)", f"Δ gen_q (MVAr)"]
-    tab = [["PP", "0.00", "0.00", "0.00"]]
+    tab = [["PP (ref)", "0.00", "0.00", "0.00"]]
     if lightsim2grid.SolverType.GaussSeidel in solver_types and no_gs is False:
         tab.append(["LS+GS",
                     f"{np.max(np.abs(aor_gs - aor_pp)):.2e}",
@@ -158,6 +178,11 @@ def main(max_ts, env_name_input, test=True,
                     f"{np.max(np.abs(aor_klu - aor_pp)):.2e}",
                     f"{np.max(np.abs(gen_p_klu - gen_p_pp)):.2e}",
                     f"{np.max(np.abs(gen_q_klu - gen_q_pp)):.2e}"])
+    if lightsim2grid.SolverType.NICSLU in solver_types:
+        tab.append(["LS+NICSLU",
+                    f"{np.max(np.abs(aor_nicslu - aor_pp)):.2e}",
+                    f"{np.max(np.abs(gen_p_nicslu - gen_p_pp)):.2e}",
+                    f"{np.max(np.abs(gen_q_nicslu - gen_q_pp)):.2e}"])
 
     if TABULATE_AVAIL:
         res_use_with_grid2op_2 = tabulate(tab, headers=hds,  tablefmt="rst")
diff --git a/benchmarks/debug_speed_improvments.py b/benchmarks/debug_speed_improvments.py
new file mode 100644
index 0000000..d2212c4
--- /dev/null
+++ b/benchmarks/debug_speed_improvments.py
@@ -0,0 +1,4 @@
+import grid2op
+from lightsim2grid import LightSimBackend
+
+env = grid2op.make("l2rpn_case14_sandbox", backend=LightSimBackend())
diff --git a/benchmarks/utils_benchmark.py b/benchmarks/utils_benchmark.py
index 5de716c..8bd191d 100644
--- a/benchmarks/utils_benchmark.py
+++ b/benchmarks/utils_benchmark.py
@@ -49,8 +49,14 @@ def get_rest(env_pp, env_KLU, env_SLU, env_GS):
 
 
 def run_env(env, max_ts, agent, chron_id=None, keep_forecast=False, with_type_solver=False, env_seed=None):
+    pbar_desc = None
     if with_type_solver:
-        print(f"begin: {env.backend._grid.get_solver_type()}")
+        try:
+            pbar_desc = f"{env.backend._grid.get_solver_type()}".split(".")[1]
+        except Exception:
+            # it just means I will not be able to print the fancy name on the progress bar...
+            pass
+
     nb_rows = min(env.chronics_handler.max_timestep(), max_ts)
     if nb_rows == -1:
         # -1 indicated "infinite data"
@@ -84,15 +90,12 @@ def run_env(env, max_ts, agent, chron_id=None, keep_forecast=False, with_type_so
     # it's not 0. because a powerflow is run when the backend is reset, and this is one more powerflow than the
     # number of steps
 
-    if with_type_solver:
-        print(f"after reset: {env.backend._grid.get_solver_type()}")
-
     done = False
     reward = env.reward_range[0]
     nb_ts = 0
     prev_act = None
     beg_ = time.time()
-    with tqdm(total=nb_rows) as pbar:
+    with tqdm(total=nb_rows, desc=pbar_desc) as pbar:
         while not done:
             act = agent.act(obs, reward, done)
             obs, reward, done, info = env.step(act)
@@ -108,8 +111,6 @@ def run_env(env, max_ts, agent, chron_id=None, keep_forecast=False, with_type_so
             prev_act = act
     end_ = time.time()
     total_time = end_ - beg_
-    if with_type_solver:
-        print(f"end: {env.backend._grid.get_solver_type()}")
     return nb_ts, total_time, aor, gen_p, gen_q
 
 
@@ -121,4 +122,37 @@ def str2bool(v):
     elif v.lower() in ('no', 'false', 'f', 'n', '0'):
         return False
     else:
-        raise argparse.ArgumentTypeError('Boolean value expected.')
\ No newline at end of file
+        raise argparse.ArgumentTypeError('Boolean value expected.')
+
+
+def print_configuration():
+    try:
+        import platform
+        print(f"- system: {platform.system()} {platform.release()}")
+    except ImportError:
+        print(f"- system: please install the `platform` to have this information")
+
+    try:
+        import distro
+        print(f"- OS: {distro.linux_distribution(full_distribution_name=False)[0]} "
+              f"{distro.linux_distribution(full_distribution_name=False)[1]}")
+    except ImportError:
+        print(f"- OS: please install the `distro` to have this information")
+
+    try:
+        import cpuinfo
+        info_ = cpuinfo.get_cpu_info()
+        print(f"- processor: {info_['brand_raw']}")
+        print(f"- python version: {info_['python_version']}")
+    except ImportError:
+        print(f"- processor: please install the `py-cpuinfo` to have this information")
+        print(f"- python version: please install the `py-cpuinfo` to have this information")
+
+    import pandas as pd
+    import pandapower as pp
+    import lightsim2grid
+    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()
diff --git a/lightsim2grid/LightSimBackend.py b/lightsim2grid/LightSimBackend.py
index 6f2fbad..3a99553 100644
--- a/lightsim2grid/LightSimBackend.py
+++ b/lightsim2grid/LightSimBackend.py
@@ -84,13 +84,6 @@ def __init__(self, detailed_infos_for_cascading_failures=False):
 
         self.thermal_limit_a = None
 
-        self._iref_slack = None
-        self._id_bus_added = None
-        self._fact_mult_gen = -1
-        self._what_object_where = None
-        self._number_true_line = -1
-        self._corresp_name_fun = {}
-        self._get_vector_inj = {}
         self.dim_topo = -1
         self._init_action_to_set = None
         self._backend_action_class = None
@@ -130,8 +123,16 @@ def get_theta(self):
         storage_theta: ``numpy.ndarray``
             Gives the voltage angle to the bus at which each storage unit is connected
         """
-        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()))
+        line_or_theta = np.empty(self.n_line)
+        line_or_theta[:self.__nb_powerline] = self._grid.get_lineor_theta()
+        line_or_theta[self.__nb_powerline:] = self._grid.get_trafohv_theta()
+
+        line_ex_theta = np.empty(self.n_line)
+        line_ex_theta[:self.__nb_powerline] = self._grid.get_lineex_theta()
+        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()))
         load_theta = self.cst_1 * self._grid.get_load_theta()
         gen_theta = self.cst_1 * self._grid.get_gen_theta()
         storage_theta = self.cst_1 * self._grid.get_storage_theta()
@@ -563,41 +564,49 @@ def runpf(self, is_dc=False):
                     self.V = np.ones(self.nb_bus_total, dtype=np.complex_) * self._grid.get_init_vm_pu()
                 V = self._grid.dc_pf(self.V, self.max_it, self.tol)
                 if V.shape[0] == 0:
-                    raise DivergingPowerFlow("divergence of powerflow (non connected grid)")
+                    raise DivergingPowerFlow("Divergence of powerflow (non connected grid)")
             else:
                 if (self.V is None) or (self.V.shape[0] == 0):
-                    # init from dc approx in this case
+                    # create the vector V as it is not created
                     self.V = np.ones(self.nb_bus_total, dtype=np.complex_) * self._grid.get_init_vm_pu()
 
                 if self.initdc:
                     self._grid.deactivate_result_computation()
-                    V = self._grid.dc_pf(copy.deepcopy(self.V), self.max_it, self.tol)
+                    # if I init with dc values, it should depends on previous state
+                    self.V[:] = self._grid.get_init_vm_pu()  # see issue 30
+                    Vdc = self._grid.dc_pf(copy.deepcopy(self.V), self.max_it, self.tol)
                     self._grid.reactivate_result_computation()
-                    if V.shape[0] == 0:
+                    if Vdc.shape[0] == 0:
                         raise DivergingPowerFlow("divergence of powerflow (non connected grid)")
-                    self.V[:] = V
-
-                V = self._grid.ac_pf(self.V, self.max_it, self.tol)
+                    V_init = Vdc
+                else:
+                    V_init = copy.deepcopy(self.V)
+                V = self._grid.ac_pf(V_init, self.max_it, self.tol)
                 if V.shape[0] == 0:
                     # V = self._grid.ac_pf(self.V, self.max_it, self.tol)
-                    raise DivergingPowerFlow("divergence of powerflow")
+                    raise DivergingPowerFlow(f"divergence of powerflow (more than {self.max_it} iterations)")
 
             self.comp_time += self._grid.get_computation_time()
             self.V[:] = V
-            lpor, lqor, lvor, laor = self._grid.get_lineor_res()
-            lpex, lqex, lvex, laex = self._grid.get_lineex_res()
-            tpor, tqor, tvor, taor = self._grid.get_trafohv_res()
-            tpex, tqex, tvex, taex = self._grid.get_trafolv_res()
-
-            self.p_or[:] = np.concatenate((lpor, tpor))
-            self.q_or[:] = np.concatenate((lqor, tqor))
-            self.v_or[:] = np.concatenate((lvor, tvor))
-            self.a_or[:] = 1000. * np.concatenate((laor, taor))
-
-            self.p_ex[:] = np.concatenate((lpex, tpex))
-            self.q_ex[:] = np.concatenate((lqex, tqex))
-            self.v_ex[:] = np.concatenate((lvex, tvex))
-            self.a_ex[:] = 1000. * np.concatenate((laex, taex))
+            (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.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.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.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_or *= 1000.  # amps in lightsim, kA expected in grid2op
+            self.a_ex *= 1000.  # amps in lightsim, kA expected in grid2op
 
             self.a_or[~np.isfinite(self.a_or)] = 0.
             self.v_or[~np.isfinite(self.v_or)] = 0.
@@ -650,7 +659,7 @@ def _fill_nans(self):
             self.storage_p[:] = np.NaN
             self.storage_q[:] = np.NaN
             self.storage_v[:] = np.NaN
-        res = False
+        self.V[:] = self._grid.get_init_vm_pu()  # reset the V to its "original" value (see issue 30)
 
     def copy(self):
         # i can perform a regular copy, everything has been initialized
@@ -665,11 +674,48 @@ def copy(self):
         self._grid = None
         self.__me_at_init = None
         self.init_pp_backend = None
-        res = copy.deepcopy(self)
 
+        ####################
+        # res = copy.deepcopy(self)  # super slow
+        res = type(self).__new__(type(self))
+
+        # copy the regular attribute
+        res.__has_storage = self.__has_storage
+        res.__current_solver_type = self.__current_solver_type
+        res.__nb_powerline = self.__nb_powerline
+        res.__nb_bus_before = self.__nb_bus_before
+        res.cst_1 = dt_float(1.0)
+        li_regular_attr = ["detailed_infos_for_cascading_failures", "comp_time", "can_output_theta", "_is_loaded",
+                           "nb_bus_total", "initdc",
+                           "_big_topo_to_obj", "max_it", "tol", "dim_topo",
+                           "_idx_hack_storage"]
+        for attr_nm in li_regular_attr:
+            setattr(res, attr_nm, copy.deepcopy(getattr(self, attr_nm)))
+
+        # copy the numpy array
+        res.__nb_bus_before = copy.deepcopy(self.__nb_bus_before)
+        li_attr_npy = ["thermal_limit_a", "_sh_vnkv", "_init_bus_load", "_init_bus_gen",
+                       "_init_bus_lor", "_init_bus_lex", "nb_obj_per_bus", "next_prod_p", "topo_vect",
+                       "shunt_topo_vect", "V", "prod_pu_to_kv", "load_pu_to_kv", "lines_or_pu_to_kv",
+                       "lines_ex_pu_to_kv",
+                       "p_or", "q_or", "v_or", "a_or",
+                       "p_ex", "q_ex", "v_ex", "a_ex",
+                       "load_p", "load_q", "load_v",
+                       "prod_p", "prod_q", "prod_v",
+                       "storage_p", "storage_q", "storage_v",
+                       ]
+        for attr_nm in li_attr_npy:
+            setattr(res, attr_nm, copy.deepcopy(getattr(self, attr_nm)))
+        ###############
+
+        # handle the most complicated
         res._grid = mygrid.copy()
-        res.__me_at_init = __me_at_init.copy()
-        res.init_pp_backend = inippbackend.copy()
+        res.__me_at_init = __me_at_init.copy()  # this is const
+        res.init_pp_backend = inippbackend  # this is const
+        res._init_action_to_set = copy.deepcopy(self._init_action_to_set)
+        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
 
         self._grid = mygrid
         self.init_pp_backend = inippbackend
@@ -735,7 +781,6 @@ def get_current_solver_type(self):
         return self.__current_solver_type
 
     def reset(self, grid_path, grid_filename=None):
-        self.V = None
         self._fill_nans()
         self._grid = self.__me_at_init.copy()
         self._grid.change_solver(self.__current_solver_type)
diff --git a/lightsim2grid/_aux_add_gen.py b/lightsim2grid/_aux_add_gen.py
index 590cf01..256f155 100644
--- a/lightsim2grid/_aux_add_gen.py
+++ b/lightsim2grid/_aux_add_gen.py
@@ -7,6 +7,7 @@
 # This file is part of LightSim2grid, LightSim2grid implements a c++ backend targeting the Grid2Op platform.
 import numpy as np
 
+
 def _aux_add_gen(model, pp_net):
     """
     Add the generators of the pp_net into the lightsim2grid "model"
diff --git a/lightsim2grid/solver/__init__.py b/lightsim2grid/solver/__init__.py
index 8c01286..ecef04c 100644
--- a/lightsim2grid/solver/__init__.py
+++ b/lightsim2grid/solver/__init__.py
@@ -7,6 +7,13 @@
     # KLU is not available
     pass
 
+try:
+    from lightsim2grid_cpp import NICSLUSolver
+    __all__.append("NICSLUSolver")
+except Exception as exc_:
+    # NICSLU is not available
+    pass
+
 from lightsim2grid_cpp import SparseLUSolver
 from lightsim2grid_cpp import GaussSeidelSolver
 from lightsim2grid_cpp import GaussSeidelSynchSolver
diff --git a/setup.py b/setup.py
index cc2a891..56464a6 100644
--- a/setup.py
+++ b/setup.py
@@ -150,6 +150,46 @@
     src_files.append("src/KLUSolver.cpp")
     extra_compile_args_tmp.append("-DKLU_SOLVER_AVAILABLE")
 
+# Try to locate the NICSLU sparse linaer solver
+if "PATH_NICSLU" in os.environ:
+    # user indicate the path for the NICSLU library (see https://github.com/chenxm1986/nicslu)
+    # eg "/home/user/Documents/nicslu/nicslu202103/"
+
+    path_nicslu = os.path.abspath(os.environ["PATH_NICSLU"])
+    include_nicslu = True
+    # check for appropriate license
+    if not os.path.exists(path_nicslu):
+        include_nicslu = False
+    if include_nicslu and not os.path.exists(os.path.join(path_nicslu, "license")):
+        # license not located at the right directory
+        include_nicslu = False
+    if include_nicslu and not os.path.exists(os.path.join(path_nicslu, "license", "nicslu.lic")):
+        # no license found
+        include_nicslu = False
+    libnicslu_path = None
+    if include_nicslu:
+        if sys.platform.startswith("linux") or sys.platform.startswith("darwin"):
+            libnicslu_path = os.path.join(path_nicslu, "linux/lib_centos6_x64_gcc482_fma/int32/libnicslu.so")
+            if not os.path.exists(libnicslu_path):
+                include_nicslu = False
+                libnicslu_path = None
+        elif sys.platform.startswith("win"):
+            libnicslu_path = os.path.join(path_nicslu, "windows/lib_win7_x64_fma/int32/nicslu.lib")
+            if not os.path.exists(libnicslu_path):
+                include_nicslu = False
+                libnicslu_path = None
+
+    if include_nicslu and libnicslu_path is not None:
+        LIBS.append(os.path.join(path_nicslu, libnicslu_path))
+        include_dirs.append(os.path.join(path_nicslu, "include"))
+        src_files.append("src/NICSLUSolver.cpp")
+        extra_compile_args.append("-DNICSLU_SOLVER_AVAILABLE")
+
+if "__COUT_TIMES" in os.environ:
+    # to add extra info in cout for the computation times, we do not recommend to use it !
+    if os.environ["__COUT_TIMES"] == "1":
+        extra_compile_args.append("-D__COUT_TIMES")
+
 ext_modules = [
     Pybind11Extension(
         'lightsim2grid_cpp',
@@ -177,7 +217,9 @@
         "benchmark": [
             "tabulate",
             "grid2op>=1.5.0",
-            "numpy"
+            "numpy",
+            "distro",
+            "py-cpuinfo"
         ],
         "recommended": [
             "grid2op>=1.5.0"
diff --git a/src/BaseNRSolver.cpp b/src/BaseNRSolver.cpp
index 3ae8f73..4465a11 100644
--- a/src/BaseNRSolver.cpp
+++ b/src/BaseNRSolver.cpp
@@ -93,18 +93,19 @@ bool BaseNRSolver::compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
         res = false;
     }
     timer_total_nr_ += timer.duration();
-//    std::cout << "Computation time: " << "\n\t timer_initialize_: " << timer_initialize_
-//              << "\n\t timer_dSbus_: " << timer_dSbus_
-//              << "\n\t timer_fillJ_: " << timer_fillJ_
-//              << "\n\t timer_Fx_: " << timer_Fx_
-//              << "\n\t timer_check_: " << timer_check_
-//              << "\n\t timer_solve_: " << timer_solve_
-//              << "\n\t timer_total_nr_: " << timer_total_nr_
-//              << "\n\n";
+    #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_fillJ_: " << timer_fillJ_
+                  << "\n\t timer_Fx_: " << timer_Fx_
+                  << "\n\t timer_check_: " << timer_check_
+                  << "\n\t timer_solve_: " << timer_solve_
+                  << "\n\t timer_total_nr_: " << timer_total_nr_
+                  << "\n\n";
+    #endif // __COUT_TIMES
     return res;
 }
 
-
 void BaseNRSolver::reset(){
     BaseSolver::reset();
     // reset specific attributes
@@ -117,53 +118,53 @@ void BaseNRSolver::reset(){
 void BaseNRSolver::_dSbus_dV(const Eigen::Ref<const Eigen::SparseMatrix<cplx_type> > & Ybus,
                              const Eigen::Ref<const CplxVect > & V){
     auto timer = CustTimer();
-    auto size_dS = V.size();
-    CplxVect Vnorm = V.array() / V.array().abs();
-    CplxVect Ibus = Ybus * V;
+    const auto size_dS = V.size();
+    const CplxVect Vnorm = V.array() / V.array().abs();
+    const CplxVect Ibus = Ybus * V;
 
-    // TODO see if i can reuse previous values, i am not sure
     dS_dVm_ = Ybus;
     dS_dVa_ = Ybus;
 
-    // i fill the buffer columns per columns
-    for (int k=0; k < size_dS; ++k){
-        for (Eigen::SparseMatrix<cplx_type>::InnerIterator it(dS_dVm_,k); it; ++it)
+    cplx_type * ds_dvm_x_ptr = dS_dVm_.valuePtr();
+    cplx_type * ds_dva_x_ptr = dS_dVa_.valuePtr();
+    unsigned int pos_el = 0;
+    for (int col_=0; col_ < size_dS; ++col_){
+        for (Eigen::Ref<const Eigen::SparseMatrix<cplx_type> >::InnerIterator it(Ybus, col_); it; ++it)
         {
-            it.valueRef() *= Vnorm(it.col());  // dS_dVm[k] *= Vnorm[Yj[k]]
-            it.valueRef() = std::conj(it.valueRef()) * V(it.row());  // dS_dVm[k] = conj(dS_dVm[k]) * V[r]
-            if(it.col() == it.row()){
-                // diagonal element
-                it.valueRef() += std::conj(Ibus(it.row())) * Vnorm(it.row()); // dS_dVm[k] += buffer[r] # buffer being conj(Ibus) * Vnorm
-            }
-        }
-    }
+            const int row_id = it.row();
+            const int col_id = it.col();
+            cplx_type & ds_dvm_el = ds_dvm_x_ptr[pos_el];
+            cplx_type & ds_dva_el = ds_dva_x_ptr[pos_el];
 
-    for (int k=0; k < size_dS; ++k){
-        for (Eigen::SparseMatrix<cplx_type>::InnerIterator it(dS_dVa_,k); it; ++it)
-        {
-            it.valueRef() *= V(it.col());  // dS_dVa[k] *= V[Yj[k]]
-            if(it.col() == it.row()){
-                // diagonal element
-                it.valueRef() -= Ibus(it.row());  // dS_dVa[k] = -Ibus[r] + dS_dVa[k]
+            ds_dvm_el *= Vnorm(col_id);  // dS_dVm[k] *= Vnorm[Yj[k]]
+            ds_dvm_el = std::conj(ds_dvm_el) * V(row_id);  // dS_dVm[k] = conj(dS_dVm[k]) * V[r]
+
+            ds_dva_el *= V(col_id);  // dS_dVa[k] *= V[Yj[k]]
+
+            if(col_id == row_id)
+            {
+                ds_dvm_el += std::conj(Ibus(row_id)) * Vnorm(row_id); // dS_dVm[k] += conj(Ibus) * Vnorm
+                ds_dva_x_ptr[pos_el] -= Ibus(row_id);  // dS_dVa[k] = -Ibus[r] + dS_dVa[k]
             }
-            cplx_type tmp = my_i * V(it.row());
-            it.valueRef() = std::conj(-it.valueRef()) * tmp;  // dS_dVa[k] = conj(-dS_dVa[k]) * (1j * V[r])
+            cplx_type tmp = my_i * V(row_id);
+            ds_dva_el = std::conj(-ds_dva_el) * tmp;  // dS_dVa[k] = conj(-dS_dVa[k]) * (1j * V[r])
+
+            // go to next element
+            ++pos_el;
         }
     }
-    dS_dVa_.makeCompressed();
-    dS_dVm_.makeCompressed();
     timer_dSbus_ += timer.duration();
 }
 
 void BaseNRSolver::_get_values_J(int & nb_obj_this_col,
-                              std::vector<int> & inner_index,
-                              std::vector<real_type> & values,
-                              const Eigen::SparseMatrix<real_type> & mat,  // ex. dS_dVa_r
-                              const std::vector<int> & index_row_inv, // ex. pvpq_inv
-                              const Eigen::VectorXi & index_col, // ex. pvpq
-                              int col_id,
-                              int row_lag  // 0 for J11 for example, n_pvpq for J12
-                             )
+                                 std::vector<int> & inner_index,
+                                 std::vector<real_type> & values,
+                                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & mat,  // ex. dS_dVa_r
+                                 const std::vector<int> & index_row_inv, // ex. pvpq_inv
+                                 const Eigen::VectorXi & index_col, // ex. pvpq
+                                 int col_id,
+                                 int row_lag  // 0 for J11 for example, n_pvpq for J12
+                                 )
 {
     /**
     This function will fill the "inner_index" and "values" with the non zero values
@@ -177,17 +178,18 @@ void BaseNRSolver::_get_values_J(int & nb_obj_this_col,
     **/
     int col_id_mat = index_col(col_id);
 
-    int start_id = mat.outerIndexPtr()[col_id_mat];
-    int end_id = mat.outerIndexPtr()[col_id_mat+1];
+    const int start_id = mat.outerIndexPtr()[col_id_mat];
+    const int end_id = mat.outerIndexPtr()[col_id_mat+1];
+    const real_type * val_prt = mat.valuePtr();
     for(int obj_id = start_id; obj_id < end_id; ++obj_id)
     {
-        int row_id_dS_dVa = mat.innerIndexPtr()[obj_id];
+        const int row_id_dS_dVa = mat.innerIndexPtr()[obj_id];
         // I add the value only if the rows was selected in the indexes
-        int row_id = index_row_inv[row_id_dS_dVa];
+        const int row_id = index_row_inv[row_id_dS_dVa];
         if(row_id >= 0)
         {
             inner_index.push_back(row_id+row_lag);
-            values.push_back(mat.valuePtr()[obj_id]);
+            values.push_back(val_prt[obj_id]);
             nb_obj_this_col++;
         }
     }
@@ -205,7 +207,7 @@ void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix<cplx_type> & Y
     J has the shape
     | J11 | J12 |               | (pvpq, pvpq) | (pvpq, pq) |
     | --------- | = dimensions: | ------------------------- |
-    | J21 | J22 |               |  (pq, pvpq)  | (pvpq, pq) |
+    | J21 | J22 |               |  (pq, pvpq)  | (pq, pq) |
     python implementation:
     J11 = dS_dVa[array([pvpq]).T, pvpq].real
     J12 = dS_dVm[array([pvpq]).T, pq].real
@@ -215,30 +217,96 @@ void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix<cplx_type> & Y
 
     auto timer = CustTimer();
     _dSbus_dV(Ybus, V);
-    Eigen::SparseMatrix<real_type> dS_dVa_r = dS_dVa_.real();
-    Eigen::SparseMatrix<real_type> dS_dVa_i = dS_dVa_.imag();
-    Eigen::SparseMatrix<real_type> dS_dVm_r = dS_dVm_.real();
-    Eigen::SparseMatrix<real_type> dS_dVm_i = dS_dVm_.imag();
+
+    const Eigen::SparseMatrix<real_type> dS_dVa_r = dS_dVa_.real();
+    const Eigen::SparseMatrix<real_type> dS_dVa_i = dS_dVa_.imag();
+    const Eigen::SparseMatrix<real_type> dS_dVm_r = dS_dVm_.real();
+    const Eigen::SparseMatrix<real_type> dS_dVm_i = dS_dVm_.imag();
 
     const int n_pvpq = pvpq.size();
     const int n_pq = pq.size();
-
     const int size_j = n_pvpq + n_pq;
 
+    // TODO to gain a bit more time below, try to compute directly, in _dSbus_dV(Ybus, V);
+    // TODO the `dS_dVa_[pvpq, pvpq]`
+    // TODO so that it's easier to retrieve in the next few lines !
+    if(J_.cols() != size_j)
+    {
+        // auto timer2 = CustTimer();
+        // first time i initialized the matrix, so i need to compute its sparsity pattern
+        fill_jacobian_matrix_unkown_sparsity_pattern(dS_dVa_r, dS_dVa_i, dS_dVm_r, dS_dVm_i,
+                                                     Ybus, V, pq, pvpq, pq_inv, pvpq_inv
+                                                     );
+        // std::cout << "fill_jacobian_matrix_unkown_sparsity_pattern : " << timer2.duration() << std::endl;
+    }else{
+        // the sparsity pattern of J_ is already known, i can reuse it to fill it
+        // properly and faster
+        // auto timer3 = CustTimer();
+        fill_jacobian_matrix_kown_sparsity_pattern(dS_dVa_r, dS_dVa_i, dS_dVm_r, dS_dVm_i,
+                                                     Ybus, V, pq, pvpq, pq_inv, pvpq_inv
+                                                     );
+        // std::cout << "fill_jacobian_matrix_kown_sparsity_pattern : " << timer3.duration() << std::endl;
+    }
+    timer_fillJ_ += timer.duration();
+}
+
+void BaseNRSolver::fill_jacobian_matrix_unkown_sparsity_pattern(
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_r,
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_i,
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_r,
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_i,
+        const Eigen::SparseMatrix<cplx_type> & Ybus,
+        const CplxVect & V,
+        const Eigen::VectorXi & pq,
+        const Eigen::VectorXi & pvpq,
+        const std::vector<int> & pq_inv,
+        const std::vector<int> & pvpq_inv
+    )
+{
+    /**
+    This functions fills the jacobian matrix when its sparsity pattern is not know in advance (typically
+    the first iteration of the Newton Raphson)
+
+    For that we need to perform relatively expensive computation from dS_dV* in order to retrieve it.
+
+    This function is NOT optimized for speed...
+
+    Remember:
+
+    J has the shape
+    | J11 | J12 |               | (pvpq, pvpq) | (pvpq, pq) |
+    | --------- | = dimensions: | ------------------------- |
+    | J21 | J22 |               |  (pq, pvpq)  | (pq, pq) |
+    python implementation:
+    J11 = dS_dVa[array([pvpq]).T, pvpq].real
+    J12 = dS_dVm[array([pvpq]).T, pq].real
+    J21 = dS_dVa[array([pq]).T, pvpq].imag
+    J22 = dS_dVm[array([pq]).T, pq].imag
+
+    **/
     bool need_insert = false;  // i optimization: i don't need to insert the coefficient in the matrix
+    const int n_pvpq = pvpq.size();
+    const int n_pq = pq.size();
+    const int size_j = n_pvpq + n_pq;
+
+    // Method (1) seems to be faster than the others
+
+    // optim : if the matrix was already computed, i don't initialize it, i instead reuse as much as i can
+    // i can do that because the matrix will ALWAYS have the same non zero coefficients.
+    // in this if, i allocate it in a "large enough" place to avoid copy when first filling it
     if(J_.cols() != size_j)
     {
-        // optim : if the matrix was already computed, i don't initalize it, i instead reuse as much as i can
-        // i can do that because the matrix will ALWAYS have the same non zero coefficients.
-        // in this if, i allocate it in a "large enough" place to avoid copy when first filling it
         need_insert = true;
         J_ = Eigen::SparseMatrix<real_type>(size_j,size_j);
         // pre allocate a large enough matrix
         J_.reserve(2*(dS_dVa_.nonZeros()+dS_dVm_.nonZeros()));
-        // from an experiment, outerIndexPtr is inialized, with the number of columns
+        // from an experiment, outerIndexPtr is initialized, with the number of columns
         // innerIndexPtr and valuePtr are not.
     }
 
+    // std::vector<Eigen::Triplet<double> >coeffs;  // HERE FOR PERF OPTIM (3)
+    // coeffs.reserve(2*(dS_dVa_.nonZeros()+dS_dVm_.nonZeros()));  // HERE FOR PERF OPTIM (3)
+
     // i fill the buffer columns per columns
     int nb_obj_this_col = 0;
     std::vector<int> inner_index;
@@ -268,23 +336,11 @@ void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix<cplx_type> & Y
         // "efficient" insert of the element in the matrix
         for(int in_ind=0; in_ind < nb_obj_this_col; ++in_ind){
             int row_id = inner_index[in_ind];
-            if(need_insert) J_.insert(row_id, col_id) = values[in_ind];
-            else J_.coeffRef(row_id, col_id) = values[in_ind];
+            if(need_insert) J_.insert(row_id, col_id) = values[in_ind];  // HERE FOR PERF OPTIM (1)
+            else J_.coeffRef(row_id, col_id) = values[in_ind];  // HERE FOR PERF OPTIM (1)
+            // J_.insert(row_id, col_id) = values[in_ind];  // HERE FOR PERF OPTIM (2)
+            // coeffs.push_back(Eigen::Triplet<double>(row_id, col_id, values[in_ind]));   // HERE FOR PERF OPTIM (3)
         }
-//                if(need_insert){
-//                    for(int in_ind=0; in_ind < nb_obj_this_col; ++in_ind){
-//                        int row_id = inner_index[in_ind];
-//                        J_.insert(row_id, col_id) = values[in_ind];
-//                        // else J_.coeffRef(row_id, col_id) = values[in_ind];
-//                    }
-//                }else{
-//                    int in_ind=0;
-//                    for (Eigen::SparseMatrix<real_type>::InnerIterator it(J_,col_id); it; ++it, ++in_ind)
-//                    {
-////                        int row_id = inner_index[it.row()];
-//                        it.valueRef() = values[it.row()];
-//                    }
-//                }
     }
 
     //TODO make same for the second part (have a funciton for previous loop)
@@ -312,24 +368,100 @@ void BaseNRSolver::fill_jacobian_matrix(const Eigen::SparseMatrix<cplx_type> & Y
         // "efficient" insert of the element in the matrix
         for(int in_ind=0; in_ind < nb_obj_this_col; ++in_ind){
             int row_id = inner_index[in_ind];
-            if(need_insert) J_.insert(row_id, col_id + n_pvpq) = values[in_ind];
-            else J_.coeffRef(row_id, col_id + n_pvpq) = values[in_ind];
+            if(need_insert) J_.insert(row_id, col_id + n_pvpq) = values[in_ind];  // HERE FOR PERF OPTIM (1)
+            else J_.coeffRef(row_id, col_id + n_pvpq) = values[in_ind];  // HERE FOR PERF OPTIM (1)
+            // J_.insert(row_id, col_id + n_pvpq) = values[in_ind];  // HERE FOR PERF OPTIM (2)
+            // coeffs.push_back(Eigen::Triplet<double>(row_id, col_id + n_pvpq, values[in_ind]));   // HERE FOR PERF OPTIM (3)
         }
-//                if(need_insert){
-//                    for(int in_ind=0; in_ind < nb_obj_this_col; ++in_ind){
-//                        int row_id = inner_index[in_ind];
-//                        J_.insert(row_id, col_id + n_pvpq) = values[in_ind];
-//                        // else J_.coeffRef(row_id, col_id) = values[in_ind];
-//                    }
-//                }else{
-//                    int in_ind=0;
-//                    for (Eigen::SparseMatrix<real_type>::InnerIterator it(J_, col_id + n_pvpq); it; ++it, ++in_ind)
-//                    {
-//                        int row_id = inner_index[it.row()];
-//                        it.valueRef() = values[row_id];
-//                    }
-//                }
     }
+    // J_.setFromTriplets(coeffs.begin(), coeffs.end());  // HERE FOR PERF OPTIM (3)
     J_.makeCompressed();
-    timer_fillJ_ += timer.duration();
 }
+
+void BaseNRSolver::fill_jacobian_matrix_kown_sparsity_pattern(
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_r,
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_i,
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_r,
+        const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_i,
+        const Eigen::SparseMatrix<cplx_type> & Ybus,
+        const CplxVect & V,
+        const Eigen::VectorXi & pq,
+        const Eigen::VectorXi & pvpq,
+        const std::vector<int> & pq_inv,
+        const std::vector<int> & pvpq_inv
+    )
+{
+    /**
+    This functions fills the jacobian matrix when its sparsity pattern is KNOWN in advance (typically
+    the second and next iterations of the Newton Raphson)
+
+    We don't need to perform heavy computation but just to read the right data from the right matrix!
+
+    This function is optimized for speed. In particular, it does not "uncompress" the J_ matrix and only
+    change the value pointer (not the inner nor outer pointer)
+
+    It is optimized only if J_ is in default Eigen format (column)
+
+    Remember:
+
+    J has the shape
+    | J11 | J12 |               | (pvpq, pvpq) | (pvpq, pq) |
+    | --------- | = dimensions: | ------------------------- |
+    | J21 | J22 |               |  (pq, pvpq)  | (pq, pq) |
+    python implementation:
+    J11 = dS_dVa[array([pvpq]).T, pvpq].real
+    J12 = dS_dVm[array([pvpq]).T, pq].real
+    J21 = dS_dVa[array([pq]).T, pvpq].imag
+    J22 = dS_dVm[array([pq]).T, pq].imag
+
+    **/
+
+    const int n_pvpq = pvpq.size();
+
+    real_type * J_x_ptr = J_.valuePtr();
+    const int n_row = J_.cols();
+    unsigned int pos_el = 0;
+    for (int col_=0; col_ < n_row; ++col_){
+        for (Eigen::SparseMatrix<real_type>::InnerIterator it(J_, col_); it; ++it)
+        {
+            const int row_id = it.row();
+            const int col_id = it.col();  // it's equal to "col_"
+            real_type & this_el = J_x_ptr[pos_el];
+            // std::cout << "filling J_[" << row_id <<" ,"<< col_id << "]" << std::endl;
+            if((col_id < n_pvpq) && (row_id < n_pvpq)){
+                // this is the J11 part (dS_dVa_r)
+                const int row_id_dS_dVa_r = pvpq[row_id];
+                const int col_id_dS_dVa_r = pvpq[col_id];
+                this_el = dS_dVa_r.coeff(row_id_dS_dVa_r, col_id_dS_dVa_r);
+
+                // I don't need to perform these checks: if they failed, the element would not be in J_ in the first place
+                // const int is_row_non_null = pq_inv[row_id_dS_dVa_r];
+                // const int is_col_non_null = pq_inv[col_id_dS_dVa_r];
+                // if(is_row_non_null >= 0 && is_col_non_null >= 0)
+                //     this_el = dS_dVa_r.coeff(row_id_dS_dVa_r, col_id_dS_dVa_r);
+                // else
+                //     std::cout << "dS_dVa_r: missed" << std::endl;
+
+            }else if((col_id < n_pvpq) && (row_id >= n_pvpq)){
+                // this is the J21 part (dS_dVa_i)
+                const int row_id_dS_dVa_i = pq[row_id - n_pvpq];
+                const int col_id_dS_dVa_i = pvpq[col_id];
+                this_el = dS_dVa_i.coeff(row_id_dS_dVa_i, col_id_dS_dVa_i);
+            }else if((col_id >= n_pvpq) && (row_id < n_pvpq)){
+                // this is the J12 part (dS_dVm_r)
+                const int row_id_dS_dVm_r = pvpq[row_id];
+                const int col_id_dS_dVm_r = pq[col_id - n_pvpq];
+                this_el = dS_dVm_r.coeff(row_id_dS_dVm_r, col_id_dS_dVm_r);
+            }else if((col_id >= n_pvpq) && (row_id >= n_pvpq)){
+                // this is the J22 part (dS_dVm_i)
+                const int row_id_dS_dVm_i = pq[row_id - n_pvpq];
+                const int col_id_dS_dVm_i = pq[col_id - n_pvpq];
+                this_el = dS_dVm_i.coeff(row_id_dS_dVm_i, col_id_dS_dVm_i);
+
+            }
+
+            // go to the next element
+            ++pos_el;
+        }
+    }
+}
\ No newline at end of file
diff --git a/src/BaseNRSolver.h b/src/BaseNRSolver.h
index 21fca9e..11217e8 100644
--- a/src/BaseNRSolver.h
+++ b/src/BaseNRSolver.h
@@ -67,7 +67,7 @@ class BaseNRSolver : public BaseSolver
         void _get_values_J(int & nb_obj_this_col,
                            std::vector<int> & inner_index,
                            std::vector<real_type> & values,
-                           const Eigen::SparseMatrix<real_type> & mat,  // ex. dS_dVa_r
+                           const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & mat,  // ex. dS_dVa_r
                            const std::vector<int> & index_row_inv, // ex. pvpq_inv
                            const Eigen::VectorXi & index_col, // ex. pvpq
                            int col_id,
@@ -81,6 +81,30 @@ class BaseNRSolver : public BaseSolver
                                   const std::vector<int> & pq_inv,
                                   const std::vector<int> & pvpq_inv
                                   );
+        void fill_jacobian_matrix_kown_sparsity_pattern(
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_r,
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_i,
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_r,
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_i,
+                 const Eigen::SparseMatrix<cplx_type> & Ybus,
+                 const CplxVect & V,
+                 const Eigen::VectorXi & pq,
+                 const Eigen::VectorXi & pvpq,
+                 const std::vector<int> & pq_inv,
+                 const std::vector<int> & pvpq_inv
+                 );
+        void fill_jacobian_matrix_unkown_sparsity_pattern(
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_r,
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVa_i,
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_r,
+                 const Eigen::Ref<const Eigen::SparseMatrix<real_type> > & dS_dVm_i,
+                 const Eigen::SparseMatrix<cplx_type> & Ybus,
+                 const CplxVect & V,
+                 const Eigen::VectorXi & pq,
+                 const Eigen::VectorXi & pvpq,
+                 const std::vector<int> & pq_inv,
+                 const std::vector<int> & pvpq_inv
+                 );
 
     protected:
 
diff --git a/src/BaseSolver.h b/src/BaseSolver.h
index b10b958..3290f76 100644
--- a/src/BaseSolver.h
+++ b/src/BaseSolver.h
@@ -42,13 +42,13 @@ class BaseSolver : public BaseConstants
 
         ~BaseSolver(){}
 
-        Eigen::Ref<RealVect> get_Va(){
+        Eigen::Ref<const RealVect> get_Va() const{
             return Va_;
         }
-        Eigen::Ref<RealVect> get_Vm(){
+        Eigen::Ref<const RealVect> get_Vm() const{
             return Vm_;
         }
-        Eigen::Ref<CplxVect> get_V(){
+        Eigen::Ref<const CplxVect> get_V() const{
             return V_;
         }
         int get_error(){
diff --git a/src/ChooseSolver.cpp b/src/ChooseSolver.cpp
index 54f41a6..3396606 100644
--- a/src/ChooseSolver.cpp
+++ b/src/ChooseSolver.cpp
@@ -10,17 +10,17 @@
 #include <iostream>
 // template specialization
 template<SolverType ST>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp()
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp() const
 {
     throw std::runtime_error("ChooseSolver::get_V_tmp: Unknown solver type.");
 }
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::SparseLU>()
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::SparseLU>() const
 {
     return _solver_lu.get_V();
 }
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::KLU>()
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::KLU>() const
 {
     #ifndef KLU_SOLVER_AVAILABLE
         // I asked result of KLU solver without the required libraries
@@ -30,17 +30,27 @@ Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::KLU>()
     #endif
 }
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidel>()
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::NICSLU>() const
+{
+    #ifndef NICSLU_SOLVER_AVAILABLE
+        // I asked result of KLU solver without the required libraries
+        throw std::runtime_error("ChooseSolver::get_V_tmp: Impossible to use the NICSLU solver, that is not available on your plaform.");
+    #else
+        return _solver_nicslu.get_V();
+    #endif
+}
+template<>
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidel>() const
 {
     return _solver_gaussseidel.get_V();
 }
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidelSynch>()
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidelSynch>() const
 {
     return _solver_gaussseidelsynch.get_V();
 }
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::DC>()
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::DC>() const
 {
     return _solver_dc.get_V();
 }
@@ -123,6 +133,23 @@ bool ChooseSolver::compute_pf_tmp<SolverType::KLU>(const Eigen::SparseMatrix<cpl
         return _solver_klu.compute_pf(Ybus, V, Sbus, pv, pq, max_iter, tol);
     #endif
 }
+template<>
+bool ChooseSolver::compute_pf_tmp<SolverType::NICSLU>(const Eigen::SparseMatrix<cplx_type> & Ybus,
+                       CplxVect & V,
+                       const CplxVect & Sbus,
+                       const Eigen::VectorXi & pv,
+                       const Eigen::VectorXi & pq,
+                       int max_iter,
+                       real_type tol
+                       )
+{
+    #ifndef NICSLU_SOLVER_AVAILABLE
+        // I asked result of KLU solver without the required libraries
+        throw std::runtime_error("ChooseSolver::compute_pf_tmp: Impossible to use the NICSLU solver, that is not available on your plaform.");
+    #else
+        return _solver_nicslu.compute_pf(Ybus, V, Sbus, pv, pq, max_iter, tol);
+    #endif
+}
 
 template<SolverType ST>
 Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp()
@@ -159,34 +186,44 @@ Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::KLU>()
         return _solver_klu.get_J();
     #endif
 }
+template<>
+Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::NICSLU>()
+{
+    #ifndef NICSLU_SOLVER_AVAILABLE
+        // I asked result of KLU solver without the required libraries
+        throw std::runtime_error("ChooseSolver::get_J_tmp: Impossible to use the NICSLU solver, that is not available on your plaform.");
+    #else
+        return _solver_nicslu.get_J();
+    #endif
+}
 
 template<SolverType ST>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp()
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp() const
 {
     throw std::runtime_error("ChooseSolver::get_Va_tmp: Unknown solver type.");
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::SparseLU>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::SparseLU>() const
 {
     return _solver_lu.get_Va();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidel>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidel>() const
 {
     return _solver_gaussseidel.get_Va();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidelSynch>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidelSynch>() const
 {
     return _solver_gaussseidelsynch.get_Va();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::DC>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::DC>() const
 {
     return _solver_dc.get_Va();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::KLU>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::KLU>() const
 {
     #ifndef KLU_SOLVER_AVAILABLE
         // I asked result of KLU solver without the required libraries
@@ -195,33 +232,43 @@ Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::KLU>()
         return _solver_klu.get_Va();
     #endif
 }
+template<>
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::NICSLU>() const
+{
+    #ifndef NICSLU_SOLVER_AVAILABLE
+        // I asked result of KLU solver without the required libraries
+        throw std::runtime_error("ChooseSolver::get_Va_tmp: Impossible to use the NICSLU solver, that is not available on your plaform.");
+    #else
+        return _solver_nicslu.get_Va();
+    #endif
+}
 template<SolverType ST>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp()
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp() const
 {
     throw std::runtime_error("ChooseSolver::get_Vm_tmp: Unknown solver type.");
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::SparseLU>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::SparseLU>() const
 {
     return _solver_lu.get_Vm();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::DC>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::DC>() const
 {
     return _solver_dc.get_Vm();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidel>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidel>() const
 {
     return _solver_gaussseidel.get_Vm();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidelSynch>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidelSynch>() const
 {
     return _solver_gaussseidelsynch.get_Vm();
 }
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::KLU>()
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::KLU>() const
 {
     #ifndef KLU_SOLVER_AVAILABLE
         // I asked result of KLU solver without the required libraries
@@ -230,6 +277,16 @@ Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::KLU>()
         return _solver_klu.get_Vm();
     #endif
 }
+template<>
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::NICSLU>() const
+{
+    #ifndef NICSLU_SOLVER_AVAILABLE
+        // I asked result of KLU solver without the required libraries
+        throw std::runtime_error("ChooseSolver::get_Vm_tmp: Impossible to use the NICSLU solver, that is not available on your plaform.");
+    #else
+        return _solver_nicslu.get_Vm();
+    #endif
+}
 
 template<SolverType ST>
 double ChooseSolver::get_computation_time_tmp()
@@ -266,6 +323,17 @@ double ChooseSolver::get_computation_time_tmp<SolverType::KLU>()
     #endif
 }
 template<>
+double ChooseSolver::get_computation_time_tmp<SolverType::NICSLU>()
+{
+    #ifndef NICSLU_SOLVER_AVAILABLE
+        // I asked result of KLU solver without the required libraries
+        throw std::runtime_error("ChooseSolver::get_computation_time_tmp: Impossible to use the NICSLU solver, that is not available on your plaform.");
+    #else
+        const auto & res =  _solver_nicslu.get_timers();
+        return std::get<3>(res);
+    #endif
+}
+template<>
 double ChooseSolver::get_computation_time_tmp<SolverType::DC>()
 {
    const auto & res =  _solver_dc.get_timers();
@@ -274,13 +342,15 @@ double ChooseSolver::get_computation_time_tmp<SolverType::DC>()
 //TODO refactor all the functions above by making a template function "get_solver"
 
 // function definition
-Eigen::Ref<CplxVect> ChooseSolver::get_V(){
+Eigen::Ref<const CplxVect> ChooseSolver::get_V() const {
     check_right_solver();
     if(_solver_type == SolverType::SparseLU)
     {
          return get_V_tmp<SolverType::SparseLU>();
     }else if(_solver_type == SolverType::KLU){
          return get_V_tmp<SolverType::KLU>();
+    }else if(_solver_type == SolverType::NICSLU){
+         return get_V_tmp<SolverType::NICSLU>();
     }else if(_solver_type == SolverType::GaussSeidel){
          return get_V_tmp<SolverType::GaussSeidel>();
     }else if(_solver_type == SolverType::GaussSeidelSynch){
@@ -292,13 +362,15 @@ Eigen::Ref<CplxVect> ChooseSolver::get_V(){
     }
 }
 
-Eigen::Ref<RealVect> ChooseSolver::get_Va(){
+Eigen::Ref<const RealVect> ChooseSolver::get_Va() const {
     check_right_solver();
     if(_solver_type == SolverType::SparseLU)
     {
          return get_Va_tmp<SolverType::SparseLU>();
     }else if(_solver_type == SolverType::KLU){
          return get_Va_tmp<SolverType::KLU>();
+    }else if(_solver_type == SolverType::NICSLU){
+         return get_Va_tmp<SolverType::NICSLU>();
     }else if(_solver_type == SolverType::GaussSeidel){
          return get_Va_tmp<SolverType::GaussSeidel>();
     }else if(_solver_type == SolverType::GaussSeidelSynch){
@@ -309,13 +381,15 @@ Eigen::Ref<RealVect> ChooseSolver::get_Va(){
         throw std::runtime_error("ChooseSolver::get_Va: Unknown solver type.");
     }
 }
-Eigen::Ref<RealVect> ChooseSolver::get_Vm(){
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm() const {
     check_right_solver();
     if(_solver_type == SolverType::SparseLU)
     {
          return get_Vm_tmp<SolverType::SparseLU>();
     }else if(_solver_type == SolverType::KLU){
          return get_Vm_tmp<SolverType::KLU>();
+    }else if(_solver_type == SolverType::NICSLU){
+         return get_Vm_tmp<SolverType::NICSLU>();
     }else if(_solver_type == SolverType::GaussSeidel){
          return get_Vm_tmp<SolverType::GaussSeidel>();
     }else if(_solver_type == SolverType::GaussSeidelSynch){
@@ -342,6 +416,8 @@ bool ChooseSolver::compute_pf(const Eigen::SparseMatrix<cplx_type> & Ybus,
         return compute_pf_tmp<SolverType::SparseLU>(Ybus, V, Sbus, pv, pq, max_iter, tol);
     }else if(_solver_type == SolverType::KLU){
         return compute_pf_tmp<SolverType::KLU>(Ybus, V, Sbus, pv, pq, max_iter, tol);
+    }else if(_solver_type == SolverType::NICSLU){
+        return compute_pf_tmp<SolverType::NICSLU>(Ybus, V, Sbus, pv, pq, max_iter, tol);
     }else if(_solver_type == SolverType::GaussSeidel){
         return compute_pf_tmp<SolverType::GaussSeidel>(Ybus, V, Sbus, pv, pq, max_iter, tol);
     }else if(_solver_type == SolverType::GaussSeidelSynch){
@@ -360,6 +436,8 @@ Eigen::SparseMatrix<real_type> ChooseSolver::get_J(){
          return get_J_tmp<SolverType::SparseLU>();
     }else if(_solver_type == SolverType::KLU){
          return get_J_tmp<SolverType::KLU>();
+    }else if(_solver_type == SolverType::NICSLU){
+         return get_J_tmp<SolverType::NICSLU>();
     }else if(_solver_type == SolverType::GaussSeidel){
          return get_J_tmp<SolverType::GaussSeidel>();
     }else if(_solver_type == SolverType::GaussSeidelSynch){
@@ -379,6 +457,8 @@ double ChooseSolver::get_computation_time()
          return get_computation_time_tmp<SolverType::SparseLU>();
     }else if(_solver_type == SolverType::KLU){
          return get_computation_time_tmp<SolverType::KLU>();
+    }else if(_solver_type == SolverType::NICSLU){
+         return get_computation_time_tmp<SolverType::NICSLU>();
     }else if(_solver_type == SolverType::GaussSeidel){
          return get_computation_time_tmp<SolverType::GaussSeidel>();
     }else if(_solver_type == SolverType::GaussSeidelSynch){
diff --git a/src/ChooseSolver.h b/src/ChooseSolver.h
index a990ad4..7a04722 100644
--- a/src/ChooseSolver.h
+++ b/src/ChooseSolver.h
@@ -12,13 +12,18 @@
 #include<vector>
 
 // import newton raphson solvers using different linear algebra solvers
-#include "KLUSolver.h"
+#ifdef KLU_SOLVER_AVAILABLE
+    #include "KLUSolver.h"
+#endif
+#ifdef NICSLU_SOLVER_AVAILABLE
+    #include "NICSLUSolver.h"
+#endif
 #include "SparseLUSolver.h"
 #include "GaussSeidelSolver.h"
 #include "GaussSeidelSynchSolver.h"
 #include "DCSolver.h"
 
-enum class SolverType { SparseLU, KLU, GaussSeidel, DC, GaussSeidelSynch};
+enum class SolverType { SparseLU, KLU, GaussSeidel, DC, GaussSeidelSynch, NICSLU};
 
 
 // NB: when adding a new solver, you need to specialize the *tmp method (eg get_Va_tmp)
@@ -36,19 +41,29 @@ class ChooseSolver
             res.push_back(SolverType::SparseLU);
             res.push_back(SolverType::GaussSeidel);
             res.push_back(SolverType::DC);
+            res.push_back(SolverType::GaussSeidelSynch);
             #ifdef KLU_SOLVER_AVAILABLE
                 res.push_back(SolverType::KLU);
             #endif
-            res.push_back(SolverType::GaussSeidelSynch);
+            #ifdef NICSLU_SOLVER_AVAILABLE
+                res.push_back(SolverType::NICSLU);
+            #endif
             return res;
         }
         SolverType get_type() const {return _solver_type;}
         void change_solver(const SolverType & type)
         {
             if(type == _solver_type) return;
+            std::string msg;
             #ifndef KLU_SOLVER_AVAILABLE
                 // TODO better handling of that :-/
-                if(type == SolverType::KLU) throw std::runtime_error("Impossible to change for the KLU solver, that is not available on your platform.");
+                msg = "Impossible to change for the KLU solver, that is not available on your platform.";
+                if(type == SolverType::KLU) throw std::runtime_error(msg);
+            #endif
+            #ifndef NICSLU_SOLVER_AVAILABLE
+                // TODO better handling of that :-/
+                msg = "Impossible to change for the NICSLU solver, that is not available on your platform.";
+                if(type == SolverType::NICSLU) throw std::runtime_error(msg);
             #endif
             _solver_type = type;
         }
@@ -61,6 +76,9 @@ class ChooseSolver
             #ifdef KLU_SOLVER_AVAILABLE
                 _solver_klu.reset();
             #endif  // KLU_SOLVER_AVAILABLE
+            #ifdef NICSLU_SOLVER_AVAILABLE
+                _solver_nicslu.reset();
+            #endif  // NICSLU_SOLVER_AVAILABLE
             _solver_gaussseidelsynch.reset();
         }
 
@@ -74,29 +92,29 @@ class ChooseSolver
                         int max_iter,
                         real_type tol
                         );
-        Eigen::Ref<CplxVect> get_V();
+        Eigen::Ref<const CplxVect> get_V() const;
+        Eigen::Ref<const RealVect> get_Va() const;
+        Eigen::Ref<const RealVect> get_Vm() const;
         Eigen::SparseMatrix<real_type> get_J();
-        Eigen::Ref<RealVect> get_Va();
-        Eigen::Ref<RealVect> get_Vm();
         double get_computation_time();
 
     private:
-        void check_right_solver()
+        void check_right_solver() const
         {
             if(_solver_type != _type_used_for_nr) throw std::runtime_error("ChooseSolver: Solver mismatch: current solver is not the last solver used to perform a powerflow");
         }
 
         template<SolverType ST>
-        Eigen::Ref<CplxVect> get_V_tmp();
+        Eigen::Ref<const CplxVect> get_V_tmp() const;
 
         template<SolverType ST>
-        Eigen::SparseMatrix<real_type> get_J_tmp();
+        Eigen::Ref<const RealVect> get_Va_tmp() const;
 
         template<SolverType ST>
-        Eigen::Ref<RealVect> get_Va_tmp();
+        Eigen::Ref<const RealVect> get_Vm_tmp() const;
 
         template<SolverType ST>
-        Eigen::Ref<RealVect> get_Vm_tmp();
+        Eigen::SparseMatrix<real_type> get_J_tmp();
 
         template<SolverType ST>
         double get_computation_time_tmp();
@@ -123,6 +141,9 @@ class ChooseSolver
         #ifdef KLU_SOLVER_AVAILABLE
             KLUSolver _solver_klu;
         #endif  // KLU_SOLVER_AVAILABLE
+        #ifdef NICSLU_SOLVER_AVAILABLE
+            KLUSolver _solver_nicslu;
+        #endif  // NICSLU_SOLVER_AVAILABLE
 
 };
 
@@ -130,15 +151,17 @@ class ChooseSolver
 // template specialization
 
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::SparseLU>();
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::SparseLU>() const;
+template<>
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::KLU>() const;
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::KLU>();
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::NICSLU>() const;
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidel>();
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidel>() const;
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidelSynch>();
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::GaussSeidelSynch>() const;
 template<>
-Eigen::Ref<CplxVect> ChooseSolver::get_V_tmp<SolverType::DC>();
+Eigen::Ref<const CplxVect> ChooseSolver::get_V_tmp<SolverType::DC>() const;
 
 template<>
 bool ChooseSolver::compute_pf_tmp<SolverType::SparseLU>(const Eigen::SparseMatrix<cplx_type> & Ybus,
@@ -159,6 +182,15 @@ bool ChooseSolver::compute_pf_tmp<SolverType::KLU>(const Eigen::SparseMatrix<cpl
                        real_type tol
                        );
 template<>
+bool ChooseSolver::compute_pf_tmp<SolverType::NICSLU>(const Eigen::SparseMatrix<cplx_type> & Ybus,
+                       CplxVect & V,
+                       const CplxVect & Sbus,
+                       const Eigen::VectorXi & pv,
+                       const Eigen::VectorXi & pq,
+                       int max_iter,
+                       real_type tol
+                       );
+template<>
 bool ChooseSolver::compute_pf_tmp<SolverType::GaussSeidel>(const Eigen::SparseMatrix<cplx_type> & Ybus,
                        CplxVect & V,
                        const CplxVect & Sbus,
@@ -191,6 +223,8 @@ Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::SparseLU>();
 template<>
 Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::KLU>();
 template<>
+Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::NICSLU>();
+template<>
 Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::GaussSeidel>();
 template<>
 Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::GaussSeidelSynch>();
@@ -198,25 +232,27 @@ template<>
 Eigen::SparseMatrix<real_type> ChooseSolver::get_J_tmp<SolverType::DC>();
 
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::SparseLU>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::SparseLU>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::KLU>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::KLU>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidel>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::NICSLU>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidelSynch>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidel>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Va_tmp<SolverType::DC>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::GaussSeidelSynch>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::SparseLU>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Va_tmp<SolverType::DC>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::KLU>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::SparseLU>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidel>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::KLU>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidelSynch>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidel>() const;
 template<>
-Eigen::Ref<RealVect> ChooseSolver::get_Vm_tmp<SolverType::DC>();
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::GaussSeidelSynch>() const;
+template<>
+Eigen::Ref<const RealVect> ChooseSolver::get_Vm_tmp<SolverType::DC>() const;
 
 // computation times
 template<>
@@ -224,6 +260,8 @@ double ChooseSolver::get_computation_time_tmp<SolverType::SparseLU>();
 template<>
 double ChooseSolver::get_computation_time_tmp<SolverType::KLU>();
 template<>
+double ChooseSolver::get_computation_time_tmp<SolverType::NICSLU>();
+template<>
 double ChooseSolver::get_computation_time_tmp<SolverType::GaussSeidel>();
 template<>
 double ChooseSolver::get_computation_time_tmp<SolverType::GaussSeidelSynch>();
diff --git a/src/DataGen.cpp b/src/DataGen.cpp
index 9bc678c..d2d11af 100644
--- a/src/DataGen.cpp
+++ b/src/DataGen.cpp
@@ -127,9 +127,9 @@ void DataGen::fillpv(std::vector<int> & bus_pv,
     }
 }
 
-void DataGen::compute_results(const Eigen::Ref<RealVect> & Va,
-                               const Eigen::Ref<RealVect> & Vm,
-                               const Eigen::Ref<CplxVect> & V,
+void DataGen::compute_results(const Eigen::Ref<const RealVect> & Va,
+                               const Eigen::Ref<const RealVect> & Vm,
+                               const Eigen::Ref<const CplxVect> & V,
                                const std::vector<int> & id_grid_to_solver,
                                const RealVect & bus_vn_kv,
                                real_type sn_mva)
diff --git a/src/DataGen.h b/src/DataGen.h
index fb43832..061957f 100644
--- a/src/DataGen.h
+++ b/src/DataGen.h
@@ -148,9 +148,9 @@ class DataGen: public DataGeneric
                         const std::vector<int> & id_grid_to_solver);
     void init_q_vector(int nb_bus); // delta_q_per_gen_
 
-    void compute_results(const Eigen::Ref<RealVect> & Va,
-                         const Eigen::Ref<RealVect> & Vm,
-                         const Eigen::Ref<CplxVect> & V,
+    void compute_results(const Eigen::Ref<const RealVect> & Va,
+                         const Eigen::Ref<const RealVect> & Vm,
+                         const Eigen::Ref<const CplxVect> & V,
                          const std::vector<int> & id_grid_to_solver,
                          const RealVect & bus_vn_kv,
                          real_type sn_mva);
@@ -167,7 +167,7 @@ class DataGen: public DataGeneric
     void set_vm(CplxVect & V, const std::vector<int> & id_grid_to_solver);
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
-    const RealVect & get_theta() const {return 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/DataGeneric.cpp b/src/DataGeneric.cpp
index 19e5d82..afdbfe8 100644
--- a/src/DataGeneric.cpp
+++ b/src/DataGeneric.cpp
@@ -92,8 +92,8 @@ int DataGeneric::_get_bus(int el_id, const std::vector<bool> & status_, const Ei
     return res;
 }
 
-void DataGeneric::v_kv_from_vpu(const Eigen::Ref<RealVect> & Va,
-                                const Eigen::Ref<RealVect> & Vm,
+void DataGeneric::v_kv_from_vpu(const Eigen::Ref<const RealVect> & Va,
+                                const Eigen::Ref<const RealVect> & Vm,
                                 const std::vector<bool> & status,
                                 int nb_element,
                                 const Eigen::VectorXi & bus_me_id,
@@ -119,8 +119,8 @@ void DataGeneric::v_kv_from_vpu(const Eigen::Ref<RealVect> & Va,
 }
 
 
-void DataGeneric::v_deg_from_va(const Eigen::Ref<RealVect> & Va,
-                                const Eigen::Ref<RealVect> & Vm,
+void DataGeneric::v_deg_from_va(const Eigen::Ref<const RealVect> & Va,
+                                const Eigen::Ref<const RealVect> & Vm,
                                 const std::vector<bool> & status,
                                 int nb_element,
                                 const Eigen::VectorXi & bus_me_id,
diff --git a/src/DataGeneric.h b/src/DataGeneric.h
index fe76544..9be5486 100644
--- a/src/DataGeneric.h
+++ b/src/DataGeneric.h
@@ -98,8 +98,8 @@ class DataGeneric : public BaseConstants
         /**
         convert v from pu to v in kv (and assign it to the right element...)
         **/
-        void v_kv_from_vpu(const Eigen::Ref<RealVect> & Va,
-                           const Eigen::Ref<RealVect> & Vm,
+        void v_kv_from_vpu(const Eigen::Ref<const RealVect> & Va,
+                           const Eigen::Ref<const RealVect> & Vm,
                            const std::vector<bool> & status,
                            int nb_element,
                            const Eigen::VectorXi & bus_me_id,
@@ -111,8 +111,8 @@ class DataGeneric : public BaseConstants
         /**
         compute va in degree from va in rad.
         **/
-        void v_deg_from_va(const Eigen::Ref<RealVect> & Va,
-                           const Eigen::Ref<RealVect> & Vm,
+        void v_deg_from_va(const Eigen::Ref<const RealVect> & Va,
+                           const Eigen::Ref<const RealVect> & Vm,
                            const std::vector<bool> & status,
                            int nb_element,
                            const Eigen::VectorXi & bus_me_id,
diff --git a/src/DataLine.cpp b/src/DataLine.cpp
index edcacef..dbf3802 100644
--- a/src/DataLine.cpp
+++ b/src/DataLine.cpp
@@ -134,7 +134,9 @@ void DataLine::fillYbus(std::vector<Eigen::Triplet<cplx_type> > & res,
         res.push_back(Eigen::Triplet<cplx_type> (bus_ex_solver_id, bus_ex_solver_id, tmp));
     }
 }
-void DataLine::fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver)
+void DataLine::fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res,
+                              bool ac,
+                              const std::vector<int> & id_grid_to_solver)
 {
 
     //TODO this is no more used!!!! see the other fillYbus
@@ -217,9 +219,9 @@ void DataLine::reset_results()
 }
 
 
-void DataLine::compute_results(const Eigen::Ref<RealVect> & Va,
-                               const Eigen::Ref<RealVect> & Vm,
-                               const Eigen::Ref<CplxVect> & V,
+void DataLine::compute_results(const Eigen::Ref<const RealVect> & Va,
+                               const Eigen::Ref<const RealVect> & Vm,
+                               const Eigen::Ref<const CplxVect> & V,
                                const std::vector<int> & id_grid_to_solver,
                                const RealVect & bus_vn_kv,
                                real_type sn_mva)
diff --git a/src/DataLine.h b/src/DataLine.h
index f78c67b..a7758c1 100644
--- a/src/DataLine.h
+++ b/src/DataLine.h
@@ -176,9 +176,9 @@ class DataLine : public DataGeneric
                           );
     virtual void fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver);
 
-    void compute_results(const Eigen::Ref<RealVect> & Va,
-                         const Eigen::Ref<RealVect> & Vm,
-                         const Eigen::Ref<CplxVect> & V,
+    void compute_results(const Eigen::Ref<const RealVect> & Va,
+                         const Eigen::Ref<const RealVect> & Vm,
+                         const Eigen::Ref<const CplxVect> & V,
                          const std::vector<int> & id_grid_to_solver,
                          const RealVect & bus_vn_kv,
                          real_type sn_mva);
@@ -189,8 +189,8 @@ class DataLine : public DataGeneric
     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_);}
 
-    const RealVect & get_theta_or() const {return res_powerline_thetaor_;}
-    const RealVect & get_theta_ex() const {return 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_;}
 
     protected:
diff --git a/src/DataLoad.cpp b/src/DataLoad.cpp
index 38ddf98..1ffc6d0 100644
--- a/src/DataLoad.cpp
+++ b/src/DataLoad.cpp
@@ -70,9 +70,9 @@ void DataLoad::fillSbus(CplxVect & Sbus, bool ac, const std::vector<int> & id_gr
     }
 }
 
-void DataLoad::compute_results(const Eigen::Ref<RealVect> & Va,
-                               const Eigen::Ref<RealVect> & Vm,
-                               const Eigen::Ref<CplxVect> & V,
+void DataLoad::compute_results(const Eigen::Ref<const RealVect> & Va,
+                               const Eigen::Ref<const RealVect> & Vm,
+                               const Eigen::Ref<const CplxVect> & V,
                                const std::vector<int> & id_grid_to_solver,
                                const RealVect & bus_vn_kv,
                                real_type sn_mva)
diff --git a/src/DataLoad.h b/src/DataLoad.h
index ebb062d..fbf9449 100644
--- a/src/DataLoad.h
+++ b/src/DataLoad.h
@@ -67,9 +67,9 @@ class DataLoad : public DataGeneric
 
     virtual void fillSbus(CplxVect & Sbus, bool ac, const std::vector<int> & id_grid_to_solver);
 
-    void compute_results(const Eigen::Ref<RealVect> & Va,
-                         const Eigen::Ref<RealVect> & Vm,
-                         const Eigen::Ref<CplxVect> & V,
+    void compute_results(const Eigen::Ref<const RealVect> & Va,
+                         const Eigen::Ref<const RealVect> & Vm,
+                         const Eigen::Ref<const CplxVect> & V,
                          const std::vector<int> & id_grid_to_solver,
                          const RealVect & bus_vn_kv,
                          real_type sn_mva);
@@ -78,7 +78,7 @@ class DataLoad : public DataGeneric
     virtual void get_q(std::vector<real_type>& q_by_bus);
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
-    const RealVect & get_theta() const {return res_theta_;}
+    Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
     const std::vector<bool>& get_status() const {return status_;}
 
     protected:
diff --git a/src/DataSGen.cpp b/src/DataSGen.cpp
index e1426ea..2320d02 100644
--- a/src/DataSGen.cpp
+++ b/src/DataSGen.cpp
@@ -105,9 +105,9 @@ void DataSGen::fillSbus(CplxVect & Sbus, bool ac, const std::vector<int> & id_gr
     }
 }
 
-void DataSGen::compute_results(const Eigen::Ref<RealVect> & Va,
-                               const Eigen::Ref<RealVect> & Vm,
-                               const Eigen::Ref<CplxVect> & V,
+void DataSGen::compute_results(const Eigen::Ref<const RealVect> & Va,
+                               const Eigen::Ref<const RealVect> & Vm,
+                               const Eigen::Ref<const CplxVect> & V,
                                const std::vector<int> & id_grid_to_solver,
                                const RealVect & bus_vn_kv,
                                real_type sn_mva)
@@ -160,7 +160,7 @@ real_type DataSGen::get_p_slack(int slack_bus_id)
     for(int sgen_id = 0; sgen_id < nb_element; ++sgen_id)
     {
         if(!status_[sgen_id]) continue;
-        if(bus_id_(sgen_id) == slack_bus_id) res += res_p_(sgen_id);
+        if(bus_id_(sgen_id) == slack_bus_id) res -= res_p_(sgen_id);
     }
     return res;
 }
@@ -172,6 +172,6 @@ void DataSGen::get_q(std::vector<real_type>& q_by_bus)
     {
         if(!status_[sgen_id]) continue;
         int bus_id = bus_id_[sgen_id];
-        q_by_bus[bus_id] += res_q_(sgen_id); //TODO weird that i need to put a + here and a - for the active!
+        q_by_bus[bus_id] -= res_q_(sgen_id); //TODO weird that i need to put a + here and a - for the active!
     }
 }
diff --git a/src/DataSGen.h b/src/DataSGen.h
index 2c28825..f0c869b 100644
--- a/src/DataSGen.h
+++ b/src/DataSGen.h
@@ -72,9 +72,9 @@ class DataSGen: public DataGeneric
 
     virtual void fillSbus(CplxVect & Sbus, bool ac, const std::vector<int> & id_grid_to_solver);
 
-    void compute_results(const Eigen::Ref<RealVect> & Va,
-                         const Eigen::Ref<RealVect> & Vm,
-                         const Eigen::Ref<CplxVect> & V,
+    void compute_results(const Eigen::Ref<const RealVect> & Va,
+                         const Eigen::Ref<const RealVect> & Vm,
+                         const Eigen::Ref<const CplxVect> & V,
                          const std::vector<int> & id_grid_to_solver,
                          const RealVect & bus_vn_kv,
                          real_type sn_mva);
@@ -83,7 +83,7 @@ class DataSGen: public DataGeneric
     virtual void get_q(std::vector<real_type>& q_by_bus);
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
-    const RealVect & get_theta() const {return res_theta_;}
+    Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
     const std::vector<bool>& get_status() const {return status_;}
 
     protected:
diff --git a/src/DataShunt.cpp b/src/DataShunt.cpp
index 0dadd35..b95026b 100644
--- a/src/DataShunt.cpp
+++ b/src/DataShunt.cpp
@@ -116,9 +116,9 @@ void DataShunt::fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, co
     //TODO this is no more used!!!! see the other fillYbus
 }
 
-void DataShunt::compute_results(const Eigen::Ref<RealVect> & Va,
-                                const Eigen::Ref<RealVect> & Vm,
-                                const Eigen::Ref<CplxVect> & V,
+void DataShunt::compute_results(const Eigen::Ref<const RealVect> & Va,
+                                const Eigen::Ref<const RealVect> & Vm,
+                                const Eigen::Ref<const CplxVect> & V,
                                 const std::vector<int> & id_grid_to_solver,
                                 const RealVect & bus_vn_kv,
                                 real_type sn_mva)
diff --git a/src/DataShunt.h b/src/DataShunt.h
index 411e1a8..cc8794a 100644
--- a/src/DataShunt.h
+++ b/src/DataShunt.h
@@ -66,9 +66,9 @@ class DataShunt : public DataGeneric
     virtual void fillYbus_spmat(Eigen::SparseMatrix<cplx_type> & res, bool ac, const std::vector<int> & id_grid_to_solver);
     virtual void fillSbus(CplxVect & Sbus, bool ac, const std::vector<int> & id_grid_to_solver);  // in DC i need that
 
-    void compute_results(const Eigen::Ref<RealVect> & Va,
-                         const Eigen::Ref<RealVect> & Vm,
-                         const Eigen::Ref<CplxVect> & V,
+    void compute_results(const Eigen::Ref<const RealVect> & Va,
+                         const Eigen::Ref<const RealVect> & Vm,
+                         const Eigen::Ref<const CplxVect> & V,
                          const std::vector<int> & id_grid_to_solver,
                          const RealVect & bus_vn_kv,
                          real_type sn_mva);
@@ -77,7 +77,7 @@ class DataShunt : public DataGeneric
     virtual void get_q(std::vector<real_type>& q_by_bus);
 
     tuple3d get_res() const {return tuple3d(res_p_, res_q_, res_v_);}
-    const RealVect & get_theta() const {return res_theta_;}
+    Eigen::Ref<const RealVect> get_theta() const {return res_theta_;}
     const std::vector<bool>& get_status() const {return status_;}
 
     protected:
diff --git a/src/DataTrafo.cpp b/src/DataTrafo.cpp
index 56808d0..327f309 100644
--- a/src/DataTrafo.cpp
+++ b/src/DataTrafo.cpp
@@ -299,9 +299,9 @@ void DataTrafo::fillSbus(CplxVect & Sbus, bool ac, const std::vector<int> & id_g
     }
 }
 
-void DataTrafo::compute_results(const Eigen::Ref<RealVect> & Va,
-                                const Eigen::Ref<RealVect> & Vm,
-                                const Eigen::Ref<CplxVect> & V,
+void DataTrafo::compute_results(const Eigen::Ref<const RealVect> & Va,
+                                const Eigen::Ref<const RealVect> & Vm,
+                                const Eigen::Ref<const CplxVect> & V,
                                 const std::vector<int> & id_grid_to_solver,
                                 const RealVect & bus_vn_kv,
                                 real_type sn_mva
diff --git a/src/DataTrafo.h b/src/DataTrafo.h
index 33219b9..511b795 100644
--- a/src/DataTrafo.h
+++ b/src/DataTrafo.h
@@ -182,9 +182,9 @@ class DataTrafo : public DataGeneric
                           real_type sn_mva);
     virtual void fillSbus(CplxVect & Sbus, bool ac, const std::vector<int> & id_grid_to_solver);  // needed for dc mode
 
-    void compute_results(const Eigen::Ref<RealVect> & Va,
-                         const Eigen::Ref<RealVect> & Vm,
-                         const Eigen::Ref<CplxVect> & V,
+    void compute_results(const Eigen::Ref<const RealVect> & Va,
+                         const Eigen::Ref<const RealVect> & Vm,
+                         const Eigen::Ref<const CplxVect> & V,
                          const std::vector<int> & id_grid_to_solver,
                          const RealVect & bus_vn_kv,
                          real_type sn_mva);
@@ -194,8 +194,8 @@ class DataTrafo : public DataGeneric
 
     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_);}
-    const RealVect & get_theta_hv() const {return res_theta_hv_;}
-    const RealVect & get_theta_lv() const {return 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_;}
 
     const std::vector<bool>& get_status() const {return status_;}
 
diff --git a/src/GridModel.cpp b/src/GridModel.cpp
index 28a442f..495aa0b 100644
--- a/src/GridModel.cpp
+++ b/src/GridModel.cpp
@@ -21,7 +21,7 @@ GridModel::GridModel(const GridModel & other)
 
     // copy the powersystem representation
     // 1. bus
-    bus_vn_kv_.array() = other.bus_vn_kv_;
+    bus_vn_kv_ = other.bus_vn_kv_;
     bus_status_ = other.bus_status_;
 
     // 2. powerline
@@ -466,6 +466,7 @@ void GridModel::compute_results(){
     const auto & Va = _solver.get_Va();
     const auto & Vm = _solver.get_Vm();
     const auto & V = _solver.get_V();
+
     // for powerlines
     powerlines_.compute_results(Va, Vm, V, id_me_to_solver_, bus_vn_kv_, sn_mva_);
     // for trafo
diff --git a/src/GridModel.h b/src/GridModel.h
index 2145513..005b2ca 100644
--- a/src/GridModel.h
+++ b/src/GridModel.h
@@ -278,43 +278,47 @@ class GridModel : public DataGeneric
         tuple3d get_storages_res() const {return storages_.get_res();}
         const std::vector<bool>& get_storages_status() const { return storages_.get_status();}
 
-        const RealVect & get_gen_theta() const  {return generators_.get_theta();}
-        const RealVect & get_load_theta() const  {return loads_.get_theta();}
-        const RealVect & get_shunt_theta() const  {return shunts_.get_theta();}
-        const RealVect & get_storage_theta() const  {return storages_.get_theta();}
-        const RealVect & get_lineor_theta() const {return powerlines_.get_theta_or();}
-        const RealVect & get_lineex_theta() const {return powerlines_.get_theta_ex();}
-        const RealVect & get_trafohv_theta() const {return trafos_.get_theta_hv();}
-        const RealVect & get_trafolv_theta() const {return trafos_.get_theta_lv();}
+        Eigen::Ref<const RealVect> get_gen_theta() const  {return generators_.get_theta();}
+        Eigen::Ref<const RealVect> get_load_theta() const  {return loads_.get_theta();}
+        Eigen::Ref<const RealVect> get_shunt_theta() const  {return shunts_.get_theta();}
+        Eigen::Ref<const RealVect> get_storage_theta() const  {return storages_.get_theta();}
+        Eigen::Ref<const RealVect> get_lineor_theta() const {return powerlines_.get_theta_or();}
+        Eigen::Ref<const RealVect> get_lineex_theta() const {return powerlines_.get_theta_ex();}
+        Eigen::Ref<const RealVect> get_trafohv_theta() const {return trafos_.get_theta_hv();}
+        Eigen::Ref<const RealVect> get_trafolv_theta() const {return trafos_.get_theta_lv();}
 
         // 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(){
-            return Ybus_;
+            return Ybus_;  // This is copied to python
         }
-        CplxVect get_Sbus(){
+        Eigen::Ref<CplxVect> get_Sbus(){
             return Sbus_;
         }
-        Eigen::VectorXi get_pv(){
+        Eigen::Ref<const Eigen::VectorXi> get_pv() const{
             return bus_pv_;
         }
-        Eigen::VectorXi get_pq(){
+        Eigen::Ref<const Eigen::VectorXi> get_pq() const{
             return bus_pq_;
         }
-        Eigen::Ref<RealVect> get_Va(){
+        Eigen::Ref<const CplxVect> get_V() const{
+            return _solver.get_V();
+        }
+        Eigen::Ref<const RealVect> get_Va() const{
             return _solver.get_Va();
         }
-        Eigen::Ref<RealVect> get_Vm(){
+        Eigen::Ref<const RealVect> get_Vm() const{
             return _solver.get_Vm();
         }
         Eigen::SparseMatrix<real_type> get_J(){
-            return _solver.get_J();
+            return _solver.get_J();  // This is copied to python
         }
         real_type get_computation_time(){ return _solver.get_computation_time();}
 
         // part dedicated to grid2op backend, optimized for grid2op data representation (for speed)
         // this is not recommended to use it outside of its intended usage within grid2op !
-        void update_bus_status(int nb_bus_before, Eigen::Ref<Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> > active_bus);
+        void update_bus_status(int nb_bus_before,
+                               Eigen::Ref<Eigen::Array<bool, Eigen::Dynamic, 2, Eigen::RowMajor> > active_bus);
         void update_gens_p(Eigen::Ref<Eigen::Array<bool, Eigen::Dynamic, Eigen::RowMajor> > has_changed,
                            Eigen::Ref<Eigen::Array<float, Eigen::Dynamic, Eigen::RowMajor> > new_values);
         void update_gens_v(Eigen::Ref<Eigen::Array<bool, Eigen::Dynamic, Eigen::RowMajor> > has_changed,
diff --git a/src/NICSLUSolver.cpp b/src/NICSLUSolver.cpp
new file mode 100644
index 0000000..7b27a3e
--- /dev/null
+++ b/src/NICSLUSolver.cpp
@@ -0,0 +1,133 @@
+// Copyright (c) 2020, 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.
+
+#include "NICSLUSolver.h"
+
+void NICSLUSolver::reset(){
+    BaseNRSolver::reset();
+
+    // free everything
+    n_ = -1;
+    solver_.Free();
+    delete [] ai_;  // created in NICSLUSolver::initialize
+    delete [] ap_;  // created in NICSLUSolver::initialize
+    ai_ = nullptr;
+    ap_ = nullptr;
+
+    solver_ = CNicsLU();
+    int ret = solver_.Initialize();
+    if (ret < 0){   // fail
+        if (ret == -100){
+            // fail to initialize because of a license issue
+            std::string msg = "Fail to initilize the NICSLU solver because we cannot find the nicslu.lic file.";
+            msg += "Please copy this file at the location you want to use the NICSLU solver.";
+            std::cout << msg << std::endl;
+        }
+        std::cout << "NICSLUSolver::reset() solver_.Initialize error: "  << ret << std::endl;
+        err_ = 1;  // TODO clarify that
+        return;
+    }
+    // solver_.SetConfiguration(0, 1.); //enable timer
+}
+
+void NICSLUSolver::initialize(){
+    // default Eigen representation: column major, which is good for klu !
+    auto timer = CustTimer();
+    int ret;
+    n_ = J_.cols(); // should be equal to J_.nrows()
+    err_ = 0; // reset error message
+    //    symbolic_ = klu_analyze(n_, J_.outerIndexPtr(), J_.innerIndexPtr(), &common_);
+    //    numeric_ = klu_factor(J_.outerIndexPtr(), J_.innerIndexPtr(), J_.valuePtr(), symbolic_, &common_);
+
+    // solver.Analyze(n, ax, ai, ap, MATRIX_ROW_REAL);
+    // convert int * (storage class for Eign) to unsigned int * (expected storage class for NICSLU...)
+
+    // std::cout << "n_ = " << n_ << std::endl;  // DEBUG
+    const unsigned int nnz = J_.nonZeros();  // DEBUG
+    // std::cout << "nnz = " << nnz << std::endl;  // DEBUG
+
+    ai_ = new unsigned int [nnz];  // deleted in destructor and NICSLUSolver::reset
+    const int * ref_ai = J_.innerIndexPtr();
+    for(unsigned int i = 0; i < nnz; ++i){
+        ai_[i] = static_cast<unsigned int>(ref_ai[i]);
+        // std::cout << "ai_["<<i<<"] = " << ai_[i] << std::endl;  // DEBUG
+    }
+    // std::cout << std::endl;  // DEBUG
+
+    ap_ = new unsigned int [n_+1];  // deleted in destructor and NICSLUSolver::reset
+    const int * ref_ap = J_.outerIndexPtr();
+    for(int i = 0; i < n_+1; ++i){
+        ap_[i] = static_cast<unsigned int>(ref_ap[i]);
+        // std::cout << "ap_["<<i<<"] = " << ap_[i] << std::endl;  // DEBUG
+    }
+    // std::cout << std::endl;  // DEBUG
+
+    // const double * ref_ax = J_.valuePtr();  // DEBUG
+    // for(unsigned int i = 0; i < nnz; ++i){  // DEBUG
+    //     std::cout << "ax["<<i<<"] = " << ref_ax[i] << std::endl;  // DEBUG
+    // }  // DEBUG
+
+    ret = solver_.Analyze(n_,
+                          J_.valuePtr(),
+                          ai_,
+                          ap_,
+                          MATRIX_COLUMN_REAL);  // MATRIX_COLUMN_REAL, MATRIX_ROW_REAL
+    if (ret < 0){
+        err_ = 1;
+        std::cout << "NICSLUSolver::initialize() solver_.Analyze error: "  << ret << std::endl;
+        // https://github.com/chenxm1986/nicslu/blob/master/nicslu202103/include/nicslu.h for error info
+        return;
+    }
+    // solver.FactorizeMatrix(ax, 0); //use all created threads
+    ret = solver_.FactorizeMatrix(J_.valuePtr(), nb_thread_);
+    if (ret < 0){
+        err_ = 1;
+        std::cout << "NICSLUSolver::initialize() solver_.FactorizeMatrix error: "  << ret << std::endl;
+        // https://github.com/chenxm1986/nicslu/blob/master/nicslu202103/include/nicslu.h for error info
+        return;
+    }
+    need_factorize_ = false;
+    timer_initialize_ += timer.duration();
+}
+
+void NICSLUSolver::solve(RealVect & b, bool has_just_been_inialized){
+    // solves (for x) the linear system J.x = b
+    // supposes that the solver has been initialized (call klu_solver.analyze() before calling that)
+    // J is const even if it does not compile if said const
+    auto timer = CustTimer();
+    int ret;
+    bool stop = false;
+    RealVect x;
+    if(!has_just_been_inialized){
+        // if the call to "klu_factor" has been made this iteration, there is no need
+        // to re factor again the matrix
+        // i'm in the case where it has not
+        //        ret = klu_refactor(J_.outerIndexPtr(), J_.innerIndexPtr(), J_.valuePtr(), symbolic_, numeric_, &common_);
+
+        // solver.FactorizeMatrix(ax, 0); //use all created threads
+        ret  = solver_.FactorizeMatrix(J_.valuePtr(), nb_thread_);
+        if (ret < 0) {
+            std::cout << "NICSLUSolver::solve solver_.FactorizeMatrix error: " << ret << std::endl;
+            // https://github.com/chenxm1986/nicslu/blob/master/nicslu202103/include/nicslu.h for error info
+            err_ = 2;
+            stop = true;
+        }
+    }
+    if(!stop){
+        // ok = klu_solve(symbolic_, numeric_, n_, 1, &b(0), &common_);
+        // solver.Solve(b, x);
+        x = RealVect(n_);
+        ret = solver_.Solve(&b(0), &x(0));
+        if (ret < 0) {
+            std::cout << "NICSLUSolver::solve solver_.Solve error: " << ret << std::endl;
+            err_ = 3;
+        }
+        b = x;
+    }
+    timer_solve_ += timer.duration();
+}
diff --git a/src/NICSLUSolver.h b/src/NICSLUSolver.h
new file mode 100644
index 0000000..6c5a057
--- /dev/null
+++ b/src/NICSLUSolver.h
@@ -0,0 +1,83 @@
+// Copyright (c) 2020, 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.
+
+#ifdef NICSLU_SOLVER_AVAILABLE
+#ifndef NICSLUSOLVER_H
+#define NICSLUSOLVER_H
+
+#include <iostream>
+#include <vector>
+#include <stdio.h>
+#include <cstdint> // for int32
+#include <chrono>
+#include <cmath>  // for PI
+
+// eigen is necessary to easily pass data from numpy to c++ without any copy.
+// and to optimize the matrix operations
+#include "Utils.h"
+#include "Eigen/Core"
+#include "Eigen/Dense"
+#include "Eigen/SparseCore"
+
+// import klu package
+#include "nicslu_cpp.inl"
+
+#include "CustTimer.h"
+#include "BaseNRSolver.h"
+/**
+class to handle the solver using newton-raphson method, using NICSLU algorithm and sparse matrices.
+
+As long as the admittance matrix of the system does not change, you can reuse the same solver.
+Reusing the same solver is possible, but "reset" method must be called.
+
+Otherwise, unexpected behaviour might follow, including "segfault".
+
+NB: the code of NICSLU is not included in this repository. This class is only compiled if the "setup.py"
+can find a version of `https://github.com/chenxm1986/nicslu`. Be careful though, this code is under some
+specific license.
+
+**/
+
+// TODO use the cpp API instead !
+class NICSLUSolver: public BaseNRSolver
+{
+    public:
+        NICSLUSolver():BaseNRSolver(), solver_(), nb_thread_(1){
+        }
+
+        ~NICSLUSolver()
+         {
+             solver_.Free();
+             delete [] ai_;
+             delete [] ap_;
+         }
+
+        virtual void reset();
+
+    protected:
+        virtual
+        void initialize();
+
+        virtual
+        void solve(RealVect & b, bool has_just_been_inialized);
+
+    private:
+        // solver initialization
+        CNicsLU solver_;
+        const unsigned int nb_thread_;
+        unsigned int * ai_;
+        unsigned int * ap_;
+
+        // no copy allowed
+        NICSLUSolver( const NICSLUSolver & ) ;
+        NICSLUSolver & operator=( const NICSLUSolver & ) ;
+
+};
+
+#endif // NICSLUSOLVER_H
+#endif  // NICSLU_SOLVER_AVAILABLE
diff --git a/src/Utils.h b/src/Utils.h
index 1d3c185..a6d9a7c 100644
--- a/src/Utils.h
+++ b/src/Utils.h
@@ -22,8 +22,13 @@ 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 std::tuple<EigenPythonNumType, EigenPythonNumType, EigenPythonNumType> tuple3d;
-typedef std::tuple<EigenPythonNumType, EigenPythonNumType, EigenPythonNumType, EigenPythonNumType> tuple4d;
+typedef std::tuple<Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType> > tuple3d;
+typedef std::tuple<Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType>,
+                   Eigen::Ref<const EigenPythonNumType> > tuple4d;
 typedef Eigen::Matrix<real_type, Eigen::Dynamic, 1> RealVect;
 typedef Eigen::Matrix<cplx_type, Eigen::Dynamic, 1> CplxVect;
 
diff --git a/src/main.cpp b/src/main.cpp
index 5d36eaf..d3ad1fe 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -24,8 +24,9 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .value("SparseLU", SolverType::SparseLU)
         .value("KLU", SolverType::KLU)
         .value("GaussSeidel", SolverType::GaussSeidel)
-        .value("DC", SolverType::DC)
         .value("GaussSeidelSynch", SolverType::GaussSeidelSynch)
+        .value("DC", SolverType::DC)
+        .value("NICSLU", SolverType::NICSLU)
         .export_values();
 
     #ifdef KLU_SOLVER_AVAILABLE
@@ -41,7 +42,22 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("compute_pf", &KLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>())  // perform the newton raphson optimization
         .def("get_timers", &KLUSolver::get_timers)  // returns the timers corresponding to times the solver spent in different part
         .def("solve", &KLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>() );  // perform the newton raphson optimization
-    #endif
+    #endif  // KLU_SOLVER_AVAILABLE
+
+    #ifdef NICSLU_SOLVER_AVAILABLE
+    py::class_<NICSLUSolver>(m, "NICSLUSolver")
+        .def(py::init<>())
+        .def("get_J", &NICSLUSolver::get_J)  // (get the jacobian matrix, sparse csc matrix)
+        .def("get_Va", &NICSLUSolver::get_Va)  // get the voltage angle vector (vector of double)
+        .def("get_Vm", &NICSLUSolver::get_Vm)  // get the voltage magnitude vector (vector of double)
+        .def("get_error", &NICSLUSolver::get_error)  // get the error message, see the definition of "err_" for more information
+        .def("get_nb_iter", &NICSLUSolver::get_nb_iter)  // return the number of iteration performed at the last optimization
+        .def("reset", &NICSLUSolver::reset)  // reset the solver to its original state
+        .def("converged", &NICSLUSolver::converged)  // whether the solver has converged
+        .def("compute_pf", &NICSLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>())  // perform the newton raphson optimization
+        .def("get_timers", &NICSLUSolver::get_timers)  // returns the timers corresponding to times the solver spent in different part
+        .def("solve", &NICSLUSolver::compute_pf, py::call_guard<py::gil_scoped_release>() );  // perform the newton raphson optimization
+    #endif  // NICSLU_SOLVER_AVAILABLE
 
     py::class_<SparseLUSolver>(m, "SparseLUSolver")
         .def(py::init<>())
@@ -285,6 +301,7 @@ PYBIND11_MODULE(lightsim2grid_cpp, m)
         .def("change_q_storage", &GridModel::change_q_storage)
 
         // get back the results
+        .def("get_V", &GridModel::get_V)
         .def("get_Va", &GridModel::get_Va)
         .def("get_Vm", &GridModel::get_Vm)
         .def("get_J", &GridModel::get_J)


Test patch:
diff --git a/lightsim2grid/tests/test_Checker.py b/lightsim2grid/tests/test_Checker.py
index 6a257e0..5a6d458 100644
--- a/lightsim2grid/tests/test_Checker.py
+++ b/lightsim2grid/tests/test_Checker.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 os
 import pandas as pd
 import unittest
diff --git a/lightsim2grid/tests/test_DataConverter.py b/lightsim2grid/tests/test_DataConverter.py
index 7e72b09..790b963 100644
--- a/lightsim2grid/tests/test_DataConverter.py
+++ b/lightsim2grid/tests/test_DataConverter.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 unittest
 import numpy as np
 import pdb
diff --git a/lightsim2grid/tests/test_Deactivatedbus.py b/lightsim2grid/tests/test_Deactivatedbus.py
index 8340dbb..59623e7 100644
--- a/lightsim2grid/tests/test_Deactivatedbus.py
+++ b/lightsim2grid/tests/test_Deactivatedbus.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 os
 import unittest
 import copy
diff --git a/lightsim2grid/tests/test_GaussSeidelSolver.py b/lightsim2grid/tests/test_GaussSeidelSolver.py
index 30018c8..e5d77e4 100644
--- a/lightsim2grid/tests/test_GaussSeidelSolver.py
+++ b/lightsim2grid/tests/test_GaussSeidelSolver.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 os
 import unittest
 import numpy as np
diff --git a/lightsim2grid/tests/test_GridModel.py b/lightsim2grid/tests/test_GridModel.py
index 206d368..e8b27a2 100644
--- a/lightsim2grid/tests/test_GridModel.py
+++ b/lightsim2grid/tests/test_GridModel.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 unittest
 import numpy as np
 import pandapower.networks as pn
diff --git a/lightsim2grid/tests/test_Iterators.py b/lightsim2grid/tests/test_Iterators.py
index 5bd5f8d..64932a2 100644
--- a/lightsim2grid/tests/test_Iterators.py
+++ b/lightsim2grid/tests/test_Iterators.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 os
 import unittest
 import grid2op
diff --git a/lightsim2grid/tests/test_KLUBackendAction.py b/lightsim2grid/tests/test_KLUBackendAction.py
index 94a4f26..0e25188 100644
--- a/lightsim2grid/tests/test_KLUBackendAction.py
+++ b/lightsim2grid/tests/test_KLUBackendAction.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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.
+
 from grid2op import make
 from grid2op.Agent import DoNothingAgent
 from grid2op.Parameters import Parameters
diff --git a/lightsim2grid/tests/test_KLUBackendDN.py b/lightsim2grid/tests/test_KLUBackendDN.py
index 414e824..ef8af27 100644
--- a/lightsim2grid/tests/test_KLUBackendDN.py
+++ b/lightsim2grid/tests/test_KLUBackendDN.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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.
+
 from grid2op import make
 from grid2op.Agent import DoNothingAgent
 from grid2op.Parameters import Parameters
diff --git a/lightsim2grid/tests/test_KLUBackendTestAgent.py b/lightsim2grid/tests/test_KLUBackendTestAgent.py
index bff5689..de0fcc6 100644
--- a/lightsim2grid/tests/test_KLUBackendTestAgent.py
+++ b/lightsim2grid/tests/test_KLUBackendTestAgent.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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.
+
 from grid2op import make
 from grid2op.Agent import AgentWithConverter
 from grid2op.Parameters import Parameters
@@ -222,6 +230,10 @@ def test_do_nothing(self):
             self.tol = 2e-4
         elif env_name == "l2rpn_case14_sandbox":
             self.tol = 2e-4
+        elif env_name == "rte_case5_example":
+            self.tol = 2e-4
+        elif env_name == "l2rpn_neurips_2020_track1":
+            self.tol = 4e-5
 
         assert nb_ts_klu == nb_ts_pp, "not same number of timesteps for {}: lightsim: {}, pp: {}" \
                                       "".format(env_name, nb_ts_klu, nb_ts_pp)
diff --git a/lightsim2grid/tests/test_KLUSolver.py b/lightsim2grid/tests/test_KLUSolver.py
index d71529e..67ec943 100644
--- a/lightsim2grid/tests/test_KLUSolver.py
+++ b/lightsim2grid/tests/test_KLUSolver.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 os
 import unittest
 import numpy as np
@@ -55,8 +63,7 @@ def load_path(self, path):
     def _load_vect(self, myzip, iter_max, nm):
         try:
             res = self.load_array(myzip, "{}_{}.npy".format(nm, iter_max))
-        except:
-            pdb.set_trace()
+        except Exception as exc_:
             raise RuntimeError("{} is not defined for iteration {} for case {}".format(nm, iter_max, self.path))
         return res
 
diff --git a/lightsim2grid/tests/test_NICSLUSolver.py b/lightsim2grid/tests/test_NICSLUSolver.py
new file mode 100644
index 0000000..88a9408
--- /dev/null
+++ b/lightsim2grid/tests/test_NICSLUSolver.py
@@ -0,0 +1,124 @@
+# Copyright (c) 2020, 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 os
+import unittest
+import numpy as np
+import pdb
+import zipfile
+from scipy import sparse
+import re
+NICSLU_AVAILBLE = False
+try:
+    from lightsim2grid_cpp import NICSLUSolver
+    NICSLU_AVAILBLE = True
+except ImportError:
+    # KLU solver is not available, these tests cannot be carried out
+    pass
+
+
+class MakeTests(unittest.TestCase):
+    def __init__(self, methodName='runTest'):
+        unittest.TestCase.__init__(self, methodName=methodName)
+        self.methodName = methodName
+
+        self.max_it = 10
+        self.tol = 1e-8  # tolerance for the solver
+        self.tol_test = 1e-4  # tolerance for the test (2 matrices are equal if the l_1 of their difference is less than this)
+        if not NICSLU_AVAILBLE:
+            return
+        self.solver = NICSLUSolver()
+
+        self.path = None
+        self.V_init = None
+        self.pq = None
+        self.pv = None
+        self.Sbus = None
+        self.Ybus = None
+
+    def load_array(self, myzip, nm):
+        arr = myzip.extract(nm)
+        res = np.load(arr)
+        os.remove(arr)
+        return res
+
+    def load_path(self, path):
+        try:
+            self.path = path
+            with zipfile.ZipFile(self.path) as myzip:
+                self.V_init = self.load_array(myzip, "V0.npy")
+                self.pq = self.load_array(myzip, "pq.npy")
+                self.pv = self.load_array(myzip, "pv.npy")
+                self.Sbus = self.load_array(myzip, "Sbus.npy")
+                self.Ybus = self.load_array(myzip, "Ybus.npy")
+                self.Ybus = sparse.csc_matrix(self.Ybus)
+            return True
+        except Exception as exc_:
+            return False
+
+    def _load_vect(self, myzip, iter_max, nm):
+        try:
+            res = self.load_array(myzip, "{}_{}.npy".format(nm, iter_max))
+        except Exception as exc_:
+            raise RuntimeError("{} is not defined for iteration {} for case {}".format(nm, iter_max, self.path))
+        return res
+
+    def load_res(self, iter_max):
+        with zipfile.ZipFile(self.path) as myzip:
+            J = self._load_vect(myzip, iter_max, "J")
+            V = self._load_vect(myzip, iter_max, "V")
+        return J, V
+
+    def compare_sparse_mat(self, J, J_pp, pv, pq):
+        """
+        Test that the matrices J and J_pp are almost equal
+        :param J:
+        :param J_pp:
+        :param pv:
+        :param pq:
+        :return:
+        """
+        pvpq = np.r_[pv, pq]
+
+        comp_val = np.abs(J - J_pp)
+        comp_val = comp_val
+        assert np.sum(np.abs(comp_val[:len(pvpq), :len(pvpq)])) <= self.tol_test, "J11 (dS_dVa_r) are not equal"
+        assert np.sum(np.abs(comp_val[len(pvpq):, :len(pvpq)])) <= self.tol_test, "J21 (dS_dVa_i) are not equal"
+        assert np.sum(np.abs(comp_val[:len(pvpq), len(pvpq):])) <= self.tol_test, "J12 (dS_dVm_r) are not equal"
+        assert np.sum(np.abs(comp_val[len(pvpq):, len(pvpq):])) <= self.tol_test, "J22 (dS_dVm_i) are not equal"
+
+    def solver_aux(self):
+        self.solver.reset()
+        has_conv = self.solver.compute_pf(self.Ybus, self.V_init, self.Sbus, self.pv, self.pq, self.max_it, self.tol)
+        assert has_conv, "the load flow has diverged for {}".format(self.path)
+        J = self.solver.get_J()
+        Va = self.solver.get_Va()
+        Vm = self.solver.get_Vm()
+        J_pp, V_pp = self.load_res(iter_max=self.solver.get_nb_iter())
+        self.compare_sparse_mat(J, J_pp, self.pv, self.pq)
+        Va_pp = np.angle(V_pp)
+        Vm_pp = np.abs(V_pp)
+        assert np.sum(np.abs(Va - Va_pp)) <= self.tol_test, "voltages angles are not the same"
+        assert np.sum(np.abs(Vm - Vm_pp)) <= self.tol_test, "voltages magnitude are not the same"
+
+    def test_dir(self):
+        if not NICSLU_AVAILBLE:
+            self.skipTest("NICSLU is not installed")
+        nb_tested = 0
+        for path in os.listdir("."):
+            _, ext = os.path.splitext(path)
+            if ext == ".zip":
+                path_ok = self.load_path(path)
+                if path_ok:
+                    self.solver_aux()
+                    nb_tested += 1
+        assert nb_tested == 5, "incorrect number of test cases found, found {} while there should be 5".format(nb_tested)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/lightsim2grid/tests/test_NewtonPF.py b/lightsim2grid/tests/test_NewtonPF.py
index 3b5dd8c..e49b9c8 100644
--- a/lightsim2grid/tests/test_NewtonPF.py
+++ b/lightsim2grid/tests/test_NewtonPF.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 os
 import unittest
 import numpy as np
diff --git a/lightsim2grid/tests/test_SameResPP.py b/lightsim2grid/tests/test_SameResPP.py
index 6313457..e81afe5 100644
--- a/lightsim2grid/tests/test_SameResPP.py
+++ b/lightsim2grid/tests/test_SameResPP.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 unittest
 import tempfile
 import pandapower as pp
@@ -25,6 +33,7 @@ def setUp(self) -> None:
 
     def test_case14(self):
         case = pn.case14()
+        self.tol = 2e-3
         self._aux_test(case)
 
     def test_case39(self):
@@ -37,6 +46,7 @@ def test_case118(self):
 
     def test_case1888rte(self):
         case = pn.case1888rte()
+        self.tol = 3e-4
         self._aux_test(case)
 
     # def test_case300(self):
@@ -266,7 +276,7 @@ def _aux_test(self, pn_net):
         v_or_me2 = np.concatenate((lvor, tvor))
         a_or_me2 = 1000. * np.concatenate((laor, taor))
         test_ok = True
-        # pdb.set_trace()
+
         max_mis = np.max(np.abs(p_or_me2 - por_pp))
         assert np.max(error_q) <= self.tol, f"\t Error: por do not match, maximum absolute error is {max_mis:.5f} MW"
         max_mis = np.max(np.abs(q_or_me2 - qor_pp))
diff --git a/lightsim2grid/tests/test_SparseLUSolver.py b/lightsim2grid/tests/test_SparseLUSolver.py
index 9243dbe..ec4b101 100644
--- a/lightsim2grid/tests/test_SparseLUSolver.py
+++ b/lightsim2grid/tests/test_SparseLUSolver.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 os
 import unittest
 import numpy as np
diff --git a/lightsim2grid/tests/test_issue_30.py b/lightsim2grid/tests/test_issue_30.py
new file mode 100644
index 0000000..9c6d178
--- /dev/null
+++ b/lightsim2grid/tests/test_issue_30.py
@@ -0,0 +1,61 @@
+# Copyright (c) 2020, 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 warnings
+import grid2op
+import unittest
+import numpy as np
+from lightsim2grid import LightSimBackend
+from grid2op.Action import CompleteAction
+import pdb
+
+
+class TestBug(unittest.TestCase):
+    def setUp(self) -> None:
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore")
+            self.env = grid2op.make("l2rpn_case14_sandbox",
+                                    test=True,
+                                    action_class=CompleteAction,
+                                    backend=LightSimBackend())
+
+    def test_diverge(self):
+        obs = self.env.reset()
+
+        # prove that it should converged
+        sim_o, sim_r, sim_d, sim_i = obs.simulate(self.env.action_space())
+        assert not sim_d, "should have converged"
+
+        # put the grid in a terrible spot
+        act10 = self.env.action_space({"set_bus": {"lines_or_id": [(2, 1), (3, 2), (4, 1)],
+                                                   "lines_ex_id": [(0, 2)],
+                                                   "loads_id": [(0, 2)],
+                                                   "generators_id": [(0, 1)]}})
+        obs, reward, done, info = self.env.step(act10)
+        assert not done, "should have converged here"
+
+        # i set the self.V (used to anyway) to a super low value with a grid
+        # that barely converges (in the bug it was with a set of actions carefully chosen)
+        obs._obs_env.backend.V[:] = np.array([9.66548436e-01-3.55087542e-01j, 9.39831081e-01-4.20738137e-01j,
+                                              1.01642368e+00+1.64881592e-01j, 9.92535287e-01+2.97637335e-02j,
+                                              7.32930199e-01-3.22013532e-01j, 1.00021220e+00-4.57794282e-01j,
+                                              6.07088708e-18-2.53076324e-18j, 1.01531177e+00-4.23251771e-01j,
+                                              1.00220145e-01-4.02978083e-02j, 9.17213244e-01-3.85631892e-01j,
+                                              9.53039556e-01-4.21909511e-01j, 9.06932626e-01-4.29640029e-01j,
+                                              8.31803137e-01-3.93133296e-01j, 3.52206227e-01-1.62093322e-01j,
+                                              1.05583333e+00+0.00000000e+00j, 8.67128154e-01-3.45062851e-01j,
+                                              8.81047620e-01-5.81841146e-01j, 6.26448055e-01-2.78619470e-01j,
+                                              1.05583333e+00+0.00000000e+00j, 1.05583333e+00+0.00000000e+00j,
+                                              1.05583333e+00+0.00000000e+00j, 1.05583333e+00+0.00000000e+00j,
+                                              1.05583333e+00+0.00000000e+00j, 1.05583333e+00+0.00000000e+00j,
+                                              1.05583333e+00+0.00000000e+00j, 1.05583333e+00+0.00000000e+00j,
+                                              1.05583333e+00+0.00000000e+00j, 1.05583333e+00+0.00000000e+00j])
+
+        # and it used to make this diverge
+        sim_o, sim_r, sim_d, sim_i = obs.simulate(self.env.action_space())
+        assert not sim_d, "should have converged"
diff --git a/lightsim2grid/tests/test_pickleable.py b/lightsim2grid/tests/test_pickleable.py
index d728925..e736ff8 100644
--- a/lightsim2grid/tests/test_pickleable.py
+++ b/lightsim2grid/tests/test_pickleable.py
@@ -1,3 +1,11 @@
+# Copyright (c) 2020, 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 tempfile
 import pickle
 import os