quantum blur demo

main.py

import pew
import quantumblur as qb
from microqiskit import QuantumCircuit
from microqiskit import pi

pew.init()
screen = pew.Pix()

fps = 6

def draw_cursor(x,y,undraw=False):
    for dx,dy in [(-1,0),(1,0),(0,-1),(0,1)]:
        if x+dx in range(8) and y+dy in range(8):
            if undraw:
                screen.pixel(x+dx,y+dy,image[x+dx,y+dy])
            else:
                screen.pixel(x+dx,y+dy,3 )
    pew.show(screen)

def draw_image(image):
    for x in range(8):
        for y in range(8):
            screen.pixel(x,y,image[x,y])
    pew.show(screen)

def height2image(height):
    image = {}
    for x,y in height:
        h = height[x,y]
        if h<=1/3:
            image[x,y] = 0
        elif h<=2/3:
            image[x,y] = 1
        else:
            image[x,y] = 2
    return image

image = {}
for x in range(8):
    for y in range(8):
        image[x,y] = 1
image[0,0] = 0
image[7,7] = 0
image[0,7] = 2
image[7,0] = 2

X,Y = 4,4
idle_frames = 0
draw_image(image)
n = 0
theta = pi/100
while True:

    (dX,dY) = (0,0)
    keys = pew.keys()
    if keys==pew.K_UP:
        dY = -1
    elif keys==pew.K_DOWN:
        dY = +1
    if keys==pew.K_LEFT:
        dX = -1
    elif keys==pew.K_RIGHT:
        dX = +1
    if X+dX in range(8):
        X += dX
    if Y+dY in range(8):
        Y += dY

    if idle_frames<5*fps:
        
        idle = (X+dX,Y+dY)==(X,Y)

        if keys==pew.K_O or keys==pew.K_X:
            if keys==pew.K_O:
                image[X,Y] = (image[X,Y]+1)%3
            elif keys==pew.K_X:
                image[X,Y] = (image[X,Y]-1)%3
            screen.pixel(X,Y,image[X,Y])
            idle = False

        draw_cursor(X,Y)
        pew.tick(1/fps)
        draw_cursor(X,Y,undraw=True)

        if idle:
            idle_frames += 1
        else:
            idle_frames = 0

    else:
        if n==0:
            qc = qb.height2circuit(image)
            draw_image(image)
        else:
            qc_rot = QuantumCircuit(6)
            for q in range(6):
                qc_rot.rx(n*theta,q)
            draw_image(height2image(qb.circuit2height(qc+qc_rot)))
        n += 1
        #pew.tick(1/16)

microqiskit

import random
from math import cos,sin,pi
r2=0.70710678118 
class QuantumCircuit:
  def __init__(self,n,m=0):
    self.num_qubits=n
    self.num_clbits=m
    self.name = ''
    self.data=[]
  def __add__(self,self2):
    self3=QuantumCircuit(max(self.num_qubits,self2.num_qubits),max(self.num_clbits,self2.num_clbits))
    self3.data=self.data+self2.data
    self3.name = self.name
    return self3
  def initialize(self,k):
    self.data[:] = [] 
    self.data.append(('init',[e for e in k])) 
  def x(self,q):
    self.data.append(('x',q))
  def rx(self,theta,q):
    self.data.append(('rx',theta,q))
  def rz(self,theta,q):
    self.data.append(('rz',theta,q))
  def h(self,q):
    self.data.append(('h',q))
  def cx(self,s,t):
    self.data.append(('cx',s,t))
  def crx(self,theta,s,t):
    self.data.append(('crx',theta,s,t))
  def measure(self,q,b):
    assert b<self.num_clbits, 'Index for output bit out of range.'
    assert q<self.num_qubits, 'Index for qubit out of range.'
    self.data.append(('m',q,b))
  def ry(self,theta,q):
    self.rx(pi/2,q)
    self.rz(theta,q)
    self.rx(-pi/2,q)
  def z(self,q):
    self.rz(pi,q)
  def y(self,q):
    self.rz(pi,q)
    self.x(q)
def simulate(qc,shots=1024,get='counts',noise_model=[]):
  def superpose(x,y):
    return [r2*(x[j]+y[j])for j in range(2)],[r2*(x[j]-y[j])for j in range(2)]
  def turn(x,y,theta):
    theta = float(theta)
    return [x[0]*cos(theta/2)+y[1]*sin(theta/2),x[1]*cos(theta/2)-y[0]*sin(theta/2)],[y[0]*cos(theta/2)+x[1]*sin(theta/2),y[1]*cos(theta/2)-x[0]*sin(theta/2)]
  def phaseturn(x,y,theta):
     theta = float(theta)
     return [[x[0]*cos(theta/2) - x[1]*sin(-theta/2),x[1]*cos(theta/2) + x[0]*sin(-theta/2)],[y[0]*cos(theta/2) - y[1]*sin(+theta/2),y[1]*cos(theta/2) + y[0]*sin(+theta/2)]]   
  k = [[0,0] for _ in range(2**qc.num_qubits)] 
  k[0] = [1.0,0.0] 
  if noise_model:
    if type(noise_model)==float:
       noise_model = [noise_model]*qc.num_qubits
  outputnum_clbitsap = {}
  for gate in qc.data:
    if gate[0]=='init': 
      if type(gate[1][0])==list:
        k = [e for e in gate[1]]
      else: 
        k = [[e,0] for e in gate[1]]
    elif gate[0]=='m': 
      outputnum_clbitsap[gate[2]] = gate[1]
    elif gate[0] in ['x','h','rx','rz']: 
      j = gate[-1] 
      for i0 in range(2**j):
        for i1 in range(2**(qc.num_qubits-j-1)):
          b0=i0+2**(j+1)*i1 
          b1=b0+2**j 
          if gate[0]=='x': 
            k[b0],k[b1]=k[b1],k[b0]
          elif gate[0]=='h': 
            k[b0],k[b1]=superpose(k[b0],k[b1])
          elif gate[0]=='rx': 
            theta = gate[1]
            k[b0],k[b1]=turn(k[b0],k[b1],theta)
          elif gate[0]=='rz': 
            theta = gate[1]
            k[b0],k[b1]=phaseturn(k[b0],k[b1],theta) 
    elif gate[0] in ['cx','crx']: 
      if gate[0]=='cx': 
        [s,t] = gate[1:]
      else:
        theta = gate[1]
        [s,t] = gate[2:]
      [l,h] = sorted([s,t])
      for i0 in range(2**l):
        for i1 in range(2**(h-l-1)):
          for i2 in range(2**(qc.num_qubits-h-1)):
            b0=i0+2**(l+1)*i1+2**(h+1)*i2+2**s 
            b1=b0+2**t  
            if gate[0]=='cx':
                k[b0],k[b1]=k[b1],k[b0] 
            else:
                k[b0],k[b1]=turn(k[b0],k[b1],theta) 
  if get=='statevector':
    return k
  else:
    probs = [e[0]**2+e[1]**2 for e in k]
    if noise_model:
      for j in range(qc.num_qubits):
        p_meas = noise_model[j]
        for i0 in range(2**j):
          for i1 in range(2**(qc.num_qubits-j-1)):
            b0=i0+2**(j+1)*i1 
            b1=b0+2**j 
            p0 = probs[b0]
            p1 = probs[b1]
            probs[b0] = (1-p_meas)*p0 + p_meas*p1
            probs[b1] = (1-p_meas)*p1 + p_meas*p0
    if get=='probabilities_dict':
      return {('{0:0'+str(qc.num_qubits)+'b}').format(j):p for j,p in enumerate(probs)}
    elif get in ['counts', 'memory']:
      m = [False for _ in range(qc.num_qubits)]
      for gate in qc.data:
        for j in range(qc.num_qubits):
          assert  not ((gate[-1]==j) and m[j]), 'Incorrect or missing measure command.'
          m[j] = (gate==('m',j,j))
      m=[]
      for _ in range(shots):
        cumu=0
        un=True
        r=random.random()
        for j,p in enumerate(probs):
          cumu += p
          if r<cumu and un:    
            raw_out=('{0:0'+str(qc.num_qubits)+'b}').format(j)
            out_list = ['0']*qc.num_clbits
            for bit in outputnum_clbitsap:
              out_list[qc.num_clbits-1-bit] = raw_out[qc.num_qubits-1-outputnum_clbitsap[bit]]
            out = ''.join(out_list)
            m.append(out)
            un=False
      if get=='memory':
        return m
      else:
        counts = {}
        for out in m:
          if out in counts:
            counts[out] += 1
          else:
            counts[out] = 1
        return counts

quantumblur.py

# -*- coding: utf-8 -*-

# This code is part of Qiskit.
#
# (C) Copyright IBM 2020s.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

import math
import random

from microqiskit import QuantumCircuit, simulate
simple_python = True

def _kron(vec0,vec1):
    new_vec = []
    for amp0 in vec0:
        for amp1 in vec1:
            new_vec.append(amp0*amp1)
    return new_vec


def _get_size(height):
    Lx = 0
    Ly = 0
    for (x,y) in height:
        Lx = max(x+1,Lx)
        Ly = max(y+1,Ly)
    return Lx,Ly


def circuit2probs(qc):

    probs = simulate(qc,get='probabilities_dict')
    
    return probs


def make_line ( length ):

    # number of bits required
    n = int(math.ceil(math.log(length)/math.log(2)))
    
    # iteratively build list
    line = ['0','1']
    for j in range(n-1):
        # first append a reverse-ordered version of the current list
        line = line + line[::-1]
        # then add a '0' onto the end of all bit strings in the first half
        for j in range(int(float(len(line))/2)):
            line[j] += '0'
        # and a '1' for the second half
        for j in range(int(float(len(line))/2),int(len(line))):
            line[j] += '1'
            
    return line


def normalize(ket):

    N = 0
    for amp in ket:
        N += amp*amp.conjugate()
    for j,amp in enumerate(ket):
        ket[j] = float(amp)/math.sqrt(N)
    return ket


def make_grid(Lx,Ly=None):

    # set Ly if not supplied
    if not Ly:
        Ly = Lx
    
    # make the lines
    line_x = make_line( Lx )
    line_y = make_line( Ly )
    
    # make the grid
    grid = {}
    for x in range(Lx):
        for y in range(Ly):
            grid[ line_x[x]+line_y[y] ] = (x,y)
            
    # determine length of the bit strings
    n = len(line_x[0]+line_y[0])
            
    return grid, n


def height2circuit(height, log=False, eps=1e-2):

    # get bit strings for the grid
    Lx,Ly = _get_size(height)
    grid, n = make_grid(Lx,Ly)
    
    # create required state vector
    state = [0]*(2**n)
    if log:
        # normalize heights
        max_h = max(height.values())
        height = {pos:float(height[pos])/max_h for pos in height}
        # find minimum (not too small) normalized height
        min_h = min([height[pos] for pos in height if height[pos] > eps])
        # this minimum value defines the base
        base = 1.0/min_h
    for bitstring in grid:
        (x,y) = grid[bitstring]
        if (x,y) in height:
            h = height[x,y]
            if log:
                state[ int(bitstring,2) ] = math.sqrt(base**(float(h)/min_h))
            else:
                state[ int(bitstring,2) ] = math.sqrt( h )
    state = normalize(state)
        
    # define and initialize quantum circuit            
    qc = QuantumCircuit(n)
    # microqiskit style
    qc.initialize(state)
    qc.name = '('+str(Lx)+','+str(Ly)+')'

    return qc


def probs2height(probs, size=None, log=False):

    
    # get grid info
    if size:
        (Lx,Ly) = size
    else:
        Lx = int(2**(len(list(probs.keys())[0])/2))
        Ly = Lx
    grid,_ = make_grid(Lx,Ly)
    
    # set height to probs value, rescaled such that the maximum is 1
    max_h = max( probs.values() )   
    height = {(x,y):0.0 for x in range(Lx) for y in range(Ly)}
    for bitstring in probs:
        if bitstring in grid:
            height[grid[bitstring]] = float(probs[bitstring])/max_h
         
    # take logs if required
    if log:
        min_h = min([height[pos] for pos in height if height[pos] !=0])
        alt_min_h = min([height[pos] for pos in height])
        base = 1/min_h
        for pos in height:
            if height[pos]>0:
                height[pos] = max(math.log(height[pos]/min_h)/math.log(base),0)
            else:
                height[pos] = 0.0
                        
    return height

    
def circuit2height(qc, log=False):
    
    probs = circuit2probs(qc)
    try:
        # get size from circuit
        size = eval(qc.name)
    except:
        # if not in circuit name, infer it from qubit number
        L = int(2**(qc.num_qubits/2))
        size = (L,L)
    return probs2height(probs, size=size, log=log)


def combine_circuits(qc0,qc1):

    warning = "Combined circuits should contain only initialization."

    # create a circuit with the combined number of qubits
    num_qubits = qc0.num_qubits + qc1.num_qubits
    combined_qc = QuantumCircuit(num_qubits)

    # extract statevectors for any initialization commands
    kets = [None,None]
    for j,qc in enumerate([qc0, qc1]):
        for gate in qc.data:
            assert gate[0]=='init', warning
            kets[j] = gate[1]

    # combine into a statevector for all the qubits
    ket = None
    if kets[0] and kets[1]:
        ket = _kron(kets[0], kets[1])
    elif kets[0]:
        ket = _kron(kets[0], [1]+[0]*(2**qc1.num_qubits-1))
    elif kets[1]:
        ket = _kron([1]+[0]*(2**qc0.num_qubits-1),kets[1])

    # use this to initialize
    if ket:
        combined_qc.initialize(ket)
    
    # prevent circuit name from being used for size determination
    combined_qc.name = 'None'
            
    return combined_qc