Generalized dot product for xarray objects. Like np.einsum, but
    provides a simpler interface based on array dimensions.

    Parameters
    ----------
    arrays: DataArray (or Variable) objects
        Arrays to compute.
    dims: str or tuple of strings, optional
        Which dimensions to sum over.
        If not speciified, then all the common dimensions are summed over.
    **kwargs: dict
        Additional keyword arguments passed to numpy.einsum or
        dask.array.einsum

    Returns
    -------
    dot: DataArray

    Examples
    --------

    >>> da_a = xr.DataArray(np.arange(3 * 4).reshape(3, 4), dims=['a', 'b'])
    >>> da_b = xr.DataArray(np.arange(3 * 4 * 5).reshape(3, 4, 5),
    >>>                     dims=['a', 'b', 'c'])
    >>> da_c = xr.DataArray(np.arange(5 * 6).reshape(5, 6), dims=['c', 'd'])
    >>>
    >>> xr.dot(da_a, da_b, dims=['a', 'b']).dims
    ('c', )
    >>> xr.dot(da_a, da_b, dims=['a']).dims
    ('b', 'c')
    >>> xr.dot(da_a, da_b, da_c, dims=['b', 'c']).dims
    ('a', 'd')

def dot(*arrays, dims=None, **kwargs):
    """Generalized dot product for xarray objects. Like np.einsum, but
    provides a simpler interface based on array dimensions.

    Parameters
    ----------
    arrays: DataArray (or Variable) objects
        Arrays to compute.
    dims: str or tuple of strings, optional
        Which dimensions to sum over.
        If not speciified, then all the common dimensions are summed over.
    **kwargs: dict
        Additional keyword arguments passed to numpy.einsum or
        dask.array.einsum

    Returns
    -------
    dot: DataArray

    Examples
    --------

    >>> da_a = xr.DataArray(np.arange(3 * 4).reshape(3, 4), dims=['a', 'b'])
    >>> da_b = xr.DataArray(np.arange(3 * 4 * 5).reshape(3, 4, 5),
    >>>                     dims=['a', 'b', 'c'])
    >>> da_c = xr.DataArray(np.arange(5 * 6).reshape(5, 6), dims=['c', 'd'])
    >>>
    >>> xr.dot(da_a, da_b, dims=['a', 'b']).dims
    ('c', )
    >>> xr.dot(da_a, da_b, dims=['a']).dims
    ('b', 'c')
    >>> xr.dot(da_a, da_b, da_c, dims=['b', 'c']).dims
    ('a', 'd')
    """
    from .dataarray import DataArray
    from .variable import Variable

    if any(not isinstance(arr, (Variable, DataArray)) for arr in arrays):
        raise TypeError('Only xr.DataArray and xr.Variable are supported.'
                        'Given {}.'.format([type(arr) for arr in arrays]))

    if len(arrays) == 0:
        raise TypeError('At least one array should be given.')

    if isinstance(dims, str):
        dims = (dims, )

    common_dims = set.intersection(*[set(arr.dims) for arr in arrays])
    all_dims = []
    for arr in arrays:
        all_dims += [d for d in arr.dims if d not in all_dims]

    einsum_axes = 'abcdefghijklmnopqrstuvwxyz'
    dim_map = {d: einsum_axes[i] for i, d in enumerate(all_dims)}

    if dims is None:
        # find dimensions that occur more than one times
        dim_counts = Counter()
        for arr in arrays:
            dim_counts.update(arr.dims)
        dims = tuple(d for d, c in dim_counts.items() if c > 1)

    dims = tuple(dims)  # make dims a tuple

    # dimensions to be parallelized
    broadcast_dims = tuple(d for d in all_dims
                           if d in common_dims and d not in dims)
    input_core_dims = [[d for d in arr.dims if d not in broadcast_dims]
                       for arr in arrays]
    output_core_dims = [tuple(d for d in all_dims if d not in
                              dims + broadcast_dims)]

    # older dask than 0.17.4, we use tensordot if possible.
    if isinstance(arr.data, dask_array_type):
        import dask
        if LooseVersion(dask.__version__) < LooseVersion('0.17.4'):
            if len(broadcast_dims) == 0 and len(arrays) == 2:
                axes = [[arr.get_axis_num(d) for d in arr.dims if d in dims]
                        for arr in arrays]
                return apply_ufunc(duck_array_ops.tensordot, *arrays,
                                   dask='allowed',
                                   input_core_dims=input_core_dims,
                                   output_core_dims=output_core_dims,
                                   kwargs={'axes': axes})

    # construct einsum subscripts, such as '...abc,...ab->...c'
    # Note: input_core_dims are always moved to the last position
    subscripts_list = ['...' + ''.join([dim_map[d] for d in ds]) for ds
                       in input_core_dims]
    subscripts = ','.join(subscripts_list)
    subscripts += '->...' + ''.join([dim_map[d] for d in output_core_dims[0]])

    # subscripts should be passed to np.einsum as arg, not as kwargs. We need
    # to construct a partial function for apply_ufunc to work.
    func = functools.partial(duck_array_ops.einsum, subscripts, **kwargs)
    result = apply_ufunc(func, *arrays,
                         input_core_dims=input_core_dims,
                         output_core_dims=output_core_dims,
                         dask='allowed')
    return result.transpose(*[d for d in all_dims if d in result.dims])