docs: def computeCovarianceOfSums(self, d_ij, K, a):
  
          """
          We wish to calculate the variance of a weighted sum of free energy differences.
          for example ``var(\sum a_i df_i)``.
  
          We explicitly lay out the calculations for four variables (where each variable
          is a logarithm of a partition function), then generalize.
  
          The uncertainty in the sum of two weighted differences is
  
          .. code-block:: none
  
              var(a1(f_i1 - f_j1) + a2(f_i2 - f_j2)) =
                  a1^2 var(f_i1 - f_j1) +
                  a2^2 var(f_i2 - f_j2) +
                  2 a1 a2 cov(f_i1 - f_j1, f_i2 - f_j2)
              cov(f_i1 - f_j1, f_i2 - f_j2) =
                  cov(f_i1,f_i2) -
                  cov(f_i1,f_j2) -
                  cov(f_j1,f_i2) +
                  cov(f_j1,f_j2)
  
          call:
  
          .. code-block:: none
  
              f_i1 = a
              f_j1 = b
              f_i2 = c
              f_j2 = d
              a1^2 var(a-b) + a2^2 var(c-d) + 2a1a2 cov(a-b,c-d)
  
          we want ``2cov(a-b,c-d) = 2cov(a,c)-2cov(a,d)-2cov(b,c)+2cov(b,d)``,
          since ``var(x-y) = var(x) + var(y) - 2cov(x,y)``,
          then, ``2cov(x,y) = -var(x-y) + var(x) + var(y)``. So, we get
  
          .. code-block:: none
  
              2cov(a,c) = -var(a-c) + var(a) + var(c)
              -2cov(a,d) = +var(a-d) - var(a) - var(d)
              -2cov(b,c) = +var(b-c) - var(b) - var(c)
              2cov(b,d) = -var(b-d) + var(b) + var(d)
  
          adding up, finally :
  
          .. code-block:: none
  
              2cov(a-b,c-d) = 2cov(a,c)-2cov(a,d)-2cov(b,c)+2cov(b,d) =
                  - var(a-c) + var(a-d) + var(b-c) - var(b-d)
  
              a1^2 var(a-b)+a2^2 var(c-d)+2a1a2cov(a-b,c-d) =
                  a1^2 var(a-b)+a2^2 var(c-d)+a1a2 [-var(a-c)+var(a-d)+var(b-c)-var(b-d)]
  
              var(a1(f_i1 - f_j1) + a2(f_i2 - f_j2)) =
                  a1^2 var(f_i1 - f_j1) + a2^2 var(f_i2 - f_j2) + 2a1 a2 cov(f_i1 - f_j1, f_i2 - f_j2)
              = a1^2 var(f_i1 - f_j1) + a2^2 var(f_i2 - f_j2) + a1 a2 [-var(f_i1 - f_i2) + var(f_i1 - f_j2) + var(f_j1-f_i2) - var(f_j1 - f_j2)]
  
          assume two arrays of free energy differences, and and array of constant vectors a.
          we want the variance ``var(\sum_k a_k (f_i,k - f_j,k))`` Each set is separated from the other by an offset K
          same process applies with the sum, with the single var terms and the pair terms
  
          Parameters
          ----------
          d_ij : a matrix of standard deviations of the quantities f_i - f_j
          K : The number of states in each 'chunk', has to be constant
          outputs : KxK variance matrix for the sums or differences ``\sum a_i df_i``
          """
  
          # todo: vectorize this.
          var_ij = np.square(d_ij)
          d2 = np.zeros([K,K],float)
          n = len(a)
          for i in range(K):
              for j in range(K):
                  for k in range(n):
                      d2[i,j] +=  a[k]**2 * var_ij[i+k*K,j+k*K]
                      for l in range(n):
                          d2[i,j] +=  a[k] * a[l] * (-var_ij[i+k*K,i+l*K] + var_ij[i+k*K,j+l*K] + var_ij[j+k*K,i+l*K] - var_ij[j+k*K,j+l*K])
  
          return np.sqrt(d2)