code: def measure_pure_state(prep_program, reference_state, quantum_resource,
                         variance_bound=1.0E-6):
      
      num_qubits = len(prep_program.get_qubits())
      amplitudes_to_measure = list(range(2 ** num_qubits))
  
      amplitudes = measure_wf_coefficients(prep_program, amplitudes_to_measure,
                                           reference_state,
                                           quantum_resource,
                                           variance_bound=variance_bound)
      wavefunction = np.asarray(amplitudes)
      return wavefunction.reshape((-1, 1))
docs: Measure the coefficients of the pure state
  
      :param prep_program: pyQuil program to prepare the state
      :param reference_state: Integer of the computational basis state to use as
                              a reference
      :param quantum_resource: An instance of a quantum abstract machine
      :param variance_bound: Default 1.0E-6.  variance of the monte carlo
                             estimator for the non-hermitian operator
      :return: an estimate of the wavefunction as a numpy.ndarray