diff --git a/.codecov.yml b/.codecov.yml
index a385d754..6cfd0cfd 100644
--- a/.codecov.yml
+++ b/.codecov.yml
@@ -12,7 +12,6 @@ flags:
   project:
     enabled: yes
     target: 95%
-    threshold: 1%
     if_no_uploads: error
     if_not_found: success
     if_ci_failed: error
@@ -21,7 +20,6 @@ flags:
     default:
       enabled: yes
       target: 95%
-      threshold: 1%
       if_no_uploads: error
       if_not_found: success
       if_ci_failed: error
diff --git a/.coveragerc b/.coveragerc
index d838d373..6371c705 100644
--- a/.coveragerc
+++ b/.coveragerc
@@ -5,3 +5,9 @@ source=
 omit=
     sparse/_version.py
     sparse/tests/*
+
+[report]
+exclude_lines =
+    pragma: no cover
+    return NotImplemented
+    raise NotImplementedError
diff --git a/docs/generated/sparse.full.rst b/docs/generated/sparse.full.rst
new file mode 100644
index 00000000..9b798912
--- /dev/null
+++ b/docs/generated/sparse.full.rst
@@ -0,0 +1,6 @@
+full
+====
+
+.. currentmodule:: sparse
+
+.. autofunction:: full
\ No newline at end of file
diff --git a/docs/generated/sparse.rst b/docs/generated/sparse.rst
index 451f65b3..9104380d 100644
--- a/docs/generated/sparse.rst
+++ b/docs/generated/sparse.rst
@@ -29,6 +29,8 @@ API
 
     elemwise
 
+    full
+
     nanmax
 
     nanmin
diff --git a/docs/operations.rst b/docs/operations.rst
index a6c1a10a..fa83e927 100644
--- a/docs/operations.rst
+++ b/docs/operations.rst
@@ -20,13 +20,12 @@ results for both Numpy arrays, COO arrays, or a mix of the two:
 
    np.log(X.dot(beta.T) + 1)
 
-However some operations are not supported, like inplace operations,
-operations that implicitly cause dense structures,
-or numpy functions that are not yet implemented for sparse arrays
+However some operations are not supported, like operations that
+implicitly cause dense structures, or numpy functions that are not
+yet implemented for sparse arrays.
 
 .. code-block:: python
 
-   x += y     # inplace operations not supported
    x + 1      # operations that produce dense results not supported
    np.svd(x)  # sparse svd not implemented
 
@@ -34,7 +33,7 @@ or numpy functions that are not yet implemented for sparse arrays
 This page describes those valid operations, and their limitations.
 
 :obj:`elemwise`
-~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~
 This function allows you to apply any arbitrary broadcasting function to any number of arguments
 where the arguments can be :obj:`SparseArray` objects or :obj:`scipy.sparse.spmatrix` objects.
 For example, the following will add two arrays:
@@ -155,7 +154,9 @@ be a :obj:`scipy.sparse.spmatrix` The following operators are supported:
    * :obj:`operator.lshift` (:code:`x << y`)
    * :obj:`operator.rshift` (:code:`x >> y`)
 
-.. note:: In-place operators are not supported at this time.
+.. warning::
+   While in-place operations are supported for compatibility with Numpy,
+   they are not truly in-place, and will effectively calculate the result separately.
 
 .. _operations-elemwise:
 
diff --git a/sparse/coo/__init__.py b/sparse/coo/__init__.py
index c47dbe5e..59cd03db 100644
--- a/sparse/coo/__init__.py
+++ b/sparse/coo/__init__.py
@@ -1,4 +1,4 @@
 from .core import COO
 from .umath import elemwise
 from .common import tensordot, dot, concatenate, stack, triu, tril, where, \
-    nansum, nanprod, nanmin, nanmax, nanreduce
+    nansum, nanprod, nanmin, nanmax, nanreduce, full
diff --git a/sparse/coo/common.py b/sparse/coo/common.py
index ab612c2c..89a40e68 100644
--- a/sparse/coo/common.py
+++ b/sparse/coo/common.py
@@ -616,3 +616,34 @@ def nanreduce(x, method, identity=None, axis=None, keepdims=False, **kwargs):
     """
     arr = _replace_nan(x, method.identity if identity is None else identity)
     return arr.reduce(method, axis, keepdims, **kwargs)
+
+
+def full(coords, fill_value, dtype=None, **kwargs):
+    """
+    Return a new COO of with a given sparsity pattern, filled with fill_value.
+
+    Parameters
+    ----------
+    coords : numpy.ndarray (COO.ndim, COO.nnz)
+        Index locations of nonzero values.
+    fill_value: scalar
+        Value of array at nonzero indices.
+    dtype : data-type (optional)
+        The data-type for the array. If None, defaults to
+        `np.array(fill_value).dtype`.
+    kwargs : dict (optional)
+        Additional arguments to pass to ``COO`` constructor.
+
+    Returns
+    -------
+    COO
+        Newly constructed sparse array.
+
+    See Also
+    --------
+    numpy.full : Analogous Numpy function.
+    """
+    from .core import COO
+
+    d = np.full(np.asarray(coords).shape[1], fill_value, dtype=dtype)
+    return COO(coords, data=d, **kwargs)
diff --git a/sparse/coo/core.py b/sparse/coo/core.py
index 331e21f8..ddecdf88 100644
--- a/sparse/coo/core.py
+++ b/sparse/coo/core.py
@@ -223,8 +223,9 @@ def __init__(self, coords, data=None, shape=None, has_duplicates=True,
                 data = coords[0]
                 coords = np.stack(coords[1], axis=0)
 
-        self.data = np.asarray(data)
         self.coords = np.asarray(coords)
+        self.data = np.broadcast_to(data, self.coords.shape[-1])
+
         if self.coords.ndim == 1:
             self.coords = self.coords[None, :]
 
@@ -243,7 +244,7 @@ def __init__(self, coords, data=None, shape=None, has_duplicates=True,
         else:
             dtype = np.uint8
         self.coords = self.coords.astype(dtype)
-        assert not self.shape or len(data) == self.coords.shape[1]
+        assert not self.shape or len(self.data) == self.coords.shape[1]
 
         if not sorted:
             self._sort_indices()
@@ -676,8 +677,7 @@ def sum(self, axis=None, keepdims=False, dtype=None, out=None):
         >>> s.sum()
         25
         """
-        assert out is None
-        return self.reduce(np.add, axis=axis, keepdims=keepdims, dtype=dtype)
+        return np.add.reduce(self, out=out, axis=axis, keepdims=keepdims, dtype=dtype)
 
     def max(self, axis=None, keepdims=False, out=None):
         """
@@ -738,8 +738,7 @@ def max(self, axis=None, keepdims=False, out=None):
         >>> s.max()
         8
         """
-        assert out is None
-        return self.reduce(np.maximum, axis=axis, keepdims=keepdims)
+        return np.maximum.reduce(self, out=out, axis=axis, keepdims=keepdims)
 
     def min(self, axis=None, keepdims=False, out=None):
         """
@@ -800,8 +799,7 @@ def min(self, axis=None, keepdims=False, out=None):
         >>> s.min()
         0
         """
-        assert out is None
-        return self.reduce(np.minimum, axis=axis, keepdims=keepdims)
+        return np.minimum.reduce(self, out=out, axis=axis, keepdims=keepdims)
 
     def prod(self, axis=None, keepdims=False, dtype=None, out=None):
         """
@@ -867,8 +865,7 @@ def prod(self, axis=None, keepdims=False, dtype=None, out=None):
         >>> s.prod()
         0
         """
-        assert out is None
-        return self.reduce(np.multiply, axis=axis, keepdims=keepdims, dtype=dtype)
+        return np.multiply.reduce(self, out=out, axis=axis, keepdims=keepdims, dtype=dtype)
 
     def transpose(self, axes=None):
         """
@@ -1039,16 +1036,27 @@ def __rmatmul__(self, other):
 
     def __array_ufunc__(self, ufunc, method, *inputs, **kwargs):
         out = kwargs.pop('out', None)
-        if out is not None:
+        if out is not None and not all(isinstance(x, COO) for x in out):
             return NotImplemented
 
         if method == '__call__':
-            return elemwise(ufunc, *inputs, **kwargs)
+            result = elemwise(ufunc, *inputs, **kwargs)
         elif method == 'reduce':
-            return COO._reduce(ufunc, *inputs, **kwargs)
+            result = COO._reduce(ufunc, *inputs, **kwargs)
         else:
             return NotImplemented
 
+        if out is not None:
+            (out,) = out
+            if out.shape != result.shape:
+                raise ValueError('non-broadcastable output operand with shape %s'
+                                 'doesn\'t match the broadcast shape %s' % (out.shape, result.shape))
+
+            out._make_shallow_copy_of(result)
+            return out
+
+        return result
+
     def __array__(self, dtype=None, **kwargs):
         x = self.todense()
         if dtype and x.dtype != dtype:
@@ -1366,10 +1374,11 @@ def round(self, decimals=0, out=None):
         The :code:`out` parameter is provided just for compatibility with Numpy and isn't
         actually supported.
         """
-        assert out is None
-        return elemwise(np.round, self, decimals=decimals)
+        if out is not None and not isinstance(out, tuple):
+            out = (out,)
+        return self.__array_ufunc__(np.round, '__call__', self, decimals=decimals, out=out)
 
-    def astype(self, dtype, out=None):
+    def astype(self, dtype):
         """
         Copy of the array, cast to a specified type.
 
@@ -1385,8 +1394,7 @@ def astype(self, dtype, out=None):
         The :code:`out` parameter is provided just for compatibility with Numpy and isn't
         actually supported.
         """
-        assert out is None
-        return elemwise(np.ndarray.astype, self, dtype=dtype)
+        return self.__array_ufunc__(np.ndarray.astype, '__call__', self, dtype=dtype)
 
     def maybe_densify(self, max_size=1000, min_density=0.25):
         """
diff --git a/sparse/coo/umath.py b/sparse/coo/umath.py
index 45de7890..d5ebf589 100644
--- a/sparse/coo/umath.py
+++ b/sparse/coo/umath.py
@@ -64,8 +64,7 @@ def elemwise(func, *args, **kwargs):
         elif isinstance(arg, SparseArray) and not isinstance(arg, COO):
             args[i] = COO(arg)
         elif not isinstance(arg, COO):
-            raise ValueError("Performing this operation would produce "
-                             "a dense result: %s" % str(func))
+            return NotImplemented
 
     # Filter out scalars as they are 'baked' into the function.
     func = PositinalArgumentPartial(func, pos, posargs)
diff --git a/sparse/tests/test_coo.py b/sparse/tests/test_coo.py
index 74bcce05..af1aca50 100644
--- a/sparse/tests/test_coo.py
+++ b/sparse/tests/test_coo.py
@@ -264,6 +264,22 @@ def test_elemwise(func):
     assert_eq(func(x), fs)
 
 
+@pytest.mark.parametrize('func', [np.expm1, np.log1p, np.sin, np.tan,
+                                  np.sinh, np.tanh, np.floor, np.ceil,
+                                  np.sqrt, np.conj,
+                                  np.round, np.rint, np.conjugate,
+                                  np.conj, lambda x, out: x.round(decimals=2, out=out)])
+def test_elemwise_inplace(func):
+    s = sparse.random((2, 3, 4), density=0.5)
+    x = s.todense()
+
+    func(s, out=s)
+    func(x, out=x)
+    assert isinstance(s, COO)
+
+    assert_eq(x, s)
+
+
 @pytest.mark.parametrize('func', [
     operator.mul, operator.add, operator.sub, operator.gt,
     operator.lt, operator.ne
@@ -279,6 +295,23 @@ def test_elemwise_binary(func, shape):
     assert_eq(func(xs, ys), func(x, y))
 
 
+@pytest.mark.parametrize('func', [
+    operator.imul, operator.iadd, operator.isub
+])
+@pytest.mark.parametrize('shape', [(2,), (2, 3), (2, 3, 4), (2, 3, 4, 5)])
+def test_elemwise_binary_inplace(func, shape):
+    xs = sparse.random(shape, density=0.5)
+    ys = sparse.random(shape, density=0.5)
+
+    x = xs.todense()
+    y = ys.todense()
+
+    xs = func(xs, ys)
+    x = func(x, y)
+
+    assert_eq(xs, x)
+
+
 @pytest.mark.parametrize('func', [
     lambda x, y, z: x + y + z,
     lambda x, y, z: x * y * z,
@@ -497,7 +530,7 @@ def test_ndarray_densification_fails():
     xs = sparse.random((3, 4), density=0.5)
     y = np.random.rand(3, 4)
 
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         xs + y
 
 
@@ -624,6 +657,30 @@ def test_bitwise_binary(func, shape):
     assert_eq(func(xs, ys), func(x, y))
 
 
+@pytest.mark.parametrize('func', [
+    operator.iand, operator.ior, operator.ixor
+])
+@pytest.mark.parametrize('shape', [
+    (2,),
+    (2, 3),
+    (2, 3, 4),
+    (2, 3, 4, 5)
+])
+def test_bitwise_binary_inplace(func, shape):
+    # Small arrays need high density to have nnz entries
+    # Casting floats to int will result in all zeros, hence the * 100
+    xs = (sparse.random(shape, density=0.5) * 100).astype(np.int_)
+    ys = (sparse.random(shape, density=0.5) * 100).astype(np.int_)
+
+    x = xs.todense()
+    y = ys.todense()
+
+    xs = func(xs, ys)
+    x = func(x, y)
+
+    assert_eq(xs, x)
+
+
 @pytest.mark.parametrize('func', [
     operator.lshift, operator.rshift
 ])
@@ -649,7 +706,7 @@ def test_bitshift_binary(func, shape):
 
 
 @pytest.mark.parametrize('func', [
-    operator.and_
+    operator.ilshift, operator.irshift
 ])
 @pytest.mark.parametrize('shape', [
     (2,),
@@ -657,13 +714,37 @@ def test_bitshift_binary(func, shape):
     (2, 3, 4),
     (2, 3, 4, 5)
 ])
-def test_bitwise_scalar(func, shape):
+def test_bitshift_binary_inplace(func, shape):
     # Small arrays need high density to have nnz entries
     # Casting floats to int will result in all zeros, hence the * 100
     xs = (sparse.random(shape, density=0.5) * 100).astype(np.int_)
 
     # Can't merge into test_bitwise_binary because left/right shifting
     # with something >= 64 isn't defined.
+    ys = (sparse.random(shape, density=0.5) * 64).astype(np.int_)
+
+    x = xs.todense()
+    y = ys.todense()
+
+    xs = func(xs, ys)
+    x = func(x, y)
+
+    assert_eq(xs, x)
+
+
+@pytest.mark.parametrize('func', [
+    operator.and_
+])
+@pytest.mark.parametrize('shape', [
+    (2,),
+    (2, 3),
+    (2, 3, 4),
+    (2, 3, 4, 5)
+])
+def test_bitwise_scalar(func, shape):
+    # Small arrays need high density to have nnz entries
+    # Casting floats to int will result in all zeros, hence the * 100
+    xs = (sparse.random(shape, density=0.5) * 100).astype(np.int_)
     y = np.random.randint(100)
 
     x = xs.todense()
@@ -1376,6 +1457,17 @@ def test_two_arg_where():
         sparse.where(cs, xs)
 
 
+@pytest.mark.parametrize('func', [
+    operator.imul, operator.iadd, operator.isub
+])
+def test_inplace_invalid_shape(func):
+    xs = sparse.random((3, 4), density=0.5)
+    ys = sparse.random((2, 3, 4), density=0.5)
+
+    with pytest.raises(ValueError):
+        func(xs, ys)
+
+
 def test_nonzero():
     s = sparse.random((2, 3, 4), density=0.5)
     x = s.todense()
@@ -1395,3 +1487,22 @@ def test_argwhere():
     x = s.todense()
 
     assert_eq(np.argwhere(s), np.argwhere(x), compare_dtype=False)
+
+
+def test_full():
+    coords = np.tile(np.arange(5), (2, 1))
+
+    a = sparse.full(coords, 1)
+    e = np.diag(np.ones(5, dtype=int))
+    assert_eq(e, a)
+    assert_eq(e, COO(coords, 1))
+
+    a = sparse.full(coords, 1.0)
+    e = np.eye(5)
+    assert_eq(e, a)
+    assert_eq(e, COO(coords, 1.0))
+
+    a = sparse.full(coords.tolist(), -1.0)
+    e = -np.eye(5)
+    assert_eq(e, a)
+    assert_eq(e, COO(coords, -1.0))