Add initialise simulator instruction to to statevector_state about qiskit-aer HOT 21 CLOSED

@chriseclectic Can you please assign @ShellyGarion to this and remove me?

from qiskit-aer.

@chriseclectic : It seems that Qiskit/qiskit#1896 is an issue and not a PR.
Is there a relevant PR in Terra? Could we develop the code in Aer independently of Terra?
At least for testing we may need the code in Terra.

from qiskit-aer.

@ajavadia Can you point shelly to your branch implementing this?

The code in Aer can be developed independently of the terra piece. I've added parsing for the qobj instruction in PR #88, so if you manually create a qobj instruction you can test.

In the Op class the initialize instruction is stored so that:

Op init;  // initialize op
init.type = OpType::initialize
init.qubits = [a, b, c] // n-qubits to initialize
init.params = psi  // n-qubit vector

To add it you can modify the statevector_state class to add a case for this new op type which can implement the reset on the necessary qubits, then you may need to add a method to the qubitvector class for updating the statevector. You can look at the lambda functions for applying n-qubit matrices to see how the indexing works.

You will need a general N-qubit version of something that for the two qubit case does:

Suppose are going to initalize qubit-0 (for example) to state [a, b] for a system in arbitrary two qubit state.

1. Reset qubit-0: leaves system in state |phi>\otimes|0>. As vector this is [n, m, 0, 0].
2. Extract |phi> from block: phi = [n, m].
3. Tensor on new psi = [a,b], which involves setting blocks to a * [n,m], b * [n,m] = [a*n, a*m, b*n, b*m]

from qiskit-aer.

@chriseclectic - regarding your 3 suggested steps:

1. Can we use the function apply_reset in statevector_state.hpp here?
2. Is there a function for this extraction/projection?
   If not - should this function be added to statevector_state.hpp or to qubitvector.hpp?
3. Is there a function for computing the tensor product of the state vectors?
   Should there be a special function for handling the indexes or is it a part of the tensor calculation?
   (e.g. if there are 3 qubits q1,q2,q3 then initializing q1,q2 has different indexing from initializing q1,q3).
   Again, should these functions be added to statevector_state.hpp or to qubitvector.hpp?

Basically, I think that the most complicated part here is to do the indexing properly, especially since the functions in qubitvector.hpp are written in an optimized/vectorized way.

from qiskit-aer.

Hi @ShellyGarion So i think you can do something like this for the update

• Assume N qubit state, and we are initializing M qubits that have been reset to |0,..,0⟩ ( you can use the existing reset instruction for this).
• The statevector is now in a state Swap( |psi⟩⊗|0,..,0⟩ ), where |psi⟩ is the post-reset state of the K=N-M non-initialized qubits (Note: the swap is just so i could write as a tensor product easily, but they could be any position qubits)
• If we get the index for the reset qubits being in zero, we can index the remaining K qubits using the index0_* functions from statevector/indexing.hpp. By looping over these indexes in the underlying statevector we get the psi[k] component from data_[j] with j = index0_state(sorted_init_qubits, k).
• We can get all the indices for the block for a given k using the indexes_* functions (note that the non 0 parts are currently zero due to the reset.
• We can use this to set those values of the statevector to the component for the new initialized state.

Something like this:

template <typename data_t>
template <size_t N>
void QubitVector<data_t>::initialize_component(const areg_t<N> &qs, const cvector_t &state) {

  // Lambda function for initializing component
  auto lambda = [&](const areg_t<1ULL << N> &inds, const cvector_t &_state)->void {
    const uint_t DIM = 1ULL << N;
    complex_t cache = inds[0];  // the k-th component of non-initialized vector
    for (size_t i = 0; i < DIM; i++) {
      data_[inds[i]] = cache * _state[i];  // set component to psi[k] * state[i]
    }
  // Use the lambda function
  apply_matrix_lambda(lambda, qs, state);
}

from qiskit-aer.

In order to be able to debug and add printouts to statevector_state.hpp code, I needed to comment the line:
#validate_qobj_against_schema(self._qobj)
in: aerjob.py - since Terra does not support yet reset of more than one qubit.
Now I can see the printouts.

from qiskit-aer.

@chriseclectic: A few questions:

• What should the above initialize_component function do?
  Should it do both the extraction of psi and the tensor product of psi and the new state ?
  Or only the tensor product?
• What should it get as input? Are these the qubits that should be initialized?
  (In #117, I changed: const areg_t<N> &qs to: const reg_t &qubits following #109)
• What is the difference between the functions: index0, indexes and the lambda function for initializing component above?
• Should we use the functions index0 and/or indexes in addition to initialize_component?

from qiskit-aer.

@ShellyGarion

1. The above code snippet should initialize the specified qubits to a desired statevector (leaving the other qubits in their current state) assuming the qubits being initialized have already been reset to the |000..> state.
2. It gets as input qubits being initialized (eg: [2,4]) and the vector they will be set to (eg [psi00, psi01, psi10, psi11], and yes should be changed to const reg_t & since the #109 merge.
3. index0 returns only the integer representation of a number of bits set to zero inserted into an arbitrary bit string. Eg: for qubits 0,2 in a state k = ba ( ba = 00 => k=0, etc).

 indexes0([1], k) -> int(b0a)
 indexes0([1,3], k) -> int(0b0a)

indexes returns the array of all bit values for the specified qubits (eg for previous example [int(b0a), int(b1a)], if it were two qubits inserted say at 1,3 it would be:

indexes([1,3], [1,3], k) -> [int(0b0a), int(0b1a), int(1b0a), (1b1a)]`

If the qubits were passed in reverse order it would swap qubit position in the list:

indexes([3,1], [1,3], k) -> `[int(0b0a), int(1b0a), int(0b1a), (1b1a)]

4. I think the initialize_component, along with explicit reset of the qubits, should be all you need (it handles the indexing in its lambda function).

from qiskit-aer.

@chriseclectic - Thank you for the detailed explanation. I added some more explanatory comments to the code.
Now the code compiles and running, and I need to verify that it does the correct calculation. I needed to make a few changes in your function above (define N, and change the types of qubits and inds):

template <typename data_t>
void QubitVector<data_t>::initialize_component(const reg_t &qubits, const cvector_t &state) {

  // Lambda function for initializing component
  const size_t N = qubits.size();
  auto lambda = [&](indexes_t inds, const cvector_t &_state)->void {
    const uint_t DIM = 1ULL << N;
    complex_t cache = inds[0];  // the k-th component of non-initialized vector
    for (size_t i = 0; i < DIM; i++) {
      data_[inds[i]] = cache * _state[i];  // set component to psi[k] * state[i]
    }    // (where psi is is the post-reset state of the non-initialized qubits)
   };
  // Use the lambda function
  apply_matrix_lambda(lambda, qubits, state);
}

So, now the new function in the Statevector class is simply:

template <class statevec_t>
void State<statevec_t>::apply_initialize(const reg_t &qubits,
                                         const cvector_t &params,
                                         RngEngine &rng) {
   // Apply reset to qubits
   apply_reset(qubits, rng);
   // Apply initialize_component
   BaseState::qreg_.initialize_component(qubits, params);
}

from qiskit-aer.

If you check out the master branch of Terra, you can test this.

from qiskit import *
qr = QuantumRegister(2)

# Method 1
circ1 = QuantumCircuit(qr)
circ1.initialize([0, 0, 0, 1], qr[:])

# Method 2
circ2 = QuantumCircuit(qr)
circ2.initialize([0, 0, 0, 1], [qr[0], qr[1]])

# Method 3
circ3 = QuantumCircuit(qr)
circ3.initialize([0, 1], [qr[0]])
circ3.initialize([0, 1], [qr[1]])


# Implementation
circ0 = QuantumCircuit(qr)
circ0.reset(qr[0])
circ0.x(qr[0])
circ0.reset(qr[1])
circ0.x(qr[1])

# assemble into a qobj
from qiskit.compiler import assemble_circuits
qobj = assemble_circuits([circ0, circ1, circ2, circ3])

This qobj contains 4 circuits, all of which should yield the same statevector.

from qiskit-aer.

Indeed, here is the output of @ajavadia code above:
statevector of circ1: [0.+0.j 0.+0.j 0.+0.j 1.+0.j]
statevector of circ2: [0.+0.j 0.+0.j 0.+0.j 1.+0.j]
statevector of circ3: [0.+0.j 0.+0.j 0.+0.j 1.+0.j]
statevector of circ0: [0.+0.j 0.+0.j 0.+0.j 1.+0.j]

from qiskit-aer.

It seems that the above API assumes that the initialization is only to a real state.
Is it possible/required to initialize an unreal/imagionary state?
I tried the following API:
circ3.initialize([0, 0+1.j], [qr[0]])
and got at error:
TypeError: can't convert complex to float

from qiskit-aer.

There is currently a difference between the API of initialize in Aer/C++ and in Terra/Pyhton.

In Terra/Python qobj:
QasmQobjInstruction(name='initialize', params=[0, 1.00000000000000], qubits=[0])
Namely, params is a list (vector) of complex numbers

and in Aer/C++ qobj:
QasmQobjInstruction(name='initialize', params=[[0, 0], [1, 0]], qubits=[0])
params has type cvector_t which is a vector of pairs of doubles:
using complex_t = std::complex<double>;
using cvector_t = std::vector<complex_t>;

from qiskit-aer.

The JSON parsing on the C++ side should be able to parse either a complex or real vector into a complex vector, so as long as it works for both cases I'm happy on our end (this true for Aer/C++ parsing of all complex vectors or complex matrices) . If you change Ali's examples to use dtype=complex in the states it should serialize as you expect, and you can test with an actual complex state too.

However the real case might have some issues for schema validation on the Terra end so might be worth opening an issue on the Terra GitHub for that @ajavadia. In that case it's because on the Terra side the way it unrolls Numpy arrays doesnt assume anything about dtype. It just uses array.tolist() and then on a second pass if values are complex it unrolls to real pairs. So you can either make the schema allow real or complex vectors, or you would have to convert the input to a complex array as part of the instruction (using params = np.array(params, dtype=complex) or something along those lines)

from qiskit-aer.

@chriseclectic - can you perhaps provide an example?

For example, I use the following code:

qr = QuantumRegister(3)

circ4 = QuantumCircuit(qr)
circ4.h(qr[0])
circ4.h(qr[1])
circ4.h(qr[2])

circ5 = QuantumCircuit(qr)
circ5.h(qr[0])
circ5.h(qr[1])
circ5.h(qr[2])
circ5.reset (qr[0])

circ6 = QuantumCircuit(qr)
circ6.h(qr[0])
circ6.h(qr[1])
circ6.h(qr[2])
circ6.initialize([0, 1], [qr[0]])

Then the output is:

statevector of circ4: [0.35355339+0.j 0.35355339+0.j 0.35355339+0.j 0.35355339+0.j
 0.35355339+0.j 0.35355339+0.j 0.35355339+0.j 0.35355339+0.j]

statevector of circ5: [0.5+0.j 0. +0.j 0.5+0.j 0. +0.j 0.5+0.j 0. +0.j 0.5+0.j 0. +0.j]

statevector of circ6: [-0.35355339+0.j  0.35355339+0.j -0.35355339+0.j  0.35355339+0.j
 -0.35355339+0.j  0.35355339+0.j -0.35355339+0.j  0.35355339+0.j]

so the last statevector is wrong - it should have been:
[0. +0.j 0.5+0.j 0. +0.j 0.5+0.j 0. +0.j 0.5+0.j 0. +0.j 0.5+0.j]

However, when I generate the qobj and explicitly provide the following instruction:
qobj.experiments[0].instructions.append(QasmQobjInstruction(name='initialize', qubits=[0], params=[[0,0],[1,0]]))
then I get the correct statevector.

from qiskit-aer.

@chriseclectic - I tried the following code as you suggested:

circ6 = QuantumCircuit(qr)
circ6.h(qr[0])
circ6.h(qr[1])
circ6.h(qr[2])
circ6.initialize(np.array([0,1], dtype=complex), [qr[0]])

But I still got the wrong statevector:
statevector of circ6: [-0.35355339+0.j 0.35355339+0.j -0.35355339+0.j 0.35355339+0.j -0.35355339+0.j 0.35355339+0.j -0.35355339+0.j 0.35355339+0.j]

from qiskit-aer.

@ShellyGarion complex should work. see the example here which runs successfully:
https://github.com/Qiskit/qiskit-terra/blob/master/examples/python/initialize.py

@chriseclectic i think terra puts the instruction name as "initialize" in the qobj, not "init". Let's keep it that way.

from qiskit-aer.

The following code works OK:

qr = QuantumRegister(3)
cr = ClassicalRegister(3)
circ = QuantumCircuit(qr, cr)
circ.h(qr[0])
circ.h(qr[1])
circ.h(qr[2])
simulator = Aer.get_backend('statevector_simulator')
qobj = qiskit.compile(circ, backend=simulator)
qobj.experiments[0].instructions.append(QasmQobjInstruction(name='initialize', qubits=[0], params=[[0,0],[1,0]]))
result = simulator.run(qobj).result()
statevector = result.get_statevector(circ)

from qiskit-aer.

@chriseclectic can we close this one?

from qiskit-aer.

@chriseclectic - where should we add the initialize gate to the aer tutorial?
Should we add it to this notebook:
https://github.com/Qiskit/qiskit-tutorials/blob/master/qiskit/aer/aer_provider.ipynb

from qiskit-aer.

@ShellyGarion For the tutorial I would suggest rewriting the section Starting simulation with a custom initial state to for the QasmSimulator and StatevectorSimulator to use the initialize function instead of the backend_options={"initial_statevector': vec} (You can leave the starting the unitary simulator in a custom initial state alone since this case isn't handled by the initialize function).

It is worth mentioning that the simulator supports this option directly for efficiency, but that it can also be unrolled to standard gates for execution on actual devices (unlike the old initial_statevector backend option).

Closing this issue now since the instruction is supported in Aer.

from qiskit-aer.