def nabla_f_loglik(self, F_estimand):
    f_matrix = np.reshape(F_estimand, (self.B, self.K), order='F')
    f_full_matrix = gp_inflate_duplicate(f_matrix, self.block_assignment_numeric, self.N, self.K)
    Z = self.current_Z
    f_gradient = np.zeros(f_matrix.shape)
    temp_matrix = (1 - np.exp((f_full_matrix - logsumexp(f_full_matrix, axis=1)[:, np.newaxis])))
    inflated_output_matrix = np.multiply(Z, temp_matrix)
    f_gradient += gp_deflate_sum(inflated_output_matrix, self.block_assignment_numeric, self.N, self.B, self.K)
    for k in range(self.K):
        f_gradient[:, k] += self.f_prior.get_nabla_f(f=f_matrix[:, k])
    return f_gradient.flatten(order='F')
def nabla_f_loglik(self, F_estimand):
    f_matrix = np.reshape(F_estimand, (self.B, self.K), order='F')
    f_full_matrix = gp_inflate_duplicate(f_matrix, self.block_assignment_numeric, self.N, self.K)
    Z = self.current_Z
    f_gradient = np.zeros(f_matrix.shape)
    temp_matrix = (1 - np.exp((f_full_matrix - logsumexp(f_full_matrix, axis=1)[:, np.newaxis])))
    inflated_output_matrix = np.multiply(Z, temp_matrix)
    f_gradient += gp_deflate_sum(inflated_output_matrix, self.block_assignment_numeric, self.N, self.B, self.K)
    for k in range(self.K):
        f_gradient[:, k] += self.f_prior.get_nabla_f(f=f_matrix[:, k])
    return f_gradient.flatten(order='F')