{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.7.6"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":4104,"databundleVersionId":46661,"sourceType":"competition"},{"sourceId":2822362,"sourceType":"datasetVersion","datasetId":1725862}],"dockerImageVersionId":30527,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install qiskit\n!pip install qiskit-aer","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:49:17.155004Z","iopub.execute_input":"2023-08-20T22:49:17.155631Z","iopub.status.idle":"2023-08-20T22:49:46.647133Z","shell.execute_reply.started":"2023-08-20T22:49:17.155563Z","shell.execute_reply":"2023-08-20T22:49:46.645759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import torch\nfrom torch.autograd import Function\nimport torch.optim as optim\nimport torch.nn as nn\n\n# if torch.cuda.is_available():\n#     device = torch.device(\"cuda:0\")\n#     print(\"Running on the GPU\")\n# else:\n#     device = torch.device(\"cpu\")\n#     print(\"Running on the CPU\")\n    \n# torch.device(\"cpu\")","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:48:19.842119Z","iopub.execute_input":"2023-08-20T22:48:19.842777Z","iopub.status.idle":"2023-08-20T22:48:23.774226Z","shell.execute_reply.started":"2023-08-20T22:48:19.842735Z","shell.execute_reply":"2023-08-20T22:48:23.772803Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from qiskit import execute\nfrom qiskit.circuit import Parameter,ControlledGate\nfrom qiskit import Aer\nimport qiskit\nimport numpy as np","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:48:26.745359Z","iopub.execute_input":"2023-08-20T22:48:26.746307Z","iopub.status.idle":"2023-08-20T22:48:26.753625Z","shell.execute_reply.started":"2023-08-20T22:48:26.746265Z","shell.execute_reply":"2023-08-20T22:48:26.751856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from tqdm import tqdm","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:48:36.316586Z","iopub.execute_input":"2023-08-20T22:48:36.31709Z","iopub.status.idle":"2023-08-20T22:48:36.327424Z","shell.execute_reply.started":"2023-08-20T22:48:36.31702Z","shell.execute_reply":"2023-08-20T22:48:36.326365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from matplotlib import pyplot as plt\n%matplotlib inline","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:48:37.907275Z","iopub.execute_input":"2023-08-20T22:48:37.907727Z","iopub.status.idle":"2023-08-20T22:48:37.914213Z","shell.execute_reply.started":"2023-08-20T22:48:37.90769Z","shell.execute_reply":"2023-08-20T22:48:37.913298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.random.seed = 42\n\nNUM_QUBITS = 4\nNUM_SHOTS = 10000\nSHIFT = np.pi/2\nLEARNING_RATE = 0.01\nMOMENTUM = 0.5\n\nSIMULATOR = Aer.get_backend('qasm_simulator')","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:49:54.550207Z","iopub.execute_input":"2023-08-20T22:49:54.550612Z","iopub.status.idle":"2023-08-20T22:49:54.558238Z","shell.execute_reply.started":"2023-08-20T22:49:54.550582Z","shell.execute_reply":"2023-08-20T22:49:54.557123Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create list of all possible outputs of quantum circuit (2**NUM_QUBITS possible)\nimport itertools\ndef create_QC_OUTPUTS():\n    # for this circuit, there are only 2 outputs: '0' and '1'\n    return ['0', '1']\n\nQC_OUTPUTS = create_QC_OUTPUTS()\nprint(QC_OUTPUTS)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:49:46.736509Z","iopub.execute_input":"2023-08-20T22:49:46.736836Z","iopub.status.idle":"2023-08-20T22:49:46.742994Z","shell.execute_reply.started":"2023-08-20T22:49:46.736806Z","shell.execute_reply":"2023-08-20T22:49:46.741698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Define function to translate Q-Circuit parameters from pytorch back to QISKIT","metadata":{}},{"cell_type":"markdown","source":"## 3. Contruct QuantumCircuit QFT Class","metadata":{}},{"cell_type":"code","source":"class QiskitCircuit():\n    \n    def __init__(self, n_qubits, backend, shots):\n        # --- Circuit definition ---\n        self.circuit = qiskit.QuantumCircuit(n_qubits, 1)\n        self.n_qubits = n_qubits\n        self.thetas ={k : Parameter('Theta'+str(k))for k in range(self.n_qubits)}\n        \n        all_qubits = [i for i in range(n_qubits)]\n        self.circuit.h(0)\n        for k in range(n_qubits-1):\n            self.circuit.cx(k, k+1)\n            \n        self.circuit.barrier()\n        for k in range(n_qubits):\n            self.circuit.ry(self.thetas[k], k)\n        self.circuit.barrier()\n\n        for i in reversed(range(n_qubits-1)):\n            self.circuit.cx(i, i+1)\n        self.circuit.h(0)\n        self.circuit.measure(0,0)\n        \n        \n        # ---------------------------\n        \n        self.backend = backend\n        self.shots = shots\n        \n#             check = perc\n#             for i in range(nr_qubits):\n#                 check *= (float(key[i])-1/2)*2\n#             expects += check   \n        \n    def N_qubit_expectation_Z(self,counts, shots, nr_qubits):\n        expects = np.zeros(len(QC_OUTPUTS))\n        for k in range(len(QC_OUTPUTS)):\n            key = QC_OUTPUTS[k]\n            perc = counts.get(key, 0) /shots\n            expects[k] = perc\n        return expects\n    \n    def run(self, i):\n        params = i\n#         print('params = {}'.format(len(params)))\n        backend = Aer.get_backend('qasm_simulator')\n    \n        job_sim = execute(self.circuit,\n                              self.backend,\n                              shots=self.shots,\n                              parameter_binds = [{self.thetas[k] : params[k].item() for k in range(NUM_QUBITS)}])\n#         \n        result_sim = job_sim.result()\n        counts = result_sim.get_counts(self.circuit)\n        return self.N_qubit_expectation_Z(counts,self.shots, 1)","metadata":{"ExecuteTime":{"end_time":"2019-10-01T16:09:30.59873Z","start_time":"2019-10-01T16:09:30.567861Z"},"execution":{"iopub.status.busy":"2023-08-20T22:57:41.53179Z","iopub.execute_input":"2023-08-20T22:57:41.53345Z","iopub.status.idle":"2023-08-20T22:57:41.548878Z","shell.execute_reply.started":"2023-08-20T22:57:41.533386Z","shell.execute_reply":"2023-08-20T22:57:41.547241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Creating a Test Circuit and Measuring Its Expectation","metadata":{}},{"cell_type":"code","source":"circuit = QiskitCircuit(NUM_QUBITS, SIMULATOR, NUM_SHOTS)\ntest_params = [np.pi / 2**k for k in range(NUM_QUBITS)]\ntest_names = ['pi/{}'.format(2**k) for k in range (NUM_QUBITS)]\njob = execute(circuit, backend)\n\nprint('Expected value for rotation {}: {}'.format(test_names, circuit.run(torch.Tensor(test_params))))\n# circuit.circuit.draw(output='mpl', filename='MNIST01-bell/Figures/{}-qubit circuit bell.jpg'.format(NUM_QUBITS))","metadata":{"execution":{"iopub.status.busy":"2023-08-20T23:02:37.234749Z","iopub.execute_input":"2023-08-20T23:02:37.235159Z","iopub.status.idle":"2023-08-20T23:02:37.244011Z","shell.execute_reply.started":"2023-08-20T23:02:37.235129Z","shell.execute_reply":"2023-08-20T23:02:37.24233Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_names","metadata":{"execution":{"iopub.status.busy":"2023-08-20T23:02:57.963166Z","iopub.execute_input":"2023-08-20T23:02:57.963608Z","iopub.status.idle":"2023-08-20T23:02:57.97136Z","shell.execute_reply.started":"2023-08-20T23:02:57.963575Z","shell.execute_reply":"2023-08-20T23:02:57.969853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### TorchCircuit()\n\nA pytorch layer always has two functions. One for the forward pass and one for the backward pass. The forward pass simply takes the Quantum Circuits variational parameters from the previous pytorch layer and runs the circuit on the defined hardware (defined in `QiskitCircuit.run()`) and returns the measurements from the quantum hardware.\nThese measurements will be the inputs of the next pytorch layer.\n\nThe backward pass returns the gradients of the quantum circuit. In this case here it is finite difference.\n\nthe `forward_tensor` is saved from the forward pass. So we just have to do one evaluation of the Q-Circuit in the backpass for the finite difference.\n\nThe `gradient` variable here is as well hard coded to 3 parameters. This should be updated in the future and made more general.\n\nThe loop `for k in range(len(input_numbers)):` goes through all the parameters (in this case 3), and shifts them by a small $\\epsilon$. Then it runs the circuit and takes the diefferences of the ouput for the parameters $\\Theta$ and $\\Theta + \\epsilon$. This is the finite difference. ","metadata":{}},{"cell_type":"code","source":"class TorchCircuit(Function):    \n\n    @staticmethod\n    def forward(ctx, i):\n        if not hasattr(ctx, 'QiskitCirc'):\n            ctx.QiskitCirc = QiskitCircuit(NUM_QUBITS, SIMULATOR, shots=NUM_SHOTS)\n            \n        exp_value = ctx.QiskitCirc.run(i)\n        \n        result = torch.tensor([exp_value])\n        \n        \n        ctx.save_for_backward(result, i)\n        \n        return result\n    \n    @staticmethod\n    def backward(ctx, grad_output):\n        \n        forward_tensor, i = ctx.saved_tensors\n#         print('forward_tensor = {}'.format(forward_tensor))\n        input_numbers = i\n#         print('input_numbers = {}'.format(input_numbers))\n        gradients = torch.Tensor()\n        \n        for k in range(NUM_QUBITS):\n            shift_right = input_numbers.detach().clone()\n            shift_right[k] = shift_right[k] + SHIFT\n            shift_left = input_numbers.detach().clone()\n            shift_left[k] = shift_left[k] - SHIFT\n            \n#             print('shift_right = {}, shift_left = {}'.format(shift_right, shift_left))\n            \n            expectation_right = ctx.QiskitCirc.run(shift_right)\n            expectation_left  = ctx.QiskitCirc.run(shift_left)\n#             print('expectation_right = {}, \\nexpectation_left = {}'.format(expectation_right, expectation_left))\n            \n            gradient = torch.tensor([expectation_right]) - torch.tensor([expectation_left])\n            # rescale gradient\n#             gradient = gradient / torch.norm(gradient)\n#             print('gradient for k={}: {}'.format(k, gradient))\n            gradients = torch.cat((gradients, gradient.float()))\n            \n        result = torch.Tensor(gradients)\n#         print('gradients = {}'.format(result))\n#         print('grad_output = {}'.format(grad_output))\n\n        return (result.float() * grad_output.float()).T","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:52:57.361195Z","iopub.execute_input":"2023-08-20T22:52:57.361726Z","iopub.status.idle":"2023-08-20T22:52:57.375395Z","shell.execute_reply.started":"2023-08-20T22:52:57.361691Z","shell.execute_reply":"2023-08-20T22:52:57.374127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x = torch.tensor([np.pi/4]*NUM_QUBITS, requires_grad=True)\n\nqc = TorchCircuit.apply\ny1 = qc(x)\nprint('y1 after quantum layer: {}'.format(y1))\ny1 = nn.Linear(2,1)(y1.float())\ny1.backward()\nprint('x.grad = {}'.format(x.grad))","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:53:03.233236Z","iopub.execute_input":"2023-08-20T22:53:03.233814Z","iopub.status.idle":"2023-08-20T22:53:03.550321Z","shell.execute_reply.started":"2023-08-20T22:53:03.23377Z","shell.execute_reply":"2023-08-20T22:53:03.548412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test the Quantum Circuit's Gradient Descent\n\nFirst, we want the \"neural net\" consisting of just the quantum circuit (with its 4 inputs and 4 outputs) and a linear layer (from 4 inputs to 1 output) that scales measurement 1 by 1, measurement 2 by 2, etc., until it converges to a target value (-1). So, we define a cost function where the cost is defined as the square distance from the target value.\n\n`x` is the initialization of the parameters. Here, every angle in the quantum circuit starts at $\\pi/4$. We should see that the loss eventually goes down.","metadata":{}},{"cell_type":"code","source":"qc = TorchCircuit.apply\n\ndef cost(x):\n    target = -1\n    expval = qc(x)[0]\n    # simple linear layer: average all outputs of quantum layer\n#     print(expval)\n    val = sum([(i+1)*expval[i] for i in range(2)]) / 2\n#     print(val)\n    return torch.abs(val - target) ** 2, expval\n\nx = torch.tensor(test_params, requires_grad=True)\nopt = torch.optim.Adam([x], lr=.05)\n\nnum_epoch = 100\n\nloss_list = []\nexpval_list = []\n\nfor i in tqdm(range(num_epoch)):\n# for i in range(num_epoch):\n    opt.zero_grad()\n    loss, expval = cost(x)\n    loss.backward()\n    opt.step()\n    loss_list.append(loss.item())\n    expval_list.append(expval)\n\nplt.plot(loss_list)","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:53:23.475352Z","iopub.execute_input":"2023-08-20T22:53:23.475837Z","iopub.status.idle":"2023-08-20T22:53:23.791955Z","shell.execute_reply.started":"2023-08-20T22:53:23.475802Z","shell.execute_reply":"2023-08-20T22:53:23.790465Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### MNIST in pytorch","metadata":{}},{"cell_type":"code","source":"import torch\nimport torch.nn as nn\nimport torch.nn.functional as F\nimport torch.optim as optim","metadata":{"execution":{"iopub.status.busy":"2023-08-20T22:53:31.565141Z","iopub.execute_input":"2023-08-20T22:53:31.565564Z","iopub.status.idle":"2023-08-20T22:53:31.572375Z","shell.execute_reply.started":"2023-08-20T22:53:31.565529Z","shell.execute_reply":"2023-08-20T22:53:31.570541Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load MNIST (0-1) Dataset\n\n**Training Data**","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport torchvision\nfrom torchvision import datasets, transforms\n\n# Concentrating on the first 100 samples\nn_samples = 150\n\nX_train = datasets.MNIST(root='./data', train=True, download=True,\n                         transform=transforms.Compose([transforms.ToTensor()]))\n\n# Leaving only labels 0 and 1 \nidx = np.append(np.where(X_train.targets == 0)[0][:n_samples], \n                np.where(X_train.targets == 1)[0][:n_samples])\n\nX_train.data = X_train.data[idx]\nX_train.targets = X_train.targets[idx]\n\n\ntrain_loader = torch.utils.data.DataLoader(X_train, batch_size=1, shuffle=True, pin_memory=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Testing Data**","metadata":{}},{"cell_type":"code","source":"n_samples = 200\n\nX_test = datasets.MNIST(root='./data', train=False, download=True,\n                        transform=transforms.Compose([transforms.ToTensor()]))\n\nidx = np.append(np.where(X_test.targets == 0)[0][n_samples:], \n                np.where(X_test.targets == 1)[0][n_samples:])\n\nX_test.data = X_test.data[idx]\nX_test.targets = X_test.targets[idx]\n\ntest_loader = torch.utils.data.DataLoader(X_test, batch_size=1, shuffle=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Define Neural Network with Q-node\n\nThis NN is  2 layers of ConvNN and a fully connected layer, with a Q-Node as a classifier.","metadata":{}},{"cell_type":"code","source":"class Net(nn.Module):\n    def __init__(self):\n        super(Net, self).__init__()\n        self.conv1 = nn.Conv2d(1, 10, kernel_size=5)\n        self.conv2 = nn.Conv2d(10, 20, kernel_size=5)\n        self.conv2_drop = nn.Dropout2d()\n        self.fc1 = nn.Linear(320, 50)\n        self.fc2 = nn.Linear(50, NUM_QUBITS)\n        self.qc = TorchCircuit.apply\n        self.qcsim = nn.Linear(NUM_QUBITS, 1)\n        self.fc3 = nn.Linear(1, 2)\n\n    def forward(self, x):\n        x = F.relu(F.max_pool2d(self.conv1(x), 2))\n        x = F.relu(F.max_pool2d(self.conv2_drop(self.conv2(x)), 2))\n        x = x.view(-1, 320)\n        x = F.relu(self.fc1(x))\n        x = F.dropout(x, training=self.training)\n        x = self.fc2(x)\n        x = np.pi*torch.tanh(x)\n        \n#         print('params to QC: {}'.format(x))\n\n        MODE = 'QC' # 'QC' or 'QC_sim'\n    \n        if MODE == 'QC': \n            x = qc(x[0]) # QUANTUM LAYER\n        \n        else:\n            x = self.qcsim(x)\n            \n#         print('output of QC = {}'.format(x))\n        \n#         # softmax rather than sigmoid\n#         x = self.fc3(x.float())\n#         print('output of Linear(1, 2): {}'.format(x))\n#         x = F.softmax(x, 1)\n\n        x = torch.sigmoid(x)\n        x = torch.cat((x, 1-x), -1)\n#         print(x)\n        return x\n    \n    \n    def predict(self, x):\n        # apply softmax\n        pred = self.forward(x)\n#         print(pred)\n        ans = torch.argmax(pred[0]).item()\n        return torch.tensor(ans)\n    \nnetwork = Net()#.to(device)\noptimizer = optim.Adam(network.parameters(), lr=0.001)\n\n# optimizer = optim.Adam(network.parameters(), lr=learning_rate)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"epochs = 20\nloss_list = []\nloss_func = nn.CrossEntropyLoss()\n\nfor epoch in range(epochs):\n    total_loss = []\n    for batch_idx, (data, target) in enumerate(train_loader):\n#         print(batch_idx)\n        optimizer.zero_grad()        \n        # Forward pass\n        output = network(data)\n        # Calculating loss\n        loss = loss_func(output, target)\n        # Backward pass\n        loss.backward()\n        # Optimize the weights\n        optimizer.step()\n        \n        total_loss.append(loss.item())\n        \n    loss_list.append(sum(total_loss)/len(total_loss))\n    print('Training [{:.0f}%]\\tLoss: {:.4f}'.format(\n        100. * (epoch + 1) / epochs, loss_list[-1]))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(loss_list)\nplt.title('Hybrid NN Training Convergence for {}-qubit'.format(NUM_QUBITS))\nplt.xlabel('Training Iterations')\nplt.ylabel('Cross Entropy Loss')\nplt.savefig('MNIST01-bell/Figures/Hybrid {}-qubit Loss Curve bell.jpg'.format(NUM_QUBITS))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Test accuracy of NN\n\nThe outcome is not always the same because the prediction is probabilistic.","metadata":{}},{"cell_type":"code","source":"accuracy = 0\nnumber = 0\nfor batch_idx, (data, target) in enumerate(test_loader):\n    number +=1\n    output = network.predict(data).item()\n    accuracy += (output == target[0].item())*1","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Performance on test data is is: {}/{} = {}%\".format(accuracy,number,100*accuracy/number))    ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"n_samples_shape = (8, 6)\ncount = 0\nfig, axes = plt.subplots(nrows=n_samples_shape[0], ncols=n_samples_shape[1], figsize=(10, 2*n_samples_shape[0]))\n\nnetwork.eval()\nwith torch.no_grad():\n    for batch_idx, (data, target) in enumerate(test_loader):\n        if count == n_samples_shape[0]*n_samples_shape[1]:\n            break\n        pred = network.predict(data).item()\n\n        axes[count//n_samples_shape[1]][count%n_samples_shape[1]].imshow(data[0].numpy().squeeze(), cmap='gray')\n\n        axes[count//n_samples_shape[1]][count%n_samples_shape[1]].set_xticks([])\n        axes[count//n_samples_shape[1]][count%n_samples_shape[1]].set_yticks([])\n        axes[count//n_samples_shape[1]][count%n_samples_shape[1]].set_title('Predicted {}'.format(pred))\n        \n        count += 1","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}