from functools import reduce


def pauli_expectation_value(pauli, meas_results):
    """ Compute and return the expectation value of a given Pauli
    based on the measurement outcomes observed in result.


    Parameters:
        pauli (string): A string indicating the Pauli whose expectation value
            we want to compute, e.g. "ZZIZZ"
        meas_results (Dict): A dictionary containing the results of an experiment run on qiskit.
            The key value pairs are computational basis states and number of times that
            state was observed, e.g. {'1001' : 24, '1000' : 36}, etc.

    Returns:
        The expectation value of pauli.        
    """
    pauli_list = list(pauli)
    n_qubits = len(pauli_list)
    n_shots = sum(meas_results.values())

    if any([op not in ['I', 'X', 'Y', 'Z'] for op in pauli_list]):
        raise ValueError("Pauli string must consist only of I, X, Y, or Z.")
    
    # Determine whether the computational basis states in meas_results are +1 or -1
    # eigenstates of the Pauli in question.
    eigenvalues = {}
    for basis_state in meas_results.keys():
        num_z_and_1 = [-1 if (basis_state[bit_idx] == '1' and pauli_list[bit_idx] != 'I') else 1 for bit_idx in range(n_qubits)]
        eigenvalues[basis_state] = reduce(lambda x, y: x*y, num_z_and_1)

    # Count number of +1 and -1 outcomes, i.e. 0 and 1
    num_0_outcomes = sum([meas_results[key] for key in eigenvalues.keys() if eigenvalues[key] == 1])    
    num_1_outcomes = sum([meas_results[key] for key in eigenvalues.keys() if eigenvalues[key] == -1])

    return (num_0_outcomes - num_1_outcomes) / n_shots