diff --git a/pycppad/adfun.py b/pycppad/adfun.py
index 179de9a..bbbe951 100644
--- a/pycppad/adfun.py
+++ b/pycppad/adfun.py
@@ -200,6 +200,49 @@ def independent(x) :
   msg = 'only implemented where x[j] is int, float, or a_float'
   raise NotImplementedError(msg)
 
+def independentdiag(x, opIndex) :
+  """
+  a_x = independent(x): create independent variable vector a_x, equal to x,
+  and start recording operations that use the class corresponding to ad( x[0] ).
+  """
+  #
+  # It would be better faster if all this type checking were done in the C++
+  #
+  if not isinstance(x, numpy.ndarray) :
+    raise NotImplementedError('independent(x): x is not of type numpy.array')
+  #
+  x0 = x[0]
+  if isinstance(x0, int) :
+    for j in range( len(x) ) :
+      if not isinstance(x[j], int) :
+        other = 'x[' + str(j) + '] is ' + type(x[j]).__name__
+        msg   = 'independent(x): mixed types x[0] is int and ' + other
+        raise NotImplementedError(msg)
+    x = numpy.array(x, dtype=int)       # incase dtype of x is object
+    return cppad_.independentdiag(x, 1, opIndex)     # level = 1
+  #
+  if isinstance(x0, float) :
+    for j in range( len(x) ) :
+      if not isinstance(x[j], float) :
+        other = 'x[' + str(j) + '] is ' + type(x[j]).__name__
+        msg   = 'independent(x): mixed types x[0] is float and ' + other
+        raise NotImplementedError(msg)
+    x = numpy.array(x, dtype=float)     # incase dtype of x is object
+    return cppad_.independentdiag(x, 1, opIndex)     # level = 1
+  #
+  if isinstance(x0, cppad_.a_float) :
+    for j in range( len(x) ) :
+      if not isinstance(x[j], cppad_.a_float) :
+        other = 'x[' + str(j) + '] is ' + type(x[j]).__name__
+        msg   = 'independent(x): mixed types x[0] is a_float and ' + other
+        raise NotImplementedError(msg)
+    return cppad_.independentdiag(x, 2, opIndex)     # level = 2
+  #
+  msg = 'independent(x): x[0] has type' + type(x0).__name__ + '\n'
+  msg = 'only implemented where x[j] is int, float, or a_float'
+  raise NotImplementedError(msg)
+
+
 class adfun_float(cppad_.adfun_float) :
   """
   Create a level zero function object (evaluates using floats).
diff --git a/pycppad/pycppad.cpp b/pycppad/pycppad.cpp
index e56f3d9..fa31fa8 100644
--- a/pycppad/pycppad.cpp
+++ b/pycppad/pycppad.cpp
@@ -646,6 +646,27 @@ namespace pycppad {
 		CppAD::Independent(a_x);
 		return vec2array(a_x);
 	}
+	array IndependentDiag(array& x_array, int level, int opIndex)
+	{	PYCPPAD_ASSERT(
+			level == 1 || level == 2,
+			"independent: level argument must be 1 or 2."
+		);
+		if( level == 1 )
+		{
+			double_vec      x(x_array);
+			AD_double_vec a_x(x.size() );
+			for(size_t j = 0; j < x.size(); j++)
+				a_x[j] = x[j];
+			CppAD::IndependentDiag(a_x, opIndex);
+			return vec2array(a_x);
+		}
+		AD_double_vec      x(x_array);
+		AD_AD_double_vec a_x(x.size() );
+		for(size_t j = 0; j < x.size(); j++)
+			a_x[j] = x[j];
+		CppAD::IndependentDiag(a_x, opIndex);
+		return vec2array(a_x);
+	}
 	// -------------------------------------------------------------
 	double double_(const AD_double& x)
 	{	return Value(x); }
@@ -694,6 +715,7 @@ BOOST_PYTHON_MODULE(cppad_)
 	array::set_module_and_type("numpy", "ndarray");
 	// --------------------------------------------------------------------
 	def("independent", pycppad::Independent);
+	def("independentdiag", pycppad::IndependentDiag);
 	def("float_",     pycppad::double_);
 	def("a_float_",   pycppad::AD_double_);
 	// documented in adfun.py