Skip to content

Commit

Permalink
Merge pull request ScQ-Cloud#173 from Zhaoyilunnn/master
Browse files Browse the repository at this point in the history
fix a requirement issue and format all files
  • Loading branch information
ScQ-Cloud authored Jun 17, 2024
2 parents 6ad74c3 + edd0083 commit d7d3642
Show file tree
Hide file tree
Showing 38 changed files with 1,346 additions and 1,061 deletions.
2 changes: 1 addition & 1 deletion doc/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@
"current_version": current_version,
"versions": [
["0.2.x", pages_root + "/0.2.x"],
["0.3.x", pages_root + "/0.3.x"]
["0.3.x", pages_root + "/0.3.x"],
],
}

Expand Down
4 changes: 2 additions & 2 deletions quafu/algorithms/hamiltonian.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,9 +13,9 @@
# limitations under the License.

from typing import Iterable
import scipy.sparse as sp
import numpy as np

import numpy as np
import scipy.sparse as sp
from quafu.exceptions.quafu_error import QuafuError

IMat = sp.coo_matrix(np.array([[1.0, 0.0], [0.0, 1.0]], dtype=complex))
Expand Down
109 changes: 78 additions & 31 deletions quafu/algorithms/optimizer.py
Original file line number Diff line number Diff line change
@@ -1,28 +1,52 @@
import numpy as np

def adam(func, x, gradsf, args=(), call_back=None, eta=0.01,beta1=0.9, beta2=0.999, epsilon=1e-8, maxiter=1000, tol=1e-8, verbose=False):

def adam(
func,
x,
gradsf,
args=(),
call_back=None,
eta=0.01,
beta1=0.9,
beta2=0.999,
epsilon=1e-8,
maxiter=1000,
tol=1e-8,
verbose=False,
):
n_para = len(x)
mt = np.zeros(n_para)
vt = np.zeros(n_para)

traj = []
if verbose:
print((" "*5).join(["step".ljust(10), "loss".ljust(10), "grad_norm".ljust(10)]))

grads_norm = 0.
print(
(" " * 5).join(["step".ljust(10), "loss".ljust(10), "grad_norm".ljust(10)])
)

grads_norm = 0.0
for j in range(maxiter):
traj.append(func(x, *args))
if verbose:
print((" "*5).join([("%d" %j).ljust(10), ("%.5f" %traj[j]).ljust(10), ("%.5f" %grads_norm).ljust(10)]))
if j > 0 and abs(traj[j]-traj[j-1]) <= tol:
print(
(" " * 5).join(
[
("%d" % j).ljust(10),
("%.5f" % traj[j]).ljust(10),
("%.5f" % grads_norm).ljust(10),
]
)
)
if j > 0 and abs(traj[j] - traj[j - 1]) <= tol:
return x, traj[-1], traj

grads = np.array(gradsf(x, *args))
grads_norm = np.linalg.norm(grads)
mt = beta1 * mt + (1.-beta1) * grads
vt = beta1 * vt + (1.-beta2) * grads**2
mtt = mt / (1-beta1**(j+2))
vtt = vt / (1-beta2**(j+2))
mt = beta1 * mt + (1.0 - beta1) * grads
vt = beta1 * vt + (1.0 - beta2) * grads**2
mtt = mt / (1 - beta1 ** (j + 2))
vtt = vt / (1 - beta2 ** (j + 2))
x = x - eta * mtt / (np.sqrt(vtt) + epsilon)
if call_back:
call_back(x, *args)
Expand All @@ -31,42 +55,65 @@ def adam(func, x, gradsf, args=(), call_back=None, eta=0.01,beta1=0.9, beta2=0.9


def spsa_grad(func, x, k, args=(), spsa_iter=10, c=0.1, gamma=0.101):

dim = len(x)
ck = c/(k)**gamma
ck = c / (k) ** gamma
gx = 0.0
for i in range(spsa_iter):
Delta = 2*np.round(np.random.rand(dim))-1
x1 = x + ck*Delta
x2 = x - ck*Delta
y1 = func(x1, *args)
Delta = 2 * np.round(np.random.rand(dim)) - 1
x1 = x + ck * Delta
x2 = x - ck * Delta
y1 = func(x1, *args)
y2 = func(x2, *args)
gx += (y1 - y2) / (2*ck*Delta)
gx += (y1 - y2) / (2 * ck * Delta)
gx = gx / spsa_iter
return gx

def spsa(func, x, args=(), call_back=None, spsa_iter=10, max_iter=1000, a=0.1, c=0.1, A=100, alpha=0.602, gamma=0.101, tol=1e-8, verbose=False):
'''SPSA minimize
c: at a level of standard deviation of func
A: <=10% of max_iter '''

def spsa(
func,
x,
args=(),
call_back=None,
spsa_iter=10,
max_iter=1000,
a=0.1,
c=0.1,
A=100,
alpha=0.602,
gamma=0.101,
tol=1e-8,
verbose=False,
):
"""SPSA minimize
c: at a level of standard deviation of func
A: <=10% of max_iter"""
traj = [func(x, *args)]
if verbose:
print((" "*5).join(["step".ljust(10), "loss".ljust(10), "grad_norm".ljust(10)]))

grads_norm = 0.
print(
(" " * 5).join(["step".ljust(10), "loss".ljust(10), "grad_norm".ljust(10)])
)

grads_norm = 0.0
for k in range(max_iter):
if verbose:
print((" "*5).join([("%d" %k).ljust(10), ("%.5f" %traj[k]).ljust(10), ("%.5f" %grads_norm).ljust(10)]))
if k > 0 and abs(traj[k]-traj[k-1]) <= tol:
print(
(" " * 5).join(
[
("%d" % k).ljust(10),
("%.5f" % traj[k]).ljust(10),
("%.5f" % grads_norm).ljust(10),
]
)
)
if k > 0 and abs(traj[k] - traj[k - 1]) <= tol:
return x, traj[-1], traj
ak = a/(k+1+A)**alpha
grads = spsa_grad(func, x, k+1, args, spsa_iter, c=c, gamma=gamma)

ak = a / (k + 1 + A) ** alpha
grads = spsa_grad(func, x, k + 1, args, spsa_iter, c=c, gamma=gamma)
grads_norm = np.linalg.norm(grads)
x = x - ak * grads
if call_back:
call_back(x, *args)
traj.append(func(x, *args))

return x, traj[-1], traj

63 changes: 37 additions & 26 deletions quafu/algorithms/templates/amplitude.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,11 +13,10 @@
# limitations under the License.

"""Amplitude Embedding by a decomposition into gates"""
import quafu.elements.element_gates as qeg
import numpy as np
import quafu.elements.element_gates as qeg
from quafu.elements import QuantumGate


# from .basic_entangle import BasicEntangleLayers


Expand All @@ -41,7 +40,7 @@ def __iter__(self):

def __getitem__(self, index):
return self.gate_list[index]

def __add__(self, gates):
"""Addition operator."""
out = []
Expand All @@ -51,7 +50,7 @@ def __add__(self, gates):
else:
raise TypeError("Contains unsupported gate")
return out

def __radd__(self, other):
out = []
if all(isinstance(gate, QuantumGate) for gate in other):
Expand All @@ -75,10 +74,12 @@ def _preprocess(self, state, num_qubits, pad_with, normalize):

# check shape
if len(shape) != 1:
raise ValueError(f"state must be a one-dimensional tensor; got shape {shape}.")
raise ValueError(
f"state must be a one-dimensional tensor; got shape {shape}."
)

n_state = shape[0]
dim = 2 ** num_qubits
dim = 2**num_qubits
if pad_with is None and n_state != dim:
raise ValueError(
f"The length of state should be {dim}; got length {n_state}.Please check num_qubits "
Expand Down Expand Up @@ -111,10 +112,13 @@ def _preprocess(self, state, num_qubits, pad_with, normalize):
)
new_state_batch.append(feature_set)

return np.stack(new_state_batch).astype(np.complex128) if batched else new_state_batch[0].astype(np.complex128)

def _build(self):
return (
np.stack(new_state_batch).astype(np.complex128)
if batched
else new_state_batch[0].astype(np.complex128)
)

def _build(self):
a = np.abs(self.state)
omega = np.angle(self.state)
# change order of qubits, since original code was written for IBM machines
Expand All @@ -126,7 +130,9 @@ def _build(self):
alpha_y_k = _get_alpha_y(a, len(qubits_reverse), k)
control = qubits_reverse[k:]
target = qubits_reverse[k - 1]
gate_list.extend(_apply_uniform_rotation_dagger(qeg.RYGate, alpha_y_k, control, target))
gate_list.extend(
_apply_uniform_rotation_dagger(qeg.RYGate, alpha_y_k, control, target)
)

# If necessary, apply inverse z rotation cascade to prepare correct phases of amplitudes
if not np.allclose(omega, 0):
Expand All @@ -136,11 +142,14 @@ def _build(self):
target = qubits_reverse[k - 1]
if len(alpha_z_k) > 0:
gate_list.extend(
_apply_uniform_rotation_dagger(qeg.RZGate, alpha_z_k, control, target)
_apply_uniform_rotation_dagger(
qeg.RZGate, alpha_z_k, control, target
)
)

return gate_list



## MottonenStatePreparation related functions.
def gray_code(rank):
"""Generates the Gray code of given rank.
Expand All @@ -167,18 +176,19 @@ def gray_code_recurse(g, rank):

return g


def _matrix_M_entry(row, col):

# (col >> 1) ^ col is the Gray code of col
b_and_g = row & ((col >> 1) ^ col)
sum_of_ones = 0
while b_and_g > 0:
if b_and_g & 0b1:
sum_of_ones += 1
# (col >> 1) ^ col is the Gray code of col
b_and_g = row & ((col >> 1) ^ col)
sum_of_ones = 0
while b_and_g > 0:
if b_and_g & 0b1:
sum_of_ones += 1

b_and_g = b_and_g >> 1

b_and_g = b_and_g >> 1
return (-1) ** sum_of_ones

return (-1) ** sum_of_ones

def compute_theta(alpha):
ln = alpha.shape[-1]
Expand All @@ -195,7 +205,6 @@ def compute_theta(alpha):


def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):

gate_list = []
theta = compute_theta(alpha)

Expand All @@ -220,8 +229,8 @@ def _apply_uniform_rotation_dagger(gate, alpha, control_wires, target_wire):
gate_list.append(qeg.CXGate(control_wires[control_index], target_wire))
return gate_list

def _get_alpha_z(omega, n, k):

def _get_alpha_z(omega, n, k):
indices1 = [
[(2 * j - 1) * 2 ** (k - 1) + l - 1 for l in range(1, 2 ** (k - 1) + 1)]
for j in range(1, 2 ** (n - k) + 1)
Expand All @@ -237,16 +246,18 @@ def _get_alpha_z(omega, n, k):

return np.sum(diff, axis=-1)

def _get_alpha_y(a, n, k):

def _get_alpha_y(a, n, k):
indices_numerator = [
[(2 * (j + 1) - 1) * 2 ** (k - 1) + l for l in range(2 ** (k - 1))]
for j in range(2 ** (n - k))
]
numerator = np.take(a, indices=indices_numerator, axis=-1)
numerator = np.sum(np.abs(numerator) ** 2, axis=-1)

indices_denominator = [[j * 2**k + l for l in range(2**k)] for j in range(2 ** (n - k))]
indices_denominator = [
[j * 2**k + l for l in range(2**k)] for j in range(2 ** (n - k))
]
denominator = np.take(a, indices=indices_denominator, axis=-1)
denominator = np.sum(np.abs(denominator) ** 2, axis=-1)

Expand All @@ -262,4 +273,4 @@ def _get_alpha_y(a, n, k):

division = np.where(denominator != 0.0, division, 0.0)

return 2 * np.arcsin(np.sqrt(division))
return 2 * np.arcsin(np.sqrt(division))
6 changes: 2 additions & 4 deletions quafu/elements/element_gates/clifford.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,3 @@
from .element_gates import HGate, SGate, CXGate

__all__ = ['HGate', 'SGate', 'CXGate']

from .element_gates import CXGate, HGate, SGate

__all__ = ["HGate", "SGate", "CXGate"]
Loading

0 comments on commit d7d3642

Please sign in to comment.