def __init__(self, u_kn, N_k, maximum_iterations=10000, relative_tolerance=1.0e-7, verbose=False, initial_f_k=None,
                 solver_protocol=None, initialize='zeros', x_kindices=None, **kwargs):
        """Initialize multistate Bennett acceptance ratio (MBAR) on a set of simulation data.

        Upon initialization, the dimensionless free energies for all states are computed.
        This may take anywhere from seconds to minutes, depending upon the quantity of data.
        After initialization, the computed free energies may be obtained by a call to :func:`getFreeEnergyDifferences`,
        or expectation at any state of interest can be computed by calls to :func:`computeExpectations`.

        Parameters
        ----------
        u_kn : np.ndarray, float, shape=(K, N_max)
            ``u_kn[k,n]`` is the reduced potential energy of uncorrelated
            configuration n evaluated at state ``k``.
        u_kln : np.ndarray, float, shape (K, L, N_max)
            If the simulation is in form ``u_kln[k,l,n]`` it is converted to ``u_kn`` format

            .. code-block:: python

                u_kn = [ u_1(x_1) u_1(x_2) u_1(x_3) . . . u_1(x_n)
                         u_2(x_1) u_2(x_2) u_2(x_3) . . . u_2(x_n)
                                                    . . .
                         u_k(x_1) u_k(x_2) u_k(x_3) . . . u_k(x_n)]

        N_k :  np.ndarray, int, shape=(K)
            ``N_k[k]`` is the number of uncorrelated snapshots sampled from state ``k``.
            Some may be zero, indicating that there are no samples from that state.

            We assume that the states are ordered such that the first ``N_k``
            are from the first state, the 2nd ``N_k`` the second state, and so
            forth. This only becomes important for BAR -- MBAR does not
            care which samples are from which state.  We should eventually
            allow this assumption to be overwritten by parameters passed
            from above, once ``u_kln`` is phased out.

        maximum_iterations : int, optional
            Set to limit the maximum number of iterations performed (default 1000)
        relative_tolerance : float, optional
            Set to determine the relative tolerance convergence criteria (default 1.0e-6)
        verbosity : bool, optional
            Set to True if verbose debug output is desired (default False)
        initial_f_k : np.ndarray, float, shape=(K), optional
            Set to the initial dimensionless free energies to use as a
            guess (default None, which sets all f_k = 0)
        solver_protocol : list(dict) or None, optional, default=None
            List of dictionaries to define a sequence of solver algorithms
            and options used to estimate the dimensionless free energies.
            See `pymbar.mbar_solvers.solve_mbar()` for details.  If None,
            use the developers best guess at an appropriate algorithm.

            The default will try to solve with an adaptive solver algorithm
            which alternates between self-consistent iteration and
            Newton-Raphson, where the method with the smallest
            gradient is chosen to improve numerical stability.

        initialize : 'zeros' or 'BAR', optional, Default: 'zeros'
            If equal to 'BAR', use BAR between the pairwise state to
            initialize the free energies.  Eventually, should specify a path;
            for now, it just does it zipping up the states.

            The 'BAR' option works best when the states are ordered such that adjacent states
            maximize the overlap between states. Its up to the user
            to arrange the states in such an order, or at least close to such an order.
            If you are uncertain what the order of states should be, or if it does not make
            sense to think of states as adjacent, then choose 'zeroes'.

            (default: 'zeros', unless specific values are passed in.)
        x_kindices
            Which state is each x from?  Usually doesn't matter, but does for BAR. We assume the samples
            are in ``K`` order (the first ``N_k[0]`` samples are from the 0th state, the next ``N_k[1]`` samples from
            the 1st state, and so forth.

        Notes
        -----
        The reduced potential energy ``u_kn[k,n] = u_k(x_{ln})``, where the reduced potential energy ``u_l(x)`` is
        defined (as in the text) by:
        ``u_k(x) = beta_k [ U_k(x) + p_k V(x) + mu_k' n(x) ]``
        where

        ``beta_k = 1/(kB T_k)`` is the inverse temperature of condition ``k``, where kB is Boltzmann's constant

        ``U_k(x)`` is the potential energy function for state ``k``

        ``p_k`` is the pressure at state ``k`` (if an isobaric ensemble is specified)

        ``V(x)`` is the volume of configuration ``x``

        ``mu_k`` is the M-vector of chemical potentials for the various species, if a (semi)grand ensemble is
        specified, and ``'`` denotes transpose

        ``n(x)`` is the M-vector of numbers of the various molecular species for configuration ``x``,
        corresponding to the chemical potential components of ``mu_m``.

        ``x_n`` indicates that the samples are from ``k`` different simulations of the ``n`` states. These simulations
        need only be a subset of the k states.

        The configurations ``x_ln`` must be uncorrelated.  This can be ensured by subsampling a correlated timeseries
        with a period larger than the statistical inefficiency, which can be estimated from the potential energy
        timeseries ``{u_k(x_ln)}_{n=1}^{N_k}`` using the provided utility
        :func:`pymbar.timeseries.statisticalInefficiency`.
        See the help for this function for more information.

        Examples
        --------

        >>> from pymbar import testsystems
        >>> (x_n, u_kn, N_k, s_n) = testsystems.HarmonicOscillatorsTestCase().sample(mode='u_kn')
        >>> mbar = MBAR(u_kn, N_k)

        """
        for key, val in kwargs.items():
            print("Warning: parameter {}={} is unrecognized and unused.".format(key, val))

        # Store local copies of necessary data.
        # N_k[k] is the number of samples from state k, some of which might be zero.
        self.N_k = np.array(N_k, dtype=np.int64)
        self.N = np.sum(self.N_k)

        # Get dimensions of reduced potential energy matrix, and convert to KxN form if needed.
        if len(np.shape(u_kn)) == 3:
            self.K = np.shape(u_kn)[1]  # need to set self.K, and it's the second index
            u_kn = kln_to_kn(u_kn, N_k=self.N_k)

        # u_kn[k,n] is the reduced potential energy of sample n evaluated at state k
        self.u_kn = np.array(u_kn, dtype=np.float64)

        K, N = np.shape(u_kn)

        if verbose:
            print("K (total states) = %d, total samples = %d" % (K, N))

        if np.sum(self.N_k) != N:
            raise ParameterError(
                'The sum of all N_k must equal the total number of samples (length of second dimension of u_kn.')

        # Store local copies of other data
        self.K = K  # number of thermodynamic states energies are evaluated at
        # N = \sum_{k=1}^K N_k is the total number of samples
        self.N = N  # maximum number of configurations

        # if not defined, identify from which state each sample comes from.
        if x_kindices is not None:
            self.x_kindices = x_kindices
        else:
            self.x_kindices = np.arange(N, dtype=np.int64)
            Nsum = 0
            for k in range(K):
                self.x_kindices[Nsum:Nsum+self.N_k[k]] = k
                Nsum += self.N_k[k]

        # verbosity level -- if True, will print extra debug information
        self.verbose = verbose

        # perform consistency checks on the data.

        # if, for any set of data, all reduced potential energies are the same,
        # they are probably the same state.  We check to within
        # relative_tolerance.

        self.samestates = []
        if self.verbose:
            for k in range(K):
                for l in range(k):
                    diffsum = 0
                    uzero = u_kn[k, :] - u_kn[l, :]
                    diffsum += np.dot(uzero, uzero)
                    if (diffsum < relative_tolerance):
                        self.samestates.append([k, l])
                        self.samestates.append([l, k])
                        print('')
                        print('Warning: states %d and %d have the same energies on the dataset.' % (l, k))
                        print('They are therefore likely to to be the same thermodynamic state.  This can occasionally cause')
                        print('numerical problems with computing the covariance of their energy difference, which must be')
                        print('identically zero in any case. Consider combining them into a single state.')
                        print('')

        # Print number of samples from each state.
        if self.verbose:
            print("N_k = ")
            print(self.N_k)

        # Determine list of k indices for which N_k != 0
        self.states_with_samples = np.where(self.N_k != 0)[0]
        self.states_with_samples = self.states_with_samples.astype(np.int64)

        # Number of states with samples.
        self.K_nonzero = self.states_with_samples.size
        if verbose:
            print("There are %d states with samples." % self.K_nonzero)

        # Initialize estimate of relative dimensionless free energy of each state to zero.
        # Note that f_k[0] will be constrained to be zero throughout.
        # this is default
        self.f_k = np.zeros([self.K], dtype=np.float64)

        # If an initial guess of the relative dimensionless free energies is
        # specified, start with that.
        if initial_f_k is not None:
            if self.verbose:
                print("Initializing f_k with provided initial guess.")
            # Cast to np array.
            initial_f_k = np.array(initial_f_k, dtype=np.float64)
            # Check shape
            if initial_f_k.shape != self.f_k.shape:
                raise ParameterError(
                    "initial_f_k must be a %d-dimensional np array." % self.K)
            # Initialize f_k with provided guess.
            self.f_k = initial_f_k
            if self.verbose:
                print(self.f_k)
            # Shift all free energies such that f_0 = 0.
            self.f_k[:] = self.f_k[:] - self.f_k[0]
        else:
            # Initialize estimate of relative dimensionless free energies.
            self._initializeFreeEnergies(verbose, method=initialize)

            if self.verbose:
                print("Initial dimensionless free energies with method %s" % (initialize))
                print("f_k = ")
                print(self.f_k)

        if solver_protocol is None:
            solver_protocol = ({'method': None},)
        for solver in solver_protocol:
            if 'options' not in solver:
                solver['options'] = dict()
            if 'verbose' not in solver['options']:
                # should add in other ways to get information out of the scipy solvers, not just adaptive,
                # which might involve passing in different combinations of options, and passing out other strings.
                solver['options']['verbose'] = self.verbose

        self.f_k = mbar_solvers.solve_mbar_for_all_states(self.u_kn, self.N_k, self.f_k, solver_protocol)
        self.Log_W_nk = mbar_solvers.mbar_log_W_nk(self.u_kn, self.N_k, self.f_k)

        # Print final dimensionless free energies.
        if self.verbose:
            print("Final dimensionless free energies")
            print("f_k = ")
            print(self.f_k)

        if self.verbose:
            print("MBAR initialization complete.")
def __init__ ( self , u_kn , N_k , maximum_iterations = 10000 , relative_tolerance = 1.0e-7 , verbose = False , initial_f_k = None , solver_protocol = None , initialize = 'zeros' , x_kindices = None , ** kwargs ) : for key , val in kwargs . items ( ) : print ( "Warning: parameter {}={} is unrecognized and unused." . format ( key , val ) ) self . N_k = np . array ( N_k , dtype = np . int64 ) self . N = np . sum ( self . N_k ) if len ( np . shape ( u_kn ) ) == 3 : self . K = np . shape ( u_kn ) [ 1 ] u_kn = kln_to_kn ( u_kn , N_k = self . N_k ) self . u_kn = np . array ( u_kn , dtype = np . float64 ) K , N = np . shape ( u_kn ) if verbose : print ( "K (total states) = %d, total samples = %d" % ( K , N ) ) if np . sum ( self . N_k ) != N : raise ParameterError ( 'The sum of all N_k must equal the total number of samples (length of second dimension of u_kn.' ) self . K = K self . N = N if x_kindices is not None : self . x_kindices = x_kindices else : self . x_kindices = np . arange ( N , dtype = np . int64 ) Nsum = 0 for k in range ( K ) : self . x_kindices [ Nsum : Nsum + self . N_k [ k ] ] = k Nsum += self . N_k [ k ] self . verbose = verbose self . samestates = [ ] if self . verbose : for k in range ( K ) : for l in range ( k ) : diffsum = 0 uzero = u_kn [ k , : ] - u_kn [ l , : ] diffsum += np . dot ( uzero , uzero ) if ( diffsum < relative_tolerance ) : self . samestates . append ( [ k , l ] ) self . samestates . append ( [ l , k ] ) print ( '' ) print ( 'Warning: states %d and %d have the same energies on the dataset.' % ( l , k ) ) print ( 'They are therefore likely to to be the same thermodynamic state.  This can occasionally cause' ) print ( 'numerical problems with computing the covariance of their energy difference, which must be' ) print ( 'identically zero in any case. Consider combining them into a single state.' ) print ( '' ) if self . verbose : print ( "N_k = " ) print ( self . N_k ) self . states_with_samples = np . where ( self . N_k != 0 ) [ 0 ] self . states_with_samples = self . states_with_samples . astype ( np . int64 ) self . K_nonzero = self . states_with_samples . size if verbose : print ( "There are %d states with samples." % self . K_nonzero ) self . f_k = np . zeros ( [ self . K ] , dtype = np . float64 ) if initial_f_k is not None : if self . verbose : print ( "Initializing f_k with provided initial guess." ) initial_f_k = np . array ( initial_f_k , dtype = np . float64 ) if initial_f_k . shape != self . f_k . shape : raise ParameterError ( "initial_f_k must be a %d-dimensional np array." % self . K ) self . f_k = initial_f_k if self . verbose : print ( self . f_k ) self . f_k [ : ] = self . f_k [ : ] - self . f_k [ 0 ] else : self . _initializeFreeEnergies ( verbose , method = initialize ) if self . verbose : print ( "Initial dimensionless free energies with method %s" % ( initialize ) ) print ( "f_k = " ) print ( self . f_k ) if solver_protocol is None : solver_protocol = ( { 'method' : None } , ) for solver in solver_protocol : if 'options' not in solver : solver [ 'options' ] = dict ( ) if 'verbose' not in solver [ 'options' ] : solver [ 'options' ] [ 'verbose' ] = self . verbose self . f_k = mbar_solvers . solve_mbar_for_all_states ( self . u_kn , self . N_k , self . f_k , solver_protocol ) self . Log_W_nk = mbar_solvers . mbar_log_W_nk ( self . u_kn , self . N_k , self . f_k ) if self . verbose : print ( "Final dimensionless free energies" ) print ( "f_k = " ) print ( self . f_k ) if self . verbose : print ( "MBAR initialization complete." )
Initialize multistate Bennett acceptance ratio (MBAR) on a set of simulation data.

        Upon initialization, the dimensionless free energies for all states are computed.
        This may take anywhere from seconds to minutes, depending upon the quantity of data.
        After initialization, the computed free energies may be obtained by a call to :func:`getFreeEnergyDifferences`,
        or expectation at any state of interest can be computed by calls to :func:`computeExpectations`.

        Parameters
        ----------
        u_kn : np.ndarray, float, shape=(K, N_max)
            ``u_kn[k,n]`` is the reduced potential energy of uncorrelated
            configuration n evaluated at state ``k``.
        u_kln : np.ndarray, float, shape (K, L, N_max)
            If the simulation is in form ``u_kln[k,l,n]`` it is converted to ``u_kn`` format

            .. code-block:: python

                u_kn = [ u_1(x_1) u_1(x_2) u_1(x_3) . . . u_1(x_n)
                         u_2(x_1) u_2(x_2) u_2(x_3) . . . u_2(x_n)
                                                    . . .
                         u_k(x_1) u_k(x_2) u_k(x_3) . . . u_k(x_n)]

        N_k :  np.ndarray, int, shape=(K)
            ``N_k[k]`` is the number of uncorrelated snapshots sampled from state ``k``.
            Some may be zero, indicating that there are no samples from that state.

            We assume that the states are ordered such that the first ``N_k``
            are from the first state, the 2nd ``N_k`` the second state, and so
            forth. This only becomes important for BAR -- MBAR does not
            care which samples are from which state.  We should eventually
            allow this assumption to be overwritten by parameters passed
            from above, once ``u_kln`` is phased out.

        maximum_iterations : int, optional
            Set to limit the maximum number of iterations performed (default 1000)
        relative_tolerance : float, optional
            Set to determine the relative tolerance convergence criteria (default 1.0e-6)
        verbosity : bool, optional
            Set to True if verbose debug output is desired (default False)
        initial_f_k : np.ndarray, float, shape=(K), optional
            Set to the initial dimensionless free energies to use as a
            guess (default None, which sets all f_k = 0)
        solver_protocol : list(dict) or None, optional, default=None
            List of dictionaries to define a sequence of solver algorithms
            and options used to estimate the dimensionless free energies.
            See `pymbar.mbar_solvers.solve_mbar()` for details.  If None,
            use the developers best guess at an appropriate algorithm.

            The default will try to solve with an adaptive solver algorithm
            which alternates between self-consistent iteration and
            Newton-Raphson, where the method with the smallest
            gradient is chosen to improve numerical stability.

        initialize : 'zeros' or 'BAR', optional, Default: 'zeros'
            If equal to 'BAR', use BAR between the pairwise state to
            initialize the free energies.  Eventually, should specify a path;
            for now, it just does it zipping up the states.

            The 'BAR' option works best when the states are ordered such that adjacent states
            maximize the overlap between states. Its up to the user
            to arrange the states in such an order, or at least close to such an order.
            If you are uncertain what the order of states should be, or if it does not make
            sense to think of states as adjacent, then choose 'zeroes'.

            (default: 'zeros', unless specific values are passed in.)
        x_kindices
            Which state is each x from?  Usually doesn't matter, but does for BAR. We assume the samples
            are in ``K`` order (the first ``N_k[0]`` samples are from the 0th state, the next ``N_k[1]`` samples from
            the 1st state, and so forth.

        Notes
        -----
        The reduced potential energy ``u_kn[k,n] = u_k(x_{ln})``, where the reduced potential energy ``u_l(x)`` is
        defined (as in the text) by:
        ``u_k(x) = beta_k [ U_k(x) + p_k V(x) + mu_k' n(x) ]``
        where

        ``beta_k = 1/(kB T_k)`` is the inverse temperature of condition ``k``, where kB is Boltzmann's constant

        ``U_k(x)`` is the potential energy function for state ``k``

        ``p_k`` is the pressure at state ``k`` (if an isobaric ensemble is specified)

        ``V(x)`` is the volume of configuration ``x``

        ``mu_k`` is the M-vector of chemical potentials for the various species, if a (semi)grand ensemble is
        specified, and ``'`` denotes transpose

        ``n(x)`` is the M-vector of numbers of the various molecular species for configuration ``x``,
        corresponding to the chemical potential components of ``mu_m``.

        ``x_n`` indicates that the samples are from ``k`` different simulations of the ``n`` states. These simulations
        need only be a subset of the k states.

        The configurations ``x_ln`` must be uncorrelated.  This can be ensured by subsampling a correlated timeseries
        with a period larger than the statistical inefficiency, which can be estimated from the potential energy
        timeseries ``{u_k(x_ln)}_{n=1}^{N_k}`` using the provided utility
        :func:`pymbar.timeseries.statisticalInefficiency`.
        See the help for this function for more information.

        Examples
        --------

        >>> from pymbar import testsystems
        >>> (x_n, u_kn, N_k, s_n) = testsystems.HarmonicOscillatorsTestCase().sample(mode='u_kn')
        >>> mbar = MBAR(u_kn, N_k)