Re-base the date-time branch back to the trunk.

18 files changed, 1716 insertions, 503 deletions

diff --git a/numpy/lib/arraysetops.py b/numpy/lib/arraysetops.py
index c5e7822f2..b8ae9a9f3 100644
--- a/numpy/lib/arraysetops.py
+++ b/numpy/lib/arraysetops.py
@@ -3,40 +3,36 @@ Set operations for 1D numeric arrays based on sorting.
 
 :Contains:
   ediff1d,
-  unique1d,
+  unique,
   intersect1d,
-  intersect1d_nu,
   setxor1d,
-  setmember1d,
-  setmember1d_nu,
+  in1d,
   union1d,
   setdiff1d
 
+:Deprecated:
+  unique1d,
+  intersect1d_nu,
+  setmember1d
+
 :Notes:
 
-All functions work best with integer numerical arrays on input (e.g. indices).
-For floating point arrays, innacurate results may appear due to usual round-off
+For floating point arrays, inaccurate results may appear due to usual round-off
 and floating point comparison issues.
 
-Except unique1d, union1d and intersect1d_nu, all functions expect inputs with
-unique elements. Speed could be gained in some operations by an implementaion of
-sort(), that can provide directly the permutation vectors, avoiding thus calls
-to argsort().
+Speed could be gained in some operations by an implementation of
+sort(), that can provide directly the permutation vectors, avoiding
+thus calls to argsort().
 
-Run _test_unique1d_speed() to compare performance of numpy.unique1d() and
-numpy.unique() - it should be the same.
-
-To do: Optionally return indices analogously to unique1d for all functions.
-
-created:       01.11.2005
-last revision: 07.01.2007
+To do: Optionally return indices analogously to unique for all functions.
 
 :Author: Robert Cimrman
 """
 __all__ = ['ediff1d', 'unique1d', 'intersect1d', 'intersect1d_nu', 'setxor1d',
-           'setmember1d', 'setmember1d_nu', 'union1d', 'setdiff1d']
+           'setmember1d', 'union1d', 'setdiff1d', 'unique', 'in1d']
 
 import numpy as np
+from numpy.lib.utils import deprecate_with_doc
 
 def ediff1d(ary, to_end=None, to_begin=None):
     """
@@ -50,7 +46,7 @@ def ediff1d(ary, to_end=None, to_begin=None):
         If provided, this number will be tacked onto the end of the returned
         differences.
     to_begin : number, optional
-        If provided, this number will be taked onto the beginning of the
+        If provided, this number will be tacked onto the beginning of the
         returned differences.
 
     Returns
@@ -73,26 +69,26 @@ def ediff1d(ary, to_end=None, to_begin=None):
         arrays.append(to_end)
 
     if len(arrays) != 1:
-        # We'll save ourselves a copy of a potentially large array in the common
-        # case where neither to_begin or to_end was given.
+        # We'll save ourselves a copy of a potentially large array in
+        # the common case where neither to_begin or to_end was given.
         ed = np.hstack(arrays)
 
     return ed
 
-def unique1d(ar1, return_index=False, return_inverse=False):
+def unique(ar, return_index=False, return_inverse=False):
     """
     Find the unique elements of an array.
 
     Parameters
     ----------
-    ar1 : array_like
+    ar : array_like
         This array will be flattened if it is not already 1-D.
     return_index : bool, optional
         If True, also return the indices against `ar1` that result in the
         unique array.
     return_inverse : bool, optional
         If True, also return the indices against the unique array that
-        result in `ar1`.
+        result in `ar`.
 
     Returns
     -------
@@ -112,17 +108,17 @@ def unique1d(ar1, return_index=False, return_inverse=False):
 
     Examples
     --------
-    >>> np.unique1d([1, 1, 2, 2, 3, 3])
+    >>> np.unique([1, 1, 2, 2, 3, 3])
     array([1, 2, 3])
     >>> a = np.array([[1, 1], [2, 3]])
-    >>> np.unique1d(a)
+    >>> np.unique(a)
     array([1, 2, 3])
 
     Reconstruct the input from unique values:
 
-    >>> np.unique1d([1,2,6,4,2,3,2], return_index=True)
+    >>> np.unique([1,2,6,4,2,3,2], return_index=True)
     >>> x = [1,2,6,4,2,3,2]
-    >>> u, i = np.unique1d(x, return_inverse=True)
+    >>> u, i = np.unique(x, return_inverse=True)
     >>> u
     array([1, 2, 3, 4, 6])
     >>> i
@@ -131,14 +127,15 @@ def unique1d(ar1, return_index=False, return_inverse=False):
     [1, 2, 6, 4, 2, 3, 2]
 
     """
-    if return_index:
-        import warnings
-        warnings.warn("The order of the output arguments for "
-                      "`return_index` has changed.  Before, "
-                      "the output was (indices, unique_arr), but "
-                      "has now been reversed to be more consistent.")
+    try:
+        ar = ar.flatten()
+    except AttributeError:
+        if not return_inverse and not return_index:
+            items = sorted(set(ar))
+            return np.asarray(items)
+        else:
+            ar = np.asanyarray(ar).flatten()
 
-    ar = np.asanyarray(ar1).flatten()
     if ar.size == 0:
         if return_inverse and return_index:
             return ar, np.empty(0, np.bool), np.empty(0, np.bool)
@@ -166,44 +163,18 @@ def unique1d(ar1, return_index=False, return_inverse=False):
         flag = np.concatenate(([True], ar[1:] != ar[:-1]))
         return ar[flag]
 
-def intersect1d(ar1, ar2):
-    """
-    Intersection returning repeated or unique elements common to both arrays.
-
-    Parameters
-    ----------
-    ar1,ar2 : array_like
-        Input arrays.
-
-    Returns
-    -------
-    out : ndarray, shape(N,)
-        Sorted 1D array of common elements with repeating elements.
-
-    See Also
-    --------
-    intersect1d_nu : Returns only unique common elements.
-    numpy.lib.arraysetops : Module with a number of other functions for
-                            performing set operations on arrays.
-
-    Examples
-    --------
-    >>> np.intersect1d([1,3,3],[3,1,1])
-    array([1, 1, 3, 3])
-
-    """
-    aux = np.concatenate((ar1,ar2))
-    aux.sort()
-    return aux[aux[1:] == aux[:-1]]
 
-def intersect1d_nu(ar1, ar2):
+def intersect1d(ar1, ar2, assume_unique=False):
     """
     Intersection returning unique elements common to both arrays.
 
     Parameters
     ----------
-    ar1,ar2 : array_like
+    ar1, ar2 : array_like
         Input arrays.
+    assume_unique : bool
+        If True, the input arrays are both assumed to be unique, which
+        can speed up the calculation.  Default is False.
 
     Returns
     -------
@@ -212,34 +183,34 @@ def intersect1d_nu(ar1, ar2):
 
     See Also
     --------
-    intersect1d : Returns repeated or unique common elements.
     numpy.lib.arraysetops : Module with a number of other functions for
                             performing set operations on arrays.
 
     Examples
     --------
-    >>> np.intersect1d_nu([1,3,3],[3,1,1])
+    >>> np.intersect1d([1,3,3], [3,1,1])
     array([1, 3])
 
     """
-    # Might be faster than unique1d( intersect1d( ar1, ar2 ) )?
-    aux = np.concatenate((unique1d(ar1), unique1d(ar2)))
+    if not assume_unique:
+        # Might be faster than unique( intersect1d( ar1, ar2 ) )?
+        ar1 = unique(ar1)
+        ar2 = unique(ar2)
+    aux = np.concatenate( (ar1, ar2) )
     aux.sort()
     return aux[aux[1:] == aux[:-1]]
 
-def setxor1d(ar1, ar2):
+def setxor1d(ar1, ar2, assume_unique=False):
     """
-    Set exclusive-or of 1D arrays with unique elements.
-
-    Use unique1d() to generate arrays with only unique elements to use as
-    inputs to this function.
+    Set exclusive-or of two 1D arrays.
 
     Parameters
     ----------
-    ar1 : array_like
-        Input array.
-    ar2 : array_like
-        Input array.
+    ar1, ar2 : array_like
+        Input arrays.
+    assume_unique : bool
+        If True, the input arrays are both assumed to be unique, which
+        can speed up the calculation.  Default is False.
 
     Returns
     -------
@@ -252,7 +223,11 @@ def setxor1d(ar1, ar2):
                             performing set operations on arrays.
 
     """
-    aux = np.concatenate((ar1, ar2))
+    if not assume_unique:
+        ar1 = unique(ar1)
+        ar2 = unique(ar2)
+
+    aux = np.concatenate( (ar1, ar2) )
     if aux.size == 0:
         return aux
 
@@ -263,98 +238,68 @@ def setxor1d(ar1, ar2):
     flag2 = flag[1:] == flag[:-1]
     return aux[flag2]
 
-def setmember1d(ar1, ar2):
+def in1d(ar1, ar2, assume_unique=False):
     """
-    Return a boolean array set True where first element is in second array.
-
-    Boolean array is the shape of `ar1` containing True where the elements
-    of `ar1` are in `ar2` and False otherwise.
+    Test whether each element of an array is also present in a second array.
 
-    Use unique1d() to generate arrays with only unique elements to use as
-    inputs to this function.
+    Returns a boolean array the same length as `ar1` that is True
+    where an element of `ar1` is in `ar2` and False otherwise.
 
     Parameters
     ----------
-    ar1 : array_like
-        Input array.
-    ar2 : array_like
-        Input array.
+    ar1, ar2 : array_like
+        Input arrays.
+    assume_unique : bool
+        If True, the input arrays are both assumed to be unique, which
+        can speed up the calculation.  Default is False.
 
     Returns
     -------
     mask : ndarray, bool
         The values `ar1[mask]` are in `ar2`.
 
-
     See Also
     --------
-    setmember1d_nu : Works for arrays with non-unique elements.
     numpy.lib.arraysetops : Module with a number of other functions for
                             performing set operations on arrays.
 
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
     Examples
     --------
     >>> test = np.arange(5)
     >>> states = [0, 2]
-    >>> mask = np.setmember1d(test,states)
+    >>> mask = np.setmember1d(test, states)
     >>> mask
     array([ True, False,  True, False, False], dtype=bool)
     >>> test[mask]
     array([0, 2])
 
     """
-    # We need this to be a stable sort, so always use 'mergesort' here. The
-    # values from the first array should always come before the values from the
-    # second array.
-    ar = np.concatenate( (ar1, ar2 ) )
+    if not assume_unique:
+        ar1, rev_idx = np.unique(ar1, return_inverse=True)
+        ar2 = np.unique(ar2)
+
+    ar = np.concatenate( (ar1, ar2) )
+    # We need this to be a stable sort, so always use 'mergesort'
+    # here. The values from the first array should always come before
+    # the values from the second array.
     order = ar.argsort(kind='mergesort')
     sar = ar[order]
     equal_adj = (sar[1:] == sar[:-1])
     flag = np.concatenate( (equal_adj, [False] ) )
-
     indx = order.argsort(kind='mergesort')[:len( ar1 )]
-    return flag[indx]
-
-def setmember1d_nu(ar1, ar2):
-    """
-    Return a boolean array set True where first element is in second array.
-
-    Boolean array is the shape of `ar1` containing True where the elements
-    of `ar1` are in `ar2` and False otherwise.
-
-    Unlike setmember1d(), this version works also for arrays with duplicate
-    values.  It uses setmember1d() internally. For arrays with unique
-    entries it is slower than calling setmember1d() directly.
-
-    Parameters
-    ----------
-    ar1 : array_like
-        Input array.
-    ar2 : array_like
-        Input array.
 
-    Returns
-    -------
-    mask : ndarray, bool
-        The values `ar1[mask]` are in `ar2`.
-
-    See Also
-    --------
-    setmember1d : Faster for arrays with unique elements.
-    numpy.lib.arraysetops : Module with a number of other functions for
-                            performing set operations on arrays.
-
-    """
-    unique_ar1, rev_idx = np.unique1d(ar1, return_inverse=True)
-    mask = np.setmember1d(unique_ar1, np.unique1d(ar2))
-    return mask[rev_idx]
+    if assume_unique:
+        return flag[indx]
+    else:
+        return flag[indx][rev_idx]
 
 def union1d(ar1, ar2):
     """
-    Union of 1D arrays with unique elements.
-
-    Use unique1d() to generate arrays with only unique elements to use as
-    inputs to this function.
+    Union of two 1D arrays.
 
     Parameters
     ----------
@@ -374,14 +319,11 @@ def union1d(ar1, ar2):
                             performing set operations on arrays.
 
     """
-    return unique1d( np.concatenate( (ar1, ar2) ) )
+    return unique( np.concatenate( (ar1, ar2) ) )
 
-def setdiff1d(ar1, ar2):
+def setdiff1d(ar1, ar2, assume_unique=False):
     """
-    Set difference of 1D arrays with unique elements.
-
-    Use unique1d() to generate arrays with only unique elements to use as
-    inputs to this function.
+    Set difference of two 1D arrays.
 
     Parameters
     ----------
@@ -389,6 +331,9 @@ def setdiff1d(ar1, ar2):
         Input array.
     ar2 : array_like
         Input comparison array.
+    assume_unique : bool
+        If True, the input arrays are both assumed to be unique, which
+        can speed up the calculation.  Default is False.
 
     Returns
     -------
@@ -401,8 +346,80 @@ def setdiff1d(ar1, ar2):
                             performing set operations on arrays.
 
     """
-    aux = setmember1d(ar1,ar2)
+    if not assume_unique:
+        ar1 = unique(ar1)
+        ar2 = unique(ar2)
+    aux = in1d(ar1, ar2, assume_unique=True)
     if aux.size == 0:
         return aux
     else:
         return np.asarray(ar1)[aux == 0]
+
+@deprecate_with_doc('')
+def unique1d(ar1, return_index=False, return_inverse=False):
+    """
+    This function is deprecated. Use unique() instead.
+    """
+    if return_index:
+        import warnings
+        warnings.warn("The order of the output arguments for "
+                      "`return_index` has changed.  Before, "
+                      "the output was (indices, unique_arr), but "
+                      "has now been reversed to be more consistent.")
+
+    ar = np.asanyarray(ar1).flatten()
+    if ar.size == 0:
+        if return_inverse and return_index:
+            return ar, np.empty(0, np.bool), np.empty(0, np.bool)
+        elif return_inverse or return_index:
+            return ar, np.empty(0, np.bool)
+        else:
+            return ar
+
+    if return_inverse or return_index:
+        perm = ar.argsort()
+        aux = ar[perm]
+        flag = np.concatenate(([True], aux[1:] != aux[:-1]))
+        if return_inverse:
+            iflag = np.cumsum(flag) - 1
+            iperm = perm.argsort()
+            if return_index:
+                return aux[flag], perm[flag], iflag[iperm]
+            else:
+                return aux[flag], iflag[iperm]
+        else:
+            return aux[flag], perm[flag]
+
+    else:
+        ar.sort()
+        flag = np.concatenate(([True], ar[1:] != ar[:-1]))
+        return ar[flag]
+
+@deprecate_with_doc('')
+def intersect1d_nu(ar1, ar2):
+    """
+    This function is deprecated.  Use intersect1d()
+    instead.
+    """
+    # Might be faster than unique1d( intersect1d( ar1, ar2 ) )?
+    aux = np.concatenate((unique1d(ar1), unique1d(ar2)))
+    aux.sort()
+    return aux[aux[1:] == aux[:-1]]
+
+@deprecate_with_doc('')
+def setmember1d(ar1, ar2):
+    """
+    This function is deprecated.  Use in1d(assume_unique=True)
+    instead.
+    """
+    # We need this to be a stable sort, so always use 'mergesort' here. The
+    # values from the first array should always come before the values from the
+    # second array.
+    ar = np.concatenate( (ar1, ar2 ) )
+    order = ar.argsort(kind='mergesort')
+    sar = ar[order]
+    equal_adj = (sar[1:] == sar[:-1])
+    flag = np.concatenate( (equal_adj, [False] ) )
+
+    indx = order.argsort(kind='mergesort')[:len( ar1 )]
+    return flag[indx]
diff --git a/numpy/lib/financial.py b/numpy/lib/financial.py
index 0cef1c4d2..503d43647 100644
--- a/numpy/lib/financial.py
+++ b/numpy/lib/financial.py
@@ -28,6 +28,18 @@ def fv(rate, nper, pmt, pv, when='end'):
     """
     Compute the future value.
 
+    Given:
+     * a present value, `pv`
+     * an interest `rate` compounded once per period, of which
+       there are
+     * `nper` total
+     * a (fixed) payment, `pmt`, paid either
+     * at the beginning (`when` = {'begin', 1}) or the end
+       (`when` = {'end', 0}) of each period
+
+    Return:
+       the value at the end of the `nper` periods
+
     Parameters
     ----------
     rate : scalar or array_like of shape(M, )
@@ -61,6 +73,17 @@ def fv(rate, nper, pmt, pv, when='end'):
 
      fv + pv + pmt * nper == 0
 
+    References
+    ----------
+    .. [WRW] Wheeler, D. A., E. Rathke, and R. Weir (Eds.) (2009, May).
+       Open Document Format for Office Applications (OpenDocument)v1.2,
+       Part 2: Recalculated Formula (OpenFormula) Format - Annotated Version,
+       Pre-Draft 12. Organization for the Advancement of Structured Information
+       Standards (OASIS). Billerica, MA, USA. [ODT Document].
+       Available:
+       http://www.oasis-open.org/committees/documents.php?wg_abbrev=office-formula
+       OpenDocument-formula-20090508.odt
+
     Examples
     --------
     What is the future value after 10 years of saving $100 now, with
@@ -94,6 +117,19 @@ def pmt(rate, nper, pv, fv=0, when='end'):
     """
     Compute the payment against loan principal plus interest.
 
+    Given:
+     * a present value, `pv` (e.g., an amount borrowed)
+     * a future value, `fv` (e.g., 0)
+     * an interest `rate` compounded once per period, of which
+       there are
+     * `nper` total
+     * and (optional) specification of whether payment is made
+       at the beginning (`when` = {'begin', 1}) or the end
+       (`when` = {'end', 0}) of each period
+
+    Return:
+       the (fixed) periodic payment.
+
     Parameters
     ----------
     rate : array_like
@@ -102,8 +138,8 @@ def pmt(rate, nper, pv, fv=0, when='end'):
         Number of compounding periods
     pv : array_like
         Present value
-    fv : array_like
-        Future value
+    fv : array_like (optional)
+        Future value (default = 0)
     when : {{'begin', 1}, {'end', 0}}, {string, int}
         When payments are due ('begin' (1) or 'end' (0))
 
@@ -117,7 +153,7 @@ def pmt(rate, nper, pv, fv=0, when='end'):
 
     Notes
     -----
-    The payment ``pmt`` is computed by solving the equation::
+    The payment is computed by solving the equation::
 
      fv +
      pv*(1 + rate)**nper +
@@ -127,16 +163,37 @@ def pmt(rate, nper, pv, fv=0, when='end'):
 
       fv + pv + pmt * nper == 0
 
+    for ``pmt``.
+
+    Note that computing a monthly mortgage payment is only
+    one use for this function.  For example, pmt returns the
+    periodic deposit one must make to achieve a specified
+    future balance given an initial deposit, a fixed,
+    periodically compounded interest rate, and the total
+    number of periods.
+
+    References
+    ----------
+    .. [WRW] Wheeler, D. A., E. Rathke, and R. Weir (Eds.) (2009, May).
+       Open Document Format for Office Applications (OpenDocument)v1.2,
+       Part 2: Recalculated Formula (OpenFormula) Format - Annotated Version,
+       Pre-Draft 12. Organization for the Advancement of Structured Information
+       Standards (OASIS). Billerica, MA, USA. [ODT Document].
+       Available:
+       http://www.oasis-open.org/committees/documents.php?wg_abbrev=office-formula
+       OpenDocument-formula-20090508.odt
+
     Examples
     --------
-    What would the monthly payment need to be to pay off a $200,000 loan in 15
+    What is the monthly payment needed to pay off a $200,000 loan in 15
     years at an annual interest rate of 7.5%?
 
     >>> np.pmt(0.075/12, 12*15, 200000)
     -1854.0247200054619
 
-    In order to pay-off (i.e. have a future-value of 0) the $200,000 obtained
-    today, a monthly payment of $1,854.02 would be required.
+    In order to pay-off (i.e., have a future-value of 0) the $200,000 obtained
+    today, a monthly payment of $1,854.02 would be required.  Note that this
+    example illustrates usage of `fv` having a default value of 0.
 
     """
     when = _convert_when(when)
@@ -282,6 +339,18 @@ def pv(rate, nper, pmt, fv=0.0, when='end'):
     """
     Compute the present value.
 
+    Given:
+     * a future value, `fv`
+     * an interest `rate` compounded once per period, of which
+       there are
+     * `nper` total
+     * a (fixed) payment, `pmt`, paid either
+     * at the beginning (`when` = {'begin', 1}) or the end
+       (`when` = {'end', 0}) of each period
+
+    Return:
+       the value now
+
     Parameters
     ----------
     rate : array_like
@@ -302,7 +371,7 @@ def pv(rate, nper, pmt, fv=0.0, when='end'):
 
     Notes
     -----
-    The present value ``pv`` is computed by solving the equation::
+    The present value is computed by solving the equation::
 
      fv +
      pv*(1 + rate)**nper +
@@ -312,6 +381,45 @@ def pv(rate, nper, pmt, fv=0.0, when='end'):
 
      fv + pv + pmt * nper = 0
 
+    for `pv`, which is then returned.
+
+    References
+    ----------
+    .. [WRW] Wheeler, D. A., E. Rathke, and R. Weir (Eds.) (2009, May).
+       Open Document Format for Office Applications (OpenDocument)v1.2,
+       Part 2: Recalculated Formula (OpenFormula) Format - Annotated Version,
+       Pre-Draft 12. Organization for the Advancement of Structured Information
+       Standards (OASIS). Billerica, MA, USA. [ODT Document].
+       Available:
+       http://www.oasis-open.org/committees/documents.php?wg_abbrev=office-formula
+       OpenDocument-formula-20090508.odt
+
+    Examples
+    --------
+    What is the present value (e.g., the initial investment)
+    of an investment that needs to total $15692.93
+    after 10 years of saving $100 every month?  Assume the
+    interest rate is 5% (annually) compounded monthly.
+
+    >>> np.pv(0.05/12, 10*12, -100, 15692.93)
+    -100.00067131625819
+
+    By convention, the negative sign represents cash flow out
+    (i.e., money not available today).  Thus, to end up with
+    $15,692.93 in 10 years saving $100 a month at 5% annual
+    interest, one's initial deposit should also be $100.
+
+    If any input is array_like, ``pv`` returns an array of equal shape.
+    Let's compare different interest rates in the example above:
+
+    >>> a = np.array((0.05, 0.04, 0.03))/12
+    >>> np.pv(a, 10*12, -100, 15692.93)
+    array([ -100.00067132,  -649.26771385, -1273.78633713])
+
+    So, to end up with the same $15692.93 under the same $100 per month
+    "savings plan," for annual interest rates of 4% and 3%, one would
+    need initial investments of $649.27 and $1273.79, respectively.
+
     """
     when = _convert_when(when)
     rate, nper, pmt, fv, when = map(np.asarray, [rate, nper, pmt, fv, when])
@@ -391,24 +499,54 @@ def irr(values):
     """
     Return the Internal Rate of Return (IRR).
 
-    This is the rate of return that gives a net present value of 0.0.
+    This is the "average" periodically compounded rate of return
+    that gives a net present value of 0.0; for a more complete explanation,
+    see Notes below.
 
     Parameters
     ----------
     values : array_like, shape(N,)
-        Input cash flows per time period.  At least the first value would be
-        negative to represent the investment in the project.
+        Input cash flows per time period.  By convention, net "deposits"
+        are negative and net "withdrawals" are positive.  Thus, for example,
+        at least the first element of `values`, which represents the initial
+        investment, will typically be negative.
 
     Returns
     -------
     out : float
         Internal Rate of Return for periodic input values.
 
+    Notes
+    -----
+    The IRR is perhaps best understood through an example (illustrated
+    using np.irr in the Examples section below).  Suppose one invests
+    100 units and then makes the following withdrawals at regular
+    (fixed) intervals: 39, 59, 55, 20.  Assuming the ending value is 0,
+    one's 100 unit investment yields 173 units; however, due to the
+    combination of compounding and the periodic withdrawals, the
+    "average" rate of return is neither simply 0.73/4 nor (1.73)^0.25-1.
+    Rather, it is the solution (for :math:`r`) of the equation:
+
+    .. math:: -100 + \\frac{39}{1+r} + \\frac{59}{(1+r)^2}
+     + \\frac{55}{(1+r)^3} + \\frac{20}{(1+r)^4} = 0
+
+    In general, for `values` :math:`= [v_0, v_1, ... v_M]`,
+    irr is the solution of the equation: [G]_
+
+    .. math:: \\sum_{t=0}^M{\\frac{v_t}{(1+irr)^{t}}} = 0
+
+    References
+    ----------
+    .. [G] L. J. Gitman, "Principles of Managerial Finance, Brief," 3rd ed.,
+       Addison-Wesley, 2003, pg. 348.
+
     Examples
     --------
     >>> np.irr([-100, 39, 59, 55, 20])
     0.2809484211599611
 
+    (Compare with the Example given for numpy.lib.financial.npv)
+
     """
     res = np.roots(values[::-1])
     # Find the root(s) between 0 and 1
@@ -430,8 +568,14 @@ def npv(rate, values):
     rate : scalar
         The discount rate.
     values : array_like, shape(M, )
-        The values of the time series of cash flows.  Must be the same
-        increment as the `rate`.
+        The values of the time series of cash flows.  The (fixed) time
+        interval between cash flow "events" must be the same as that
+        for which `rate` is given (i.e., if `rate` is per year, then
+        precisely a year is understood to elapse between each cash flow
+        event).  By convention, investments or "deposits" are negative,
+        income or "withdrawals" are positive; `values` must begin with
+        the initial investment, thus `values[0]` will typically be
+        negative.
 
     Returns
     -------
@@ -440,9 +584,21 @@ def npv(rate, values):
 
     Notes
     -----
-    Returns the result of:
+    Returns the result of: [G]_
+
+    .. math :: \\sum_{t=0}^M{\\frac{values_t}{(1+rate)^{t}}}
+
+    References
+    ----------
+    .. [G] L. J. Gitman, "Principles of Managerial Finance, Brief," 3rd ed.,
+       Addison-Wesley, 2003, pg. 346.
 
-    .. math :: \\sum_{t=1}^M{\\frac{values_t}{(1+rate)^{t}}}
+    Examples
+    --------
+    >>> np.npv(0.281,[-100, 39, 59, 55, 20])
+    -0.0066187288356340801
+
+    (Compare with the Example given for numpy.lib.financial.irr)
 
     """
     values = np.asarray(values)
@@ -456,7 +612,7 @@ def mirr(values, finance_rate, reinvest_rate):
     ----------
     values : array_like
         Cash flows (must contain at least one positive and one negative value)
-        or nan is returned.
+        or nan is returned.  The first value is considered a sunk cost at time zero.
     finance_rate : scalar
         Interest rate paid on the cash flows
     reinvest_rate : scalar
@@ -469,13 +625,13 @@ def mirr(values, finance_rate, reinvest_rate):
 
     """
 
-    values = np.asarray(values)
+    values = np.asarray(values, dtype=np.double)
+    n = values.size
     pos = values > 0
     neg = values < 0
-    if not (pos.size > 0 and neg.size > 0):
+    if not (pos.any() and neg.any()):
         return np.nan
+    numer = np.abs(npv(reinvest_rate, values*pos))*(1 + reinvest_rate)
+    denom = np.abs(npv(finance_rate, values*neg))*(1 + finance_rate)
+    return (numer/denom)**(1.0/(n - 1))*(1 + reinvest_rate) - 1
 
-    n = pos.size + neg.size
-    numer = -npv(reinvest_rate, values[pos])*((1+reinvest_rate)**n)
-    denom = npv(finance_rate, values[neg])*(1+finance_rate)
-    return (numer / denom)**(1.0/(n-1)) - 1
diff --git a/numpy/lib/function_base.py b/numpy/lib/function_base.py
index b493801df..663c3d2ef 100644
--- a/numpy/lib/function_base.py
+++ b/numpy/lib/function_base.py
@@ -3,7 +3,7 @@ __all__ = ['logspace', 'linspace',
            'select', 'piecewise', 'trim_zeros',
            'copy', 'iterable',
            'diff', 'gradient', 'angle', 'unwrap', 'sort_complex', 'disp',
-           'unique', 'extract', 'place', 'nansum', 'nanmax', 'nanargmax',
+           'extract', 'place', 'nansum', 'nanmax', 'nanargmax',
            'nanargmin', 'nanmin', 'vectorize', 'asarray_chkfinite', 'average',
            'histogram', 'histogramdd', 'bincount', 'digitize', 'cov',
            'corrcoef', 'msort', 'median', 'sinc', 'hamming', 'hanning',
@@ -28,6 +28,7 @@ from numpy.lib.twodim_base import diag
 from _compiled_base import _insert, add_docstring
 from _compiled_base import digitize, bincount, interp as compiled_interp
 from arraysetops import setdiff1d
+from utils import deprecate_with_doc
 import numpy as np
 
 #end Fernando's utilities
@@ -377,11 +378,11 @@ def histogram(a, bins=10, range=None, normed=False, weights=None, new=None):
 
         n = np.diff(n)
 
-        if normed is False:
-            return n, bins
-        elif normed is True:
+        if normed:
             db = array(np.diff(bins), float)
             return n/(n*db).sum(), bins
+        else:
+            return n, bins
 
 
 def histogramdd(sample, bins=10, range=None, normed=False, weights=None):
@@ -717,9 +718,9 @@ def piecewise(x, condlist, funclist, *args, **kw):
 
     Parameters
     ----------
-    x : (N,) ndarray
+    x : ndarray
         The input domain.
-    condlist : list of M (N,)-shaped boolean arrays
+    condlist : list of bool arrays
         Each boolean array corresponds to a function in `funclist`.  Wherever
         `condlist[i]` is True, `funclist[i](x)` is used as the output value.
 
@@ -727,24 +728,24 @@ def piecewise(x, condlist, funclist, *args, **kw):
         and should therefore be of the same shape as `x`.
 
         The length of `condlist` must correspond to that of `funclist`.
-        If one extra function is given, i.e. if the length of `funclist` is
-        M+1, then that extra function is the default value, used wherever
-        all conditions are false.
-    funclist : list of M or M+1 callables, f(x,*args,**kw), or values
+        If one extra function is given, i.e. if
+        ``len(funclist) - len(condlist) == 1``, then that extra function
+        is the default value, used wherever all conditions are false.
+    funclist : list of callables, f(x,*args,**kw), or scalars
         Each function is evaluated over `x` wherever its corresponding
         condition is True.  It should take an array as input and give an array
         or a scalar value as output.  If, instead of a callable,
-        a value is provided then a constant function (``lambda x: value``) is
+        a scalar is provided then a constant function (``lambda x: scalar``) is
         assumed.
     args : tuple, optional
         Any further arguments given to `piecewise` are passed to the functions
-        upon execution, i.e., if called ``piecewise(...,...,1,'a')``, then
-        each function is called as ``f(x,1,'a')``.
-    kw : dictionary, optional
+        upon execution, i.e., if called ``piecewise(..., ..., 1, 'a')``, then
+        each function is called as ``f(x, 1, 'a')``.
+    kw : dict, optional
         Keyword arguments used in calling `piecewise` are passed to the
         functions upon execution, i.e., if called
-        ``piecewise(...,...,lambda=1)``, then each function is called as
-        ``f(x,lambda=1)``.
+        ``piecewise(..., ..., lambda=1)``, then each function is called as
+        ``f(x, lambda=1)``.
 
     Returns
     -------
@@ -754,6 +755,11 @@ def piecewise(x, condlist, funclist, *args, **kw):
         as defined by the boolean arrays in `condlist`.  Portions not covered
         by any condition have undefined values.
 
+
+    See Also
+    --------
+    choose, select, where
+
     Notes
     -----
     This is similar to choose or select, except that functions are
@@ -773,8 +779,8 @@ def piecewise(x, condlist, funclist, *args, **kw):
     --------
     Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``.
 
-    >>> x = np.arange(6) - 2.5 # x runs from -2.5 to 2.5 in steps of 1
-    >>> np.piecewise(x, [x < 0, x >= 0.5], [-1,1])
+    >>> x = np.arange(6) - 2.5
+    >>> np.piecewise(x, [x < 0, x >= 0], [-1, 1])
     array([-1., -1., -1.,  1.,  1.,  1.])
 
     Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for
@@ -836,39 +842,35 @@ def select(condlist, choicelist, default=0):
 
     Parameters
     ----------
-    condlist : list of N boolean arrays of length M
-            The conditions C_0 through C_(N-1) which determine
-            from which vector the output elements are taken.
-    choicelist : list of N arrays of length M
-            Th vectors V_0 through V_(N-1), from which the output
-            elements are chosen.
+    condlist : list of bool ndarrays
+        The list of conditions which determine from which array in `choicelist`
+        the output elements are taken. When multiple conditions are satisfied,
+        the first one encountered in `condlist` is used.
+    choicelist : list of ndarrays
+        The list of arrays from which the output elements are taken. It has
+        to be of the same length as `condlist`.
+    default : scalar, optional
+        The element inserted in `output` when all conditions evaluate to False.
 
     Returns
     -------
-    output : 1-dimensional array of length M
-            The output at position m is the m-th element of the first
-            vector V_n for which C_n[m] is non-zero.  Note that the
-            output depends on the order of conditions, since the
-            first satisfied condition is used.
-
-    Notes
-    -----
-    Equivalent to:
-    ::
+    output : ndarray
+        The output at position m is the m-th element of the array in
+        `choicelist` where the m-th element of the corresponding array in
+        `condlist` is True.
 
-                output = []
-                for m in range(M):
-                    output += [V[m] for V,C in zip(values,cond) if C[m]]
-                              or [default]
+    See Also
+    --------
+    where : Return elements from one of two arrays depending on condition.
+    take, choose, compress, diag, diagonal
 
     Examples
     --------
-    >>> t = np.arange(10)
-    >>> s = np.arange(10)*100
-    >>> condlist = [t == 4, t > 5]
-    >>> choicelist = [s, t]
+    >>> x = np.arange(10)
+    >>> condlist = [x<3, x>5]
+    >>> choicelist = [x, x**2]
     >>> np.select(condlist, choicelist)
-    array([  0,   0,   0,   0, 400,   0,   6,   7,   8,   9])
+    array([ 0,  1,  2,  0,  0,  0, 36, 49, 64, 81])
 
     """
     n = len(condlist)
@@ -960,11 +962,17 @@ def gradient(f, *varargs):
 
     Examples
     --------
-    >>> np.gradient(np.array([[1,1],[3,4]]))
-    [array([[ 2.,  3.],
-           [ 2.,  3.]]),
-     array([[ 0.,  0.],
-           [ 1.,  1.]])]
+    >>> x = np.array([1, 2, 4, 7, 11, 16], dtype=np.float)
+    >>> np.gradient(x)
+    array([ 1. ,  1.5,  2.5,  3.5,  4.5,  5. ])
+    >>> np.gradient(x, 2)
+    array([ 0.5 ,  0.75,  1.25,  1.75,  2.25,  2.5 ])
+
+    >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float))
+    [array([[ 2.,  2., -1.],
+           [ 2.,  2., -1.]]),
+    array([[ 1. ,  2.5,  4. ],
+           [ 1. ,  1. ,  1. ]])]
 
     """
     N = len(f.shape)  # number of dimensions
@@ -1026,7 +1034,11 @@ def gradient(f, *varargs):
 
 def diff(a, n=1, axis=-1):
     """
-    Calculate the nth order discrete difference along given axis.
+    Calculate the n-th order discrete difference along given axis.
+
+    The first order difference is given by ``out[n] = a[n+1] - a[n]`` along
+    the given axis, higher order differences are calculated by using `diff`
+    recursively.
 
     Parameters
     ----------
@@ -1035,26 +1047,31 @@ def diff(a, n=1, axis=-1):
     n : int, optional
         The number of times values are differenced.
     axis : int, optional
-        The axis along which the difference is taken.
+        The axis along which the difference is taken, default is the last axis.
 
     Returns
     -------
     out : ndarray
-        The `n` order differences.  The shape of the output is the same as `a`
-        except along `axis` where the dimension is `n` less.
+        The `n` order differences. The shape of the output is the same as `a`
+        except along `axis` where the dimension is smaller by `n`.
+
+    See Also
+    --------
+    gradient, ediff1d
 
     Examples
     --------
-    >>> x = np.array([0,1,3,9,5,10])
+    >>> x = np.array([1, 2, 4, 7, 0])
     >>> np.diff(x)
-    array([ 1,  2,  6, -4,  5])
-    >>> np.diff(x,n=2)
-    array([  1,   4, -10,   9])
-    >>> x = np.array([[1,3,6,10],[0,5,6,8]])
+    array([ 1,  2,  3, -7])
+    >>> np.diff(x, n=2)
+    array([  1,   1, -10])
+
+    >>> x = np.array([[1, 3, 6, 10], [0, 5, 6, 8]])
     >>> np.diff(x)
     array([[2, 3, 4],
-    [5, 1, 2]])
-    >>> np.diff(x,axis=0)
+           [5, 1, 2]])
+    >>> np.diff(x, axis=0)
     array([[-1,  2,  0, -2]])
 
     """
@@ -1201,15 +1218,34 @@ def unwrap(p, discont=pi, axis=-1):
     ----------
     p : array_like
         Input array.
-    discont : float
-        Maximum discontinuity between values.
-    axis : integer
-        Axis along which unwrap will operate.
+    discont : float, optional
+        Maximum discontinuity between values, default is ``pi``.
+    axis : int, optional
+        Axis along which unwrap will operate, default is the last axis.
 
     Returns
     -------
     out : ndarray
-        Output array
+        Output array.
+
+    See Also
+    --------
+    rad2deg, deg2rad
+
+    Notes
+    -----
+    If the discontinuity in `p` is smaller than ``pi``, but larger than
+    `discont`, no unwrapping is done because taking the 2*pi complement
+    would only make the discontinuity larger.
+
+    Examples
+    --------
+    >>> phase = np.linspace(0, np.pi, num=5)
+    >>> phase[3:] += np.pi
+    >>> phase
+    array([ 0.        ,  0.78539816,  1.57079633,  5.49778714,  6.28318531])
+    >>> np.unwrap(phase)
+    array([ 0.        ,  0.78539816,  1.57079633, -0.78539816,  0.        ])
 
     """
     p = asarray(p)
@@ -1310,31 +1346,11 @@ import sys
 if sys.hexversion < 0x2040000:
     from sets import Set as set
 
+@deprecate_with_doc('')
 def unique(x):
     """
-    Return the sorted, unique elements of an array or sequence.
-
-    Parameters
-    ----------
-    x : ndarray or sequence
-        Input array.
-
-    Returns
-    -------
-    y : ndarray
-        The sorted, unique elements are returned in a 1-D array.
-
-    Examples
-    --------
-    >>> np.unique([1, 1, 2, 2, 3, 3])
-    array([1, 2, 3])
-    >>> a = np.array([[1, 1], [2, 3]])
-    >>> np.unique(a)
-    array([1, 2, 3])
-
-    >>> np.unique([True, True, False])
-    array([False,  True], dtype=bool)
-
+    This function is deprecated.  Use numpy.lib.arraysetops.unique()
+    instead.
     """
     try:
         tmp = x.flatten()
@@ -1365,53 +1381,64 @@ def extract(condition, arr):
 
     See Also
     --------
-    take, put, putmask
+    take, put, putmask, compress
 
     Examples
     --------
-    >>> arr = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12]])
+    >>> arr = np.arange(12).reshape((3, 4))
     >>> arr
-    array([[ 1,  2,  3,  4],
-           [ 5,  6,  7,  8],
-           [ 9, 10, 11, 12]])
+    array([[ 0,  1,  2,  3],
+           [ 4,  5,  6,  7],
+           [ 8,  9, 10, 11]])
     >>> condition = np.mod(arr, 3)==0
     >>> condition
-    array([[False, False,  True, False],
-           [False,  True, False, False],
-           [ True, False, False,  True]], dtype=bool)
+    array([[ True, False, False,  True],
+           [False, False,  True, False],
+           [False,  True, False, False]], dtype=bool)
     >>> np.extract(condition, arr)
-    array([ 3,  6,  9, 12])
+    array([0, 3, 6, 9])
+
 
     If `condition` is boolean:
 
     >>> arr[condition]
-    array([ 3,  6,  9, 12])
+    array([0, 3, 6, 9])
 
     """
     return _nx.take(ravel(arr), nonzero(ravel(condition))[0])
 
 def place(arr, mask, vals):
     """
-    Changes elements of an array based on conditional and input values.
+    Change elements of an array based on conditional and input values.
 
-    Similar to ``putmask(a, mask, vals)`` but the 1D array `vals` has the
-    same number of elements as the non-zero values of `mask`. Inverse of
-    ``extract``.
+    Similar to ``np.putmask(a, mask, vals)``, the difference is that `place`
+    uses the first N elements of `vals`, where N is the number of True values
+    in `mask`, while `putmask` uses the elements where `mask` is True.
 
-    Sets `a`.flat[n] = `values`\\[n] for each n where `mask`.flat[n] is true.
+    Note that `extract` does the exact opposite of `place`.
 
     Parameters
     ----------
     a : array_like
         Array to put data into.
     mask : array_like
-        Boolean mask array.
-    values : array_like, shape(number of non-zero `mask`, )
-        Values to put into `a`.
+        Boolean mask array. Must have the same size as `a`.
+    vals : 1-D sequence
+        Values to put into `a`. Only the first N elements are used, where
+        N is the number of True values in `mask`. If `vals` is smaller
+        than N it will be repeated.
 
     See Also
     --------
-    putmask, put, take
+    putmask, put, take, extract
+
+    Examples
+    --------
+    >>> x = np.arange(6).reshape(2, 3)
+    >>> np.place(x, x>2, [44, 55])
+    >>> x
+    array([[ 0,  1,  2],
+           [44, 55, 44]])
 
     """
     return _insert(arr, mask, vals)
@@ -2841,6 +2868,25 @@ def trapz(y, x=None, dx=1.0, axis=-1):
     axis : int, optional
         Specify the axis.
 
+    Returns
+    -------
+    out : float
+        Definite integral as approximated by trapezoidal rule.
+
+    Notes
+    -----
+    Image [2]_ illustrates trapezoidal rule -- y-axis locations of points will
+    be taken from `y` array, by default x-axis distances between points will be
+    1.0, alternatively they can be provided with `x` array or with `dx` scalar.
+    Return value will be equal to combined area under the red lines.
+
+
+    References
+    ----------
+    .. [1] Wikipedia page: http://en.wikipedia.org/wiki/Trapezoidal_rule
+
+    .. [2] Illustration image: http://en.wikipedia.org/wiki/File:Composite_trapezoidal_rule_illustration.png
+
     Examples
     --------
     >>> np.trapz([1,2,3])
diff --git a/numpy/lib/index_tricks.py b/numpy/lib/index_tricks.py
index b8add9ed7..eeb1d37aa 100644
--- a/numpy/lib/index_tricks.py
+++ b/numpy/lib/index_tricks.py
@@ -3,29 +3,38 @@ __all__ = ['unravel_index',
            'ogrid',
            'r_', 'c_', 's_',
            'index_exp', 'ix_',
-           'ndenumerate','ndindex']
+           'ndenumerate','ndindex',
+           'fill_diagonal','diag_indices','diag_indices_from']
 
 import sys
 import numpy.core.numeric as _nx
-from numpy.core.numeric import asarray, ScalarType, array
+from numpy.core.numeric import ( asarray, ScalarType, array, alltrue, cumprod,
+                                 arange )
 from numpy.core.numerictypes import find_common_type
 import math
 
 import function_base
 import numpy.core.defmatrix as matrix
+from function_base import diff
 makemat = matrix.matrix
 
 # contributed by Stefan van der Walt
 def unravel_index(x,dims):
     """
-    Convert a flat index into an index tuple for an array of given shape.
+    Convert a flat index to an index tuple for an array of given shape.
 
     Parameters
     ----------
     x : int
         Flattened index.
-    dims : shape tuple
-        Input shape.
+    dims : tuple of ints
+        Input shape, the shape of an array into which indexing is
+        required.
+
+    Returns
+    -------
+    idx : tuple of ints
+        Tuple of the same shape as `dims`, containing the unraveled index.
 
     Notes
     -----
@@ -34,7 +43,7 @@ def unravel_index(x,dims):
 
     Examples
     --------
-    >>> arr = np.arange(20).reshape(5,4)
+    >>> arr = np.arange(20).reshape(5, 4)
     >>> arr
     array([[ 0,  1,  2,  3],
            [ 4,  5,  6,  7],
@@ -72,21 +81,45 @@ def unravel_index(x,dims):
     return tuple(x/dim_prod % dims)
 
 def ix_(*args):
-    """ Construct an open mesh from multiple sequences.
+    """
+    Construct an open mesh from multiple sequences.
+
+    This function takes N 1-D sequences and returns N outputs with N
+    dimensions each, such that the shape is 1 in all but one dimension
+    and the dimension with the non-unit shape value cycles through all
+    N dimensions.
 
-    This function takes n 1-d sequences and returns n outputs with n
-    dimensions each such that the shape is 1 in all but one dimension and
-    the dimension with the non-unit shape value cycles through all n
-    dimensions.
+    Using `ix_` one can quickly construct index arrays that will index
+    the cross product. ``a[np.ix_([1,3],[2,5])]`` returns the array
+    ``[a[1,2] a[1,5] a[3,2] a[3,5]]``.
 
-    Using ix_() one can quickly construct index arrays that will index
-    the cross product.
+    Parameters
+    ----------
+    args : 1-D sequences
 
-    a[ix_([1,3,7],[2,5,8])]  returns the array
+    Returns
+    -------
+    out : ndarrays
+        N arrays with N dimensions each, with N the number of input
+        sequences. Together these arrays form an open mesh.
+
+    See Also
+    --------
+    ogrid, mgrid, meshgrid
+
+    Examples
+    --------
+    >>> a = np.arange(10).reshape(2, 5)
+    >>> ixgrid = np.ix_([0,1], [2,4])
+    >>> ixgrid
+    (array([[0],
+           [1]]), array([[2, 4]]))
+    >>> print ixgrid[0].shape, ixgrid[1].shape
+    (2, 1) (1, 2)
+    >>> a[ixgrid]
+    array([[2, 4],
+           [7, 9]])
 
-    a[1,2]  a[1,5]  a[1,8]
-    a[3,2]  a[3,5]  a[3,8]
-    a[7,2]  a[7,5]  a[7,8]
     """
     out = []
     nd = len(args)
@@ -215,7 +248,11 @@ mgrid.__doc__ = None # set in numpy.add_newdocs
 ogrid.__doc__ = None # set in numpy.add_newdocs
 
 class AxisConcatenator(object):
-    """Translates slice objects to concatenation along an axis.
+    """
+    Translates slice objects to concatenation along an axis.
+
+    For detailed documentation on usage, see `r_`.
+
     """
     def _retval(self, res):
         if self.matrix:
@@ -338,11 +375,96 @@ class AxisConcatenator(object):
 # in help(r_)
 
 class RClass(AxisConcatenator):
-    """Translates slice objects to concatenation along the first axis.
+    """
+    Translates slice objects to concatenation along the first axis.
+
+    This is a simple way to build up arrays quickly. There are two use cases.
+
+    1. If the index expression contains comma separated arrays, then stack
+       them along their first axis.
+    2. If the index expression contains slice notation or scalars then create
+       a 1-D array with a range indicated by the slice notation.
+
+    If slice notation is used, the syntax ``start:stop:step`` is equivalent
+    to ``np.arange(start, stop, step)`` inside of the brackets. However, if
+    ``step`` is an imaginary number (i.e. 100j) then its integer portion is
+    interpreted as a number-of-points desired and the start and stop are
+    inclusive. In other words ``start:stop:stepj`` is interpreted as
+    ``np.linspace(start, stop, step, endpoint=1)`` inside of the brackets.
+    After expansion of slice notation, all comma separated sequences are
+    concatenated together.
+
+    Optional character strings placed as the first element of the index
+    expression can be used to change the output. The strings 'r' or 'c' result
+    in matrix output. If the result is 1-D and 'r' is specified a 1 x N (row)
+    matrix is produced. If the result is 1-D and 'c' is specified, then a N x 1
+    (column) matrix is produced. If the result is 2-D then both provide the
+    same matrix result.
+
+    A string integer specifies which axis to stack multiple comma separated
+    arrays along. A string of two comma-separated integers allows indication
+    of the minimum number of dimensions to force each entry into as the
+    second integer (the axis to concatenate along is still the first integer).
+
+    A string with three comma-separated integers allows specification of the
+    axis to concatenate along, the minimum number of dimensions to force the
+    entries to, and which axis should contain the start of the arrays which
+    are less than the specified number of dimensions. In other words the third
+    integer allows you to specify where the 1's should be placed in the shape
+    of the arrays that have their shapes upgraded. By default, they are placed
+    in the front of the shape tuple. The third argument allows you to specify
+    where the start of the array should be instead. Thus, a third argument of
+    '0' would place the 1's at the end of the array shape. Negative integers
+    specify where in the new shape tuple the last dimension of upgraded arrays
+    should be placed, so the default is '-1'.
+
+    Parameters
+    ----------
+    Not a function, so takes no parameters
 
-    For example:
+
+    Returns
+    -------
+    A concatenated ndarray or matrix.
+
+    See Also
+    --------
+    concatenate : Join a sequence of arrays together.
+    c_ : Translates slice objects to concatenation along the second axis.
+
+    Examples
+    --------
     >>> np.r_[np.array([1,2,3]), 0, 0, np.array([4,5,6])]
     array([1, 2, 3, 0, 0, 4, 5, 6])
+    >>> np.r_[-1:1:6j, [0]*3, 5, 6]
+    array([-1. , -0.6, -0.2,  0.2,  0.6,  1. ,  0. ,  0. ,  0. ,  5. ,  6. ])
+
+    String integers specify the axis to concatenate along or the minimum
+    number of dimensions to force entries into.
+
+    >>> np.r_['-1', a, a] # concatenate along last axis
+    array([[0, 1, 2, 0, 1, 2],
+           [3, 4, 5, 3, 4, 5]])
+    >>> np.r_['0,2', [1,2,3], [4,5,6]] # concatenate along first axis, dim>=2
+    array([[1, 2, 3],
+           [4, 5, 6]])
+
+    >>> np.r_['0,2,0', [1,2,3], [4,5,6]]
+    array([[1],
+           [2],
+           [3],
+           [4],
+           [5],
+           [6]])
+    >>> np.r_['1,2,0', [1,2,3], [4,5,6]]
+    array([[1, 4],
+           [2, 5],
+           [3, 6]])
+
+    Using 'r' or 'c' as a first string argument creates a matrix.
+
+    >>> np.r_['r',[1,2,3], [4,5,6]]
+    matrix([[1, 2, 3, 4, 5, 6]])
 
     """
     def __init__(self):
@@ -351,11 +473,21 @@ class RClass(AxisConcatenator):
 r_ = RClass()
 
 class CClass(AxisConcatenator):
-    """Translates slice objects to concatenation along the second axis.
+    """
+    Translates slice objects to concatenation along the second axis.
+
+    This is short-hand for ``np.r_['-1,2,0', index expression]``, which is
+    useful because of its common occurrence. In particular, arrays will be
+    stacked along their last axis after being upgraded to at least 2-D with
+    1's post-pended to the shape (column vectors made out of 1-D arrays).
 
-    For example:
+    For detailed documentation, see `r_`.
+
+    Examples
+    --------
     >>> np.c_[np.array([[1,2,3]]), 0, 0, np.array([[4,5,6]])]
-    array([1, 2, 3, 0, 0, 4, 5, 6])
+    array([[1, 2, 3, 0, 0, 4, 5, 6]])
+
     """
     def __init__(self):
         AxisConcatenator.__init__(self, -1, ndmin=2, trans1d=0)
@@ -373,9 +505,13 @@ class ndenumerate(object):
     a : ndarray
       Input array.
 
+    See Also
+    --------
+    ndindex, flatiter
+
     Examples
     --------
-    >>> a = np.array([[1,2],[3,4]])
+    >>> a = np.array([[1, 2], [3, 4]])
     >>> for index, x in np.ndenumerate(a):
     ...     print index, x
     (0, 0) 1
@@ -388,6 +524,17 @@ class ndenumerate(object):
         self.iter = asarray(arr).flat
 
     def next(self):
+        """
+        Standard iterator method, returns the index tuple and array value.
+
+        Returns
+        -------
+        coords : tuple of ints
+            The indices of the current iteration.
+        val : scalar
+            The array element of the current iteration.
+
+        """
         return self.iter.coords, self.iter.next()
 
     def __iter__(self):
@@ -399,17 +546,21 @@ class ndindex(object):
     An N-dimensional iterator object to index arrays.
 
     Given the shape of an array, an `ndindex` instance iterates over
-    the N-dimensional index of the array. At each iteration, the index of the
-    last dimension is incremented by one.
+    the N-dimensional index of the array. At each iteration a tuple
+    of indices is returned, the last dimension is iterated over first.
 
     Parameters
     ----------
-    `*args` : integers
-      The size of each dimension in the counter.
+    `*args` : ints
+      The size of each dimension of the array.
+
+    See Also
+    --------
+    ndenumerate, flatiter
 
     Examples
     --------
-    >>> for index in np.ndindex(3,2,1):
+    >>> for index in np.ndindex(3, 2, 1):
     ...     print index
     (0, 0, 0)
     (0, 1, 0)
@@ -442,9 +593,25 @@ class ndindex(object):
             self._incrementone(axis-1)
 
     def ndincr(self):
+        """
+        Increment the multi-dimensional index by one.
+
+        `ndincr` takes care of the "wrapping around" of the axes.
+        It is called by `ndindex.next` and not normally used directly.
+
+        """
         self._incrementone(self.nd-1)
 
     def next(self):
+        """
+        Standard iterator method, updates the index and returns the index tuple.
+
+        Returns
+        -------
+        val : tuple of ints
+            Returns a tuple containing the indices of the current iteration.
+
+        """
         if (self.index >= self.total):
             raise StopIteration
         val = tuple(self.ind)
@@ -501,3 +668,167 @@ index_exp = IndexExpression(maketuple=True)
 s_ = IndexExpression(maketuple=False)
 
 # End contribution from Konrad.
+
+
+# The following functions complement those in twodim_base, but are
+# applicable to N-dimensions.
+
+def fill_diagonal(a, val):
+    """Fill the main diagonal of the given array of any dimensionality.
+
+    For an array with ndim > 2, the diagonal is the list of locations with
+    indices a[i,i,...,i], all identical.
+
+    This function modifies the input array in-place, it does not return a
+    value.
+
+    This functionality can be obtained via diag_indices(), but internally this
+    version uses a much faster implementation that never constructs the indices
+    and uses simple slicing.
+
+    Parameters
+    ----------
+    a : array, at least 2-dimensional.
+      Array whose diagonal is to be filled, it gets modified in-place.
+
+    val : scalar
+      Value to be written on the diagonal, its type must be compatible with
+      that of the array a.
+
+    See also
+    --------
+    diag_indices, diag_indices_from
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    Examples
+    --------
+    >>> a = zeros((3,3),int)
+    >>> fill_diagonal(a,5)
+    >>> a
+    array([[5, 0, 0],
+           [0, 5, 0],
+           [0, 0, 5]])
+
+    The same function can operate on a 4-d array:
+    >>> a = zeros((3,3,3,3),int)
+    >>> fill_diagonal(a,4)
+
+    We only show a few blocks for clarity:
+    >>> a[0,0]
+    array([[4, 0, 0],
+           [0, 0, 0],
+           [0, 0, 0]])
+    >>> a[1,1]
+    array([[0, 0, 0],
+           [0, 4, 0],
+           [0, 0, 0]])
+    >>> a[2,2]
+    array([[0, 0, 0],
+           [0, 0, 0],
+           [0, 0, 4]])
+
+    """
+    if a.ndim < 2:
+        raise ValueError("array must be at least 2-d")
+    if a.ndim == 2:
+        # Explicit, fast formula for the common case.  For 2-d arrays, we
+        # accept rectangular ones.
+        step = a.shape[1] + 1
+    else:
+        # For more than d=2, the strided formula is only valid for arrays with
+        # all dimensions equal, so we check first.
+        if not alltrue(diff(a.shape)==0):
+            raise ValueError("All dimensions of input must be of equal length")
+        step = 1 + (cumprod(a.shape[:-1])).sum()
+
+    # Write the value out into the diagonal.
+    a.flat[::step] = val
+
+
+def diag_indices(n, ndim=2):
+    """Return the indices to access the main diagonal of an array.
+
+    This returns a tuple of indices that can be used to access the main
+    diagonal of an array with ndim (>=2) dimensions and shape (n,n,...,n).  For
+    ndim=2 this is the usual diagonal, for ndim>2 this is the set of indices
+    to access A[i,i,...,i] for i=[0..n-1].
+
+    Parameters
+    ----------
+    n : int
+      The size, along each dimension, of the arrays for which the returned
+      indices can be used.
+
+    ndim : int, optional
+      The number of dimensions.
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    See also
+    --------
+    diag_indices_from
+
+    Examples
+    --------
+    Create a set of indices to access the diagonal of a (4,4) array:
+    >>> di = diag_indices(4)
+
+    >>> a = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
+    >>> a
+    array([[ 1,  2,  3,  4],
+           [ 5,  6,  7,  8],
+           [ 9, 10, 11, 12],
+           [13, 14, 15, 16]])
+    >>> a[di] = 100
+    >>> a
+    array([[100,   2,   3,   4],
+           [  5, 100,   7,   8],
+           [  9,  10, 100,  12],
+           [ 13,  14,  15, 100]])
+
+    Now, we create indices to manipulate a 3-d array:
+    >>> d3 = diag_indices(2,3)
+
+    And use it to set the diagonal of a zeros array to 1:
+    >>> a = zeros((2,2,2),int)
+    >>> a[d3] = 1
+    >>> a
+    array([[[1, 0],
+            [0, 0]],
+
+           [[0, 0],
+            [0, 1]]])
+
+    """
+    idx = arange(n)
+    return (idx,) * ndim
+
+
+def diag_indices_from(arr):
+    """Return the indices to access the main diagonal of an n-dimensional array.
+
+    See diag_indices() for full details.
+
+    Parameters
+    ----------
+    arr : array, at least 2-d
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    """
+
+    if not arr.ndim >= 2:
+        raise ValueError("input array must be at least 2-d")
+    # For more than d=2, the strided formula is only valid for arrays with
+    # all dimensions equal, so we check first.
+    if not alltrue(diff(arr.shape) == 0):
+        raise ValueError("All dimensions of input must be of equal length")
+
+    return diag_indices(arr.shape[0], arr.ndim)
diff --git a/numpy/lib/info.py b/numpy/lib/info.py
index f93234d57..4a781a2ca 100644
--- a/numpy/lib/info.py
+++ b/numpy/lib/info.py
@@ -135,12 +135,11 @@ Set operations for 1D numeric arrays based on sort() function.
 
 ================ ===================
 ediff1d          Array difference (auxiliary function).
-unique1d         Unique elements of 1D array.
+unique           Unique elements of an array.
 intersect1d      Intersection of 1D arrays with unique elements.
-intersect1d_nu   Intersection of 1D arrays with any elements.
 setxor1d         Set exclusive-or of 1D arrays with unique elements.
-setmember1d      Return an array of shape of ar1 containing 1 where
-                 the elements of ar1 are in ar2 and 0 otherwise.
+in1d             Test whether elements in a 1D array are also present in
+                 another array.
 union1d          Union of 1D arrays with unique elements.
 setdiff1d        Set difference of 1D arrays with unique elements.
 ================ ===================
diff --git a/numpy/lib/io.py b/numpy/lib/io.py
index 98d071fab..3a962c7e1 100644
--- a/numpy/lib/io.py
+++ b/numpy/lib/io.py
@@ -118,6 +118,27 @@ class NpzFile(object):
         else:
             raise KeyError, "%s is not a file in the archive" % key
 
+
+    def __iter__(self):
+        return iter(self.files)
+
+    def items(self):
+        return [(f, self[f]) for f in self.files]
+
+    def iteritems(self):
+        for f in self.files:
+            yield (f, self[f])
+
+    def keys(self):
+        return self.files
+
+    def iterkeys(self):
+        return self.__iter__()
+
+    def __contains__(self, key):
+        return self.files.__contains__(key)
+
+
 def load(file, mmap_mode=None):
     """
     Load a pickled, ``.npy``, or ``.npz`` binary file.
@@ -126,6 +147,7 @@ def load(file, mmap_mode=None):
     ----------
     file : file-like object or string
         The file to read.  It must support ``seek()`` and ``read()`` methods.
+        If the filename extension is ``.gz``, the file is first decompressed.
     mmap_mode: {None, 'r+', 'r', 'w+', 'c'}, optional
         If not None, then memory-map the file, using the given mode
         (see `numpy.memmap`).  The mode has no effect for pickled or
@@ -146,6 +168,11 @@ def load(file, mmap_mode=None):
     IOError
         If the input file does not exist or cannot be read.
 
+    See Also
+    --------
+    save, savez, loadtxt
+    memmap : Create a memory-map to an array stored in a file on disk.
+
     Notes
     -----
     - If the file contains pickle data, then whatever is stored in the
@@ -202,20 +229,20 @@ def load(file, mmap_mode=None):
 
 def save(file, arr):
     """
-    Save an array to a binary file in NumPy format.
+    Save an array to a binary file in NumPy ``.npy`` format.
 
     Parameters
     ----------
-    f : file or string
+    file : file or string
         File or filename to which the data is saved.  If the filename
         does not already have a ``.npy`` extension, it is added.
-    x : array_like
-        Array data.
+    arr : array_like
+        Array data to be saved.
 
     See Also
     --------
-    savez : Save several arrays into an .npz compressed archive
-    savetxt : Save an array to a file as plain text
+    savez : Save several arrays into a .npz compressed archive
+    savetxt, load
 
     Examples
     --------
@@ -225,7 +252,7 @@ def save(file, arr):
     >>> x = np.arange(10)
     >>> np.save(outfile, x)
 
-    >>> outfile.seek(0)
+    >>> outfile.seek(0) # only necessary in this example (with tempfile)
     >>> np.load(outfile)
     array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
 
@@ -273,6 +300,12 @@ def savez(file, *args, **kwds):
     The .npz file format is a zipped archive of files named after the variables
     they contain.  Each file contains one variable in .npy format.
 
+    Examples
+    --------
+    >>> x = np.random.random((3, 3))
+    >>> y = np.zeros((3, 2))
+    >>> np.savez('data', x=x, y=y)
+
     """
 
     # Import is postponed to here since zipfile depends on gzip, an optional
@@ -523,20 +556,20 @@ def loadtxt(fname, dtype=float, comments='#', delimiter=None, converters=None,
 
 def savetxt(fname, X, fmt='%.18e',delimiter=' '):
     """
-    Save an array to file.
+    Save an array to a text file.
 
     Parameters
     ----------
-    fname : filename or a file handle
-        If the filename ends in .gz, the file is automatically saved in
-        compressed gzip format.  The load() command understands gzipped
-        files transparently.
+    fname : filename or file handle
+        If the filename ends in ``.gz``, the file is automatically saved in
+        compressed gzip format.  `loadtxt` understands gzipped files
+        transparently.
     X : array_like
-        Data.
-    fmt : string or sequence of strings
+        Data to be saved to a text file.
+    fmt : str or sequence of strs
         A single format (%10.5f), a sequence of formats, or a
         multi-format string, e.g. 'Iteration %d -- %10.5f', in which
-        case delimiter is ignored.
+        case `delimiter` is ignored.
     delimiter : str
         Character separating columns.
 
@@ -588,15 +621,20 @@ def savetxt(fname, X, fmt='%.18e',delimiter=' '):
 
         ``x,X`` : unsigned hexadecimal integer
 
-    This is not an exhaustive specification.
-
+    This explanation of ``fmt`` is not complete, for an exhaustive
+    specification see [1]_.
 
+    References
+    ----------
+    .. [1] `Format Specification Mini-Language
+           <http://docs.python.org/library/string.html#
+           format-specification-mini-language>`_, Python Documentation.
 
     Examples
     --------
-    >>> savetxt('test.out', x, delimiter=',') # X is an array
-    >>> savetxt('test.out', (x,y,z)) # x,y,z equal sized 1D arrays
-    >>> savetxt('test.out', x, fmt='%1.4e') # use exponential notation
+    >>> savetxt('test.out', x, delimiter=',')   # X is an array
+    >>> savetxt('test.out', (x,y,z))   # x,y,z equal sized 1D arrays
+    >>> savetxt('test.out', x, fmt='%1.4e')   # use exponential notation
 
     """
 
@@ -712,15 +750,13 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, skiprows=0,
 
     Each line past the first `skiprows` ones is split at the `delimiter`
     character, and characters following the `comments` character are discarded.
-    
-
 
     Parameters
     ----------
-    fname : file or string
-        File or filename to read.  If the filename extension is `.gz` or `.bz2`,
-        the file is first decompressed.
-    dtype : data-type
+    fname : {file, string}
+        File or filename to read.  If the filename extension is `.gz` or
+        `.bz2`, the file is first decompressed.
+    dtype : dtype
         Data type of the resulting array.  If this is a flexible data-type,
         the resulting array will be 1-dimensional, and each row will be
         interpreted as an element of the array. In this case, the number
@@ -729,20 +765,20 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, skiprows=0,
         of the dtype.
         If None, the dtypes will be determined by the contents of each
         column, individually.
-    comments : {string}, optional
+    comments : string, optional
         The character used to indicate the start of a comment.
         All the characters occurring on a line after a comment are discarded
-    delimiter : {string}, optional
+    delimiter : string, optional
         The string used to separate values.  By default, any consecutive
         whitespace act as delimiter.
-    skiprows : {int}, optional
+    skiprows : int, optional
         Numbers of lines to skip at the beginning of the file.
     converters : {None, dictionary}, optional
         A dictionary mapping column number to a function that will convert
         values in the column to a number. Converters can also be used to
         provide a default value for missing data:
         ``converters = {3: lambda s: float(s or 0)}``.
-    missing : {string}, optional
+    missing : string, optional
         A string representing a missing value, irrespective of the column where
         it appears (e.g., `'missing'` or `'unused'`).
     missing_values : {None, dictionary}, optional
@@ -757,20 +793,21 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, skiprows=0,
         If `names` is a sequence or a single-string of comma-separated names,
         the names will be used to define the field names in a flexible dtype.
         If `names` is None, the names of the dtype fields will be used, if any.
-    excludelist : {sequence}, optional
+    excludelist : sequence, optional
         A list of names to exclude. This list is appended to the default list
         ['return','file','print']. Excluded names are appended an underscore:
         for example, `file` would become `file_`.
-    deletechars : {string}, optional
-        A string combining invalid characters that must be deleted from the names.
+    deletechars : string, optional
+        A string combining invalid characters that must be deleted from the
+        names.
     case_sensitive : {True, False, 'upper', 'lower'}, optional
         If True, field names are case_sensitive.
         If False or 'upper', field names are converted to upper case.
         If 'lower', field names are converted to lower case.
-    unpack : {bool}, optional
+    unpack : bool, optional
         If True, the returned array is transposed, so that arguments may be
         unpacked using ``x, y, z = loadtxt(...)``
-    usemask : {bool}, optional
+    usemask : bool, optional
         If True, returns a masked array.
         If False, return a regular standard array.
 
@@ -779,23 +816,20 @@ def genfromtxt(fname, dtype=float, comments='#', delimiter=None, skiprows=0,
     out : MaskedArray
         Data read from the text file.
 
-    Notes
+    See Also
     --------
+    numpy.loadtxt : equivalent function when no data is missing.
+
+    Notes
+    -----
     * When spaces are used as delimiters, or when no delimiter has been given
       as input, there should not be any missing data between two fields.
     * When the variable are named (either by a flexible dtype or with `names`,
-      there must not be any header in the file (else a :exc:ValueError exception
-      is raised).
-
-    Warnings
-    --------
+      there must not be any header in the file (else a :exc:ValueError
+      exception is raised).
     * Individual values are not stripped of spaces by default.
       When using a custom converter, make sure the function does remove spaces.
 
-    See Also
-    --------
-    numpy.loadtxt : equivalent function when no data is missing.
-
     """
     #
     if usemask:
@@ -1128,20 +1162,21 @@ def recfromtxt(fname, dtype=None, comments='#', delimiter=None, skiprows=0,
                excludelist=None, deletechars=None, case_sensitive=True,
                usemask=False):
     """
-    Load ASCII data stored in fname and returns a standard recarray (if 
+    Load ASCII data stored in fname and returns a standard recarray (if
     `usemask=False`) or a MaskedRecords (if `usemask=True`).
-    
+
     Complete description of all the optional input parameters is available in
     the docstring of the `genfromtxt` function.
-    
+
     See Also
     --------
     numpy.genfromtxt : generic function
 
-    Warnings
-    --------
+    Notes
+    -----
     * by default, `dtype=None`, which means that the dtype of the output array
       will be determined from the data.
+
     """
     kwargs = dict(dtype=dtype, comments=comments, delimiter=delimiter, 
                   skiprows=skiprows, converters=converters,
diff --git a/numpy/lib/scimath.py b/numpy/lib/scimath.py
index 269d332bf..0e1bafa91 100644
--- a/numpy/lib/scimath.py
+++ b/numpy/lib/scimath.py
@@ -166,7 +166,8 @@ def _fix_real_abs_gt_1(x):
     return x
 
 def sqrt(x):
-    """Return the square root of x.
+    """
+    Return the square root of x.
 
     Parameters
     ----------
@@ -174,12 +175,29 @@ def sqrt(x):
 
     Returns
     -------
-    array_like output.
+    out : array_like
+
+    Notes
+    -----
+
+    As the numpy.sqrt, this returns the principal square root of x, which is
+    what most people mean when they use square root; the principal square root
+    of x is not any number z such as z^2 = x.
+
+    For positive numbers, the principal square root is defined as the positive
+    number z such as z^2 = x.
+
+    The principal square root of -1 is i, the principal square root of any
+    negative number -x is defined a i * sqrt(x). For any non zero complex
+    number, it is defined by using the following branch cut: x = r e^(i t) with
+    r > 0 and -pi < t <= pi. The principal square root is then
+    sqrt(r) e^(i t/2).
 
     Examples
     --------
 
     For real, non-negative inputs this works just like numpy.sqrt():
+
     >>> np.lib.scimath.sqrt(1)
     1.0
 
@@ -187,33 +205,20 @@ def sqrt(x):
     array([ 1.,  2.])
 
     But it automatically handles negative inputs:
+
     >>> np.lib.scimath.sqrt(-1)
     (0.0+1.0j)
 
     >>> np.lib.scimath.sqrt([-1,4])
     array([ 0.+1.j,  2.+0.j])
 
-    Notes
-    -----
-
-    As the numpy.sqrt, this returns the principal square root of x, which is
-    what most people mean when they use square root; the principal square root
-    of x is not any number z such as z^2 = x.
-
-    For positive numbers, the principal square root is defined as the positive
-    number z such as z^2 = x.
-
-    The principal square root of -1 is i, the principal square root of any
-    negative number -x is defined a i * sqrt(x). For any non zero complex
-    number, it is defined by using the following branch cut: x = r e^(i t) with
-    r > 0 and -pi < t <= pi. The principal square root is then
-    sqrt(r) e^(i t/2).
     """
     x = _fix_real_lt_zero(x)
     return nx.sqrt(x)
 
 def log(x):
-    """Return the natural logarithm of x.
+    """
+    Return the natural logarithm of x.
 
     If x contains negative inputs, the answer is computed and returned in the
     complex domain.
@@ -224,7 +229,7 @@ def log(x):
 
     Returns
     -------
-    array_like
+    out : array_like
 
     Examples
     --------
@@ -237,12 +242,14 @@ def log(x):
 
     >>> np.lib.scimath.log(-math.exp(1)) == (1+1j*math.pi)
     True
+
     """
     x = _fix_real_lt_zero(x)
     return nx.log(x)
 
 def log10(x):
-    """Return the base 10 logarithm of x.
+    """
+    Return the base 10 logarithm of x.
 
     If x contains negative inputs, the answer is computed and returned in the
     complex domain.
@@ -253,12 +260,13 @@ def log10(x):
 
     Returns
     -------
-    array_like
+    out : array_like
 
     Examples
     --------
 
     (We set the printing precision so the example can be auto-tested)
+
     >>> np.set_printoptions(precision=4)
 
     >>> np.lib.scimath.log10([10**1,10**2])
@@ -267,12 +275,14 @@ def log10(x):
 
     >>> np.lib.scimath.log10([-10**1,-10**2,10**2])
     array([ 1.+1.3644j,  2.+1.3644j,  2.+0.j    ])
+
     """
     x = _fix_real_lt_zero(x)
     return nx.log10(x)
 
 def logn(n, x):
-    """Take log base n of x.
+    """
+    Take log base n of x.
 
     If x contains negative inputs, the answer is computed and returned in the
     complex domain.
@@ -283,12 +293,13 @@ def logn(n, x):
 
     Returns
     -------
-    array_like
+    out : array_like
 
     Examples
     --------
 
     (We set the printing precision so the example can be auto-tested)
+
     >>> np.set_printoptions(precision=4)
 
     >>> np.lib.scimath.logn(2,[4,8])
@@ -296,13 +307,15 @@ def logn(n, x):
 
     >>> np.lib.scimath.logn(2,[-4,-8,8])
     array([ 2.+4.5324j,  3.+4.5324j,  3.+0.j    ])
+
     """
     x = _fix_real_lt_zero(x)
     n = _fix_real_lt_zero(n)
     return nx.log(x)/nx.log(n)
 
 def log2(x):
-    """ Take log base 2 of x.
+    """
+    Take log base 2 of x.
 
     If x contains negative inputs, the answer is computed and returned in the
     complex domain.
@@ -313,12 +326,13 @@ def log2(x):
 
     Returns
     -------
-    array_like
+    out : array_like
 
     Examples
     --------
 
     (We set the printing precision so the example can be auto-tested)
+
     >>> np.set_printoptions(precision=4)
 
     >>> np.lib.scimath.log2([4,8])
@@ -326,12 +340,14 @@ def log2(x):
 
     >>> np.lib.scimath.log2([-4,-8,8])
     array([ 2.+4.5324j,  3.+4.5324j,  3.+0.j    ])
+
     """
     x = _fix_real_lt_zero(x)
     return nx.log2(x)
 
 def power(x, p):
-    """Return x**p.
+    """
+    Return x**p.
 
     If x contains negative values, it is converted to the complex domain.
 
@@ -344,11 +360,12 @@ def power(x, p):
 
     Returns
     -------
-    array_like
+    out : array_like
 
     Examples
     --------
     (We set the printing precision so the example can be auto-tested)
+
     >>> np.set_printoptions(precision=4)
 
     >>> np.lib.scimath.power([2,4],2)
@@ -359,6 +376,7 @@ def power(x, p):
 
     >>> np.lib.scimath.power([-2,4],2)
     array([  4.+0.j,  16.+0.j])
+
     """
     x = _fix_real_lt_zero(x)
     p = _fix_int_lt_zero(p)
@@ -393,7 +411,8 @@ def arccos(x):
     return nx.arccos(x)
 
 def arcsin(x):
-    """Compute the inverse sine of x.
+    """
+    Compute the inverse sine of x.
 
     For real x with abs(x)<=1, this returns the principal value.
 
@@ -410,6 +429,7 @@ def arcsin(x):
     Examples
     --------
     (We set the printing precision so the example can be auto-tested)
+
     >>> np.set_printoptions(precision=4)
 
     >>> np.lib.scimath.arcsin(0)
@@ -417,12 +437,14 @@ def arcsin(x):
 
     >>> np.lib.scimath.arcsin([0,1])
     array([ 0.    ,  1.5708])
+
     """
     x = _fix_real_abs_gt_1(x)
     return nx.arcsin(x)
 
 def arctanh(x):
-    """Compute the inverse hyperbolic tangent of x.
+    """
+    Compute the inverse hyperbolic tangent of x.
 
     For real x with abs(x)<=1, this returns the principal value.
 
@@ -434,7 +456,7 @@ def arctanh(x):
 
     Returns
     -------
-    array_like
+    out : array_like
 
     Examples
     --------
@@ -446,6 +468,7 @@ def arctanh(x):
 
     >>> np.lib.scimath.arctanh([0,2])
     array([ 0.0000+0.j    ,  0.5493-1.5708j])
+
     """
     x = _fix_real_abs_gt_1(x)
     return nx.arctanh(x)
diff --git a/numpy/lib/shape_base.py b/numpy/lib/shape_base.py
index 19dd54f7a..a5bf4d0ea 100644
--- a/numpy/lib/shape_base.py
+++ b/numpy/lib/shape_base.py
@@ -892,6 +892,19 @@ def dsplit(ary,indices_or_sections):
         raise ValueError, 'vsplit only works on arrays of 3 or more dimensions'
     return split(ary,indices_or_sections,2)
 
+def get_array_prepare(*args):
+    """Find the wrapper for the array with the highest priority.
+
+    In case of ties, leftmost wins. If no wrapper is found, return None
+    """
+    wrappers = [(getattr(x, '__array_priority__', 0), -i,
+                 x.__array_prepare__) for i, x in enumerate(args)
+                                   if hasattr(x, '__array_prepare__')]
+    wrappers.sort()
+    if wrappers:
+        return wrappers[-1][-1]
+    return None
+
 def get_array_wrap(*args):
     """Find the wrapper for the array with the highest priority.
 
@@ -975,7 +988,6 @@ def kron(a,b):
     True
 
     """
-    wrapper = get_array_wrap(a, b)
     b = asanyarray(b)
     a = array(a,copy=False,subok=True,ndmin=b.ndim)
     ndb, nda = b.ndim, a.ndim
@@ -998,6 +1010,10 @@ def kron(a,b):
     axis = nd-1
     for _ in xrange(nd):
         result = concatenate(result, axis=axis)
+    wrapper = get_array_prepare(a, b)
+    if wrapper is not None:
+        result = wrapper(result)
+    wrapper = get_array_wrap(a, b)
     if wrapper is not None:
         result = wrapper(result)
     return result
@@ -1007,6 +1023,19 @@ def tile(A, reps):
     """
     Construct an array by repeating A the number of times given by reps.
 
+    If `reps` has length ``d``, the result will have dimension of
+    ``max(d, A.ndim)``.
+
+    If ``A.ndim < d``, `A` is promoted to be d-dimensional by prepending new
+    axes. So a shape (3,) array is promoted to (1, 3) for 2-D replication,
+    or shape (1, 1, 3) for 3-D replication. If this is not the desired
+    behavior, promote `A` to d-dimensions manually before calling this
+    function.
+
+    If ``A.ndim > d``, `reps` is promoted to `A`.ndim by pre-pending 1's to it.
+    Thus for an `A` of shape (2, 3, 4, 5), a `reps` of (2, 2) is treated as
+    (1, 1, 2, 2).
+
     Parameters
     ----------
     A : array_like
@@ -1017,24 +1046,11 @@ def tile(A, reps):
     Returns
     -------
     c : ndarray
-        The output array.
+        The tiled output array.
 
     See Also
     --------
-    repeat
-
-    Notes
-    -----
-    If `reps` has length d, the result will have dimension of max(d, `A`.ndim).
-
-    If `A`.ndim < d, `A` is promoted to be d-dimensional by prepending new
-    axes. So a shape (3,) array is promoted to (1,3) for 2-D replication,
-    or shape (1,1,3) for 3-D replication. If this is not the desired behavior,
-    promote `A` to d-dimensions manually before calling this function.
-
-    If `A`.ndim > d, `reps` is promoted to `A`.ndim by pre-pending 1's to it.
-    Thus for an `A` of shape (2,3,4,5), a `reps` of (2,2) is treated as
-    (1,1,2,2).
+    repeat : Repeat elements of an array.
 
     Examples
     --------
@@ -1046,7 +1062,6 @@ def tile(A, reps):
            [0, 1, 2, 0, 1, 2]])
     >>> np.tile(a, (2, 1, 2))
     array([[[0, 1, 2, 0, 1, 2]],
-    <BLANKLINE>
            [[0, 1, 2, 0, 1, 2]]])
 
     >>> b = np.array([[1, 2], [3, 4]])
diff --git a/numpy/lib/tests/test_arraysetops.py b/numpy/lib/tests/test_arraysetops.py
index 40bc11f6e..92305129a 100644
--- a/numpy/lib/tests/test_arraysetops.py
+++ b/numpy/lib/tests/test_arraysetops.py
@@ -9,39 +9,61 @@ from numpy.lib.arraysetops import *
 import warnings
 
 class TestAso(TestCase):
-    def test_unique1d( self ):
+    def test_unique( self ):
         a = np.array( [5, 7, 1, 2, 1, 5, 7] )
 
         ec = np.array( [1, 2, 5, 7] )
-        c = unique1d( a )
+        c = unique( a )
         assert_array_equal( c, ec )
 
         warnings.simplefilter('ignore', Warning)
-        unique, indices = unique1d( a, return_index=True )
+        vals, indices = unique( a, return_index=True )
         warnings.resetwarnings()
 
         ed = np.array( [2, 3, 0, 1] )
-        assert_array_equal(unique, ec)
+        assert_array_equal(vals, ec)
         assert_array_equal(indices, ed)
 
-        assert_array_equal([], unique1d([]))
+        warnings.simplefilter('ignore', Warning)
+        vals, ind0, ind1 = unique( a, return_index=True,
+                                     return_inverse=True )
+        warnings.resetwarnings()
+
+        ee = np.array( [2, 3, 0, 1, 0, 2, 3] )
+        assert_array_equal(vals, ec)
+        assert_array_equal(ind0, ed)
+        assert_array_equal(ind1, ee)
+
+        assert_array_equal([], unique([]))
 
     def test_intersect1d( self ):
+        # unique inputs
         a = np.array( [5, 7, 1, 2] )
         b = np.array( [2, 4, 3, 1, 5] )
 
         ec = np.array( [1, 2, 5] )
-        c = intersect1d( a, b )
+        c = intersect1d( a, b, assume_unique=True )
         assert_array_equal( c, ec )
 
+        # non-unique inputs
+        a = np.array( [5, 5, 7, 1, 2] )
+        b = np.array( [2, 1, 4, 3, 3, 1, 5] )
+
+        ed = np.array( [1, 2, 5] )
+        c = intersect1d( a, b )
+        assert_array_equal( c, ed )
+
         assert_array_equal([], intersect1d([],[]))
 
     def test_intersect1d_nu( self ):
+        # This should be removed when intersect1d_nu is removed.
         a = np.array( [5, 5, 7, 1, 2] )
         b = np.array( [2, 1, 4, 3, 3, 1, 5] )
 
         ec = np.array( [1, 2, 5] )
+        warnings.simplefilter('ignore', Warning)
         c = intersect1d_nu( a, b )
+        warnings.resetwarnings()
         assert_array_equal( c, ec )
 
         assert_array_equal([], intersect1d_nu([],[]))
@@ -83,11 +105,14 @@ class TestAso(TestCase):
         assert_array_equal([1],ediff1d(two_elem))
 
     def test_setmember1d( self ):
+        # This should be removed when setmember1d is removed.
         a = np.array( [5, 7, 1, 2] )
         b = np.array( [2, 4, 3, 1, 5] )
 
         ec = np.array( [True, False, True, True] )
+        warnings.simplefilter('ignore', Warning)
         c = setmember1d( a, b )
+        warnings.resetwarnings()
         assert_array_equal( c, ec )
 
         a[0] = 8
@@ -102,51 +127,77 @@ class TestAso(TestCase):
 
         assert_array_equal([], setmember1d([],[]))
 
-    def test_setmember1d_nu(self):
+    def test_in1d(self):
+        a = np.array( [5, 7, 1, 2] )
+        b = np.array( [2, 4, 3, 1, 5] )
+
+        ec = np.array( [True, False, True, True] )
+        c = in1d( a, b, assume_unique=True )
+        assert_array_equal( c, ec )
+
+        a[0] = 8
+        ec = np.array( [False, False, True, True] )
+        c = in1d( a, b, assume_unique=True )
+        assert_array_equal( c, ec )
+
+        a[0], a[3] = 4, 8
+        ec = np.array( [True, False, True, False] )
+        c = in1d( a, b, assume_unique=True )
+        assert_array_equal( c, ec )
+        
         a = np.array([5,4,5,3,4,4,3,4,3,5,2,1,5,5])
         b = [2,3,4]
 
         ec = [False, True, False, True, True, True, True, True, True, False,
               True, False, False, False]
-        c = setmember1d_nu(a, b)
+        c = in1d(a, b)
         assert_array_equal(c, ec)
 
         b = b + [5, 5, 4]
 
         ec = [True, True, True, True, True, True, True, True, True, True,
               True, False, True, True]
-        c = setmember1d_nu(a, b)
+        c = in1d(a, b)
         assert_array_equal(c, ec)
 
         a = np.array([5, 7, 1, 2])
         b = np.array([2, 4, 3, 1, 5])
 
         ec = np.array([True, False, True, True])
-        c = setmember1d_nu(a, b)
+        c = in1d(a, b)
         assert_array_equal(c, ec)
 
         a = np.array([5, 7, 1, 1, 2])
         b = np.array([2, 4, 3, 3, 1, 5])
 
         ec = np.array([True, False, True, True, True])
-        c = setmember1d_nu(a, b)
+        c = in1d(a, b)
         assert_array_equal(c, ec)
 
         a = np.array([5])
         b = np.array([2])
 
         ec = np.array([False])
-        c = setmember1d_nu(a, b)
+        c = in1d(a, b)
         assert_array_equal(c, ec)
 
         a = np.array([5, 5])
         b = np.array([2, 2])
 
         ec = np.array([False, False])
-        c = setmember1d_nu(a, b)
+        c = in1d(a, b)
         assert_array_equal(c, ec)
 
-        assert_array_equal(setmember1d_nu([], []), [])
+        assert_array_equal(in1d([], []), [])
+
+    def test_in1d_char_array( self ):
+        a = np.array(['a', 'b', 'c','d','e','c','e','b'])
+        b = np.array(['a','c'])
+
+        ec = np.array([True, False, True, False, False, True, False, False])
+        c = in1d(a, b)
+
+        assert_array_equal(c, ec)
 
     def test_union1d( self ):
         a = np.array( [5, 4, 7, 1, 2] )
@@ -159,7 +210,7 @@ class TestAso(TestCase):
         assert_array_equal([], union1d([],[]))
 
     def test_setdiff1d( self ):
-        a = np.array( [6, 5, 4, 7, 1, 2] )
+        a = np.array( [6, 5, 4, 7, 1, 2, 7, 4] )
         b = np.array( [2, 4, 3, 3, 2, 1, 5] )
 
         ec = np.array( [6, 7] )
@@ -180,14 +231,6 @@ class TestAso(TestCase):
         assert_array_equal(setdiff1d(a,b),np.array(['c']))
 
     def test_manyways( self ):
-        nItem = 100
-        a = np.fix( nItem / 10 * np.random.random( nItem ) )
-        b = np.fix( nItem / 10 * np.random.random( nItem ) )
-
-        c1 = intersect1d_nu( a, b )
-        c2 = unique1d( intersect1d( a, b ) )
-        assert_array_equal( c1, c2 )
-
         a = np.array( [5, 7, 1, 2, 8] )
         b = np.array( [9, 8, 2, 4, 3, 1, 5] )
 
diff --git a/numpy/lib/tests/test_financial.py b/numpy/lib/tests/test_financial.py
index 1ac14b561..c1d77c517 100644
--- a/numpy/lib/tests/test_financial.py
+++ b/numpy/lib/tests/test_financial.py
@@ -36,13 +36,18 @@ class TestFinancial(TestCase):
                             117.04, 2)
 
     def test_mirr(self):
-        v1 = [-4500,-800,800,800,600,600,800,800,700,3000]
-        assert_almost_equal(np.mirr(v1,0.08,0.055),
-                            0.0665, 4)
+        val = [-4500,-800,800,800,600,600,800,800,700,3000]
+        assert_almost_equal(np.mirr(val, 0.08, 0.055), 0.0666, 4)
+
+        val = [-120000,39000,30000,21000,37000,46000]
+        assert_almost_equal(np.mirr(val, 0.10, 0.12), 0.126094, 6)
+
+        val = [100,200,-50,300,-200]
+        assert_almost_equal(np.mirr(val, 0.05, 0.06), 0.3428, 4)
+
+        val = [39000,30000,21000,37000,46000]
+        assert_(np.isnan(np.mirr(val, 0.10, 0.12)))
 
-        v2 = [-120000,39000,30000,21000,37000,46000]
-        assert_almost_equal(np.mirr(v2,0.10,0.12),
-                            0.1344, 4)
 
 
 def test_unimplemented():
diff --git a/numpy/lib/tests/test_index_tricks.py b/numpy/lib/tests/test_index_tricks.py
index 47529502d..d7e61799a 100644
--- a/numpy/lib/tests/test_index_tricks.py
+++ b/numpy/lib/tests/test_index_tricks.py
@@ -1,5 +1,8 @@
 from numpy.testing import *
-from numpy import array, ones, r_, mgrid, unravel_index
+import numpy as np
+from numpy import ( array, ones, r_, mgrid, unravel_index, zeros, where,
+                    ndenumerate, fill_diagonal, diag_indices,
+                    diag_indices_from )
 
 class TestUnravelIndex(TestCase):
     def test_basic(self):
@@ -62,5 +65,60 @@ class TestConcatenator(TestCase):
         assert_array_equal(d[5:,:],c)
 
 
+class TestNdenumerate(TestCase):
+    def test_basic(self):
+        a = array([[1,2], [3,4]])
+        assert_equal(list(ndenumerate(a)),
+                     [((0,0), 1), ((0,1), 2), ((1,0), 3), ((1,1), 4)])
+
+
+def test_fill_diagonal():
+    a = zeros((3, 3),int)
+    fill_diagonal(a, 5)
+    yield (assert_array_equal, a,
+           array([[5, 0, 0],
+                  [0, 5, 0],
+                  [0, 0, 5]]))
+
+    # The same function can operate on a 4-d array:
+    a = zeros((3, 3, 3, 3), int)
+    fill_diagonal(a, 4)
+    i = array([0, 1, 2])
+    yield (assert_equal, where(a != 0), (i, i, i, i))
+
+
+def test_diag_indices():
+    di = diag_indices(4)
+    a = array([[1, 2, 3, 4],
+               [5, 6, 7, 8],
+               [9, 10, 11, 12],
+               [13, 14, 15, 16]])
+    a[di] = 100
+    yield (assert_array_equal, a,
+           array([[100,   2,   3,   4],
+                  [  5, 100,   7,   8],
+                  [  9,  10, 100,  12],
+                  [ 13,  14,  15, 100]]))
+
+    # Now, we create indices to manipulate a 3-d array:
+    d3 = diag_indices(2, 3)
+
+    # And use it to set the diagonal of a zeros array to 1:
+    a = zeros((2, 2, 2),int)
+    a[d3] = 1
+    yield (assert_array_equal, a,
+           array([[[1, 0],
+                   [0, 0]],
+
+                  [[0, 0],
+                   [0, 1]]]) )
+
+def test_diag_indices_from():
+    x = np.random.random((4, 4))
+    r, c = diag_indices_from(x)
+    assert_array_equal(r, np.arange(4))
+    assert_array_equal(c, np.arange(4))
+
+
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/lib/tests/test_io.py b/numpy/lib/tests/test_io.py
index e5a73a86a..185ceef7c 100644
--- a/numpy/lib/tests/test_io.py
+++ b/numpy/lib/tests/test_io.py
@@ -916,5 +916,31 @@ def test_gzip_loadtxt_from_string():
     f = gzip.GzipFile(fileobj=s, mode="r")
     assert_array_equal(np.loadtxt(f), [1, 2, 3])
 
+def test_npzfile_dict():
+    s = StringIO.StringIO()
+    x = np.zeros((3, 3))
+    y = np.zeros((3, 3))
+
+    np.savez(s, x=x, y=y)
+    s.seek(0)
+
+    z = np.load(s)
+
+    assert 'x' in z
+    assert 'y' in z
+    assert 'x' in z.keys()
+    assert 'y' in z.keys()
+
+    for f, a in z.iteritems():
+        assert f in ['x', 'y']
+        assert_equal(a.shape, (3, 3))
+
+    assert len(z.items()) == 2
+
+    for f in z:
+        assert f in ['x', 'y']
+
+    assert 'x' in list(z.iterkeys())
+
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/lib/tests/test_regression.py b/numpy/lib/tests/test_regression.py
index b8c487962..5abf9aefe 100644
--- a/numpy/lib/tests/test_regression.py
+++ b/numpy/lib/tests/test_regression.py
@@ -48,6 +48,10 @@ class TestRegression(object):
         """Ticket 928."""
         assert_raises(ValueError, np.histogramdd, np.ones((1,10)), bins=2**10)
 
+    def test_ndenumerate_crash(self):
+        """Ticket 1140"""
+        # Shouldn't crash:
+        list(np.ndenumerate(np.array([[]])))
 
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/lib/tests/test_twodim_base.py b/numpy/lib/tests/test_twodim_base.py
index 32c4ca58e..5d850f9fd 100644
--- a/numpy/lib/tests/test_twodim_base.py
+++ b/numpy/lib/tests/test_twodim_base.py
@@ -3,8 +3,11 @@
 """
 
 from numpy.testing import *
-from numpy import arange, rot90, add, fliplr, flipud, zeros, ones, eye, \
-     array, diag, histogram2d, tri
+
+from numpy import ( arange, rot90, add, fliplr, flipud, zeros, ones, eye,
+                    array, diag, histogram2d, tri, mask_indices, triu_indices,
+                    triu_indices_from, tril_indices, tril_indices_from )
+
 import numpy as np
 
 def get_mat(n):
@@ -50,34 +53,68 @@ class TestEye(TestCase):
                                           [1,0,0],
                                           [0,1,0]]))
 
+    def test_eye_bounds(self):
+        assert_equal(eye(2, 2,  1), [[0, 1], [0, 0]])
+        assert_equal(eye(2, 2, -1), [[0, 0], [1, 0]])
+        assert_equal(eye(2, 2,  2), [[0, 0], [0, 0]])
+        assert_equal(eye(2, 2, -2), [[0, 0], [0, 0]])
+        assert_equal(eye(3, 2,  2), [[0, 0], [0, 0], [0, 0]])
+        assert_equal(eye(3, 2,  1), [[0, 1], [0, 0], [0, 0]])
+        assert_equal(eye(3, 2, -1), [[0, 0], [1, 0], [0, 1]])
+        assert_equal(eye(3, 2, -2), [[0, 0], [0, 0], [1, 0]])
+        assert_equal(eye(3, 2, -3), [[0, 0], [0, 0], [0, 0]])
+
+    def test_strings(self):
+        assert_equal(eye(2, 2, dtype='S3'), [['1', ''], ['', '1']])
+
+    def test_bool(self):
+        assert_equal(eye(2, 2, dtype=bool), [[True, False], [False, True]])
+
 class TestDiag(TestCase):
     def test_vector(self):
-        vals = (100*arange(5)).astype('l')
-        b = zeros((5,5))
+        vals = (100 * arange(5)).astype('l')
+        b = zeros((5, 5))
         for k in range(5):
-            b[k,k] = vals[k]
-        assert_equal(diag(vals),b)
-        b = zeros((7,7))
+            b[k, k] = vals[k]
+        assert_equal(diag(vals), b)
+        b = zeros((7, 7))
         c = b.copy()
         for k in range(5):
-            b[k,k+2] = vals[k]
-            c[k+2,k] = vals[k]
-        assert_equal(diag(vals,k=2), b)
-        assert_equal(diag(vals,k=-2), c)
+            b[k, k + 2] = vals[k]
+            c[k + 2, k] = vals[k]
+        assert_equal(diag(vals, k=2), b)
+        assert_equal(diag(vals, k=-2), c)
 
-    def test_matrix(self):
-        vals = (100*get_mat(5)+1).astype('l')
+    def test_matrix(self, vals=None):
+        if vals is None:
+            vals = (100 * get_mat(5) + 1).astype('l')
         b = zeros((5,))
         for k in range(5):
             b[k] = vals[k,k]
-        assert_equal(diag(vals),b)
-        b = b*0
+        assert_equal(diag(vals), b)
+        b = b * 0
         for k in range(3):
-            b[k] = vals[k,k+2]
-        assert_equal(diag(vals,2),b[:3])
+            b[k] = vals[k, k + 2]
+        assert_equal(diag(vals, 2), b[:3])
         for k in range(3):
-            b[k] = vals[k+2,k]
-        assert_equal(diag(vals,-2),b[:3])
+            b[k] = vals[k + 2, k]
+        assert_equal(diag(vals, -2), b[:3])
+
+    def test_fortran_order(self):
+        vals = array((100 * get_mat(5) + 1), order='F', dtype='l')
+        self.test_matrix(vals)
+
+    def test_diag_bounds(self):
+        A = [[1, 2], [3, 4], [5, 6]]
+        assert_equal(diag(A, k=2), [])
+        assert_equal(diag(A, k=1), [2])
+        assert_equal(diag(A, k=0), [1, 4])
+        assert_equal(diag(A, k=-1), [3, 6])
+        assert_equal(diag(A, k=-2), [5])
+        assert_equal(diag(A, k=-3), [])
+
+    def test_failure(self):
+        self.failUnlessRaises(ValueError, diag, [[[1]]])
 
 class TestFliplr(TestCase):
     def test_basic(self):
@@ -193,5 +230,76 @@ class TestTri(TestCase):
         assert_array_equal(tri(3,dtype=bool),out.astype(bool))
 
 
+def test_mask_indices():
+    # simple test without offset
+    iu = mask_indices(3, np.triu)
+    a = np.arange(9).reshape(3, 3)
+    yield (assert_array_equal, a[iu], array([0, 1, 2, 4, 5, 8]))
+    # Now with an offset
+    iu1 = mask_indices(3, np.triu, 1)
+    yield (assert_array_equal, a[iu1], array([1, 2, 5]))
+
+
+def test_tril_indices():
+    # indices without and with offset
+    il1 = tril_indices(4)
+    il2 = tril_indices(4, 2)
+
+    a = np.array([[1, 2, 3, 4],
+                  [5, 6, 7, 8],
+                  [9, 10, 11, 12],
+                  [13, 14, 15, 16]])
+
+    # indexing:
+    yield (assert_array_equal, a[il1],
+           array([ 1,  5,  6,  9, 10, 11, 13, 14, 15, 16]) )
+
+    # And for assigning values:
+    a[il1] = -1
+    yield (assert_array_equal, a,
+    array([[-1,  2,  3,  4],
+           [-1, -1,  7,  8],
+           [-1, -1, -1, 12],
+           [-1, -1, -1, -1]]) )
+
+    # These cover almost the whole array (two diagonals right of the main one):
+    a[il2] = -10
+    yield (assert_array_equal, a,
+    array([[-10, -10, -10,   4],
+           [-10, -10, -10, -10],
+           [-10, -10, -10, -10],
+           [-10, -10, -10, -10]]) )
+
+
+def test_triu_indices():
+    iu1 = triu_indices(4)
+    iu2 = triu_indices(4, 2)
+
+    a = np.array([[1, 2, 3, 4],
+                  [5, 6, 7, 8],
+                  [9, 10, 11, 12],
+                  [13, 14, 15, 16]])
+
+    # Both for indexing:
+    yield (assert_array_equal, a[iu1],
+           array([1, 2,  3,  4,  6, 7, 8, 11, 12, 16]))
+
+    # And for assigning values:
+    a[iu1] = -1
+    yield (assert_array_equal, a,
+    array([[-1, -1, -1, -1],
+           [ 5, -1, -1, -1],
+           [ 9, 10, -1, -1],
+           [13, 14, 15, -1]])  )
+
+    # These cover almost the whole array (two diagonals right of the main one):
+    a[iu2] = -10
+    yield ( assert_array_equal, a,
+    array([[ -1,  -1, -10, -10],
+           [  5,  -1,  -1, -10],
+           [  9,  10,  -1,  -1],
+           [ 13,  14,  15,  -1]]) )
+
+
 if __name__ == "__main__":
     run_module_suite()
diff --git a/numpy/lib/twodim_base.py b/numpy/lib/twodim_base.py
index f0abf3122..e794d4144 100644
--- a/numpy/lib/twodim_base.py
+++ b/numpy/lib/twodim_base.py
@@ -3,10 +3,13 @@
 """
 
 __all__ = ['diag','diagflat','eye','fliplr','flipud','rot90','tri','triu',
-           'tril','vander','histogram2d']
+           'tril','vander','histogram2d','mask_indices',
+           'tril_indices','tril_indices_from','triu_indices','triu_indices_from',
+           ]
 
 from numpy.core.numeric import asanyarray, equal, subtract, arange, \
-     zeros, greater_equal, multiply, ones, asarray
+     zeros, greater_equal, multiply, ones, asarray, alltrue, where, \
+     empty
 
 def fliplr(m):
     """
@@ -195,10 +198,16 @@ def eye(N, M=None, k=0, dtype=float):
            [ 0.,  0.,  0.]])
 
     """
-    if M is None: M = N
-    m = equal(subtract.outer(arange(N), arange(M)),-k)
-    if m.dtype != dtype:
-        m = m.astype(dtype)
+    if M is None:
+        M = N
+    m = zeros((N, M), dtype=dtype)
+    if k >= M:
+        return m
+    if k >= 0:
+        i = k
+    else:
+        i = (-k) * M
+    m[:M-k].flat[i::M+1] = 1
     return m
 
 def diag(v, k=0):
@@ -244,28 +253,26 @@ def diag(v, k=0):
     """
     v = asarray(v)
     s = v.shape
-    if len(s)==1:
+    if len(s) == 1:
         n = s[0]+abs(k)
         res = zeros((n,n), v.dtype)
-        if (k>=0):
-            i = arange(0,n-k)
-            fi = i+k+i*n
+        if k >= 0:
+            i = k
         else:
-            i = arange(0,n+k)
-            fi = i+(i-k)*n
-        res.flat[fi] = v
+            i = (-k) * n
+        res[:n-k].flat[i::n+1] = v
         return res
-    elif len(s)==2:
-        N1,N2 = s
+    elif len(s) == 2:
+        if k >= s[1]:
+            return empty(0, dtype=v.dtype)
+        if v.flags.f_contiguous:
+            # faster slicing
+            v, k, s = v.T, -k, s[::-1]
         if k >= 0:
-            M = min(N1,N2-k)
-            i = arange(0,M)
-            fi = i+k+i*N2
+            i = k
         else:
-            M = min(N1+k,N2)
-            i = arange(0,M)
-            fi = i + (i-k)*N2
-        return v.flat[fi]
+            i = (-k) * s[1]
+        return v[:s[1]-k].flat[i::s[1]+1]
     else:
         raise ValueError, "Input must be 1- or 2-d."
 
@@ -559,3 +566,233 @@ def histogram2d(x,y, bins=10, range=None, normed=False, weights=None):
         bins = [xedges, yedges]
     hist, edges = histogramdd([x,y], bins, range, normed, weights)
     return hist, edges[0], edges[1]
+
+    
+def mask_indices(n,mask_func,k=0):
+    """Return the indices to access (n,n) arrays, given a masking function.
+
+    Assume mask_func() is a function that, for a square array a of size (n,n)
+    with a possible offset argument k, when called as mask_func(a,k) returns a
+    new array with zeros in certain locations (functions like triu() or tril()
+    do precisely this).  Then this function returns the indices where the
+    non-zero values would be located.
+
+    Parameters
+    ----------
+    n : int
+      The returned indices will be valid to access arrays of shape (n,n).
+
+    mask_func : callable
+      A function whose api is similar to that of numpy.tri{u,l}.  That is,
+      mask_func(x,k) returns a boolean array, shaped like x.  k is an optional
+      argument to the function.
+
+    k : scalar
+      An optional argument which is passed through to mask_func().  Functions
+      like tri{u,l} take a second argument that is interpreted as an offset.
+
+    Returns
+    -------
+    indices : an n-tuple of index arrays.
+      The indices corresponding to the locations where mask_func(ones((n,n)),k)
+      is True.
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    Examples
+    --------
+    These are the indices that would allow you to access the upper triangular
+    part of any 3x3 array:
+    >>> iu = mask_indices(3,np.triu)
+
+    For example, if `a` is a 3x3 array:
+    >>> a = np.arange(9).reshape(3,3)
+    >>> a
+    array([[0, 1, 2],
+           [3, 4, 5],
+           [6, 7, 8]])
+
+    Then:
+    >>> a[iu]
+    array([0, 1, 2, 4, 5, 8])
+
+    An offset can be passed also to the masking function.  This gets us the
+    indices starting on the first diagonal right of the main one:
+    >>> iu1 = mask_indices(3,np.triu,1)
+
+    with which we now extract only three elements:
+    >>> a[iu1]
+    array([1, 2, 5])
+    """ 
+    m = ones((n,n),int)
+    a = mask_func(m,k)
+    return where(a != 0)
+
+
+def tril_indices(n,k=0):
+    """Return the indices for the lower-triangle of an (n,n) array.
+
+    Parameters
+    ----------
+    n : int
+      Sets the size of the arrays for which the returned indices will be valid.
+
+    k : int, optional
+      Diagonal offset (see tril() for details).
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    Examples
+    --------
+    Commpute two different sets of indices to access 4x4 arrays, one for the
+    lower triangular part starting at the main diagonal, and one starting two
+    diagonals further right:
+    
+    >>> il1 = tril_indices(4)
+    >>> il2 = tril_indices(4,2)
+
+    Here is how they can be used with a sample array:
+    >>> a = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
+    >>> a
+    array([[ 1,  2,  3,  4],
+           [ 5,  6,  7,  8],
+           [ 9, 10, 11, 12],
+           [13, 14, 15, 16]])
+
+    Both for indexing:
+    >>> a[il1]
+    array([ 1,  5,  6,  9, 10, 11, 13, 14, 15, 16])
+
+    And for assigning values:
+    >>> a[il1] = -1
+    >>> a
+    array([[-1,  2,  3,  4],
+           [-1, -1,  7,  8],
+           [-1, -1, -1, 12],
+           [-1, -1, -1, -1]])
+
+    These cover almost the whole array (two diagonals right of the main one):
+    >>> a[il2] = -10 
+    >>> a
+    array([[-10, -10, -10,   4],
+           [-10, -10, -10, -10],
+           [-10, -10, -10, -10],
+           [-10, -10, -10, -10]])
+
+    See also
+    --------
+    - triu_indices : similar function, for upper-triangular.
+    - mask_indices : generic function accepting an arbitrary mask function.
+    """
+    return mask_indices(n,tril,k)
+
+
+def tril_indices_from(arr,k=0):
+    """Return the indices for the lower-triangle of an (n,n) array.
+
+    See tril_indices() for full details.
+    
+    Parameters
+    ----------
+    n : int
+      Sets the size of the arrays for which the returned indices will be valid.
+
+    k : int, optional
+      Diagonal offset (see tril() for details).
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    """
+    if not arr.ndim==2 and arr.shape[0] == arr.shape[1]:
+        raise ValueError("input array must be 2-d and square")
+    return tril_indices(arr.shape[0],k)
+
+    
+def triu_indices(n,k=0):
+    """Return the indices for the upper-triangle of an (n,n) array.
+
+    Parameters
+    ----------
+    n : int
+      Sets the size of the arrays for which the returned indices will be valid.
+
+    k : int, optional
+      Diagonal offset (see triu() for details).
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    Examples
+    --------
+    Commpute two different sets of indices to access 4x4 arrays, one for the
+    lower triangular part starting at the main diagonal, and one starting two
+    diagonals further right:
+
+    >>> iu1 = triu_indices(4)
+    >>> iu2 = triu_indices(4,2)
+
+    Here is how they can be used with a sample array:
+    >>> a = np.array([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16]])
+    >>> a
+    array([[ 1,  2,  3,  4],
+           [ 5,  6,  7,  8],
+           [ 9, 10, 11, 12],
+           [13, 14, 15, 16]])
+
+    Both for indexing:
+    >>> a[il1]
+    array([ 1,  5,  6,  9, 10, 11, 13, 14, 15, 16])
+
+    And for assigning values:       
+    >>> a[iu] = -1
+    >>> a
+    array([[-1, -1, -1, -1],
+           [ 5, -1, -1, -1],
+           [ 9, 10, -1, -1],
+           [13, 14, 15, -1]])
+
+    These cover almost the whole array (two diagonals right of the main one):
+    >>> a[iu2] = -10
+    >>> a
+    array([[ -1,  -1, -10, -10],
+           [  5,  -1,  -1, -10],
+           [  9,  10,  -1,  -1],
+           [ 13,  14,  15,  -1]])
+
+    See also
+    --------
+    - tril_indices : similar function, for lower-triangular.
+    - mask_indices : generic function accepting an arbitrary mask function.
+    """
+    return mask_indices(n,triu,k)
+
+
+def triu_indices_from(arr,k=0):
+    """Return the indices for the lower-triangle of an (n,n) array.
+
+    See triu_indices() for full details.
+    
+    Parameters
+    ----------
+    n : int
+      Sets the size of the arrays for which the returned indices will be valid.
+
+    k : int, optional
+      Diagonal offset (see triu() for details).
+
+    Notes
+    -----
+    .. versionadded:: 1.4.0
+
+    """
+    if not arr.ndim==2 and arr.shape[0] == arr.shape[1]:
+        raise ValueError("input array must be 2-d and square")
+    return triu_indices(arr.shape[0],k)
+
diff --git a/numpy/lib/type_check.py b/numpy/lib/type_check.py
index 113cec682..69f4f2193 100644
--- a/numpy/lib/type_check.py
+++ b/numpy/lib/type_check.py
@@ -85,8 +85,8 @@ def real(val):
     Returns
     -------
     out : ndarray
-        If `val` is real, the type of `val` is used for the output.  If `val`
-        has complex elements, the returned type is float.
+        Output array. If `val` is real, the type of `val` is used for the
+        output.  If `val` has complex elements, the returned type is float.
 
     See Also
     --------
@@ -94,13 +94,13 @@ def real(val):
 
     Examples
     --------
-    >>> a = np.array([1+2j,3+4j,5+6j])
+    >>> a = np.array([1+2j, 3+4j, 5+6j])
     >>> a.real
     array([ 1.,  3.,  5.])
     >>> a.real = 9
     >>> a
     array([ 9.+2.j,  9.+4.j,  9.+6.j])
-    >>> a.real = np.array([9,8,7])
+    >>> a.real = np.array([9, 8, 7])
     >>> a
     array([ 9.+2.j,  8.+4.j,  7.+6.j])
 
@@ -109,7 +109,7 @@ def real(val):
 
 def imag(val):
     """
-    Return the imaginary part of array.
+    Return the imaginary part of the elements of the array.
 
     Parameters
     ----------
@@ -118,8 +118,22 @@ def imag(val):
 
     Returns
     -------
-    out : ndarray, real or int
-        Real part of each element, same shape as `val`.
+    out : ndarray
+        Output array. If `val` is real, the type of `val` is used for the
+        output.  If `val` has complex elements, the returned type is float.
+
+    See Also
+    --------
+    real, angle, real_if_close
+
+    Examples
+    --------
+    >>> a = np.array([1+2j, 3+4j, 5+6j])
+    >>> a.imag
+    array([ 2.,  4.,  6.])
+    >>> a.imag = np.array([8, 10, 12])
+    >>> a
+    array([ 1. +8.j,  3.+10.j,  5.+12.j])
 
     """
     return asanyarray(val).imag
diff --git a/numpy/lib/ufunclike.py b/numpy/lib/ufunclike.py
index 5dbc3f225..5e89b0930 100644
--- a/numpy/lib/ufunclike.py
+++ b/numpy/lib/ufunclike.py
@@ -176,7 +176,7 @@ def isneginf(x, y=None):
 _log2 = nx.log(2)
 def log2(x, y=None):
     """
-    Return the base 2 logarithm.
+    Return the base 2 logarithm of the input array, element-wise.
 
     Parameters
     ----------
@@ -188,7 +188,7 @@ def log2(x, y=None):
     Returns
     -------
     y : ndarray
-      The logarithm to the base 2 of `x` elementwise.
+      The logarithm to the base 2 of `x` element-wise.
       NaNs are returned where `x` is negative.
 
     See Also
@@ -197,7 +197,7 @@ def log2(x, y=None):
 
     Examples
     --------
-    >>> np.log2([-1,2,4])
+    >>> np.log2([-1, 2, 4])
     array([ NaN,   1.,   2.])
 
     """
diff --git a/numpy/lib/utils.py b/numpy/lib/utils.py
index 3de0579df..908c4995d 100644
--- a/numpy/lib/utils.py
+++ b/numpy/lib/utils.py
@@ -81,12 +81,34 @@ else:
         return func
 
 def deprecate(func, oldname=None, newname=None):
-    """Deprecate old functions.
+    """
+    Deprecate old functions.
+
     Issues a DeprecationWarning, adds warning to oldname's docstring,
     rebinds oldname.__name__ and returns new function object.
 
-    Example:
-    oldfunc = deprecate(newfunc, 'oldfunc', 'newfunc')
+    Parameters
+    ----------
+    func : function
+
+    oldname : string
+
+    newname : string
+
+    Returns
+    -------
+    old_func : function
+
+    Examples
+    --------
+    Note that olduint returns a value after printing Deprecation Warning.
+
+    >>> olduint = np.deprecate(np.uint)
+    >>> olduint(6)
+    /usr/lib/python2.5/site-packages/numpy/lib/utils.py:114:
+    DeprecationWarning: uint32 is deprecated
+      warnings.warn(str1, DeprecationWarning)
+    6
 
     """
 
@@ -186,13 +208,28 @@ def byte_bounds(a):
 
 
 def may_share_memory(a, b):
-    """Determine if two arrays can share memory
+    """
+    Determine if two arrays can share memory
 
     The memory-bounds of a and b are computed.  If they overlap then
     this function returns True.  Otherwise, it returns False.
 
     A return of True does not necessarily mean that the two arrays
     share any element.  It just means that they *might*.
+
+    Parameters
+    ----------
+    a, b : ndarray
+
+    Returns
+    -------
+    out : bool
+
+    Examples
+    --------
+    >>> np.may_share_memory(np.array([1,2]), np.array([5,8,9]))
+    False
+
     """
     a_low, a_high = byte_bounds(a)
     b_low, b_high = byte_bounds(b)
@@ -349,24 +386,46 @@ def info(object=None,maxwidth=76,output=sys.stdout,toplevel='numpy'):
 
     Parameters
     ----------
-    object : optional
-        Input object to get information about.
+    object : object or str, optional
+        Input object or name to get information about. If `object` is a
+        numpy object, its docstring is given. If it is a string, available
+        modules are searched for matching objects.
+        If None, information about `info` itself is returned.
     maxwidth : int, optional
         Printing width.
-    output : file like object open for writing, optional
-        Write into file like object.
-    toplevel : string, optional
+    output : file like object, optional
+        File like object that the output is written to, default is ``stdout``.
+        The object has to be opened in 'w' or 'a' mode.
+    toplevel : str, optional
         Start search at this level.
 
+    See Also
+    --------
+    source, lookfor
+
+    Notes
+    -----
+    When used interactively with an object, ``np.info(obj)`` is equivalent to
+    ``help(obj)`` on the Python prompt or ``obj?`` on the IPython prompt.
+
     Examples
     --------
     >>> np.info(np.polyval) # doctest: +SKIP
-
        polyval(p, x)
+         Evaluate the polynomial p at x.
+         ...
 
-         Evaluate the polymnomial p at x.
+    When using a string for `object` it is possible to get multiple results.
 
-         ...
+    >>> np.info('fft') # doctest: +SKIP
+         *** Found in numpy ***
+    Core FFT routines
+    ...
+         *** Found in numpy.fft ***
+     fft(a, n=None, axis=-1)
+    ...
+         *** Repeat reference found in numpy.fft.fftpack ***
+         *** Total of 3 references found. ***
 
     """
     global _namedict, _dictlist
@@ -512,15 +571,39 @@ def source(object, output=sys.stdout):
     """
     Print or write to a file the source code for a Numpy object.
 
+    The source code is only returned for objects written in Python. Many
+    functions and classes are defined in C and will therefore not return
+    useful information.
+
     Parameters
     ----------
     object : numpy object
-        Input object.
+        Input object. This can be any object (function, class, module, ...).
     output : file object, optional
         If `output` not supplied then source code is printed to screen
         (sys.stdout).  File object must be created with either write 'w' or
         append 'a' modes.
 
+    See Also
+    --------
+    lookfor, info
+
+    Examples
+    --------
+    >>> np.source(np.interp)
+    In file: /usr/lib/python2.6/dist-packages/numpy/lib/function_base.py
+    def interp(x, xp, fp, left=None, right=None):
+        \"\"\".... (full docstring printed)\"\"\"
+        if isinstance(x, (float, int, number)):
+            return compiled_interp([x], xp, fp, left, right).item()
+        else:
+            return compiled_interp(x, xp, fp, left, right)
+
+    The source code is only returned for objects written in Python.
+
+    >>> np.source(np.array)
+    Not available for this object.
+
     """
     # Local import to speed up numpy's import time.
     import inspect
@@ -544,28 +627,41 @@ def lookfor(what, module=None, import_modules=True, regenerate=False):
     Do a keyword search on docstrings.
 
     A list of of objects that matched the search is displayed,
-    sorted by relevance.
+    sorted by relevance. All given keywords need to be found in the
+    docstring for it to be returned as a result, but the order does
+    not matter.
 
     Parameters
     ----------
     what : str
         String containing words to look for.
-    module : str, module
-        Module whose docstrings to go through.
-    import_modules : bool
+    module : str, optional
+        Name of module whose docstrings to go through.
+    import_modules : bool, optional
         Whether to import sub-modules in packages.
-        Will import only modules in ``__all__``.
-    regenerate : bool
-        Whether to re-generate the docstring cache.
+        Will import only modules in ``__all__``. Default is True.
+    regenerate : bool, optional
+        Whether to re-generate the docstring cache. Default is False.
 
-    Examples
+    See Also
     --------
+    source, info
+
+    Notes
+    -----
+    Relevance is determined only roughly, by checking if the keywords occur
+    in the function name, at the start of a docstring, etc.
 
+    Examples
+    --------
     >>> np.lookfor('binary representation')
     Search results for 'binary representation'
     ------------------------------------------
     numpy.binary_repr
         Return the binary representation of the input number as a string.
+    numpy.base_repr
+        Return a string representation of a number in the given base system.
+    ...
 
     """
     import pydoc