diff --git a/nl-writer2/include/api/c/nl-solver-c.h b/nl-writer2/include/api/c/nl-solver-c.h
index 79022724d..b4cd7f562 100644
--- a/nl-writer2/include/api/c/nl-solver-c.h
+++ b/nl-writer2/include/api/c/nl-solver-c.h
@@ -93,7 +93,7 @@ const char* NLW2_GetErrorMessage_C(NLW2_NLSolver_C* );
 
 /// Load and solve NLW2_NLModel_C and return result.
 ///
-/// @return Solution object with computed obj_value_.
+/// @return NLW2_Solution_C object with computed obj_value_.
 ///   Valid as long as the NLW2_NLSolver_C object lives,
 ///   and until the next Solve() or ReadSolution().
 ///
@@ -149,7 +149,7 @@ int NLW2_RunSolver_C(NLW2_NLSolver_C* ,
 /// Read solution after solving a model
 /// loaded with NLW2_LoadNLModel_C().
 ///
-/// @return Solution object.
+/// @return NLW2_Solution_C object.
 ///   Valid as long as the NLW2_NLSolver_C object lives,
 ///   and until the next Solve() or ReadSolution().
 ///
diff --git a/nl-writer2/include/mp/nl-solver.h b/nl-writer2/include/mp/nl-solver.h
index b99cb2880..a4d85bc95 100644
--- a/nl-writer2/include/mp/nl-solver.h
+++ b/nl-writer2/include/mp/nl-solver.h
@@ -143,7 +143,7 @@ class NLSolver {
 
   /// Load and solve an NLModel instance.
   ///
-  /// @return Solution object with computed obj_value_
+  /// @return NLSolution object with computed obj_value_
   ///    (has operator bool() for checking
   ///     if any result was obtained.)
   ///
@@ -209,7 +209,7 @@ class NLSolver {
   /// Read solution after Solve()
   /// when the NL file was written from NLModel.
   ///
-  /// @return Solution object
+  /// @return NLSolution object
   ///    (has operator bool() for checking
   ///     if any result was obtained.)
   ///
diff --git a/nl-writer2/nlwpy/src/nlw_bindings.cc b/nl-writer2/nlwpy/src/nlw_bindings.cc
index 5d48784dc..84afc045d 100644
--- a/nl-writer2/nlwpy/src/nlw_bindings.cc
+++ b/nl-writer2/nlwpy/src/nlw_bindings.cc
@@ -33,10 +33,8 @@
 namespace py = pybind11;
 
 
-/// Variables' data by pointers
+/// Variables' data
 struct NLWPY_ColData {
-  /// Num vars
-  int num_col_;
   /// lower bounds
   std::vector<double> lower_;
   /// upper bounds
@@ -47,16 +45,14 @@ struct NLWPY_ColData {
 };
 
 /// Sparse matrix.
+///
+/// Size of the start_ array:
+/// N cols (for colwise) / N rows (for rowwise),
+/// depending on format_.
 struct NLWPY_SparseMatrix {
-  /// Size of the start_ array:
-  /// N cols (for colwise) / N rows (for rowwise),
-  /// depending on format_.
-  int num_colrow_;
   /// Format (NLW2_MatrixFormat...).
   /// Only rowwise supported.
   NLW2_MatrixFormat format_;
-  /// Nonzeros
-  size_t num_nz_;
   /// Row / col starts
   std::vector<size_t> start_;
   /// Entry index
@@ -79,16 +75,17 @@ class NLWPY_NLModel {
   { }
 
   /// Add variables (all at once.).
-  /// @todo ty can be None.
-  void SetCols(int n,
-               std::vector<double> lb,
+  /// @todo \a ty can be None.
+  void SetCols(std::vector<double> lb,
                std::vector<double> ub,
                std::vector<int> ty) {
-    vars_.num_col_ = n;
+    assert(lb.size()==ub.size());
+    assert(ty.empty() || ty.size()==lb.size());
+    num_col_ = lb.size();
     vars_.lower_ = std::move(lb);
     vars_.upper_ = std::move(ub);
     vars_.type_  = std::move(ty);
-    nlme_.SetCols({n,
+    nlme_.SetCols({num_col_,
                    vars_.lower_.data(),
                    vars_.upper_.data(),
                    vars_.type_.data()
@@ -97,6 +94,7 @@ class NLWPY_NLModel {
 
   /// Add variable names
   void SetColNames(std::vector<std::string> nm) {
+    assert(nm.size()==num_col_);
     var_names_=std::move(nm);
     var_names_c_.resize(var_names_.size());
     for (auto i=var_names_.size(); i--; )
@@ -107,24 +105,25 @@ class NLWPY_NLModel {
   /// Add linear constraints (all at once).
   /// Only rowwise matrix supported.
   void SetRows(
-      int nr,
       std::vector<double> rlb, std::vector<double> rub,
       NLW2_MatrixFormat format,
-      size_t nnz,
       std::vector<size_t> st,
       /// Entry index
       std::vector<int> ind,
       /// Entry value
       std::vector<double> val
       ) {
-    num_row_=nr; row_lb_=std::move(rlb); row_ub_=std::move(rub);
+    assert(rlb.size()==rub.size());
+    num_row_=rlb.size();
+    row_lb_=std::move(rlb); row_ub_=std::move(rub);
     A_={
-      nr, format,
-      nnz, std::move(st), std::move(ind), std::move(val)
+      format,
+      std::move(st), std::move(ind), std::move(val)
     };
-    nlme_.SetRows(nr, row_lb_.data(), row_ub_.data(),
+    nlme_.SetRows(num_row_,
+                  row_lb_.data(), row_ub_.data(),
                   {
-                    nr, format, nnz,
+                    num_row_, format, A_.value_.size(),
                     A_.start_.data(), A_.index_.data(),
                     A_.value_.data()
                   });
@@ -132,6 +131,7 @@ class NLWPY_NLModel {
 
   /// Add constraint names
   void SetRowNames(std::vector<std::string> nm) {
+    assert(nm.size()==num_row_);
     row_names_=std::move(nm);
     row_names_c_.resize(row_names_.size());
     for (auto i=row_names_.size(); i--; )
@@ -144,30 +144,33 @@ class NLWPY_NLModel {
   /// Coefficients: dense vector.
   void SetLinearObjective(NLW2_ObjSense sense, double c0,
                           std::vector<double> c) {
+    assert(c.size()==num_col_);
     obj_sense_=sense; obj_c0_=c0; obj_c_=std::move(c);
     nlme_.SetLinearObjective(sense, c0, obj_c_.data());
   }
 
   /// Add Q for the objective quadratic part 0.5 @ x.T @ Q @ x.
   /// Format: NLW2_HessianFormat...
-  void SetHessian(int nr,
-                  NLW2_HessianFormat format,     // TODO enum
-                  size_t nnz,
+  void SetHessian(NLW2_HessianFormat format,
                   std::vector<size_t> st,
                   /// Entry index
                   std::vector<int> ind,
                   /// Entry value
                   std::vector<double> val
                   ) {
+    assert(st.size()-1==num_col_);
     Q_format_ = format;
     Q_={
-      nr, NLW2_MatrixFormatIrrelevant,
-      nnz, std::move(st), std::move(ind), std::move(val)
+      NLW2_MatrixFormatIrrelevant,
+      std::move(st), std::move(ind), std::move(val)
     };
-    nlme_.SetHessian(format, {
-      nr, NLW2_MatrixFormatIrrelevant, nnz,
-      Q_.start_.data(), Q_.index_.data(),
-      Q_.value_.data()
+    nlme_.SetHessian(
+          format, {
+            (int)Q_.start_.size()-1,    // don't need the last element
+            NLW2_MatrixFormatIrrelevant,
+            Q_.value_.size(),
+            Q_.start_.data(), Q_.index_.data(),
+            Q_.value_.data()
     });
   }
 
@@ -177,6 +180,39 @@ class NLWPY_NLModel {
     nlme_.SetObjName(obj_name_.c_str());
   }
 
+  /// Set initial solution.
+  void SetWarmstart(
+      std::vector<int> i, std::vector<double> v) {
+    assert(i.size()==v.size());
+    ini_x_i_ = std::move(i);
+    ini_x_v_ = std::move(v);
+    nlme_.SetWarmstart(
+          {
+            (int)ini_x_i_.size(),
+            ini_x_i_.data(), ini_x_v_.data()
+          });
+  }
+
+  /// Set dual initial solution.
+  void SetDualWarmstart(
+      std::vector<int> i, std::vector<double> v) {
+    assert(i.size()==v.size());
+    ini_y_i_ = std::move(i);
+    ini_y_v_ = std::move(v);
+    nlme_.SetDualWarmstart(
+          {
+            (int)ini_y_i_.size(),
+            ini_y_i_.data(), ini_y_v_.data()
+          });
+  }
+
+  /// Add suffix.
+  /// @return true iff new suffix added (vs replaced.)
+  /// @note SOS constraints can be modeled as suffixes
+  ///   for some AMPL solvers.
+  bool AddSuffix(mp::NLSuffix suf)
+  { return nlme_.AddSuffix(std::move(suf)); }
+
   /// Get the model
   const mp::NLModel& GetModel() const { return nlme_; }
 
@@ -184,6 +220,7 @@ class NLWPY_NLModel {
   /// Store the strings/arrays to keep the memory
   std::string prob_name_ {"NLWPY_Model"};
   mp::NLModel nlme_;
+  int num_col_ {};
   NLWPY_ColData vars_ {};
   std::vector<std::string> var_names_ {};
   std::vector<const char*> var_names_c_ {};
@@ -199,6 +236,11 @@ class NLWPY_NLModel {
   NLW2_HessianFormat Q_format_ {};
   NLWPY_SparseMatrix Q_ {};
   std::string obj_name_ {"obj[1]"};
+
+  std::vector<int> ini_x_i_;
+  std::vector<double> ini_x_v_;
+  std::vector<int> ini_y_i_;
+  std::vector<double> ini_y_v_;
 };
 
 mp::NLSolution NLW2_Solve(mp::NLSolver& nls,
@@ -226,6 +268,7 @@ NLW2_HessianFormat
 
 NLW2_NLOptionsBasic
            NLW2_MakeNLOptionsBasic_Default
+NLW2_NLSuffix
 NLW2_NLModel
 NLW2_NLSolution
 NLW2_NLSolver
@@ -259,6 +302,31 @@ NLW2_NLSolver
         Use this to create default options for NLModel.
     )pbdoc");
 
+    /// NLSuffix
+    py::class_<mp::NLSuffix>(m, "NLW2_NLSuffix")
+        .def(py::init<std::string, int, std::vector<double>>())
+        .def(py::init<std::string, std::string, int, std::vector<double>>())
+        .def_readwrite("name_", &mp::NLSuffix::name_)
+        .def_readwrite("table_", &mp::NLSuffix::table_)
+        .def_readwrite("kind_", &mp::NLSuffix::kind_)
+        .def_readwrite("values_", &mp::NLSuffix::values_)
+        ;
+
+    /// NLSuffixSet
+    py::class_<mp::NLSuffixSet>(m, "NLW2_NLSuffixSet")
+        .def("Find",           // Find(): return None if not found
+             [=](mp::NLSuffixSet const& ss,
+             std::string const& name, int kind) -> py::object {
+      auto pelem = ss.Find(name, kind);
+      if (pelem) {
+        return py::cast(*pelem);
+      }
+      return py::object(py::cast(nullptr));
+    })
+        .def("begin", &mp::NLSuffixSet::begin)
+        .def("end", &mp::NLSuffixSet::end)
+        ;
+
     /// NLModel
     py::class_<NLWPY_NLModel>(m, "NLW2_NLModel")
         .def(py::init<const char*>())
@@ -279,6 +347,7 @@ NLW2_NLSolver
         .def_readwrite("obj_val_", &mp::NLSolution::obj_val_)
         .def_readwrite("x_", &mp::NLSolution::x_)
         .def_readwrite("y_", &mp::NLSolution::y_)
+        .def_readwrite("suffixes_", &mp::NLSolution::suffixes_)
         ;
 
     /// NLSolver
diff --git a/nl-writer2/nlwpy/tests/test.py b/nl-writer2/nlwpy/tests/test.py
index 5d99d6f79..ff2b04b45 100644
--- a/nl-writer2/nlwpy/tests/test.py
+++ b/nl-writer2/nlwpy/tests/test.py
@@ -44,15 +44,13 @@ class ModelBuilder:
   def GetModel(self):
     nlme = m.NLW2_NLModel(self.prob_name_)
 
-    nlme.SetCols(len(self.var_lb_),
-                 self.var_lb_, self.var_ub_, self.var_type_)
+    nlme.SetCols(self.var_lb_, self.var_ub_, self.var_type_)
     nlme.SetColNames(self.var_names_)
 
     if self.A_ is not None:
       self.A_ = csr_matrix(self.A_)
-      nlme.SetRows(len(self.row_lb_), self.row_lb_, self.row_ub_,
+      nlme.SetRows(self.row_lb_, self.row_ub_,
                    self.A_format_,
-                   self.A_.nnz,
                    self.A_.indptr, self.A_.indices, self.A_.data)
     nlme.SetRowNames(self.row_names_)
 
@@ -61,9 +59,7 @@ def GetModel(self):
 
     if self.Q_ is not None:
       self.Q_ = csr_matrix(self.Q_)
-      nlme.SetHessian(self.Q_.shape[0],
-                      self.Q_format_,
-                      self.Q_.nnz,
+      nlme.SetHessian(self.Q_format_,
                       self.Q_.indptr, self.Q_.indices, self.Q_.data)
     nlme.SetObjName(self.obj_name_)