diff --git a/src/specmicp_common/io/safe_config.hpp b/src/specmicp_common/io/safe_config.hpp
index 0189037..028f5f3 100644
--- a/src/specmicp_common/io/safe_config.hpp
+++ b/src/specmicp_common/io/safe_config.hpp
@@ -1,572 +1,572 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #ifndef SPECMICP_UTILS_SAFECONFIG_HPP
 #define SPECMICP_UTILS_SAFECONFIG_HPP
 
 //! \file safe_config.hpp
 //! \brief YAML configuration file reader
 //!
 //! This is a set of wrappers over the yaml-cpp API
 
 #include "specmicp_common/types.hpp"
 
 #include "specmicp_common/pimpl_ptr.hpp"
 
 #include <string>
 #include <memory>
 
 #include <yaml-cpp/node/detail/iterator_fwd.h>
 
 namespace YAML {
     class Node;
 } //end namespace YAML
 
 namespace specmicp {
 namespace io {
 
 //! \brief Type of error that can be reported
 enum class YAMLConfigError
 {
     UnknownError,
     UnknownVariable,
     MissingRequiredAttribute,
     MissingRequiredSection,
     InvalidArgument,
     ConversionError,
     ListExpected,
     MapExpected
 };
 
 // internal store required information
 //! \internal
 struct SPECMICP_DLL_LOCAL YAMLConfigFileHandle;
 
 //! \brief Access to the configuration
 class SPECMICP_DLL_PUBLIC YAMLConfigHandle
 {
 public:
     //! \brief Report an error
     void report_error(
             YAMLConfigError error_type,
             const std::string& error_msg
             );
     //! \brief Return true if it is a sequence
     bool is_sequence();
     //! \brief Return true if sub-node exists and is a sequence
     bool is_sequence(const std::string& node);
     //! \brief Return true if sub-node exists and is a map
     bool is_map(const std::string& node);
     //! \brief Return true if it is a map
     bool is_map();
 
 
     //! \brief Return true if sub-node exist
     bool has_node(const std::string& node);
     //! \brief Return true if sub-node exist
     bool has_node(uindex_t value);
 
     //! \brief Return true if the node has the given attribute
     //!
     //! An attribute is a pair key/value, where the value is
     //! neither a map or a sequence
     bool has_attribute(const std::string& attribute);
 
     //! \brief Return true if the node has the given section
     //!
     //! A section is a pair key/value, where the value is a
     //! map or a sequence
     bool has_section(const std::string& section);
     //! \brief Return true if the node has the given section
     //!
     //! A section is a pair key/value, where the value is a
     //! map or a sequence
     bool has_section(uindex_t value);
 
     //! \brief Return the number of element in the section
     uindex_t size();
 
     //! \brief Return a section from a map
     YAMLConfigHandle get_section(const std::string& section);
     //! \brief Return a section from a list
     YAMLConfigHandle get_section(uindex_t value);
 
     //! \brief Return an attribute
     template <typename T>
     T get_attribute(const std::string& attribute);
 
     //! \brief Return a lower bounded scalar attribute
     template <typename T>
     T get_attribute(const std::string& attribute,
                     T lower_bound
                     );
     //! \brief Return a bounded scalar attribute
     template <typename T>
     T get_attribute(const std::string& attribute,
                     T lower_bound,
                     T upper_bound
                     );
     //! \brief Return a value in a list
     template <typename T>
     T get_value(uindex_t ind);
 
 
     //! \brief Return a required attribute
     template <typename T>
     T get_required_attribute(const std::string& attribute);
     //! \brief Return a required attribute
     template <typename T>
     T get_required_attribute(
             const std::string& attribute,
             T lower_bound
             );
     //! \brief Return a required attribute
     template <typename T>
     T get_required_attribute(
             const std::string& attribute,
             T lower_bound,
             T upper_bound
             );
 
     //! \brief Return an optional attribute
     template <typename T>
     T get_optional_attribute(
             const std::string& attribute,
             const T& default_value
             );
     //! \brief Return an optional attribute
     template <typename T>
     T get_optional_attribute(
             const std::string& attribute,
             const T& default_value,
             T lower_bound
             );
     //! \brief Return an optional attribute
     template <typename T>
     T get_optional_attribute(
             const std::string& attribute,
             const T& default_value,
             T lower_bound,
             T upper_bound
             );
 
 
     //! \brief Set a variable if the attribute exists
     template <typename T>
     void set_if_attribute_exists(
             T& var_to_set,
             const std::string& attribute
             );
     //! \brief Set a variable if the attribut exists, with bound checking
     template <typename T>
     void set_if_attribute_exists(
             T& var_to_set,
             const std::string& attribute,
             T lower_bound
             );
     //! \brief Set a variable if the attribut exists, with bound checking
     template <typename T>
     void set_if_attribute_exists(
             T& var_to_set,
             const std::string& attribute,
             T lower_bound,
             T upper_bound
             );
     //! \brief Set a variable if the attribut exists
     template <typename T=index_t>
     void set_if_attribute_exists(
             int& var_to_set,
             const std::string& attribute
             );
     //! \brief Set a variable if the attribut exists, with bound checking
     template <typename T=index_t>
     void set_if_attribute_exists(
             int& var_to_set,
             const std::string& attribute,
             T lower_bound
             );
 
 
     //! \brief Return a list of values for a list attribute
     template <typename T>
     std::vector<T> list_to_vector(const std::string& list);
 
     ~YAMLConfigHandle();
 
     //! \brief Copy operator
     YAMLConfigHandle& operator= (const YAMLConfigHandle& other);
     //! \brief Copy constructor
     YAMLConfigHandle(const YAMLConfigHandle& other);
     //! \brief Mobe constructor
     YAMLConfigHandle(YAMLConfigHandle&& other);
 
 
 
 protected:
     //! \brief Protectect constructor
     // can only be created from a file or another handle
     YAMLConfigHandle(
             const YAML::Node& node,
             std::shared_ptr<YAMLConfigFileHandle> const file_handle,
             const std::string& section,
             const std::string& path);
 
     //! \brief Set the file handle to manage the configuration
     void set_file_handle(std::shared_ptr<YAMLConfigFileHandle> file_handle);
 
 private:
     struct SPECMICP_DLL_LOCAL YAMLConfigHandleImpl;
     utils::pimpl_ptr<YAMLConfigHandleImpl> m_impl; //!< implementation
 
 public:
     //! \brief Iterator over a YAML map
     class MapIterator
     {
     friend YAMLConfigHandle;
 
     public:
         //! \brief move constructor
         MapIterator(MapIterator&& other);
         ~MapIterator();
 
         //! \brief Return <key, value>
         std::pair<std::string, std::string> operator* ();
 
         template <typename T>
         std::pair<std::string, T> get_as();
 
         //! \brief increment the operator
         MapIterator& operator++ ();
 
         //! \brief Equal operator
         bool operator==(const MapIterator& other);
         //! \brief Not equal operator
         bool operator!=(const MapIterator& other);
 
     private:
         //! \brief Private constructor
         MapIterator(
                 YAMLConfigHandleImpl* handle,
                 std::unique_ptr<YAML::iterator> it
                 );
 
         YAMLConfigHandleImpl* m_handle; //!< ptr to the YAML node
         std::unique_ptr<YAML::iterator> m_true_it; //!< The "true" yaml operator
 
     };
 
     MapIterator map_begin(); //!< Return an iterator to the beginning of the map
     MapIterator map_end();   //!< Return an iterator to the end of the maps
 };
 
 // ======================= //
 //      YAMLConfigFile     //
 // ======================= //
 
 //! \brief This class represent a file
 //!
 //! It must be alive while the config is accessed
 class SPECMICP_DLL_PUBLIC YAMLConfigFile:
         public YAMLConfigHandle
 {
 public:
     //! \brief Parse a file
     static YAMLConfigFile load(
             const std::string& file
             );
     //! \brief Parse a string
     static YAMLConfigFile load_from_string(
             const std::string& str,
             std::string name=""
             );
 
     //! \brief Pase a file, return a smart pointer
     static std::unique_ptr<YAMLConfigFile> make(
             const std::string& file
             );
     //! \brief Parse a file, return a smart pointer
     static std::unique_ptr<YAMLConfigFile> make_from_string(
             const std::string& str,
             std::string name=""
             );
 
     ~YAMLConfigFile();
 private:
     //! \brief Private Constructor
     YAMLConfigFile(
             const YAML::Node& node,
             const std::string name
             );
     std::shared_ptr<YAMLConfigFileHandle> m_handle; //!< Manager
 };
 
 
 // ======================= //
 // Template implementation //
 // ======================= //
 
 #ifndef SPC_DOXYGEN_SHOULD_SKIP_THIS
 // explicit specialization
 template <>
-scalar_t YAMLConfigHandle::get_attribute(const std::string& attribute);
+scalar_t YAMLConfigHandle::get_attribute(const std::string& attribute) SPECMICP_DLL_PUBLIC;
 template <>
-index_t YAMLConfigHandle::get_attribute(const std::string& attribute);
+index_t YAMLConfigHandle::get_attribute(const std::string& attribute) SPECMICP_DLL_PUBLIC;
 template <>
-uindex_t YAMLConfigHandle::get_attribute(const std::string& attribute);
+uindex_t YAMLConfigHandle::get_attribute(const std::string& attribute) SPECMICP_DLL_PUBLIC;
 template <>
-bool YAMLConfigHandle::get_attribute(const std::string& attribute);
+bool YAMLConfigHandle::get_attribute(const std::string& attribute) SPECMICP_DLL_PUBLIC;
 template <>
-std::string YAMLConfigHandle::get_attribute(const std::string& attribute);
+std::string YAMLConfigHandle::get_attribute(const std::string& attribute) SPECMICP_DLL_PUBLIC;
 #endif // SPC_SHOULD_SKIP_THIS
 
 template <typename T>
 T YAMLConfigHandle::get_attribute(
         const std::string& attribute,
         T lower_bound,
         T upper_bound)
 {
     T value = get_attribute<T>(attribute);
     if (value < lower_bound or value > upper_bound) {
         report_error(
                     YAMLConfigError::InvalidArgument,
                     "Attribute '" + attribute + "' is required to"
                     " be between " + std::to_string(lower_bound)
                     + " and " + std::to_string(upper_bound)
                     + " (Got : " + std::to_string(value) + ")."
 
                     );
     }
     return value;
 }
 
 template <typename T>
 T YAMLConfigHandle::get_attribute(
         const std::string& attribute,
         T lower_bound
         )
 {
     T value = get_attribute<T>(attribute);
     if (value < lower_bound ) {
         report_error(
                     YAMLConfigError::InvalidArgument,
                     "Attribute '" + attribute + "' is required to"
                     " be greater than " + std::to_string(lower_bound)
                     + " (Got : " + std::to_string(value) + ")."
                     );
     }
     return value;
 }
 
 #ifndef SPC_DOXYGEN_SHOULD_SKIP_THIS
 template <>
 scalar_t YAMLConfigHandle::get_value(uindex_t ind);
 template <>
 index_t YAMLConfigHandle::get_value(uindex_t ind);
 template <>
 uindex_t YAMLConfigHandle::get_value(uindex_t ind);
 template <>
 bool YAMLConfigHandle::get_value(uindex_t ind);
 template <>
 std::string YAMLConfigHandle::get_value(uindex_t ind);
 
 
 template <>
 std::vector<scalar_t> YAMLConfigHandle::list_to_vector(const std::string& list);
 template <>
 std::vector<index_t> YAMLConfigHandle::list_to_vector(const std::string& list);
 template <>
 std::vector<uindex_t> YAMLConfigHandle::list_to_vector(const std::string& list);
 template <>
 std::vector<bool> YAMLConfigHandle::list_to_vector(const std::string& list);
 template <>
 std::vector<std::string> YAMLConfigHandle::list_to_vector(const std::string& list);
 #endif // SPC_DOXYGEN_SHOULD_SKIP_THIS
 
 // Required attribute
 // ------------------
 template <typename T>
 T YAMLConfigHandle::get_required_attribute(const std::string& attribute)
 {
     if (not has_attribute(attribute)) {
         report_error(YAMLConfigError::MissingRequiredAttribute,
                      attribute);
     }
     return get_attribute<T>(attribute);
 }
 template <typename T>
 T YAMLConfigHandle::get_required_attribute(
         const std::string& attribute,
         T lower_bound
         )
 {
     if (not has_attribute(attribute)) {
         report_error(YAMLConfigError::MissingRequiredAttribute,
                      attribute);
     }
     return get_attribute<T>(attribute, lower_bound);
 }
 
 template <typename T>
 T YAMLConfigHandle::get_required_attribute(
         const std::string& attribute,
         T lower_bound,
         T upper_bound
         )
 {
     if (not has_attribute(attribute)) {
         report_error(YAMLConfigError::MissingRequiredAttribute,
                      attribute);
     }
     return get_attribute<T>(attribute, lower_bound, upper_bound);
 }
 
 // Optional attribute
 // ------------------
 template <typename T>
 T YAMLConfigHandle::get_optional_attribute(
         const std::string& attribute,
         const T& default_value)
 {
     if (not has_attribute(attribute)) {
         return default_value;
     }
     return get_attribute<T>(attribute);
 }
 template <typename T>
 T YAMLConfigHandle::get_optional_attribute(
         const std::string& attribute,
         const T& default_value,
         T lower_bound
         )
 {
     if (not has_attribute(attribute)) {
         return default_value;
     }
     return get_attribute<T>(attribute, lower_bound);
 
 }
 
 //! \brief Return an optional attribute
 template <typename T>
 T YAMLConfigHandle::get_optional_attribute(
         const std::string& attribute,
         const T& default_value,
         T lower_bound,
         T upper_bound
         )
 {
     if (not has_attribute(attribute)) {
         return default_value;
     }
     return get_attribute<T>(attribute, lower_bound, upper_bound);
 }
 
 // Set if attribute exists
 // -----------------------
 template <typename T>
 void YAMLConfigHandle::set_if_attribute_exists(
         T& var_to_set,
         const std::string& attribute
         )
 {
     if (has_attribute(attribute)) {
         var_to_set = get_attribute<T>(attribute);
     }
 }
 
 template <typename T>
 void YAMLConfigHandle::set_if_attribute_exists(
         int& var_to_set,
         const std::string& attribute
         )
 {
     if (has_attribute(attribute)) {
         var_to_set = static_cast<int>(get_attribute<T>(attribute));
     }
 }
 
 
 template <typename T>
 void YAMLConfigHandle::set_if_attribute_exists(
         T& var_to_set,
         const std::string& attribute,
         T lower_bound,
         T upper_bound
         )
 {
     if (has_attribute(attribute)) {
         var_to_set = get_attribute<T>(attribute, lower_bound, upper_bound);
     }
 }
 
 template <typename T>
 void YAMLConfigHandle::set_if_attribute_exists(
         T& var_to_set,
         const std::string& attribute,
         T lower_bound
         )
 {
     if (has_attribute(attribute)) {
         var_to_set = get_attribute<T>(attribute, lower_bound);
     }
 }
 
 template <typename T>
 void YAMLConfigHandle::set_if_attribute_exists(
         int& var_to_set,
         const std::string& attribute,
         T lower_bound
         )
 {
     if (has_attribute(attribute)) {
         var_to_set = get_attribute<T>(attribute, lower_bound);
     }
 }
 
 
 #ifndef SPC_DOXYGEN_SHOULD_SKIP_THIS
 template <>
 std::pair<std::string, scalar_t> YAMLConfigHandle::MapIterator::get_as<scalar_t>();
 template <>
 std::pair<std::string, bool> YAMLConfigHandle::MapIterator::get_as<bool>();
 template <>
 std::pair<std::string, std::string> YAMLConfigHandle::MapIterator::get_as<std::string>();
 template <>
 std::pair<std::string, index_t> YAMLConfigHandle::MapIterator::get_as<index_t>();
 #endif // SPC_DOXYGEN_SHOULD_SKIP_THIS
 
 } //end namespace io
 } //end namespace specmicp
 
 #endif // SPECMICP_UTILS_SAFECONFIG_HPP
diff --git a/src/specmicp_common/micpsolver/micpsolver.inl b/src/specmicp_common/micpsolver/micpsolver.inl
index da78ea0..af73fd6 100644
--- a/src/specmicp_common/micpsolver/micpsolver.inl
+++ b/src/specmicp_common/micpsolver/micpsolver.inl
@@ -1,450 +1,450 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #ifndef SPECMIC_MICPSOLVER_MICPSOLVER_HPP
 #include "micpsolver.hpp" // for syntaxic coloration...
 #endif // SPECMIC_MICPSOLVER_MICPSOLVER_HPP
 
 #include <Eigen/QR>
 #include "estimate_cond_number.hpp"
 
 #include "specmicp_common/log.hpp"
 
 
 //! \file micpsolver.inl implementation of the MiCP solver
 
 namespace specmicp {
 namespace micpsolver {
 
 // Main algorithm
 // ##############
 
 template <class program_t, ReformulationF reform_f>
 inline void MiCPSolver<program_t, reform_f>::setup_residuals(
         const Eigen::VectorXd& x) {
     compute_residuals(x);
     Reformulation<reform_f>::reformulate_residuals(this, x, m_residuals, m_phi_residuals);
 }
 
 template <class program_t, ReformulationF reform_f>
 inline void MiCPSolver<program_t, reform_f>::setup_jacobian(
         Eigen::VectorXd& x) {
     compute_jacobian(x);
     Reformulation<reform_f>::reformulate_jacobian(this, x);
     m_grad_phi = m_jacobian.transpose() * m_phi_residuals;
 }
 
 template <class program_t, ReformulationF reform_f>
 MiCPSolverReturnCode MiCPSolver<program_t, reform_f>::solve(Vector &x)
 {
     int cnt = 0;
     if (get_options().use_crashing) crashing(x);
     else setup_residuals(x);
     MiCPSolverReturnCode retcode = MiCPSolverReturnCode::NotConvergedYet;
     Eigen::VectorXd update(get_neq());
     while (retcode == MiCPSolverReturnCode::NotConvergedYet)
     {
         DEBUG << "Iteration : " << cnt;
         // (S.0) Initialization of the iteration
         // this is a hook for simultaneous fixed-point iterations
         // implemented to solve the non-ideality model
         bool may_have_converged = get_program()->hook_start_iteration(x, m_phi_residuals.norm());
         // (S.1) Check the convergence
         setup_residuals(x);
         get_perfs().current_residual = m_phi_residuals.norm();
         DEBUG << "Norm residuals : " << get_perfs().current_residual;
         retcode = base::check_convergence(cnt, update, x, m_phi_residuals, may_have_converged);
         get_perfs().return_code = retcode;
         if (retcode != MiCPSolverReturnCode::NotConvergedYet) break;
         ++cnt;
         // (S.2) Compute the jacobian and solve the linear problem
         setup_jacobian(x);
         if(get_options().use_scaling)
             search_direction_calculation(update);
         else
             search_direction_calculation_no_scaling(update);
         // (S.3) Linesearch
         int termcode = linesearch(update, x);
         if (termcode != 0) retcode = MiCPSolverReturnCode::LinesearchFailure;
         get_perfs().return_code = retcode;
         get_perfs().current_update = update.norm();
         DEBUG << "Return LineSearch : " << termcode;
         base::projection(x); // project x to the positive quadrant
         get_perfs().nb_iterations = cnt;
     }
     return retcode;
 }
 
 template <class program_t, ReformulationF reform_f>
 MiCPSolverReturnCode MiCPSolver<program_t, reform_f>::search_direction_calculation(Eigen::VectorXd& update)
 {
     using SolverT = Eigen::ColPivHouseholderQR<Eigen::MatrixXd>;
     Vector rscaler(Vector::Ones(m_jacobian.cols()));
     Vector cscaler(Vector::Ones(m_jacobian.rows()));
     base::scaling_jacobian(m_jacobian, m_phi_residuals, rscaler, cscaler);
     m_jacobian = rscaler.asDiagonal() * (m_jacobian) * cscaler.asDiagonal();
     SolverT solver(m_jacobian.rows(), m_jacobian.cols()); // it needs to be correctly initialized
     m_gradient_step_taken = false;
     int m;
     for (m=0; m<get_options().max_factorization_step; ++m)
     {
         const scalar_t lambda = get_options().factor_gradient_search_direction;
         solver.compute(m_jacobian);
         get_perfs().nb_factorization += 1;
         if (solver.info() != Eigen::Success or not solver.isInvertible())
         {
             DEBUG << "Solver.info : " << solver.info() << " - is invertible : " << solver.isInvertible();
             ERROR << "System cannot be solved, we try a perturbation";
             m_gradient_step_taken = true;
             m_jacobian.diagonal() += lambda*rscaler.cwiseProduct(cscaler);
             continue;
         }
         if (get_options().condition_limit > 0)
         {
             scalar_t cond = estimate_condition_number<Eigen::Upper>(solver.matrixR());
             if (cond > get_options().condition_limit)
             {
                 m_gradient_step_taken = true;
                 m_jacobian.diagonal() += lambda*rscaler.cwiseProduct(cscaler);
                 continue;
             }
         }
         update = solver.solve(-rscaler.cwiseProduct(m_phi_residuals + m*lambda*m_grad_phi));
         update  = cscaler.cwiseProduct(update);
         if (get_options().factor_descent_condition < 0 ) break; // we have a solution
         // else we compute the descent condition
         scalar_t descent_cond = m_grad_phi.dot(update);
         scalar_t norm_grad = m_grad_phi.norm();
         scalar_t norm_update = update.norm();
 
         if ((descent_cond <= -get_options().factor_descent_condition*std::min(
                      std::pow(norm_update,2),std::pow(norm_update,3)))
                and (descent_cond <= -get_options().factor_descent_condition*std::min(
                         std::pow(norm_grad,2),std::pow(norm_grad,3)))
              )
             break; // we have a solution !
         m_gradient_step_taken = true;
 
         m_jacobian.diagonal() += lambda*rscaler.cwiseProduct(cscaler);
     }
     DEBUG << "Gradient step : m = " << m;
     if (m == get_options().max_factorization_step) {
         INFO << "Full gradient step taken !";
         update = -m_grad_phi;
     }
     return MiCPSolverReturnCode::NotConvergedYet;
 }
 
 template <class program_t, ReformulationF reform_f>
 MiCPSolverReturnCode MiCPSolver<program_t, reform_f>::search_direction_calculation_no_scaling(
         Vector& update
         )
 {
     using SolverT = Eigen::ColPivHouseholderQR<Eigen::MatrixXd>;
     SolverT solver(m_jacobian.rows(), m_jacobian.cols()); // it needs to be correctly initialized
     // if everything is ok, this is the only decomposition
     m_gradient_step_taken = false;
     int m; // this is an index that increase the contribution of the gradient in the solution
     // the pure Newton solution is when m=0
     for (m=0; m<get_options().max_factorization_step; ++m)
     {
         const scalar_t lambda = get_options().factor_gradient_search_direction;
         solver.compute(m_jacobian); // we need to recompute it (?)
         get_perfs().nb_factorization += 1;
         // Check that the decomposition was ok
         if (solver.info() != Eigen::Success or not solver.isInvertible())
         {
             DEBUG << "Solver.info : " << solver.info() << " (should be " << Eigen::Success
                   << ") - is invertible : " << solver.isInvertible() << " (sould be " << true << ")";
             ERROR << "System cannot be solved, we try a perturbation";
             m_gradient_step_taken = true;
             m_jacobian.diagonal() += Eigen::VectorXd::Constant(m_jacobian.rows(), lambda);
             continue;
         }
         // Estimate the condition number
         if (get_options().condition_limit > 0)
         {
             scalar_t cond = estimate_condition_number<Eigen::Upper>(solver.matrixR());
 
             if (cond > get_options().condition_limit)
             {
                 // add the perturbation
                 m_gradient_step_taken = true;
                 m_jacobian.diagonal() += Eigen::VectorXd::Constant(m_jacobian.rows(), lambda);
             }
         }
         update = solver.solve(-(m_phi_residuals + m*lambda*m_grad_phi));
 
         if (get_options().factor_descent_condition < 0 ) break; // we have a solution
 
         scalar_t descent_cond = m_grad_phi.dot(update);
         scalar_t norm_grad = m_grad_phi.norm();
         scalar_t norm_update = update.norm();
 
         if ( (descent_cond <= -get_options().factor_descent_condition*std::min(
                          std::pow(norm_update,2),std::pow(norm_update,3)))
              and (descent_cond <= -get_options().factor_descent_condition*std::min(
                              std::pow(norm_grad,2),std::pow(norm_grad,3)))
             )
             break; // we have a solution !
 
         // add the perturbation
         m_gradient_step_taken = true;
         m_jacobian.diagonal() += Eigen::VectorXd::Constant(m_jacobian.rows(), lambda);
     }
     DEBUG << "Gradient step : m = " << m;
     if (m_gradient_step_taken && m == get_options().max_factorization_step) {
         INFO << "Full gradient step taken !";
         update = -m_grad_phi;
     }
     return MiCPSolverReturnCode::NotConvergedYet;
 }
 
 template <class program_t, ReformulationF reform_f>
 void MiCPSolver<program_t, reform_f>::crashing(Vector& x)
 {
     // steepest descent direction algorithm
     DEBUG << "Crashing ";
     const scalar_t beta = 0.5;
     const scalar_t sigma = 1e-5;
     int cnt = 0;
     while (cnt < 10)
     {
         setup_residuals(x);
         setup_jacobian(x);
         m_grad_phi = m_jacobian.transpose()*m_phi_residuals;
         Vector xp(get_neq());
         int l=0;
         const int maxl = 10;
         while (l<maxl)
         {
             xp = x - std::pow(beta, l)*m_grad_phi;
             for (int i=get_neq_free(); i<get_neq(); ++i)
             {
                 if (xp(i) < 0) xp(i) = 0;
             }
             Vector new_res(get_neq());
             base::compute_residuals(x, new_res);
             Reformulation<reform_f>::reformulate_residuals_inplace(this, x, new_res);
             scalar_t test = 0.5*(new_res.squaredNorm()
                                  - m_phi_residuals.squaredNorm()
                                  )  + sigma*m_grad_phi.dot(x - xp);
             if (test <= 0) break;
             ++l;
         }
         if (l == maxl) break;
         x = xp;
 
         ++cnt;
     }
     get_perfs().nb_crashing_iterations = cnt;
     DEBUG << "Crashing iterations : " << cnt;
     m_max_merits.reserve(4);
     SPAM << "Solution after crashing \n ------ \n " << x << "\n ----- \n";
 }
 
 // Others
 // ######
 
 template <class program_t, ReformulationF reform_f>
 int MiCPSolver<program_t, reform_f>::linesearch(
         Eigen::VectorXd& p, Eigen::VectorXd& x
         )
 {
     // References
     // ----------
     //    - Algo A6.3.1 : Dennis and Schnabel (1983)
     //    - Munson et al. (2001)
     //    - Facchinei (2003)
     //    - Nocedal & Wrigth (2006)
     DEBUG << "Linesearch";
     Eigen::VectorXd xp(get_neq());
     Eigen::VectorXd new_res(get_neq());
     double fcp;
 
     get_perfs().max_taken = false;
     int retcode = 2;
     const double alpha = 1e-4;
     double newtlen = base::is_step_too_long(p);
     double init_slope = m_grad_phi.dot(p);
     double rellength = std::abs(p(0));
     for (int i=1; i<get_neq(); ++i)
     {
         rellength = std::max(rellength, std::abs(p(i)));
     }
     double minlambda = get_options().steptol / rellength;
     double lambda = get_program()->max_lambda(x, p);
     if (not std::isfinite(lambda)) {
         ERROR << "Lambda is non finite at start of micpsolver linesearch ("
               << lambda << ")- we stop";
         return 5;
     }
     double lambda_prev = lambda;
 
     // non monotone linesearch
     // =======================
 
     double merit_value = 0.5*m_phi_residuals.squaredNorm();
     // new residual
     xp = x + lambda*p;
     base::compute_residuals(xp, new_res);
     Reformulation<reform_f>::reformulate_residuals_inplace(this, xp, new_res);
     fcp = 0.5*new_res.squaredNorm();
 
     // Skip linesearch if enough progress is done
     // ------------------------------------------
     if (fcp < get_options().coeff_accept_newton_step *merit_value)
     {
         if (m_max_merits.size() > 0) m_max_merits[m_max_merits.size()-1] = merit_value;
         else m_max_merits.push_back(merit_value);
         x = xp;
         return 0;
     }
     // Select the merit value of reference
     // -----------------------------------
     if (get_options().non_monotone_linesearch)
     {
         double mmax = merit_value;
         if (m_max_merits.size() > 0)
         {
             mmax = m_max_merits[m_max_merits.size()-1];
             // check for cycling
             // - - - - - - - - -
             if (    m_max_merits.size() ==4
                     &&  std::fabs(mmax - merit_value) < get_options().threshold_cycling_linesearch*get_options().fvectol)
             {
                 //std::cout << merit_value << std::endl;
                 WARNING << "Cycling has been detected by the linesearch - Taking the full Newton step";
                 x = xp;
                 p = lambda*p;
                 return 3;
             }
         }
         if (m_max_merits.size() < 4)
         {
             m_max_merits.push_back(merit_value);
             if (merit_value < mmax) merit_value = (3*merit_value + mmax)/4;
         }
         else if (merit_value < mmax)
         {
             m_max_merits[3] = merit_value;
             merit_value = mmax;
         }
         if (m_gradient_step_taken)
         {
             merit_value *= 100;
         }
     }
     // The linesearch
     // --------------
     double fc = merit_value;
     double fcp_prev;
     int cnt = 0;
     do
     {
         fcp = 0.5*new_res.squaredNorm();
         if (fcp <= fc - std::min(-alpha*lambda*init_slope,(1-alpha)*fc)) //pg760 Fachinei2003
         {
             retcode = 0;
-            if (lambda ==1 and (newtlen > 0.99 * get_options().maxstep)) {
+            if (lambda >=1.0 and (newtlen > 0.99 * get_options().maxstep)) {
                 get_perfs().max_taken = true;
             }
             break;
         }
         else if (lambda < minlambda)
         {
             ERROR << "Linesearch cannot find a solution bigger than the minimal step";
             retcode = 1;
             break;
         }
         else
         {
             // Select a new step length
             // - - - - - - - - - - - -
             double lambdatmp;
             if (cnt == 0) { // only a quadratic at the first
                 lambdatmp = - init_slope / (2*(fcp - fc -init_slope));
             }
             else
             {
                 const double factor = 1 /(lambda - lambda_prev);
                 const double x1 = fcp - fc - lambda*init_slope;
                 const double x2 = fcp_prev - fc - lambda_prev*init_slope;
 
                 const double a = factor * ( x1/(lambda*lambda) - x2/(lambda_prev*lambda_prev));
                 const double b = factor * ( -x1*lambda_prev/(lambda*lambda) + x2*lambda/(lambda_prev*lambda_prev));
 
-                if (a == 0)
+                if (std::abs(a)  < std::numeric_limits<scalar_t>::epsilon())
                 { // cubic interpolation is in fact a quadratic interpolation
                     lambdatmp = - init_slope/(2*b);
                 }
                 else
                 {
                     const double disc = b*b-3*a*init_slope;
                     lambdatmp = (-b+std::sqrt(disc))/(3*a);
                 }
                 if (lambdatmp > 0.5*lambda ) lambdatmp = 0.5*lambda;
             }
             lambda_prev = lambda;
             fcp_prev = fcp;
             if (lambdatmp < 0.1*lambda) {
                 lambda = 0.1 * lambda;
             } else {
                 lambda = lambdatmp;
             }
             if (not std::isfinite(lambda))
             {
                 ERROR << "Lambda is non finite in micpsolver linesearch (" << lambda << ")- we stop";
                 retcode = 3;
                 break;
             }
         }
         xp = x + lambda*p;
         base::compute_residuals(xp, new_res);
         Reformulation<reform_f>::reformulate_residuals_inplace(this, xp, new_res);
         ++cnt;
     } while(retcode == 2 and cnt < 100);
     DEBUG << "Lambda : " << lambda;
     if (cnt == 100)
     {
         ERROR << "Too much linesearch iterations ! We stop";
     }
     x = xp;
     p = lambda*p;
     return retcode;
 }
 
 
 } // end namespace micpsolver
 } // end namespace specmicp
diff --git a/src/specmicp_common/string_algorithms.cpp b/src/specmicp_common/string_algorithms.cpp
index 0862edc..a236455 100644
--- a/src/specmicp_common/string_algorithms.cpp
+++ b/src/specmicp_common/string_algorithms.cpp
@@ -1,422 +1,429 @@
 /* =============================================================================
 
  Copyright (c) 2014 - 2016
  F. Georget <fabieng@princeton.edu> Princeton University
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #include "string_algorithms.hpp"
 
 #include <algorithm>
 
 
 #define RANGE_SEP ':'
 
 
 namespace specmicp {
 namespace utils {
 
 //! \brief Split a string
 std::vector<std::string> split(const std::string& to_split, char separator)
 {
     std::vector<std::string> splitted;
     auto nb_elem = std::count(to_split.begin(), to_split.end(), separator);
     splitted.reserve(nb_elem+1);
 
     auto start = to_split.begin();
     for (auto current=to_split.cbegin(); current!= to_split.cend(); ++current)
     {
         if (*current == separator) {
             splitted.emplace_back(to_split.substr(
                                       start - to_split.cbegin(), current-start));
             start = current + 1;
         }
     }
     // last value
     if (start < to_split.cend()) {
         splitted.emplace_back(to_split.substr(start - to_split.cbegin()));
     }
     return splitted;
 }
 
 
 std::string strip(const std::string& to_trim)
 {
     auto start = to_trim.begin();
     auto end = to_trim.end();
 
     // at the beginning
     auto current=to_trim.begin();
     for (; current!=to_trim.cend(); ++current) {
         if (*current != ' ')
         {
             start = current;
             break;
         }
     }
     if (current == to_trim.cend()) {return "";} // empty string
 
     // at the end
     for (auto current=to_trim.cend()-1;current>=to_trim.begin(); --current) {
         if (*current != ' ')
         {
             end = ++current;
             break;
         }
     }
 
     return to_trim.substr(start - to_trim.begin(), end-start);
 }
 
 
 static void parse_one_term(std::vector<index_t>& numbers, const std::string& term);
 static void parse_one_term(std::vector<uindex_t>& numbers, const std::string& term);
 
 //! \brief Parse a string into a range of indices
 //!
 //! \internal
 template <typename T>
 std::vector<T> range_indices_impl(const std::string& range_str)
 {
     std::vector<T> numbers;
     const uindex_t nb_sep = std::count(range_str.cbegin(), range_str.cend(), ',');
     numbers.reserve(nb_sep+1);
     if (nb_sep == 0) { // if only one term
         parse_one_term(numbers, range_str);
     }
     else {
         const std::vector<std::string> terms = split(range_str, ',');
         for (const auto& term: terms) {
             parse_one_term(numbers, term);
         }
     }
     // sort the indices
     std::sort(numbers.begin(), numbers.end());
     return numbers;
 }
 
 template <>
 std::vector<index_t> range_indices(const std::string& range_str)
 {
     return range_indices_impl<index_t>(range_str);
 }
 
 template <>
 std::vector<uindex_t> range_indices(const std::string& range_str)
 {
     return range_indices_impl<uindex_t>(range_str);
 }
 
 void parse_one_term(std::vector<index_t>& numbers, const std::string &term)
 {
     if (term == "") return; // check for empty terms
 
     const auto to_anal = strip(term);
     const auto is_range = std::find(to_anal.cbegin(), to_anal.cend(), RANGE_SEP);
     if (is_range == to_anal.cend()) {
         // if just a number
         numbers.push_back(std::stol(to_anal));
     }
     else {
         // if a range
         const std::string first_str = to_anal.substr(0, is_range-to_anal.cbegin());
         const std::string last_str = to_anal.substr(is_range-to_anal.cbegin()+1);
         const index_t first = std::stol(first_str);
         const index_t last = std::stol(last_str);
         for (auto ind=first; ind<=last; ++ind) {
             numbers.push_back(ind);
         }
     }
 }
 
 void parse_one_term(std::vector<uindex_t>& numbers, const std::string &term)
 {
     if (term == "") return; // check for empty terms
 
     const auto to_anal = strip(term);
     const auto is_range = std::find(to_anal.cbegin(), to_anal.cend(), RANGE_SEP);
     if (is_range == to_anal.cend()) {
         // if just a number
         auto val = std::stol(to_anal);
         if (val < 0) {
             throw std::invalid_argument("Negative argurment found in "
                                         "positive ranges : '"+std::to_string(val)
                                         + "'.");
         }
         numbers.push_back(val);
     }
     else {
         // if a range
         const std::string first_str = to_anal.substr(0, is_range-to_anal.cbegin());
         const std::string last_str = to_anal.substr(is_range-to_anal.cbegin()+1);
         const index_t first = std::stol(first_str);
         if (first < 0 ) {
             throw std::invalid_argument("Negative argurment found in "
                                         "positive ranges : '"+std::to_string(first)
                                         + "'.");
         }
         const index_t last = std::stol(last_str);
         for (auto ind=first; ind<=last; ++ind) {
             numbers.push_back(ind);
         }
     }
 }
 
 template <typename T>
 class ExpressionParser
 {
 public:
     ExpressionParser(
             const std::string& expr,
             const std::unordered_map<std::string, T>& vars
             ):
         complete_expr(expr),
         variables(vars)
     {
 
     }
 
     T parse() {return parse_term(complete_expr);}
 
     T parse_value(const std::string& value);
     T parse_factor(const std::string& factor);
     T parse_term(const std::string& term);
 
 private:
     const std::string& complete_expr;
     const std::unordered_map<std::string, T>& variables;
 };
 
 template <>
 scalar_t ExpressionParser<scalar_t>::parse_value(const std::string& expr)
 {
     auto trimmed_expr = strip(expr);
     if (trimmed_expr == "") {
         throw std::invalid_argument("Error while parsing : " + complete_expr);
     }
     scalar_t val = std::nan("");
     try {
         std::size_t pos;
         val = std::stod(trimmed_expr, &pos);
         if (pos != trimmed_expr.size()) {
             throw std::logic_error("Error while processing '" + expr +
                                    "'. Just a number was expected."
                                    "Error occured while processing '"
                                    + complete_expr + "'.");
         }
     }
     catch (const std::invalid_argument& e) {
         auto it = variables.find(trimmed_expr);
         if (it == variables.cend()) {
             throw std::out_of_range("Unknown variable '"+
                                         trimmed_expr + "' in parsing of '"
                                         + complete_expr +"'.");
         }
         val = it->second;
     }
     return val;
 }
 
 template <>
 index_t ExpressionParser<index_t>::parse_value(const std::string& expr)
 {
     auto trimmed_expr = strip(expr);
     index_t val = INT_MAX;
     try {
         std::size_t pos;
         val = std::stol(trimmed_expr, &pos);
         if (pos != trimmed_expr.size()) {
             throw std::logic_error("Error while processing" + expr +
                                    ". Just a number was expected."
                                    "Error occured while processing "
                                    + complete_expr + ".");
         }
     }
     catch (const std::invalid_argument& e) {
         auto it = variables.find(trimmed_expr);
         if (it == variables.cend()) {
             throw std::out_of_range("Unknown variable '"+
                                         trimmed_expr + "' in parsing of '"
                                         + complete_expr +"'.");
         }
         val = it->second;
     }
     return val;
 }
 
 template <typename T>
 T ExpressionParser<T>::parse_factor(const std::string& expr)
 {
     auto it = expr.find_last_of("*/");
     if (it == expr.npos) {
         return parse_value(expr);
     }
     auto val1 = parse_factor(expr.substr(0, it));
     auto val2 = parse_value(expr.substr(it+1));
     if (expr[it] == '*') {
         return val1*val2;
     } else {
         return val1/val2;
     }
     return 0;
 }
 
 template <typename T>
 T ExpressionParser<T>::parse_term(const std::string& expr)
 {
     auto it = expr.find_first_of("+-");
     if (it == expr.npos) {
         return parse_factor(expr);
     }
     else if (it != 0) {
         auto val1 = parse_factor(expr.substr(0, it));
         auto val2 = parse_factor(expr.substr(it+1));
         if (expr[it] == '+') {
             return val1+val2;
         } else {
              return val1-val2;
         }
     }
     else {
         if (expr[it] == '-') {
             return -parse_term(expr.substr(1));
         } else {
             return parse_term(expr.substr(1));
         }
     }
     return 0;
 }
 
 template <>
 scalar_t ExpressionParser<scalar_t>::parse_term(const std::string& expr)
 {
     auto pred = [](const char& t) -> bool {return ((t == '+') or( t == '-'));};
     std::size_t pos = 0;
     char op;
     auto it = expr.begin();
     while (pos < expr.size())
     {
         it = std::find_if(it, expr.end(), pred);
         if (it == expr.cend()) {
             pos = expr.size();
             break;
         }
         op = *it;
         pos = it-expr.cbegin();
         // xe-y notation => it's a number !
-        if (op == '-' and expr[pos-1] == 'e' and std::isdigit(expr[pos-2]))
-        {
-            ++it; // jump to next char
-            continue;
-        }
-        else
+        if (op == '-')
         {
+            if (pos>1) {
+                if (expr[pos-1] == 'e' and std::isdigit(expr[pos-2])) {
+                    ++it;
+                    continue;
+                }
+                else {
+                    break;
+                }
+            } else {
+                break;
+            }
+        } else {
             break;
         }
     }
     if (pos == expr.size()) {
         return parse_factor(expr);
     }
 
     if (pos > 0) {
         auto val1 = parse_factor(expr.substr(0, pos));
         auto val2 = parse_factor(expr.substr(pos+1));
         if (op == '+') {
             return val1+val2;
         } else {
              return val1-val2;
         }
     }
     else {
         if (op == '-') {
             return -parse_term(expr.substr(1));
         } else {
             return parse_term(expr.substr(1));
         }
     }
     return 0;
 }
 
 template <>
 scalar_t parse_expression(
         const std::string& expr,
         const std::unordered_map<std::string, scalar_t>& variables
         )
 {
     return ExpressionParser<scalar_t>(expr, variables).parse();
 }
 
 template <>
 index_t parse_expression(
         const std::string& expr,
         const std::unordered_map<std::string, index_t>& variables
         )
 {
     return ExpressionParser<index_t>(expr, variables).parse();
 
 }
 
 
 bool string_to_bool(const std::string& to_bool_str)
 {
     // Contain the strings that matches true or false
     static const char* true_str_test[]  = {"True", "true", "1", "yes", "Yes", "Y"};
     static const char* false_str_test[] = {"False", "false", "0", "no", "No", "N"};
 
     auto trimmed = strip(to_bool_str);
     bool val = false;
     // test for true
     for (auto test: true_str_test) {
         if (test == trimmed) {
             val = true;
             goto exit;
         }
     }
     // test for false
     for (auto test: false_str_test) {
         if (test == trimmed) {
             val = false;
             goto exit;
         }
     }
     // no match => error
     throw std::invalid_argument("Unrecognized value when converting to bool '"
                                 + trimmed + "'. Recognized value : true/false."
                                 );
 
     exit:
     return val;
 }
 
 
 
 } //end namespace utils
 } //end namespace specmicp