def estimate_umbrella_sampling(
    us_trajs, us_dtrajs, us_centers, us_force_constants, md_trajs=None, md_dtrajs=None, kT=None,
    maxiter=10000, maxerr=1.0E-15, save_convergence_info=0,
    estimator='wham', lag=1, dt_traj='1 step', init=None, init_maxiter=10000, init_maxerr=1.0E-8,
    width=None, **kwargs):
    r"""
    This function acts as a wrapper for ``tram()``, ``dtram()``, ``mbar()``, and ``wham()`` and
    handles the calculation of bias energies (``bias``) and thermodynamic state trajectories
    (``ttrajs``) when the data comes from umbrella sampling and (optional) unbiased simulations.

    Parameters
    ----------
    us_trajs : list of N arrays, each of shape (T_i, d)
        List of arrays, each having T_i rows, one for each time step, and d columns where d is the
        dimensionality of the subspace in which umbrella sampling was applied. Often d=1, and thus
        us_trajs will be a list of 1d-arrays.
    us_dtrajs : list of N int arrays, each of shape (T_i,)
        The integers are indexes in 0,...,n-1 enumerating the n discrete states or the bins the
        umbrella sampling trajectory is in at any time.
    us_centers : list of N floats or d-dimensional arrays of floats
        List or array of N center positions. Each position must be a d-dimensional vector. For 1d
        umbrella sampling, one can simply pass a list of centers, e.g. [-5.0, -4.0, -3.0, ... ].
    us_force_constants : list of N floats or d- or dxd-dimensional arrays of floats
        The force constants used in the umbrellas, unit-less (e.g. kT per squared length unit). For
        multidimensional umbrella sampling, the force matrix must be used.
    md_trajs : list of M arrays, each of shape (T_i, d), optional, default=None
        Unbiased molecular dynamics simulations; format like us_trajs.
    md_dtrajs : list of M int arrays, each of shape (T_i,)
        The integers are indexes in 0,...,n-1 enumerating the n discrete states or the bins the
        unbiased trajectory is in at any time.
    kT : float or None, optional, default=None
        Use this attribute if the supplied force constants are NOT unit-less; kT must have the same
        energy unit as the force constants.
    maxiter : int, optional, default=10000
        The maximum number of self-consistent iterations before the estimator exits unsuccessfully.
    maxerr : float, optional, default=1.0E-15
        Convergence criterion based on the maximal free energy change in a self-consistent
        iteration step.
    save_convergence_info : int, optional, default=0
        Every save_convergence_info iteration steps, store the actual increment
        and the actual loglikelihood; 0 means no storage.
    estimator : str, optional, default='wham'
        Specify one of the available estimators

        | 'wham':   use WHAM
        | 'mbar':   use MBAR
        | 'dtram':  use the discrete version of TRAM
        | 'tram':  use TRAM
    lag : int or list of int, optional, default=1
        Integer lag time at which transitions are counted. Providing a list of lag times will
        trigger one estimation per lag time.
    dt_traj : str, optional, default='1 step'
        Description of the physical time corresponding to the lag. May be used by analysis
        algorithms such as plotting tools to pretty-print the axes. By default '1 step', i.e.
        there is no physical time unit.  Specify by a number, whitespace and unit. Permitted
        units are (* is an arbitrary string):

        |  'fs',   'femtosecond*'
        |  'ps',   'picosecond*'
        |  'ns',   'nanosecond*'
        |  'us',   'microsecond*'
        |  'ms',   'millisecond*'
        |  's',    'second*'
    init : str, optional, default=None
        Use a specific initialization for the self-consistent iteration:

        | None:    use a hard-coded guess for free energies and Lagrangian multipliers
        | 'wham':  perform a short WHAM estimate to initialize the free energies (only with dtram)
        | 'mbar':  perform a short MBAR estimate to initialize the free energies (only with tram)
    init_maxiter : int, optional, default=10000
        The maximum number of self-consistent iterations during the initialization.
    init_maxerr : float, optional, default=1.0E-8
        Convergence criterion for the initialization.
    width : array-like of float, optional, default=None
        Specify periodicity for individual us_traj dimensions. Each positive entry will make the
        corresponding feature periodic and use the given value as width. None/zero values will be
        treated as non-periodic.
    **kwargs : dict, optional
        You can use this to pass estimator-specific named parameters to the chosen estimator, which
        are not already coverd by ``estimate_umbrella_sampling()``.

    Returns
    -------
    A :class:`MultiThermModel <pyemma.thermo.models.multi_therm.MultiThermModel>` or :class:`MEMM <pyemma.thermo.models.memm.MEMM>` object or list thereof
        The requested estimator/model object, i.e., WHAM, MBAR, DTRAM or TRAM. If multiple lag times
        are given, a list of objects is returned (one MEMM per lag time).

    Example
    -------
    We look at a 1D umbrella sampling simulation with two umbrellas at 1.1 and 1.3 on the reaction
    coordinate with spring constant of 1.0; additionally, we have two unbiased simulations.

    We start with a joint clustering and use TRAM for the estimation:

    >>> from pyemma.coordinates import cluster_regspace as regspace
    >>> from pyemma.thermo import estimate_umbrella_sampling as estimate_us
    >>> import numpy as np
    >>> us_centers = [1.1, 1.3]
    >>> us_force_constants = [1.0, 1.0]
    >>> us_trajs = [np.array([1.0, 1.1, 1.2, 1.1, 1.0, 1.1]).reshape((-1, 1)), np.array([1.3, 1.2, 1.3, 1.4, 1.4, 1.3]).reshape((-1, 1))]
    >>> md_trajs = [np.array([0.9, 1.0, 1.1, 1.2, 1.3, 1.4]).reshape((-1, 1)), np.array([1.5, 1.4, 1.3, 1.4, 1.4, 1.5]).reshape((-1, 1))]
    >>> cluster = regspace(data=us_trajs+md_trajs, max_centers=10, dmin=0.15)
    >>> us_dtrajs = cluster.dtrajs[:2]
    >>> md_dtrajs = cluster.dtrajs[2:]
    >>> centers = cluster.clustercenters
    >>> tram = estimate_us(us_trajs, us_dtrajs, us_centers, us_force_constants, md_trajs=md_trajs, md_dtrajs=md_dtrajs, estimator='tram', lag=1)
    >>> tram.f # doctest: +ELLIPSIS
    array([ 0.63...,  1.60...,  1.31...])

    See :class:`MultiThermModel <pyemma.thermo.models.multi_therm.MultiThermModel>`
    or :class:`MEMM <pyemma.thermo.models.memm.MEMM>` for a full documentation.

    .. autoclass:: pyemma.thermo.models.multi_therm.MultiThermModel
        :members:
        :undoc-members:

        .. rubric:: Methods

        .. autoautosummary:: pyemma.thermo.models.multi_therm.MultiThermModel
           :methods:

        .. rubric:: Attributes

        .. autoautosummary:: pyemma.thermo.models.multi_therm.MultiThermModel
            :attributes:

    .. autoclass:: pyemma.thermo.models.memm.MEMM
        :members:
        :undoc-members:

        .. rubric:: Methods

        .. autoautosummary:: pyemma.thermo.models.memm.MEMM
           :methods:

        .. rubric:: Attributes

        .. autoautosummary:: pyemma.thermo.models.memm.MEMM
            :attributes:

    """
    from .util import get_umbrella_sampling_data as _get_umbrella_sampling_data
    # sanity checks
    if estimator not in ['wham', 'mbar', 'dtram', 'tram']:
        raise ValueError("unsupported estimator: %s" % estimator)
    if not isinstance(us_trajs, (list, tuple)):
        raise ValueError("The parameter us_trajs must be a list of numpy.ndarray objects")
    if not isinstance(us_centers, (list, tuple)):
        raise ValueError(
            "The parameter us_centers must be a list of floats or numpy.ndarray objects")
    if not isinstance(us_force_constants, (list, tuple)):
        raise ValueError(
            "The parameter us_force_constants must be a list of floats or numpy.ndarray objects")
    if len(us_trajs) != len(us_centers):
        raise ValueError("Unmatching number of umbrella sampling trajectories and centers: %d!=%d" \
            % (len(us_trajs), len(us_centers)))
    if len(us_trajs) != len(us_force_constants):
        raise ValueError(
            "Unmatching number of umbrella sampling trajectories and force constants: %d!=%d" \
                % (len(us_trajs), len(us_force_constants)))
    if len(us_trajs) != len(us_dtrajs):
            raise ValueError(
                "Number of continuous and discrete umbrella sampling trajectories does not " + \
                "match: %d!=%d" % (len(us_trajs), len(us_dtrajs)))
    i = 0
    for traj, dtraj in zip(us_trajs, us_dtrajs):
        if traj.shape[0] != dtraj.shape[0]:
            raise ValueError(
                "Lengths of continuous and discrete umbrella sampling trajectories with " + \
                "index %d does not match: %d!=%d" % (i, len(us_trajs), len(us_dtrajs)))
        i += 1
    if md_trajs is not None:
        if not isinstance(md_trajs, (list, tuple)):
            raise ValueError("The parameter md_trajs must be a list of numpy.ndarray objects")
        if md_dtrajs is None:
            raise ValueError("You have provided md_trajs, but md_dtrajs is None")
    if md_dtrajs is None:
        md_dtrajs = []
    else:
        if md_trajs is None:
            raise ValueError("You have provided md_dtrajs, but md_trajs is None")
        if len(md_trajs) != len(md_dtrajs):
            raise ValueError(
                "Number of continuous and discrete unbiased trajectories does not " + \
                "match: %d!=%d" % (len(md_trajs), len(md_dtrajs)))
        i = 0
        for traj, dtraj in zip(md_trajs, md_dtrajs):
            if traj.shape[0] != dtraj.shape[0]:
                raise ValueError(
                    "Lengths of continuous and discrete unbiased trajectories with " + \
                    "index %d does not match: %d!=%d" % (i, len(md_trajs), len(md_dtrajs)))
            i += 1
    # data preparation
    ttrajs, btrajs, umbrella_centers, force_constants, unbiased_state = _get_umbrella_sampling_data(
        us_trajs, us_centers, us_force_constants, md_trajs=md_trajs, kT=kT, width=width)
    estimator_obj = None
    # estimation
    if estimator == 'wham':
        estimator_obj = wham(
            ttrajs, us_dtrajs + md_dtrajs,
            _get_averaged_bias_matrix(btrajs, us_dtrajs + md_dtrajs),
            maxiter=maxiter, maxerr=maxerr,
            save_convergence_info=save_convergence_info, dt_traj=dt_traj)
    elif estimator == 'mbar':
        allowed_keys = ['direct_space']
        parsed_kwargs = dict([(i, kwargs[i]) for i in allowed_keys if i in kwargs])
        estimator_obj = mbar(
            ttrajs, us_dtrajs + md_dtrajs, btrajs,
            maxiter=maxiter, maxerr=maxerr, save_convergence_info=save_convergence_info,
            dt_traj=dt_traj, **parsed_kwargs)
    elif estimator == 'dtram':
        allowed_keys = ['count_mode', 'connectivity']
        parsed_kwargs = dict([(i, kwargs[i]) for i in allowed_keys if i in kwargs])
        estimator_obj = dtram(
            ttrajs, us_dtrajs + md_dtrajs,
            _get_averaged_bias_matrix(btrajs, us_dtrajs + md_dtrajs),
            lag, unbiased_state=unbiased_state,
            maxiter=maxiter, maxerr=maxerr, save_convergence_info=save_convergence_info,
            dt_traj=dt_traj, init=init, init_maxiter=init_maxiter, init_maxerr=init_maxerr,
            **parsed_kwargs)
    elif estimator == 'tram':
        allowed_keys = [
            'count_mode', 'connectivity', 'connectivity_factor','nn',
            'direct_space', 'N_dtram_accelerations', 'equilibrium', 'overcounting_factor', 'callback']
        parsed_kwargs = dict([(i, kwargs[i]) for i in allowed_keys if i in kwargs])
        estimator_obj = tram(
            ttrajs, us_dtrajs + md_dtrajs, btrajs, lag, unbiased_state=unbiased_state,
            maxiter=maxiter, maxerr=maxerr, save_convergence_info=save_convergence_info,
            dt_traj=dt_traj, init=init, init_maxiter=init_maxiter, init_maxerr=init_maxerr,
            **parsed_kwargs)
    # adding thermodynamic state information and return results
    try:
        estimator_obj.umbrella_centers = umbrella_centers
        estimator_obj.force_constants = force_constants
    except AttributeError:
        for obj in estimator_obj:
            obj.umbrella_centers = umbrella_centers
            obj.force_constants = force_constants
    return estimator_obj