diff --git a/lib/poek/poek/poek_pybind11.py b/lib/poek/poek/poek_pybind11.py
index 6fb387b9..1ea61fab 100644
--- a/lib/poek/poek/poek_pybind11.py
+++ b/lib/poek/poek/poek_pybind11.py
@@ -16,6 +16,20 @@ class _DummyConcreteSet(object):
     ConcreteSet = _DummyConcreteSet
 
 
+class data(object):
+    def __new__(cls, *args, **kwds):
+        #p = parameter_single(*args)
+        #if value is not None:
+        #    p.value = value
+        #return p
+        if len(args) == 0 or args[0] == 1 or type(args[0]) == str:
+            return data_(**kwds)
+        if len(args) == 1:
+            return data_(args[0], **kwds)
+        else:
+            raise RuntimeError("Data values only have one argument")
+
+
 class parameter(object):
     def __new__(cls, *args, **kwds):
         #p = parameter_single(*args)
@@ -27,7 +41,7 @@ def __new__(cls, *args, **kwds):
         if len(args) == 1:
             return parameter_(args[0], **kwds)
         else:
-            raise RuntimeError("Variables only have one argument")
+            raise RuntimeError("Parameters only have one argument")
 
 
 class variable(object):
diff --git a/lib/poek/poek/tests/test_expr.py b/lib/poek/poek/tests/test_expr.py
index 8fbc150b..33ab075f 100644
--- a/lib/poek/poek/tests/test_expr.py
+++ b/lib/poek/poek/tests/test_expr.py
@@ -43,17 +43,31 @@ def var_v():
 
 @pytest.fixture
 def param_p():
-    return variable(name="p", value=0)
+    return parameter(name="p", value=0)
 
 
 @pytest.fixture
 def param_q():
-    return variable(name="q", value=0)
+    return parameter(name="q", value=0)
 
 
 @pytest.fixture
 def param_r():
-    return variable(name="r", value=1)
+    return parameter(name="r", value=1)
+
+@pytest.fixture
+def data_d1():
+    return data(name="d1", value=0)
+
+
+@pytest.fixture
+def param_d2():
+    return data(name="d2", value=0)
+
+
+@pytest.fixture
+def param_d3():
+    return data(name="d3", value=1)
 
 
 #
@@ -96,6 +110,31 @@ def test_constraint_value():
         e = p < z
 
 
+def test_data1_value():
+    p = data(value=-1)
+    assert p.value == -1
+    p.value = 3
+    assert p.value == 3
+
+
+def test_data_float():
+    p = data(value=-1)
+    with pytest.raises(TypeError) as einfo:
+        float(p)
+    if numpy_available:
+        z = np.float32(-1)
+        p = data(value=z)
+
+
+def test_data_int():
+    p = data(value=-1)
+    with pytest.raises(TypeError) as einfo:
+        int(p)
+    if numpy_available:
+        z = np.int32(-1)
+        p = data(value=z)
+
+
 def test_param1_value():
     p = parameter(value=-1)
     assert p.value == -1
@@ -644,7 +683,7 @@ def test_paramDiff(var_a, param_p):
 
     e = a - p
     #
-    assert e.to_list() == ["+", "a", ["*", "-1.000000", "p"]]
+    assert e.to_list() == ["+", "a", ["-", "p"]]
 
     # p - a
     e = p - a
@@ -652,6 +691,20 @@ def test_paramDiff(var_a, param_p):
     assert e.to_list() == ["+", "p", ["*", "-1.000000", "a"]]
 
 
+def test_dataDiff(var_a, data_d1):
+    # a - p
+    a, d1 = var_a, data_d1
+
+    e = a - d1
+    #
+    assert e.to_list() == ["+", "a", ["-", "d1"]]
+
+    # d1 - a
+    e = d1 - a
+    #
+    assert e.to_list() == ["+", "d1", ["*", "-1.000000", "a"]]
+
+
 def test_termDiff(var_a):
     #
     # Check the structure of a simple difference with a term
@@ -764,13 +817,13 @@ def test_negation_param(param_p):
 
     e = -p
     #
-    assert e.to_list() == ["*", "-1.000000", "p"]
+    assert e.to_list() == ["-", "p"]
 
     e = -e
     #
-    # TODO: Can we detect negations of negations?
+    # TODO: Identify negations of negations
     #
-    assert e.to_list() == ["*", "1.000000", "p"]
+    assert e.to_list() == ["-", ["-", "p"]]
 
 
 def test_negation_terms(param_p, var_v):
@@ -780,10 +833,10 @@ def test_negation_terms(param_p, var_v):
     p, v = param_p, var_v
 
     e = -p * v
-    assert e.to_list() == ["*", ["*", "-1.000000", "p"], "v"]
+    assert e.to_list() == ["*", ["-", "p"], "v"]
 
     e = -e
-    assert e.to_list() == ["-", ["*", ["*", "-1.000000", "p"], "v"]]
+    assert e.to_list() == ["-", ["*", ["-", "p"], "v"]]
     #
     if True:
         e = -5 * v
@@ -833,7 +886,7 @@ def test_trivialDiff(var_a, param_p):
 
         # 0 - p
         e = 0 - p
-        assert e.to_list() == ["*", "-1.000000", "p"]
+        assert e.to_list() == ["-", "p"]
 
         # 0 - 5*a
         e = 0 - 5 * a
@@ -878,7 +931,7 @@ def test_trivialDiff(var_a, param_p):
 
         # z - p
         e = 0 - p
-        assert e.to_list() == ["*", "-1.000000", "p"]
+        assert e.to_list() == ["-", "p"]
 
         z = np.float32(0.0)
         e = a - z
diff --git a/lib/pycoek/pybind11/pycoek_pybind11.cpp b/lib/pycoek/pybind11/pycoek_pybind11.cpp
index 1cabff08..9901e12b 100644
--- a/lib/pycoek/pybind11/pycoek_pybind11.cpp
+++ b/lib/pycoek/pybind11/pycoek_pybind11.cpp
@@ -632,6 +632,13 @@ ParameterMap parameter_fn(coek::ConcreteSet& index_set, py::kwargs kwargs)
     set_kwargs_parammap(tmp, kwargs);
     return tmp;
 }
+
+DataMap data_fn(coek::ConcreteSet& index_set, py::kwargs kwargs)
+{
+    DataMap tmp(index_set);
+    set_kwargs_parammap(tmp, kwargs);
+    return tmp;
+}
 #endif
 
 VariableArray variable_fn(std::vector<int>& dimen, py::kwargs kwargs)
@@ -692,6 +699,35 @@ Parameter parameter_fn(py::kwargs kwargs)
     return tmp;
 }
 
+ParameterArray data_fn(std::vector<int>& dimen, py::kwargs kwargs)
+{
+    std::vector<size_t> _dimen(dimen.size());
+    for (size_t i = 0; i < dimen.size(); ++i) {
+        assert(dimen[i] >= 0);
+        _dimen[i] = static_cast<size_t>(dimen[i]);
+    }
+    ParameterArray tmp(_dimen);
+    set_param_kwargs(tmp, kwargs);
+    return tmp;
+}
+
+ParameterArray data_fn(int n, py::kwargs kwargs)
+{
+    if (n < 0)
+        throw std::invalid_argument("Cannot initialize data array with negative length");
+
+    ParameterArray tmp(static_cast<size_t>(n));
+    set_param_kwargs(tmp, kwargs);
+    return tmp;
+}
+
+Parameter data_fn(py::kwargs kwargs)
+{
+    Parameter tmp;
+    set_param_kwargs(tmp, kwargs);
+    return tmp;
+}
+
 #ifdef COEK_WITH_COMPACT_MODEL
 template <class T, class DTYPE>
 // coek::Expression
@@ -830,6 +866,11 @@ PYBIND11_MODULE(pycoek_pybind11, m)
           [](std::vector<int>& dimen, py::kwargs kw) { return coek::parameter_fn(dimen, kw); });
     m.def("parameter_", [](py::kwargs kw) { return coek::parameter_fn(kw); });
 
+    m.def("data_", [](int n, py::kwargs kw) { return coek::data_fn(n, kw); });
+    m.def("data_",
+          [](std::vector<int>& dimen, py::kwargs kw) { return coek::data_fn(dimen, kw); });
+    m.def("data_", [](py::kwargs kw) { return coek::data_fn(kw); });
+
 #ifdef COEK_WITH_COMPACT_MODEL
     m.def("variable_", [](coek::ConcreteSet& index_set, py::kwargs kw) {
         return coek::variable_fn(index_set, kw);
@@ -837,6 +878,9 @@ PYBIND11_MODULE(pycoek_pybind11, m)
     m.def("parameter_", [](coek::ConcreteSet& index_set, py::kwargs kw) {
         return coek::parameter_fn(index_set, kw);
     });
+    m.def("data_", [](coek::ConcreteSet& index_set, py::kwargs kw) {
+        return coek::data_fn(index_set, kw);
+    });
 #endif
     m.def("affine_expression", [](std::vector<double>& coef, std::vector<coek::Variable>& var,
                                   double offset) { return affine_expression(coef, var, offset); });
@@ -1358,6 +1402,33 @@ PYBIND11_MODULE(pycoek_pybind11, m)
         //.def("__float__", [](coek::Variable& x){return x.get_value();})
         ;
 
+    //
+    // DataArray
+    //
+    py::class_<coek::DataArray>(m, "data_array")
+        .def(py::init<>())
+        .def("__len__", [](coek::DataArray& va) { return va.size(); })
+        .def("__getitem__", [](coek::DataArray& va, int i) { return va(i); })
+        .def("__getitem__",
+             [](coek::DataArray& va, std::vector<int>& index) { return va.index(index); })
+        .def_property_readonly("name",
+                               [](coek::DataArray& x) -> py::object {
+                                   if (x.name().size() == 0)
+                                       return py::cast<std::string>("D");
+                                   else
+                                       return py::cast<std::string>(x.name());
+                               })
+        .def("is_constraint", [](const coek::DataArray&) { return false; })
+        .def("is_expression_type", [](const coek::DataArray&) { return false; })
+        .def("is_potentially_variable", [](const coek::DataArray&) { return false; })
+        .def(
+            "__iter__",
+            [](const coek::DataArray& va) {
+                typedef coek::VecKeyIterator vec_key_t;
+                return py::make_iterator(vec_key_t(), vec_key_t(va.size()));
+            },
+            py::keep_alive<0, 1>());
+
     //
     // ParameterArray
     //
@@ -1370,7 +1441,7 @@ PYBIND11_MODULE(pycoek_pybind11, m)
         .def_property_readonly("name",
                                [](coek::ParameterArray& x) -> py::object {
                                    if (x.name().size() == 0)
-                                       return py::cast<std::string>("X");
+                                       return py::cast<std::string>("P");
                                    else
                                        return py::cast<std::string>(x.name());
                                })
@@ -1420,6 +1491,18 @@ PYBIND11_MODULE(pycoek_pybind11, m)
             },
             py::keep_alive<0, 1>());
 
+    //
+    // DataMap
+    //
+#ifdef COEK_WITH_COMPACT_MODEL
+    py::class_<coek::DataMap>(m, "data_map")
+        .def("__len__", [](coek::DataMap& x) { return x.size(); })
+        .def("__getitem__",
+             [](coek::DataMap& x, py::args args) {
+                 return coek::Array_getitem<coek::DataMap, coek::Expression>(x, args);
+             });
+#endif
+
     //
     // ParameterMap
     //
@@ -2001,9 +2084,11 @@ PYBIND11_MODULE(pycoek_pybind11, m)
     //
     py::class_<coek::Model>(m, "model")
         .def(py::init<>())
+        .def("add_data", [](coek::Model& m, coek::DataArray& v) { m.add(v); })
         .def("add_parameter", [](coek::Model& m, coek::Parameter& v) { return m.add_parameter(v); })
         .def("add_parameter", [](coek::Model& m, coek::ParameterArray& v) { m.add(v); })
 #ifdef COEK_WITH_COMPACT_MODEL
+        .def("add_data", [](coek::Model& m, coek::DataMap& v) { m.add(v); })
         .def("add_parameter", [](coek::Model& m, coek::ParameterMap& v) { m.add(v); })
 #endif
         .def("add_variable_", [](coek::Model& m, coek::Variable& v) { return m.add_variable(v); })
@@ -2121,8 +2206,9 @@ PYBIND11_MODULE(pycoek_pybind11, m)
     //
     py::class_<coek::CompactModel>(m, "compact_model")
         .def(py::init<>())
-        .def("add_parameter",
-             [](coek::CompactModel& m, coek::Parameter& v) { return m.add_parameter(v); })
+        .def("add_data", [](coek::CompactModel& m, coek::DataArray& v) { m.add(v); })
+        .def("add_data", [](coek::CompactModel& m, coek::DataMap& v) { m.add(v); })
+        .def("add_parameter", [](coek::CompactModel& m, coek::Parameter& v) { return m.add_parameter(v); })
         .def("add_parameter", [](coek::CompactModel& m, coek::ParameterArray& v) { m.add(v); })
         .def("add_parameter", [](coek::CompactModel& m, coek::ParameterMap& v) { m.add(v); })
         .def("add_variable_", [](coek::CompactModel& m, coek::Variable& v) { m.add_variable(v); })