Delete `theta` output in MBAR.getFreeEnergyDifferences()?

[hovo1990]

def getFreeEnergyDifferences(self, compute_uncertainty=True, uncertainty_method=None, warning_cutoff=1.0e-10, return_theta=False): """Get the dimensionless free energy differences and uncertainties among all thermodynamic states. Parameters

    compute_uncertainty : bool, optional
        If False, the uncertainties will not be computed (default: True)
    uncertainty_method : string, optional
        Choice of method used to compute asymptotic covariance method,
        or None to use default.  See help for computeAsymptoticCovarianceMatrix()
        for more information on various methods. (default: svd)
    warning_cutoff : float, optional
        Warn if squared-uncertainty is negative and larger in magnitude
        than this number (default: 1.0e-10)
    return_theta : bool, optional
        Whether or not to return the theta matrix.  Can be useful for complicated differences.
    Returns
    -------
    Deltaf_ij : np.ndarray, float, shape=(K, K)
        Deltaf_ij[i,j] is the estimated free energy difference
    dDeltaf_ij : np.ndarray, float, shape=(K, K)
        If compute_uncertainty==True, 
        dDeltaf_ij[i,j] is the estimated statistical uncertainty 
        (one standard deviation) in Deltaf_ij[i,j].  Otherwise None.
    Theta_ij : np.ndarray, float, shape=(K, K)
        The theta_matrix if return_theta==True, otherwise None.

    Notes
    -----
    Computation of the covariance matrix may take some time for large K.
    The reported statistical uncertainty should, in the asymptotic limit, reflect one standard deviation for the normal distribution of the estimate.
    The true free energy difference should fall within the interval [-df, +df] centered on the estimate 68% of the time, and within
    the interval [-2 df, +2 df] centered on the estimate 95% of the time.
    This will break down in cases where the number of samples is not large enough to reach the asymptotic normal limit.
    See Section III of Reference [1].
    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)
    >>> Deltaf_ij, dDeltaf_ij, Theta_ij = mbar.getFreeEnergyDifferences()
    """
    Deltaf_ij, dDeltaf_ij, Theta_ij = None, None, None  # By default, returns None for dDelta and Theta

    # Compute free energy differences.
    f_i = np.matrix(self.f_k)
    Deltaf_ij = f_i - f_i.transpose()

    # zero out numerical error for thermodynamically identical states
    self._zerosamestates(Deltaf_ij)

    Deltaf_ij = np.array(Deltaf_ij)  # Convert from np.matrix to np.array

    if compute_uncertainty or return_theta:
        # Compute asymptotic covariance matrix.
        Theta_ij = self._computeAsymptoticCovarianceMatrix(
            np.exp(self.Log_W_nk), self.N_k, method=uncertainty_method)

    if compute_uncertainty:
        dDeltaf_ij = self._ErrorOfDifferences(Theta_ij, warning_cutoff=warning_cutoff)
        # zero out numerical error for thermodynamically identical states
        self._zerosamestates(dDeltaf_ij)
        # Return matrix of free energy differences and uncertainties.
        dDeltaf_ij = np.array(dDeltaf_ij)

    if return_theta == True:
         return Deltaf_ij, dDeltaf_ij, Theta_ij
    else:
         return Deltaf_ij, dDeltaf_ij

?

[kyleabeauchamp]

You aren't looking at the most recent code.

https://github.com/choderalab/pymbar/blob/master/pymbar/mbar.py#L479

[mrshirts]

Reviving this to avoid having multiple issues. Probably the best thing to do is to SAVE the theta matrix when calculating the uncertainties, and then write a method to retrieve it separately. It can be a big matrix, so one might not want to have it around in some cases, but we can deal with that a bit later (perhaps have saving be optional, but default on, then it can be listed as one of the things that can be done to save memory)

[jchodera]

Reviving this to avoid having multiple issues. Probably the best thing to do is to SAVE the theta matrix when calculating the uncertainties, and then write a method to retrieve it separately.

Sounds great!

[Lnaden]

I think we can close this as Theta is still an optional return and we'll be changing the return style in 4.0 anyways.