diff --git a/examples/Gravity_test.cpp b/examples/Gravity_test.cpp
index 7feec051..f8031bca 100755
--- a/examples/Gravity_test.cpp
+++ b/examples/Gravity_test.cpp
@@ -23,18 +23,20 @@ using namespace std;
 using namespace gravity;
 
 TEST_CASE("testing Gurobi's compound NL") {
+    double pi = 3.123456789;
+    cout << to_string_with_precision(pi, 4) << endl;
+    cout << to_string_with_precision(pi, 8) << endl;
+    cout << to_string_with_precision(pi, 10) << endl;
     Model<> M;
-    M.set_name("M");
-    var<> x("x", -2,2);
-    var<> y("y", 0,0);
-    var<> z("z", pi/2.,pi/2.);
+    M.set_name("M0");
+    var<> x("x", -1,1);
+    var<> y("y", -1,1);
     var<> resvar("resvar", -10,10);
     M.add(resvar);
     M.add(x);
     M.add(y);
-    M.add(z);
     Constraint<> C("C");
-    C = cos(cos(y) - sin(z + cos(y))) - x - 2 - resvar;
+    C = 6*sin(x)/cos(y) - resvar;
     M.add(C==0);
     M.max(resvar);
     solver<> GRB(M,gurobi);
@@ -45,7 +47,7 @@ TEST_CASE("testing Gurobi's compound NL") {
 #ifdef USE_MP
 TEST_CASE("testing readNL() function") {
     Model<> M;
-    string NL_file = string(prj_dir)+"/data_sets/NL/camshape100.nl";
+    string NL_file = "/Users/hassan.hijazi/Downloads/polygon25.nl.txt";
     int status = M.readNL(NL_file);
     solver<> GRB(M,gurobi);
     double solver_time_start = get_wall_time();
diff --git a/examples/Optimization/NonLinear/Power/ACOPF/ACOPF_main.cpp b/examples/Optimization/NonLinear/Power/ACOPF/ACOPF_main.cpp
index b35c670c..a6803386 100644
--- a/examples/Optimization/NonLinear/Power/ACOPF/ACOPF_main.cpp
+++ b/examples/Optimization/NonLinear/Power/ACOPF/ACOPF_main.cpp
@@ -60,6 +60,17 @@ int main (int argc, char * argv[])
         mtype=argv[2];
     }
 #endif
+    Model<> M;
+    int status = M.readNL(fname);
+    solver<> GRB(M,gurobi);
+    GRB.run();
+    M.print();
+    M.read_solution(fname);
+    M.write_solution();
+//    M.reset();
+//    M.is_feasible(1e-6);
+    return 0;
+    
     double total_time_start = get_wall_time();
     PowerNet grid;
     grid.readgrid(fname);
diff --git a/include/gravity/constant.h b/include/gravity/constant.h
index a1c6a96e..282ac6bf 100755
--- a/include/gravity/constant.h
+++ b/include/gravity/constant.h
@@ -52,7 +52,10 @@ namespace gravity {
         if(std::numeric_limits<T>::is_specialized && a_value==numeric_limits<T>::max()){
             return "+∞";
         }
-        out << std::setprecision(n) << a_value;
+        out.setf(ios::fixed);
+        out.precision(n);
+        out << std::scientific << a_value;
+//        out << std::setprecision(n) << a_value;
         return out.str();
     }
     /**
diff --git a/include/gravity/func.h b/include/gravity/func.h
index 07b3a890..754e2f39 100755
--- a/include/gravity/func.h
+++ b/include/gravity/func.h
@@ -353,12 +353,12 @@ namespace gravity {
     template<class T,typename enable_if<is_arithmetic<T>::value>::type* = nullptr>
     T extended_plus(T x, T y){
         
-        if(x==numeric_limits<T>::max() && y==numeric_limits<T>::lowest()){
-            throw invalid_argument("In function extended_plus cannot add +inf to -inf");
-        }
-        if(x==numeric_limits<T>::lowest() && y==numeric_limits<T>::max()){
-            throw invalid_argument("In function extended_plus cannot add -inf to +inf");
-        }
+//        if(x==numeric_limits<T>::max() && y==numeric_limits<T>::lowest()){
+//            throw invalid_argument("In function extended_plus cannot add +inf to -inf");
+//        }
+//        if(x==numeric_limits<T>::lowest() && y==numeric_limits<T>::max()){
+//            throw invalid_argument("In function extended_plus cannot add -inf to +inf");
+//        }
         //+INF
         if(x==numeric_limits<T>::max() || y==numeric_limits<T>::max()){
             return numeric_limits<T>::max();
@@ -1031,7 +1031,7 @@ namespace gravity {
                         return 0.5*bexp->_coef*(lson.get_derivative(v) + rson.get_derivative(v) + df_abs(lson-rson)*(lson.get_derivative(v)-rson.get_derivative(v)));
                         break;
                     case power_:
-                        if (!rson.is_number()) {
+                        if (!rson.is_constant()) {
                             throw invalid_argument("Function in exponent not supported yet.\n");
                         }
                         //                auto exponent = poly_eval(_rson);
@@ -4492,7 +4492,7 @@ namespace gravity {
             if(ex->is_var()){
                 auto vv = static_pointer_cast<param_>(ex);
                 parent_id++;
-                return " ( VARIABLE     , " + vv->get_name(i) + "  , +"+to_string(root_id)+" )\n";
+                return " ( VARIABLE     , " + vv->get_name(i) + "  , "+to_string(root_id)+" )\n";
             }
             else if(ex->is_function()){
                 auto f = static_pointer_cast<func>(ex);
@@ -4503,26 +4503,26 @@ namespace gravity {
                 auto uexp = static_pointer_cast<uexpr<type>>(ex);
                 switch (uexp->_otype) {
                     case cos_:
-                        str = " ( COS      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( COS      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case sin_:
-                        str = " ( SIN      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( SIN      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case sqrt_:
-                        str = " ( SQRT      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( SQRT      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case exp_:{
-                        str = " ( EXP      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( EXP      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     }
                     case log_:
-                        str = " ( LOG      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( LOG      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case abs_:
-                        str = " ( ABS      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( ABS      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case relu_:
-                        str = " ( RELU      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( RELU      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     default:
                         throw invalid_argument("Unsupported unary operation");
@@ -4539,13 +4539,13 @@ namespace gravity {
                 }
                 else if(uexp->_son->is_var()) {
                     auto vv = static_pointer_cast<param_>(uexp->_son);
-                    str += " ( VARIABLE     , " + vv->get_name(i) + "  , +"+to_string(new_root_id)+" )\n";
+                    str += " ( VARIABLE     , " + vv->get_name(i) + "  , "+to_string(new_root_id)+" )\n";
                     parent_id++;
                     return str;
                 }
                 else if(uexp->_son->is_param()) {
                     auto vv = static_pointer_cast<param_>(uexp->_son);
-                    str += " ( CONSTANT     , " + to_string(eval(vv,i)) + "  , +"+to_string(new_root_id)+" )\n";
+                    str += " ( CONSTANT     , " + to_string_with_precision(eval(vv,i), 15) + "  , "+to_string(new_root_id)+" )\n";
                     parent_id++;
                     return str;
                 }
@@ -4555,26 +4555,26 @@ namespace gravity {
                 auto bexp = static_pointer_cast<bexpr<type>>(ex);
                 switch (bexp->_otype) {
                     case plus_:
-                        str = " ( PLUS      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( PLUS      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case minus_:
-                        str = " ( MINUS      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( MINUS      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case product_:{
-                        str = " ( MULTIPLY      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( MULTIPLY      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     }
                     case div_:
-                        str = " ( DIVIDE      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( DIVIDE      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case min_:
-                        str = " ( MIN      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( MIN      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case max_:
-                        str = " ( MAX      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( MAX      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     case power_:
-                        str = " ( POW      , -1 , +"+to_string(root_id)+" )\n";
+                        str = " ( POW      , -1 , "+to_string(root_id)+" )\n";
                         break;
                     default:
                         throw invalid_argument("unsupported operation");
@@ -4590,12 +4590,12 @@ namespace gravity {
                 }
                 else if(bexp->_lson->is_var()) {
                     auto vv = static_pointer_cast<param_>(bexp->_lson);
-                    str += " ( VARIABLE     , " + vv->get_name(i) + "  , +"+to_string(new_root_id)+" )\n";
+                    str += " ( VARIABLE     , " + vv->get_name(i) + "  , "+to_string(new_root_id)+" )\n";
                     parent_id++;
                 }
                 else if(bexp->_lson->is_param()) {
                     auto vv = static_pointer_cast<param_>(bexp->_lson);
-                    str += " ( CONSTANT     , " + to_string(eval(vv,i)) + "  , +"+to_string(new_root_id)+" )\n";
+                    str += " ( CONSTANT     , " + to_string_with_precision(eval(vv,i), 15) + "  , "+to_string(new_root_id)+" )\n";
                     parent_id++;
                 }
                 if (bexp->_rson->is_function()) {
@@ -4605,12 +4605,12 @@ namespace gravity {
                 }
                 else if(bexp->_rson->is_var()) {
                     auto vv = static_pointer_cast<param_>(bexp->_rson);
-                    str += " ( VARIABLE     , " + vv->get_name(i) + "  , +"+to_string(new_root_id)+" )\n";
+                    str += " ( VARIABLE     , " + vv->get_name(i) + "  , "+to_string(new_root_id)+" )\n";
                     parent_id++;
                 }
                 else if(bexp->_rson->is_param()) {
                     auto vv = static_pointer_cast<param_>(bexp->_rson);
-                    str += " ( CONSTANT     , " + to_string(eval(vv,i)) + "  , +"+to_string(new_root_id)+" )\n";
+                    str += " ( CONSTANT     , " + to_string_with_precision(eval(vv,i), 15) + "  , "+to_string(new_root_id)+" )\n";
                     parent_id++;
                 }
                 return str;
@@ -4620,15 +4620,15 @@ namespace gravity {
         string getNLexpr(int root_id, int& parent_id, size_t i){
             if(func_is_number()){
                 parent_id++;
-                return " ( CONSTANT     , " + to_string(eval(i)) + "  , +"+to_string(root_id)+" )\n";
+                return " ( CONSTANT     , " + to_string_with_precision(eval(i),15) + "  , "+to_string(root_id)+" )\n";
             }
             if(is_linear() && _vars->size()==1 && _params->size()==0 && _cst->is_zero() && _lterms->begin()->second._sign && _lterms->begin()->second._coef->is_unit()){
                 parent_id++;
-                return " ( VARIABLE     , " + _vars->begin()->second.first->get_name(i) + "  , +"+to_string(root_id)+" )\n";
+                return " ( VARIABLE     , " + _vars->begin()->second.first->get_name(i) + "  , "+to_string(root_id)+" )\n";
             }
             if(_qterms->size()==0 && _pterms->size()==0 && !_expr && _vars->size()==0 && _params->size()==1 && _cst->is_zero() && _lterms->begin()->second._sign && _lterms->begin()->second._coef->is_unit()){
                 parent_id++;
-                return " ( CONSTANT     , " + to_string(eval(i)) + "  , +"+to_string(root_id)+" )\n";
+                return " ( CONSTANT     , " + to_string_with_precision(eval(i),15) + "  , "+to_string(root_id)+" )\n";
             }
             string str;
             int nb_elem = 0;
@@ -4647,17 +4647,17 @@ namespace gravity {
                 if(root_id==-1)
                     str += " ( PLUS     , -1  , -1 )\n";
                 else
-                    str += " ( PLUS     , -1  , +"+to_string(root_id)+" )\n";
+                    str += " ( PLUS     , -1  , "+to_string(root_id)+" )\n";
                 parent_id++;
                 current_root_id = parent_id;
                 if(!_cst->is_zero()){
-                    str += " ( CONSTANT     , " + to_string(eval(this->get_cst(), i)) + "  , +"+to_string(parent_id)+" )\n";
+                    str += " ( CONSTANT     , " + to_string_with_precision(eval(this->get_cst(), i), 15) + "  , "+to_string(parent_id)+" )\n";
                     parent_id++;
                 }
             }
             int current_parent_id = current_root_id;
             if(_lterms->size()>1){
-                str += " ( PLUS     , -1  , +"+to_string(current_root_id)+" )\n";
+                str += " ( PLUS     , -1  , "+to_string(current_root_id)+" )\n";
                 parent_id++;
                 current_parent_id = parent_id;
             }
@@ -4669,14 +4669,14 @@ namespace gravity {
                         if (!it1.second._sign) {
                             coeff *= -1;
                         }
-                        str += " ( MULTIPLY     , -1  , +"+to_string(current_parent_id)+" )\n";
+                        str += " ( MULTIPLY     , -1  , "+to_string(current_parent_id)+" )\n";
                         parent_id++;
-                        str += " ( CONSTANT     , " + to_string(coeff) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( CONSTANT     , " + to_string_with_precision(coeff, 15) + "  , "+to_string(parent_id)+" )\n";
                         if(it1.second._p->is_var()){
-                            str += " ( VARIABLE     , " + it1.second._p->get_name(i,j) + "  , +"+to_string(parent_id)+" )\n";
+                            str += " ( VARIABLE     , " + it1.second._p->get_name(i,j) + "  , "+to_string(parent_id)+" )\n";
                         }
                         else{
-                            str += " ( CONSTANT     , " + to_string(eval(it1.second._p,i,j)) + "  , +"+to_string(parent_id)+" )\n";
+                            str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p,i,j), 15) + "  , "+to_string(parent_id)+" )\n";
                         }
                         parent_id+=2;
                     }
@@ -4686,21 +4686,21 @@ namespace gravity {
                     if (!it1.second._sign) {
                         coeff *= -1;
                     }
-                    str += " ( MULTIPLY     , -1  , +"+to_string(current_parent_id)+" )\n";
+                    str += " ( MULTIPLY     , -1  , "+to_string(current_parent_id)+" )\n";
                     parent_id++;
-                    str += " ( CONSTANT     , " + to_string(coeff) + "  , +"+to_string(parent_id)+" )\n";
+                    str += " ( CONSTANT     , " + to_string_with_precision(coeff, 15) + "  , "+to_string(parent_id)+" )\n";
                     if(it1.second._p->is_var()){
-                        str += " ( VARIABLE     , " + it1.second._p->get_name(i) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( VARIABLE     , " + it1.second._p->get_name(i) + "  , "+to_string(parent_id)+" )\n";
                     }
                     else{
-                        str += " ( CONSTANT     , " + to_string(eval(it1.second._p,i)) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p,i), 15) + "  , "+to_string(parent_id)+" )\n";
                     }
                     parent_id+=2;
                 }
             }
             current_parent_id = current_root_id;
             if(_qterms->size()>1){
-                str += " ( PLUS     , -1  , +"+to_string(current_root_id)+" )\n";
+                str += " ( PLUS     , -1  , "+to_string(current_root_id)+" )\n";
                 parent_id++;
                 current_parent_id = parent_id;
             }
@@ -4712,20 +4712,20 @@ namespace gravity {
                             if (!it1.second._sign) {
                                 coeff *= -1;
                             }
-                            str += " ( MULTIPLY     , -1  , +"+to_string(current_parent_id)+" )\n";
+                            str += " ( MULTIPLY     , -1  , "+to_string(current_parent_id)+" )\n";
                             parent_id++;
-                            str += " ( CONSTANT     , " + to_string(coeff) + "  , +"+to_string(parent_id)+" )\n";
+                            str += " ( CONSTANT     , " + to_string_with_precision(coeff, 15) + "  , "+to_string(parent_id)+" )\n";
                             if(it1.second._p->first->is_var()){
-                                str += " ( VARIABLE     , " + it1.second._p->first->get_name(j) + "  , +"+to_string(parent_id)+" )\n";
+                                str += " ( VARIABLE     , " + it1.second._p->first->get_name(j) + "  , "+to_string(parent_id)+" )\n";
                             }
                             else{
-                                str += " ( CONSTANT     , " + to_string(eval(it1.second._p->first,j)) + "  , +"+to_string(parent_id)+" )\n";
+                                str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p->first,j), 15) + "  , "+to_string(parent_id)+" )\n";
                             }
                             if(it1.second._p->second->is_var()){
-                                str += " ( VARIABLE     , " + it1.second._p->second->get_name(i) + "  , +"+to_string(parent_id)+" )\n";
+                                str += " ( VARIABLE     , " + it1.second._p->second->get_name(i) + "  , "+to_string(parent_id)+" )\n";
                             }
                             else{
-                                str += " ( CONSTANT     , " + to_string(eval(it1.second._p->second,i)) + "  , +"+to_string(parent_id)+" )\n";
+                                str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p->second,i), 15) + "  , "+to_string(parent_id)+" )\n";
                             }
                             parent_id+=3;
                         }
@@ -4738,20 +4738,20 @@ namespace gravity {
                         if (!it1.second._sign) {
                             coeff *= -1;
                         }
-                        str += " ( MULTIPLY     , -1  , +"+to_string(current_parent_id)+" )\n";
+                        str += " ( MULTIPLY     , -1  , "+to_string(current_parent_id)+" )\n";
                         parent_id++;
-                        str += " ( CONSTANT     , " + to_string(coeff) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( CONSTANT     , " + to_string_with_precision(coeff, 15) + "  , "+to_string(parent_id)+" )\n";
                         if(it1.second._p->first->is_var()){
-                            str += " ( VARIABLE     , " + it1.second._p->first->get_name(i,j) + "  , +"+to_string(parent_id)+" )\n";
+                            str += " ( VARIABLE     , " + it1.second._p->first->get_name(i,j) + "  , "+to_string(parent_id)+" )\n";
                         }
                         else{
-                            str += " ( CONSTANT     , " + to_string(eval(it1.second._p->first,i,j)) + "  , +"+to_string(parent_id)+" )\n";
+                            str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p->first,i,j), 15) + "  , "+to_string(parent_id)+" )\n";
                         }
                         if(it1.second._p->second->is_var()){
-                            str += " ( VARIABLE     , " + it1.second._p->second->get_name(i,j) + "  , +"+to_string(parent_id)+" )\n";
+                            str += " ( VARIABLE     , " + it1.second._p->second->get_name(i,j) + "  , "+to_string(parent_id)+" )\n";
                         }
                         else{
-                            str += " ( CONSTANT     , " + to_string(eval(it1.second._p->second,i,j)) + "  , +"+to_string(parent_id)+" )\n";
+                            str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p->second,i,j), 15) + "  , "+to_string(parent_id)+" )\n";
                         }
                         parent_id+=3;
                     }
@@ -4761,62 +4761,76 @@ namespace gravity {
                     if (!it1.second._sign) {
                         coeff *= -1;
                     }
-                    str += " ( MULTIPLY     , -1  , +"+to_string(current_parent_id)+" )\n";
+                    str += " ( MULTIPLY     , -1  , "+to_string(current_parent_id)+" )\n";
                     parent_id++;
-                    str += " ( CONSTANT     , " + to_string(coeff) + "  , +"+to_string(parent_id)+" )\n";
+                    str += " ( CONSTANT     , " + to_string_with_precision(coeff, 15) + "  , "+to_string(parent_id)+" )\n";
                     if(it1.second._p->first->is_var()){
-                        str += " ( VARIABLE     , " + it1.second._p->first->get_name(i) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( VARIABLE     , " + it1.second._p->first->get_name(i) + "  , "+to_string(parent_id)+" )\n";
                     }
                     else{
-                        str += " ( CONSTANT     , " + to_string(eval(it1.second._p->first,i)) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p->first,i), 15) + "  , "+to_string(parent_id)+" )\n";
                     }
                     if(it1.second._p->second->is_var()){
-                        str += " ( VARIABLE     , " + it1.second._p->second->get_name(i) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( VARIABLE     , " + it1.second._p->second->get_name(i) + "  , "+to_string(parent_id)+" )\n";
                     }
                     else{
-                        str += " ( CONSTANT     , " + to_string(eval(it1.second._p->second,i)) + "  , +"+to_string(parent_id)+" )\n";
+                        str += " ( CONSTANT     , " + to_string_with_precision(eval(it1.second._p->second,i), 15) + "  , "+to_string(parent_id)+" )\n";
                     }
                     parent_id+=3;
                 }
             }
             current_parent_id = current_root_id;
             if(_pterms->size()>1){
-                str += " ( PLUS     , -1  , +"+to_string(current_root_id)+" )\n";
+                str += " ( PLUS     , -1  , "+to_string(current_root_id)+" )\n";
                 parent_id++;
                 current_parent_id = parent_id;
             }
             for (auto &pair:*this->_pterms) {
-                str += " ( MULTIPLY     , -1  , +"+to_string(current_parent_id)+" )\n";
+                str += " ( MULTIPLY     , -1  , "+to_string(current_parent_id)+" )\n";
                 parent_id++;
                 int local_current_parent_id = parent_id;
                 for (auto &vpair: *pair.second._l) {
-                    str += " ( POW     , -1  , +"+to_string(local_current_parent_id)+" )\n";
-                    parent_id++;
-                    if(vpair.first->is_var()){
-                        str += " ( VARIABLE     , " + vpair.first->get_name(i) + "  , +"+to_string(parent_id)+" )\n";
+                    bool has_power = true;
+                    if(vpair.second == unit<type>().eval())
+                        has_power = false;
+                    if(has_power){
+                        str += " ( POW     , -1  , "+to_string(local_current_parent_id)+" )\n";
+                        parent_id++;
+                        if(vpair.first->is_var()){
+                            str += " ( VARIABLE     , " + vpair.first->get_name(i) + "  , "+to_string(parent_id)+" )\n";
+                        }
+                        else{
+                            str += " ( CONSTANT     , " + to_string_with_precision(eval(vpair.first,i), 15) + "  , "+to_string(parent_id)+" )\n";
+                        }
+                        str += " ( CONSTANT     , " + to_string_with_precision(vpair.second, 15) + "  , "+to_string(parent_id)+" )\n";
+                        parent_id+=2;
                     }
                     else{
-                        str += " ( CONSTANT     , " + to_string(eval(vpair.first,i)) + "  , +"+to_string(parent_id)+" )\n";
+                        if(vpair.first->is_var()){
+                            str += " ( VARIABLE     , " + vpair.first->get_name(i) + "  , "+to_string(local_current_parent_id)+" )\n";
+                        }
+                        else{
+                            str += " ( CONSTANT     , " + to_string_with_precision(eval(vpair.first,i), 15) + "  , "+to_string(local_current_parent_id)+" )\n";
+                        }
+                        parent_id+=1;
                     }
-                    str += " ( CONSTANT     , " + to_string(vpair.second) + "  , +"+to_string(parent_id)+" )\n";
-                    parent_id+=2;
                 }
                 auto coeff = eval_coef(pair.second._coef,i);
                 if (!pair.second._sign) {
                     coeff *= -1.;
                 }
-                str += " ( CONSTANT     , " + to_string(coeff) + "  , +"+to_string(local_current_parent_id)+" )\n";
+                str += " ( CONSTANT     , " + to_string_with_precision(coeff, 15) + "  , "+to_string(local_current_parent_id)+" )\n";
                 parent_id++;
                 
             }
             if(_expr){
                 if(_expr->_coef != unit<type>().eval()){
-                    str += " ( PLUS     , -1  , +"+to_string(current_root_id)+" )\n";
+                    str += " ( PLUS     , -1  , "+to_string(current_root_id)+" )\n";
                     parent_id++;
-                    str += " ( MULTIPLY     , -1  , +"+to_string(parent_id)+" )\n";
+                    str += " ( MULTIPLY     , -1  , "+to_string(parent_id)+" )\n";
                     parent_id++;
                     current_root_id = parent_id;
-                    str += " ( CONSTANT     , " + to_string(_expr->_coef) + "  , +"+to_string(parent_id)+" )\n";
+                    str += " ( CONSTANT     , " + to_string_with_precision(_expr->_coef, 15) + "  , "+to_string(parent_id)+" )\n";
                     parent_id++;
                 }
                 str += getNLexpr(_expr, current_root_id, parent_id, i);
@@ -9684,9 +9698,9 @@ namespace gravity {
     
     template<class T1>
     func<T1> log(const func<T1>& f){
-        if(!f.is_positive()){
-            throw invalid_argument("Calling log() with a potentially negative/zero argument");
-        }
+//        if(!f.is_positive()){
+//            throw invalid_argument("Calling log() with a potentially negative/zero argument");
+//        }
         func<T1> res(uexpr<T1>(log_, f.copy()));
         if(f._range->first>=1){
             res._all_sign = non_neg_;
@@ -9749,9 +9763,9 @@ namespace gravity {
     
     template<class T1>
     func<T1> sqrt(const func<T1>& f){
-        if(!f.is_non_negative()){
-            throw invalid_argument("Calling sqrt() with a potentially negative argument");
-        }
+//        if(!f.is_non_negative()){
+//            throw invalid_argument("Calling sqrt() with a potentially negative argument");
+//        }
         func<T1> res(uexpr<T1>(sqrt_, f.copy()));
         res._all_sign = non_neg_;
         if(f.is_positive()){
@@ -10190,63 +10204,61 @@ namespace gravity {
     }
     
     template<class T, typename enable_if<is_arithmetic<T>::value>::type* = nullptr>
-    func<T> pow(const func<T>& f, int exp){
-        if(exp<0){
-            return func<T>(bexpr<T>(power_, f.copy(), make_shared<constant<int>>(exp)));
-        }
+    func<T> pow(const func<T>& f, double exp){
         if(exp==0){
             return func<T>();
         }
         if(exp==1){
             return f;
         }
-        //        if(exp==2){
-        //            return f*f;
-        //        }
+//        if(exp==2){
+//            return f*f;
+//        }
         else {
-            func<T> res(f);
-            for (int i = 1; i < exp; i++) {
-                res *= f;
-            }
-            if(f._range->first==numeric_limits<T>::lowest() || f._range->second==numeric_limits<T>::max()){
-                res._range->first = numeric_limits<T>::lowest();
-                res._range->second = numeric_limits<T>::max();
-            }
-            else {
-                res._range->first = gravity::min(std::pow(f._range->first,exp),std::pow(f._range->second,exp));
-                res._range->second = gravity::max(std::pow(f._range->first,exp),std::pow(f._range->second,exp));
-            }
-            if(exp%2==0) {
-                res._all_sign = non_neg_;
-                if(f.is_positive()){
-                    res._all_sign = pos_;
-                }
-                if(f._range->first <0 && f._range->second >0){
-                    res._range->first = 0;
-                }
-            }
-            else {
-                res._all_sign = f.get_all_sign();
-            }
-            if (f.is_linear()) {
-                if(exp%2==0) {
-                    res._all_convexity = convex_;
-                }
-                else if(f.is_non_negative()){
-                    res._all_convexity = convex_;
-                }
-                else if(f.is_non_positive()){
-                    res._all_convexity = concave_;
-                }
-                else {
-                    res._all_convexity = undet_;
-                }
-            }
-            else if(!f.is_constant()){
-                res._all_convexity = undet_;
-            }
-            res._indices = f._indices;
-            return res;
+//            func<T> res;
+            return func<T>(bexpr<T>(power_, f.copy(), func<T>(exp).copy()));
+//            for (int i = 1; i < exp; i++) {
+//                res *= f;
+//            }
+//            if(f._range->first==numeric_limits<T>::lowest() || f._range->second==numeric_limits<T>::max()){
+//                res._range->first = numeric_limits<T>::lowest();
+//                res._range->second = numeric_limits<T>::max();
+//            }
+//            else {
+//                res._range->first = gravity::min(std::pow(f._range->first,exp),std::pow(f._range->second,exp));
+//                res._range->second = gravity::max(std::pow(f._range->first,exp),std::pow(f._range->second,exp));
+//            }
+//            if(exp%2==0) {
+//                res._all_sign = non_neg_;
+//                if(f.is_positive()){
+//                    res._all_sign = pos_;
+//                }
+//                if(f._range->first <0 && f._range->second >0){
+//                    res._range->first = 0;
+//                }
+//            }
+//            else {
+//                res._all_sign = f.get_all_sign();
+//            }
+//            if (f.is_linear()) {
+//                if(exp%2==0) {
+//                    res._all_convexity = convex_;
+//                }
+//                else if(f.is_non_negative()){
+//                    res._all_convexity = convex_;
+//                }
+//                else if(f.is_non_positive()){
+//                    res._all_convexity = concave_;
+//                }
+//                else {
+//                    res._all_convexity = undet_;
+//                }
+//            }
+//            else if(!f.is_constant()){
+//                res._all_convexity = undet_;
+//            }
+//            res._indices = f._indices;
+//            return res;
         }
     }
     
diff --git a/include/gravity/model.h b/include/gravity/model.h
index 1c4b280b..7f7c8311 100755
--- a/include/gravity/model.h
+++ b/include/gravity/model.h
@@ -2951,9 +2951,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
                                 c->_violated[inst] = false;
                                 diff = std::abs(c->eval(inst));
                                 if(diff > tol) {
-                                    DebugOff("Violated equation: ");
-                                    //                        c->print(inst);
-                                    DebugOff(", violation = "<< diff << endl);
+                                    DebugOn("Violated equation: ");
+                                                            c->print(inst);
+                                    DebugOn(", violation = "<< diff << endl);
                                     nb_viol++;
                                     //                        violated = true;
                                     if (*c->_all_lazy) {
@@ -2979,9 +2979,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
                                 c->_violated[inst] = false;
                                 diff = c->eval(inst);
                                 if(diff > tol) {
-                                    DebugOff("Violated inequality: ");
-                                    //                                c->print(inst);
-                                    DebugOff(", violation = "<< diff << endl);
+                                    DebugOn("Violated inequality: ");
+                                                                    c->print(inst);
+                                    DebugOn(", violation = "<< diff << endl);
                                     nb_viol++;
                                     //                        violated = true;
                                     if (*c->_all_lazy) {
@@ -3013,9 +3013,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
                                 c->_violated[inst] = false;
                                 diff = c->eval(inst);
                                 if(diff < -tol) {
-                                    DebugOff("Violated inequality: ");
-                                    //                        c->print(inst);
-                                    DebugOff(", violation = "<< diff << endl);
+                                    DebugOn("Violated inequality: ");
+                                                            c->print(inst);
+                                    DebugOn(", violation = "<< diff << endl);
                                     nb_viol++;
                                     //                        violated = true;
                                     if (*c->_all_lazy) {
@@ -5077,11 +5077,14 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
 
             for(auto& c_p :_cons)
             {
-//                c_p.second->uneval();
                 c_p.second->_new = true;
+                c_p.second->uneval();
+                c_p.second->eval_all();
                 c_p.second->_dfdx->clear();
             }
             _obj->_new = true;
+            _obj->uneval();
+            _obj->eval_all();
             _obj->_dfdx->clear();
             for(auto &v_p: _vars)
             {
@@ -5859,6 +5862,21 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
             return size_header;
         }
         
+    /** Read solution point to file */
+        void read_solution(const string& fname);
+        
+    /** Write solution point to file */
+        void write_solution(int precision = 10){
+            ofstream myfile;
+            string fname = _name+".sol";
+            myfile.open(fname);
+            std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf
+            std::cout.rdbuf(myfile.rdbuf());
+            print_solution(precision);
+            std::cout.rdbuf(coutbuf);
+            myfile.close();
+        }
+        
         /** Write Model to file */
         void write(int precision = 10){
             ofstream myfile;
@@ -5897,6 +5915,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
         
         
         
+        bool has_int() const{
+            return !_int_vars.empty();
+        }
         
         void replace_integers(){/*< Replace internal type of integer variables so that continuous relaxations can be computed */
             bool has_int = false;
@@ -8195,11 +8216,15 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
         }
         
         var<> get_cont_int_var(int index){
-            auto x = _model->get_var<double>("x");
-            auto y = _model->get_var<double>("y");
-            if(index >= y.get_id())/* Integer variable */
-                return y(index);
+            if(_model->has_int()){
+                auto y = _model->get_var<double>("y");
+                if(y._indices->has_key(to_string(index)))/* Integer variable */
+                {
+                    return y(index);
+                }
+            }
             /* Continuous variable */
+            auto x = _model->get_var<double>("x");
             return x(index);
         }
         
@@ -8271,6 +8296,9 @@ const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<str
         func<> VisitLog(UnaryExpr e) {
             return log(Visit(e.arg()));
         }
+        func<> VisitLog10(UnaryExpr e) {
+            return log(Visit(e.arg()))/std::log(10);
+        }
 //        func<> VisitFloor(UnaryExpr e) {
 //          return IloFloor(Visit(e.arg()));
 //        }
diff --git a/include/gravity/param.h b/include/gravity/param.h
index 582fb2fb..5a84c0be 100755
--- a/include/gravity/param.h
+++ b/include/gravity/param.h
@@ -2351,10 +2351,10 @@ namespace gravity {
                     }
                     else {
                         for (size_t i = 0; i < _dim[0]; i++) {
-                            if((!_is_relaxed && !is_integer()) || std::abs(eval(i)) > 1e-4){
+//                            if((!_is_relaxed && !is_integer()) || std::abs(eval(i)) > 1e-4){
                                 str += name+"(" + _indices->_keys->at(i) + ") " + to_string_with_precision(eval(i), prec);
                                 str += " \n";
-                            }
+//                            }
                         }
                     }
                 }
diff --git a/include/gravity/solver.h b/include/gravity/solver.h
index 3bd32d99..f958cc5e 100755
--- a/include/gravity/solver.h
+++ b/include/gravity/solver.h
@@ -434,6 +434,7 @@ class solver {
                     grb_prog->_output = output;
                     //            prog.grb_prog->reset_model();
                         grb_prog->prepare_model();
+                    return 0;
                         //DebugOn("calling prep after init"<<endl);
                     optimal = grb_prog->solve(relax, tol, time_limit);
                     return_status = optimal ? 0 : -1;
diff --git a/src/GurobiProgram.cpp b/src/GurobiProgram.cpp
index aeba4301..92401fff 100755
--- a/src/GurobiProgram.cpp
+++ b/src/GurobiProgram.cpp
@@ -430,7 +430,7 @@ GurobiProgram::GurobiProgram(Model<>* m) {
         }
     }
     _model = m;
-    m->fill_in_maps();
+//    m->fill_in_maps();
     m->compute_funcs();
 }
 
@@ -814,10 +814,14 @@ bool GurobiProgram::solve(bool relax, double mipgap, double time_limit){
     return true;
 }
 
+#include <fstream>
+#include <iostream>
+#include <string>
+
 bool copyWithoutLastLine(const std::string& inputFile) {
     std::string tempFile = inputFile + ".temp"; // Temporary file
 
-    // First, perform the operation with inputFile as input and tempFile as output
+    // Open input and output files
     std::ifstream in(inputFile);
     std::ofstream out(tempFile, std::ios::trunc);
     if (!in.is_open() || !out.is_open()) {
@@ -827,17 +831,22 @@ bool copyWithoutLastLine(const std::string& inputFile) {
 
     std::string line;
     std::string prevLine;
-    while (std::getline(in, line)) {
+    bool foundGenerals = false;
+    while (!foundGenerals && std::getline(in, line)) {
+        if (line.find("Generals") != std::string::npos) {
+            foundGenerals = true; // Flag that we found the "Generals" line
+        }
         if (!prevLine.empty()) {
-            out << prevLine << '\n'; // Write previous line, ensuring not to write the last one
+            out << prevLine << '\n'; // Write previous line
         }
-        prevLine = line; // Keep track of the previous line
+        prevLine = line; // Update previous line
     }
-    
+
+    // Close files
     in.close();
     out.close();
 
-    // Now, replace the original file with the temporary file
+    // Replace original file with temporary file
     if (std::remove(inputFile.c_str()) != 0) {
         std::cerr << "Error removing original file." << std::endl;
         return false;
@@ -851,8 +860,9 @@ bool copyWithoutLastLine(const std::string& inputFile) {
 }
 
 
+
 void GurobiProgram::prepare_model(){
-    _model->fill_in_maps();
+//    _model->fill_in_maps();
     _model->compute_funcs();
     fill_in_grb_vmap();
     create_grb_constraints();
@@ -868,7 +878,7 @@ void GurobiProgram::prepare_model(){
     write_NLCstr(filename);
     
     
-    exit(0);
+    return;
     //    print_constraints();
 }
 void GurobiProgram::update_model(){
@@ -1017,40 +1027,81 @@ void GurobiProgram::write_NLCstr(const string &fname){
         cerr << "Failed to open file for appending nonlinear constraints.\n";
         return;
     }
-//    file << " zero = 0\n";
+    size_t idx = 0, idx_inst = 0, idx1 = 0, idx2 = 0, idx_inst1 = 0, idx_inst2 = 0, nb_inst = 0, inst = 0;
+    int pos_i = 0, neg_i = 0;
+    for(auto& p: _model->_cons){
+        auto c = p.second;
+        if (c->is_polynomial() || c->is_nonlinear()) {
+            nb_inst = c->get_nb_inst();
+            for (size_t i = 0; i< nb_inst; i++){
+                switch(c->get_ctype()) {
+                    case geq:
+                        pos_i++;
+                        break;
+                    case leq:
+                        neg_i++;
+                        break;
+                    default:
+                        break;
+                }
+            }
+        }
+    }
+    file << " zero(0) = 0\n";
+    for(int i = 0; i < pos_i; i++){
+        file << " 0 <= pos(" + to_string(i) + ")\n";
+    }
+    for(int i = 0; i < neg_i; i++){
+        file << " -Inf <= neg(" + to_string(i) + ") <= 0\n";
+    }
+    if(_model->has_int()){
+        auto y = _model->get_var_ptr("y");
+        file << " Generals\n";
+        file << "  ";
+        for (auto i = 0; i<y->get_dim(); i++){
+            file << y->get_name(i) << " " ;
+        }
+        file << "\n";
+    }
     file << " General Constraints\n";
 
     char sense;
-    size_t idx = 0, idx_inst = 0, idx1 = 0, idx2 = 0, idx_inst1 = 0, idx_inst2 = 0, nb_inst = 0, inst = 0;
     GRBLinExpr lterm, linlhs;
     GRBQuadExpr quadlhs;
     GRBVar gvar1, gvar2;
     double coeff;
+    pos_i = 0;
+    neg_i = 0;
     for(auto& p: _model->_cons){
         auto c = p.second;
-        if (c->is_nonlinear() && (!(c->_expr->is_uexpr() && c->get_nb_vars()==2) && !c->_expr->is_mexpr())) {
-            DebugOff("Founf a nonlinear constraint with more than two variables, writing the expression tree.\n");
-            switch(c->get_ctype()) {
-                case geq:
-                    sense = GRB_GREATER_EQUAL;
-                    break;
-                case leq:
-                    sense = GRB_LESS_EQUAL;
-                    break;
-                case eq:
-                    sense = GRB_EQUAL;
-                    break;
-                default:
-                    break;
-            }
-            assert(sense==GRB_EQUAL);/* Check what to do with inequalities */
+//        if(c->get_ctype()==eq)
+//            DebugOn(c->get_name() << endl);
+//        if(c->get_name()=="NL_C_eq_73")
+//            DebugOn(c->get_name() << endl);
+        if (c->is_polynomial() || c->is_nonlinear()) {
+            DebugOff("Found a nonlinear constraint with more than two variables, writing the expression tree.\n");
             nb_inst = c->get_nb_inst();
             for (size_t i = 0; i< nb_inst; i++){
                 if(c->_indices)
                     file <<  " " << c->get_name()+"("+c->_indices->_keys->at(i)+"):\n";
                 else
                     file <<  " " << c->get_name()+"("+to_string(i)+"):\n";
-                file << " zero(0) = NL :\n";
+                switch(c->get_ctype()) {
+                    case geq:
+                        sense = GRB_GREATER_EQUAL;
+                        file << " pos("+ to_string(pos_i++) +") = NL :\n";
+                        break;
+                    case leq:
+                        sense = GRB_LESS_EQUAL;
+                        file << " neg("+ to_string(neg_i++) +") = NL :\n";
+                        break;
+                    case eq:
+                        sense = GRB_EQUAL;
+                        file << " zero(0) = NL :\n";
+                        break;
+                    default:
+                        break;
+                }            
                 int parent_id = -1;
                 file << c->getNLexpr(-1,parent_id,i);
             }
@@ -1086,7 +1137,7 @@ void GurobiProgram::create_grb_constraints(){
         }
         c->_new = false;
         
-        if (c->is_nonlinear() && (!(c->_expr->is_uexpr() && c->get_nb_vars()==2) && !c->_expr->is_mexpr())) {
+        if (c->is_polynomial() || c->is_nonlinear()) {
             DebugOff("Gurobi cannot handle nonlinear constraints with more than two variables, ignoring it.\n");
             continue;
         }
@@ -1228,7 +1279,7 @@ void GurobiProgram::create_grb_constraints(){
                 }
                 else {
                     coeff = c->eval(lt._coef,i);
-                    if (!((coeff==1 && lt._sign) || (coeff==-1 && !lt._sign))) {
+                    if (!((coeff==1 && !lt._sign) || (coeff==-1 && lt._sign))) {
                         throw invalid_argument("Gurobi does not support this type of nonlinear constraints");
                     }
                     gvar1 = _grb_vars[lt._p->get_id() + lt._p->get_id_inst(i)];
diff --git a/src/model.cpp b/src/model.cpp
index ef2705e1..434ae74a 100755
--- a/src/model.cpp
+++ b/src/model.cpp
@@ -9,10 +9,143 @@
 #include <gravity/model.h>
 #include <gravity/solver.h>
 #include <math.h> //for setting the rounding direction
+//#include <algorithm>
+//#include <string>
+//#include <iostream>
 
 using namespace std;
 namespace gravity {
 
+static int max_line_len;
+static char* line = nullptr;
+
+char* readLine(FILE *input)
+{
+    size_t len;
+    if(fgets(line,max_line_len,input) == NULL)
+        return NULL;
+    
+    while(strrchr(line,'\n') == NULL)
+    {
+        max_line_len *= 2;
+        line = (char *) realloc(line,max_line_len);
+        len = strlen(line);
+        if(fgets(line+len,max_line_len-len,input) == NULL)
+            break;
+    }
+    return line;
+}
+
+void skip_lines(FILE *input, size_t nb_lines){
+    for (size_t i = 0; i<nb_lines; i++) {
+        readLine(input);
+    }
+}
+
+char *mystrtok(char **m,char *s,char c)
+{
+    char *p1=s?s:*m;
+    if( !*p1 )
+        return 0;
+    *m=strchr(p1,c);
+    if(*m && c==' '){
+        while(*m[0]==c){
+            *(*m)++=0;
+        }
+        if(! *m )
+            *m=p1+strlen(p1);
+        return p1;
+    }
+    if( *m )
+        *(*m)++=0;
+    else
+        *m=p1+strlen(p1);
+    return p1;
+}
+
+/** Read solution point to file */
+template <typename type>
+void Model<type>::read_solution(const string& fname){
+    std::string base_filename = fname.substr(0, fname.find_last_of("."));
+
+    string fname_sol = base_filename+".p1.sol";
+    FILE *fp = fopen(fname_sol.c_str(),"r");
+    if(fp == NULL)
+    {
+            cout << "Can’t open input file " << fname;
+            exit(1);
+    }
+    max_line_len = 1024;
+    line = new char[max_line_len];
+    auto n = get_nb_vars();
+    map<string,double> sol_map;
+    double val;
+    string vname;
+    char* p;
+    while(readLine(fp))
+    {
+        vname = mystrtok(&p,line,' ');
+        val = atof(mystrtok(&p,NULL,' '));
+        sol_map[vname] = val;
+    }
+    delete[] line;
+    fclose(fp);
+    /*read col file*/
+    string fname_col = base_filename+".col";
+    fp = fopen(fname_col.c_str(),"r");
+    if(fp == NULL)
+    {
+            cout << "Can’t open input file " << fname;
+            exit(1);
+    }
+    line = new char[max_line_len];
+    map<string,string> v_map;
+    int nb_cont = 0, nb_int = 0;
+    shared_ptr<param_> x = this->get_var_ptr("x");
+    shared_ptr<var<double>> x_cont = nullptr, y_cont = nullptr;
+    shared_ptr<var<int>> y_int = nullptr;
+    if(x){
+        x_cont = static_pointer_cast<var<double>>(x);
+        nb_cont = x_cont->get_dim();
+    }
+    shared_ptr<param_> y = this->get_var_ptr("y");
+    if(y){
+        y_int = static_pointer_cast<var<int>>(_int_vars.begin()->second);
+        y_cont = static_pointer_cast<var<double>>(y);
+        nb_int = y_int->get_dim();
+    }
+    
+    for (int i = 0; i<nb_cont; i++) {
+        readLine(fp);
+        vname = mystrtok(&p,line,'\n');
+//        vname.erase(remove(vname.begin(), vname.end(), '\n'), vname.end());
+        v_map[vname] = x_cont->get_name(i);
+        if(sol_map.count(vname)>0){
+            x_cont->_val->at(i) = sol_map[vname];
+        }
+        else{
+            x_cont->_val->at(i) = 0;
+        }
+    }
+    for (int i = 0; i<nb_int; i++) {
+        readLine(fp);
+        vname = mystrtok(&p,line,'\n');
+        v_map[vname] = y_int->get_name(i);
+        if(sol_map.count(vname)>0){
+            y_int->_val->at(i) = sol_map[vname];
+            y_cont->_val->at(i) = sol_map[vname];
+        }
+        else{
+            y_int->_val->at(i) = 0;
+            y_cont->_val->at(i) = 0;
+        }
+    }
+    auto objvar = static_pointer_cast<var<double>>(this->get_var_ptr("objvar"));
+    objvar->_val->at(0) = sol_map.at("objvar");
+    delete[] line;
+    fclose(fp);
+    this->print_solution();
+}
 
 const bool var_compare(const pair<string,shared_ptr<param_>>& v1, const pair<string,shared_ptr<param_>>& v2) {
     return v1.second->get_nb_rows() > v2.second->get_nb_rows();
@@ -6735,7 +6868,11 @@ void Model<type>::update_upper_bound(shared_ptr<Model<type>>& obbt_model, vector
     
 template <typename type>
 template<typename T,typename std::enable_if<is_arithmetic<T>::value>::type*>
-int Model<type>::readNL(const string& fname){    
+int Model<type>::readNL(const string& fname){
+    std::string base_filename = fname.substr(fname.find_last_of("/\\") + 1);
+    string::size_type const file_with_ext(base_filename.find_last_of('.'));
+    string file_without_extension = base_filename.substr(0, file_with_ext);
+    this->set_name(file_without_extension);
     mp::Problem p;
     mp::ReadNLFile(fname, p);
     auto nb_vars = p.num_vars();
@@ -6765,8 +6902,8 @@ int Model<type>::readNL(const string& fname){
     DebugOn("Number of continuous variables = " << nb_cont << endl);
     DebugOn("Number of integer variables = " << nb_int << endl);
     
-    param<> x_ub("x_ub"), x_lb("x_lb");
-    param<int> y_ub("y_ub"), y_lb("y_lb");
+    param<> x_ub("x-ub"), x_lb("x-lb");
+    param<int> y_ub("y-ub"), y_lb("y-lb");
     x_ub.in(C);x_lb.in(C);
     y_ub.in(I);y_lb.in(I);
     for (int i = 0; i<C.size(); i++) {
@@ -6784,29 +6921,75 @@ int Model<type>::readNL(const string& fname){
     int nb_lin = 0;
     int nb_nonlin = 0;
     int index = 0;
-    _name = fname;
-
-    add(x.in(C));
-    add(y.in(I));
+    
 
-    replace_integers();
+    if(!C.empty())
+        add(x.in(C));
+    if(!I.empty()){
+        add(y.in(I));
+        replace_integers();
+    }
 
     MPConverter converter(*this);
     map<int,vector<int>> constr_sparsity;
     vector<int> C_lin, C_nonlin, C_quad;
+    
+    int num_objs = p.num_objs();
+    if(num_objs>=1){
+        if(num_objs>=2){
+            DebugOn("Gravity currently supports only one objective, will only add the first one");
+        }
+        mp::Problem::Objective obj = p.obj(0);
+        mp::obj::Type main_obj_type = p.obj(0).type();
+        func<> objective;
+        auto lexpr = obj.linear_expr();
+        for (const auto term: lexpr){
+            auto coef = term.coef();
+            auto var_id = term.var_index();
+            if(coef!=0)
+                objective += coef*converter.get_cont_int_var(var_id);
+        }
+        auto nl_expr = obj.nonlinear_expr();
+        if (nl_expr){
+            auto expr = converter.Visit(nl_expr);
+            objective += expr;
+            var<> objvar("objvar");
+            add(objvar);
+            Constraint<> c("Objvar_def");
+            c += objective - objvar;
+            add(c.in(range(1,1)) == 0);
+            auto sense = main_obj_type == mp::obj::MIN;
+            if(sense){
+                min(objvar);
+            }
+            else {
+                max(objvar);
+            }
+        }
+        else{
+            auto sense = main_obj_type == mp::obj::MIN;
+            if(sense){
+                min(objective);
+            }
+            else {
+                max(objective);
+            }
+        }
+    }
+    
     for (const auto con: p.algebraic_cons()) {
         auto lexpr = con.linear_expr();
         auto nl_expr = con.nonlinear_expr();
         if (nl_expr){
             auto expr = converter.Visit(nl_expr);
-            expr.print();
             nb_nonlin++;
             C_nonlin.push_back(index);
             NonLinConstr.insert(to_string(index));
             for (const auto term: lexpr){
                 auto coef = term.coef();
                 auto var_id = term.var_index();
-                expr += coef*converter.get_cont_int_var(var_id);
+                if(coef!=0)
+                    expr += coef*converter.get_cont_int_var(var_id);
             }
             
             auto c_lb = con.lb();
@@ -6814,62 +6997,87 @@ int Model<type>::readNL(const string& fname){
             if(c_lb==c_ub){
                 Constraint<> c("NL_C_eq_"+to_string(index));
                 c += expr;
-                add(c == c_lb);
+                add(c.in(range(1,1)) == c_lb);
             }
             else {
                 if(c_lb>numeric_limits<double>::lowest()){
-                    Constraint<> c("NL_C_lb_"+to_string(index));
+                    Constraint<> c("NL_C_geq_"+to_string(index));
                     c += expr;
-                    add(c >= c_lb);
+                    add(c.in(range(1,1)) >= c_lb);
                 }
                 if(c_ub<numeric_limits<double>::max()){
-                    Constraint<> c("NL_C_ub_"+to_string(index));
+                    Constraint<> c("NL_C_leq_"+to_string(index));
                     c += expr;
-                    add(c <= c_ub);
+                    add(c.in(range(1,1)) <= c_ub);
                 }
             }
         }
         else{
-            int nb_terms = lexpr.num_terms();
-            constr_sparsity[nb_terms].push_back(index);
+            int nb_terms = 0;
             func<> expr;
             for (const auto term: lexpr){
                 auto coef = term.coef();
                 auto var_id = term.var_index();
-                expr += coef*converter.get_cont_int_var(var_id);
+                if(coef!=0){
+                    expr += coef*converter.get_cont_int_var(var_id);
+                    nb_terms++;
+                }
             }
-            
             auto c_lb = con.lb();
             auto c_ub = con.ub();
-            if(c_lb==c_ub){
-                Constraint<> c("Lin_C_eq_"+to_string(index));
-                c += expr;
-                add(c == c_lb);
+            if(false && nb_terms==1 && c_lb!=c_ub){/* this is just a bound constraint */
+                const auto term = *lexpr.begin();
+                auto coef = term.coef();
+                auto var_id = term.var_index();
+                auto v = converter.get_cont_int_var(var_id);
+                if(coef != 0 && c_lb>numeric_limits<double>::lowest()){
+                    if(v._is_relaxed){/* an integer */
+                        v._lb->_val->at(v.get_id_inst()) = c_lb/coef;
+                    }
+                    else {
+                        v._lb->_val->at(v.get_id_inst()) = c_lb/coef;
+                    }
+                }
+                if(coef != 0 && c_ub<numeric_limits<double>::max()){
+                    if(v._is_relaxed){/* an integer */
+                        v._ub->_val->at(v.get_id_inst()) = c_ub/coef;
+                    }
+                    else {
+                        v._ub->_val->at(v.get_id_inst()) = c_ub/coef;
+                    }
+                }
             }
-            else {
-                if(c_lb>numeric_limits<double>::lowest()){
-                    Constraint<> c("Lin_C_lb_"+to_string(index));
+            else{
+                constr_sparsity[nb_terms].push_back(index);
+                if(c_lb==c_ub){
+                    Constraint<> c("Lin_C_eq_"+to_string(index));
                     c += expr;
-                    add(c >= c_lb);
+                    add(c.in(range(1,1)) == c_lb);
                 }
-                if(c_ub<numeric_limits<double>::max()){
-                    Constraint<> c("Lin_C_ub_"+to_string(index));
-                    c += expr;
-                    add(c <= c_ub);
+                else {
+                    if(c_lb>numeric_limits<double>::lowest()){
+                        Constraint<> c("Lin_C_geq_"+to_string(index));
+                        c += expr;
+                        add(c.in(range(1,1)) >= c_lb);
+                    }
+                    if(c_ub<numeric_limits<double>::max()){
+                        Constraint<> c("Lin_C_leq_"+to_string(index));
+                        c += expr;
+                        add(c.in(range(1,1)) <= c_ub);
+                    }
                 }
+                LinConstr.insert(to_string(index));
+                C_lin.push_back(index);
+                nb_lin++;
             }
-            LinConstr.insert(to_string(index));
-            C_lin.push_back(index);
-            nb_lin++;
         }
-        
         index++;
     }
     DebugOn("Number of linear constraints = " << nb_lin << endl);
     DebugOn("Number of non linear constraints = " << nb_nonlin << endl);
     DebugOn("Number of sparsity degrees for linear constraints = " << constr_sparsity.size() << endl);
-
-    print();
+//    this->print();
+    
     return 0;
 }
     
@@ -6882,7 +7090,7 @@ int Model<type>::readNL(const string& fname){
     template Constraint<Cpx> Model<Cpx>::lift(Constraint<Cpx>& c, string model_type);
     template Constraint<> Model<>::lift(Constraint<>& c, string model_type);
     template int gravity::Model<double>::readNL<double, (void*)0>(const string&);
-    
+    template void gravity::Model<double>::read_solution(const string& fname);
     
     //    template void Model<double>::run_obbt(double max_time, unsigned max_iter);
     //    template func<double> constant<double>::get_real() const;