def svd(g, svdcut=1e-12, wgts=False, add_svdnoise=False):
    """ Apply SVD cuts to collection of |GVar|\s in ``g``.

    Standard usage is, for example, ::

        svdcut = ...
        gmod = svd(g, svdcut=svdcut)

    where ``g`` is an array of |GVar|\s or a dictionary containing |GVar|\s
    and/or arrays of |GVar|\s. When ``svdcut>0``, ``gmod`` is
    a copy of ``g`` whose |GVar|\s have been modified to make
    their correlation matrix less singular than that of the
    original ``g``: each eigenvalue ``eig`` of the correlation matrix is
    replaced by ``max(eig, svdcut * max_eig)`` where ``max_eig`` is
    the largest eigenvalue. This SVD cut, which is applied separately
    to each block-diagonal sub-matrix of the correlation matrix,
    increases the variance of the eigenmodes with eigenvalues smaller
    than ``svdcut * max_eig``.

    The modification of ``g``'s covariance matrix is implemented by adding
    (to ``g``) a set of |GVar|\s with zero means::

        gmod = g + gmod.svdcorrection

    where ``gmod.svdcorrection`` is an array/dictionary
    containing the |GVar|\s. If
    parameter ``add_svdnoise=True``,
    noise is included in ``gmod.svdcorrection``,
    ::

        gmod.svdcorrection += gv.sample(gmod.svdcorrection),

    before it is added to ``g``. The noise can be useful for testing fits
    and other applications.

    When ``svdcut`` is negative, eigenmodes of the correlation matrix
    whose eigenvalues are smaller than ``|svdcut| * max_eig`` are dropped
    from the new matrix and the corresponding components of ``g`` are
    zeroed out (that is, replaced by 0(0)) in ``gmod``.

    There is an additional parameter ``wgts`` in :func:`gvar.svd` whose
    default value is ``False``. Setting ``wgts=1`` or ``wgts=-1`` instead
    causes :func:`gvar.svd` to return a tuple ``(gmod, i_wgts)`` where
    ``gmod``  is the modified copy of ``g``, and ``i_wgts`` contains a
    spectral  decomposition of the covariance matrix corresponding to
    the  modified correlation matrix if ``wgts=1``, or a decomposition of its
    inverse if ``wgts=-1``. The first entry ``i, wgts = i_wgts[0]``  specifies
    the diagonal part of the matrix: ``i`` is a list of the indices in
    ``gmod.flat`` corresponding to diagonal elements, and ``wgts ** 2``
    gives the corresponding matrix elements. The second and subsequent
    entries, ``i, wgts = i_wgts[n]`` for ``n > 0``, each correspond
    to block-diagonal sub-matrices, where ``i`` is the list of
    indices corresponding to the block, and ``wgts[j]`` are eigenvectors of
    the sub-matrix rescaled so that ::

        numpy.sum(numpy.outer(wi, wi) for wi in wgts[j]

    is the sub-matrix (``wgts=1``) or its inverse (``wgts=-1``).

    To compute the inverse of the covariance matrix from ``i_wgts``,
    for example, one could use code like::

        gmod, i_wgts = svd(g, svdcut=svdcut, wgts=-1)

        inv_cov = numpy.zeros((n, n), float)
        i, wgts = i_wgts[0]                       # 1x1 sub-matrices
        if len(i) > 0:
            inv_cov[i, i] = numpy.array(wgts) ** 2
        for i, wgts in i_wgts[1:]:                # nxn sub-matrices (n>1)
            for w in wgts:
                inv_cov[i[:, None], i] += numpy.outer(w, w)

    This sets ``inv_cov`` equal to the inverse of the covariance matrix of
    the ``gmod``\s. Similarly, we can  compute the expectation value,
    ``u.dot(inv_cov.dot(v))``, between two vectors (:mod:`numpy` arrays)
    using::

        result = 0.0
        i, wgts = i_wgts[0]                       # 1x1 sub-matrices
        if len(i) > 0:
            result += numpy.sum((u[i] * wgts) * (v[i] * wgts))
        for i, wgts in i_wgts[1:]:                # nxn sub-matrices (n>1)
            result += numpy.sum(wgts.dot(u[i]) * wgts.dot(v[i]))

    where ``result`` is the desired expectation value.

    Args:

        g: An array of |GVar|\s or a dicitionary whose values are
            |GVar|\s and/or arrays of |GVar|\s.
        svdcut (None or float): If positive, replace eigenvalues ``eig``
            of the correlation matrix with ``max(eig, svdcut * max_eig)``
            where ``max_eig`` is the largest eigenvalue; if negative, discard
            eigenmodes with eigenvalues smaller than ``|svdcut| * max_eig``.
            Note ``|svdcut| < 1``. Default is 1e-12.
        wgts: Setting ``wgts=1`` causes :func:`gvar.svd` to compute
            and return a spectral decomposition of the covariance matrix of
            the modified |GVar|\s, ``gmod``. Setting ``wgts=-1`` results in
            a decomposition of the inverse of the covariance matrix. The
            default value is ``False``, in which case only ``gmod`` is returned.
        add_svdnoise: If ``True``, noise is added to the SVD correction (see
            above).

    Returns:
        A copy ``gmod`` of ``g`` whose correlation matrix is modified by
        SVD cuts. If ``wgts`` is not ``False``,
        a tuple ``(g, i_wgts)`` is returned where ``i_wgts``
        contains a spectral decomposition of ``gmod``'s
        covariance matrix or its inverse.

    Data from the SVD analysis is stored in ``gmod``:

    .. attribute:: gmod.svdcut

        SVD cut used to create ``gmod``.

    .. attribute:: gmod.dof

        Number of independent degrees of freedom left after the
        SVD cut. This is the same as the number initially unless
        ``svdcut < 0`` in which case it may be smaller.

    .. attribute:: gmod.nmod

        Number of modes whose eignevalue was modified by the
        SVD cut.

    .. attribute:: gmod.nblocks

        A dictionary where ``gmod.nblocks[s]`` contains the number of
        block-diagonal ``s``-by-``s`` sub-matrices in the correlation
        matrix.

    .. attribute:: gmod.eigen_range

        Ratio of the smallest to largest eigenvalue before SVD cuts are
        applied (but after rescaling).

    .. attribute:: gmod.logdet

        Logarithm of the determinant of the covariance matrix after SVD
        cuts are applied (excluding any omitted modes when
        ``svdcut < 0`` and any diagonal zero modes).

    .. attribute:: gmod.svdcorrection

        Array or dictionary containing the SVD corrections added to ``g``
        to create ``gmod``: ``gmod = g + gmod.svdcorrection``.
    """
    # replace g by a copy of g
    if hasattr(g,'keys'):
        is_dict = True
        g = BufferDict(g)
    else:
        is_dict = False
        class svdarray(numpy.ndarray):
            def __new__(cls, inputarray):
                obj = numpy.array(g).view(cls)
                return obj
        g = svdarray(g)
    idx_bcov = evalcov_blocks(g.flat)
    g.logdet = 0.0
    svdcorrection = numpy.zeros(len(g.flat), object)
    svdcorrection[:] = gvar(0, 0)
    g.eigen_range = 1.
    g.nmod = 0
    if wgts is not False:
        i_wgts = [([], [])] # 1st entry for all 1x1 blocks
    lost_modes = 0
    g.nblocks = {}
    for idx, block_cov in idx_bcov:
        g.nblocks[len(idx)] = g.nblocks.get(len(idx), 0) + 1
        if len(idx) == 1:
            i = idx[0]
            if block_cov[0, 0] == 0:
                g.logdet = numpy.inf
            else:
                g.logdet += numpy.log(block_cov[0, 0])
            if wgts is not False:
                i_wgts[0][0].append(i)
                i_wgts[0][1].append(block_cov[0, 0] ** (wgts * 0.5))
        else:
            s = SVD(block_cov, svdcut=svdcut, rescale=True, compute_delta=True)
            if s.D is not None:
                g.logdet -= 2 * sum(numpy.log(di) for di in s.D)
            g.logdet += sum(numpy.log(vali) for vali in s.val)
            g.nmod += s.nmod
            if s.delta is not None:
                if add_svdnoise:
                    for vali, valorigi, veci in zip(s.val, s.valorig, s.vec):
                        if vali > valorigi:
                            # add next(raniter(s.delta)) to s.delta in svdcorrection
                            s.delta += (veci / s.D) * (
                                numpy.random.normal(0.0, (vali - valorigi) ** 0.5)
                                )
                svdcorrection[idx] = s.delta
                g.flat[idx] += s.delta
            elif svdcut is not None and svdcut < 0:
                newg = numpy.zeros(len(idx), object)
                for veci in s.vec:
                    veci_D = veci / s.D
                    newg += veci_D * (veci.dot(s.D * g.flat[idx]))
                lost_modes += len(idx) - len(s.vec)
                g.flat[idx] = newg
            if wgts is not False:
                i_wgts.append(
                    (idx, [w for w in s.decomp(wgts)[::-1]])
                    )
            if s.eigen_range < g.eigen_range:
                g.eigen_range = s.eigen_range
    g.nmod += lost_modes
    g.dof = len(g.flat) - lost_modes
    g.svdcut = svdcut

    # repackage svdcorrection
    if is_dict:
        g.svdcorrection = BufferDict(g, buf=svdcorrection)
    else:
        g.svdcorrection = svdcorrection.reshape(g.shape)

    ##### for legacy code (don't use)
    svd.dof = g.dof
    svd.nmod = g.nmod
    svd.eigen_range = g.eigen_range
    svd.logdet = g.logdet
    svd.correction = g.svdcorrection.flat[:]
    svd.nblocks = g.nblocks
    ##### end of legacy code

    # repack into numpy arrays
    if wgts is not False:
        tmp = []
        for iw, wgts in i_wgts:
            tmp.append(
                (numpy.array(iw, numpy.intp), numpy.array(wgts, numpy.double))
                )
        i_wgts = tmp
        return (g, i_wgts)
    else:
        return g
Apply SVD cuts to collection of |GVar|\s in ``g``.

    Standard usage is, for example, ::

        svdcut = ...
        gmod = svd(g, svdcut=svdcut)

    where ``g`` is an array of |GVar|\s or a dictionary containing |GVar|\s
    and/or arrays of |GVar|\s. When ``svdcut>0``, ``gmod`` is
    a copy of ``g`` whose |GVar|\s have been modified to make
    their correlation matrix less singular than that of the
    original ``g``: each eigenvalue ``eig`` of the correlation matrix is
    replaced by ``max(eig, svdcut * max_eig)`` where ``max_eig`` is
    the largest eigenvalue. This SVD cut, which is applied separately
    to each block-diagonal sub-matrix of the correlation matrix,
    increases the variance of the eigenmodes with eigenvalues smaller
    than ``svdcut * max_eig``.

    The modification of ``g``'s covariance matrix is implemented by adding
    (to ``g``) a set of |GVar|\s with zero means::

        gmod = g + gmod.svdcorrection

    where ``gmod.svdcorrection`` is an array/dictionary
    containing the |GVar|\s. If
    parameter ``add_svdnoise=True``,
    noise is included in ``gmod.svdcorrection``,
    ::

        gmod.svdcorrection += gv.sample(gmod.svdcorrection),

    before it is added to ``g``. The noise can be useful for testing fits
    and other applications.

    When ``svdcut`` is negative, eigenmodes of the correlation matrix
    whose eigenvalues are smaller than ``|svdcut| * max_eig`` are dropped
    from the new matrix and the corresponding components of ``g`` are
    zeroed out (that is, replaced by 0(0)) in ``gmod``.

    There is an additional parameter ``wgts`` in :func:`gvar.svd` whose
    default value is ``False``. Setting ``wgts=1`` or ``wgts=-1`` instead
    causes :func:`gvar.svd` to return a tuple ``(gmod, i_wgts)`` where
    ``gmod``  is the modified copy of ``g``, and ``i_wgts`` contains a
    spectral  decomposition of the covariance matrix corresponding to
    the  modified correlation matrix if ``wgts=1``, or a decomposition of its
    inverse if ``wgts=-1``. The first entry ``i, wgts = i_wgts[0]``  specifies
    the diagonal part of the matrix: ``i`` is a list of the indices in
    ``gmod.flat`` corresponding to diagonal elements, and ``wgts ** 2``
    gives the corresponding matrix elements. The second and subsequent
    entries, ``i, wgts = i_wgts[n]`` for ``n > 0``, each correspond
    to block-diagonal sub-matrices, where ``i`` is the list of
    indices corresponding to the block, and ``wgts[j]`` are eigenvectors of
    the sub-matrix rescaled so that ::

        numpy.sum(numpy.outer(wi, wi) for wi in wgts[j]

    is the sub-matrix (``wgts=1``) or its inverse (``wgts=-1``).

    To compute the inverse of the covariance matrix from ``i_wgts``,
    for example, one could use code like::

        gmod, i_wgts = svd(g, svdcut=svdcut, wgts=-1)

        inv_cov = numpy.zeros((n, n), float)
        i, wgts = i_wgts[0]                       # 1x1 sub-matrices
        if len(i) > 0:
            inv_cov[i, i] = numpy.array(wgts) ** 2
        for i, wgts in i_wgts[1:]:                # nxn sub-matrices (n>1)
            for w in wgts:
                inv_cov[i[:, None], i] += numpy.outer(w, w)

    This sets ``inv_cov`` equal to the inverse of the covariance matrix of
    the ``gmod``\s. Similarly, we can  compute the expectation value,
    ``u.dot(inv_cov.dot(v))``, between two vectors (:mod:`numpy` arrays)
    using::

        result = 0.0
        i, wgts = i_wgts[0]                       # 1x1 sub-matrices
        if len(i) > 0:
            result += numpy.sum((u[i] * wgts) * (v[i] * wgts))
        for i, wgts in i_wgts[1:]:                # nxn sub-matrices (n>1)
            result += numpy.sum(wgts.dot(u[i]) * wgts.dot(v[i]))

    where ``result`` is the desired expectation value.

    Args:

        g: An array of |GVar|\s or a dicitionary whose values are
            |GVar|\s and/or arrays of |GVar|\s.
        svdcut (None or float): If positive, replace eigenvalues ``eig``
            of the correlation matrix with ``max(eig, svdcut * max_eig)``
            where ``max_eig`` is the largest eigenvalue; if negative, discard
            eigenmodes with eigenvalues smaller than ``|svdcut| * max_eig``.
            Note ``|svdcut| < 1``. Default is 1e-12.
        wgts: Setting ``wgts=1`` causes :func:`gvar.svd` to compute
            and return a spectral decomposition of the covariance matrix of
            the modified |GVar|\s, ``gmod``. Setting ``wgts=-1`` results in
            a decomposition of the inverse of the covariance matrix. The
            default value is ``False``, in which case only ``gmod`` is returned.
        add_svdnoise: If ``True``, noise is added to the SVD correction (see
            above).

    Returns:
        A copy ``gmod`` of ``g`` whose correlation matrix is modified by
        SVD cuts. If ``wgts`` is not ``False``,
        a tuple ``(g, i_wgts)`` is returned where ``i_wgts``
        contains a spectral decomposition of ``gmod``'s
        covariance matrix or its inverse.

    Data from the SVD analysis is stored in ``gmod``:

    .. attribute:: gmod.svdcut

        SVD cut used to create ``gmod``.

    .. attribute:: gmod.dof

        Number of independent degrees of freedom left after the
        SVD cut. This is the same as the number initially unless
        ``svdcut < 0`` in which case it may be smaller.

    .. attribute:: gmod.nmod

        Number of modes whose eignevalue was modified by the
        SVD cut.

    .. attribute:: gmod.nblocks

        A dictionary where ``gmod.nblocks[s]`` contains the number of
        block-diagonal ``s``-by-``s`` sub-matrices in the correlation
        matrix.

    .. attribute:: gmod.eigen_range

        Ratio of the smallest to largest eigenvalue before SVD cuts are
        applied (but after rescaling).

    .. attribute:: gmod.logdet

        Logarithm of the determinant of the covariance matrix after SVD
        cuts are applied (excluding any omitted modes when
        ``svdcut < 0`` and any diagonal zero modes).

    .. attribute:: gmod.svdcorrection

        Array or dictionary containing the SVD corrections added to ``g``
        to create ``gmod``: ``gmod = g + gmod.svdcorrection``.