From 8d88a666f4974359d97a81889fca28708c72cfa6 Mon Sep 17 00:00:00 2001 From: Sam Bonkowsky Date: Sun, 4 Aug 2024 12:51:57 -0700 Subject: [PATCH 1/4] Undo automatic type conversion in coupling_op --- SQcircuit/circuit.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index 5934aef..e34951f 100755 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -2226,8 +2226,6 @@ def coupling_op( op = self._squeeze_op(op) - if get_optim_mode(): - return sqf.qobj_to_tensor(op) return op def matrix_elements( From f6251e93baf73783746d73ee50467b6ca8454d6c Mon Sep 17 00:00:00 2001 From: Alex McKeehan Boulton Date: Tue, 6 Aug 2024 10:53:21 -0700 Subject: [PATCH 2/4] Make minor updates to circuit.py to ensure even trunc product <= K and charge/harmonic modes are balanced --- SQcircuit/circuit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index c6c0592..6d1e211 100755 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -1912,16 +1912,16 @@ def truncate_circuit(self, K: int) -> List[int]: List of truncation numbers for each mode of the circuit """ trunc_num_average = K ** (1 / len(self.omega)) - charge_cutoff = (1 / 2) * (trunc_num_average - 1) + charge_cutoff = (1 / 2) * (trunc_num_average + 1) trunc_nums = [] for mode_idx in range(len(self.omega)): # Harmonic mode if self.omega[mode_idx] != 0: - trunc_nums.append(int(np.ceil(trunc_num_average))) + trunc_nums.append(int(np.floor(trunc_num_average))) # Charge mode else: - trunc_nums.append(int(np.ceil(charge_cutoff))) + trunc_nums.append(int(np.floor(charge_cutoff))) self.set_trunc_nums(trunc_nums) return trunc_nums From ffdb161a379c0d66eb59f5578e06b88e78561ef8 Mon Sep 17 00:00:00 2001 From: taha1373 Date: Tue, 6 Aug 2024 14:03:50 -0700 Subject: [PATCH 3/4] tests for the bug reported in GitHub --- SQcircuit/tests/test_qubits.py | 81 ++++++++++++++++++++++------------ 1 file changed, 54 insertions(+), 27 deletions(-) diff --git a/SQcircuit/tests/test_qubits.py b/SQcircuit/tests/test_qubits.py index 5cc55d7..2d45c56 100755 --- a/SQcircuit/tests/test_qubits.py +++ b/SQcircuit/tests/test_qubits.py @@ -247,43 +247,70 @@ def test_coupled_fluxonium_transmon(): def test_loop_issue(): - loop1 = sq.Loop(id_str="loop1") - loop2 = sq.Loop(id_str="loop2") - ec = 1 ecj = 3 ej = 5 el = 0.5 circuit_dict = {} - # circuit_dict[(3, 0)] = [sq.Capacitor(ec, "GHz", id_str="C01")] - # circuit_dict[(3, 1)] = [sq.Capacitor(ec, "GHz", id_str="C02")] - # circuit_dict[(3, 2)] = [sq.Capacitor(ec, "GHz", id_str="C03")] - # circuit_dict[(0, 1)] = [sq.Junction(ej, "GHz", id_str="J12", loops=[loop1], - # cap=sq.Capacitor(ecj, "GHz", - # id_str="JC12"))] - # circuit_dict[(0, 2)] = [sq.Junction(ej, "GHz", id_str="J13", loops=[loop1], - # cap=sq.Capacitor(ecj, "GHz", - # id_str="JC13"))] - # circuit_dict[(1, 2)] = [ - # sq.Junction(ej, "GHz", id_str="J23", loops=[loop1, loop2], - # cap=sq.Capacitor(ecj, "GHz", id_str="JC23")), - # sq.Inductor(el, "GHz", id_str="J23", loops=[loop2])] + + loop1 = sq.Loop(id_str="loop1") + loop2 = sq.Loop(0.5, id_str="loop2") + + circuit_dict[(3, 0)] = [sq.Capacitor(ec, "GHz", id_str="C01")] + circuit_dict[(3, 1)] = [sq.Capacitor(ec, "GHz", id_str="C02")] + circuit_dict[(3, 2)] = [sq.Capacitor(ec, "GHz", id_str="C03")] + circuit_dict[(0, 1)] = [sq.Junction( + ej, "GHz", id_str="J12", loops=[loop1], + cap=sq.Capacitor(ecj, "GHz", id_str="JC12") + )] + circuit_dict[(0, 2)] = [sq.Junction( + ej, "GHz", id_str="J13", loops=[loop1], + cap=sq.Capacitor(ecj, "GHz", id_str="JC13") + )] + circuit_dict[(1, 2)] = [ + sq.Junction( + ej, "GHz", id_str="J23", loops=[loop1, loop2], + cap=sq.Capacitor(ecj, "GHz", id_str="JC23") + ), + sq.Inductor(el, "GHz", id_str="J23", loops=[loop2]) + ] + + circ_1 = sq.Circuit(circuit_dict) + circ_1.set_trunc_nums([20, 15, 1]) + efreqs, _ = circ_1.diag(2) + q_freq_1 = efreqs[1] - efreqs[0] + + circuit_dict = {} + + loop1 = sq.Loop(id_str="loop1") + loop2 = sq.Loop(0.5, id_str="loop2") circuit_dict[(0, 1)] = [sq.Capacitor(ec, "GHz", id_str="C01")] circuit_dict[(0, 2)] = [sq.Capacitor(ec, "GHz", id_str="C02")] circuit_dict[(0, 3)] = [sq.Capacitor(ec, "GHz", id_str="C03")] - circuit_dict[(1, 2)] = [sq.Junction(ej, "GHz", id_str="J12", loops=[loop1], - cap=sq.Capacitor(ecj, "GHz", - id_str="JC12"))] - circuit_dict[(1, 3)] = [sq.Junction(ej, "GHz", id_str="J13", loops=[loop1], - cap=sq.Capacitor(ecj, "GHz", - id_str="JC13"))] + circuit_dict[(1, 2)] = [sq.Junction( + ej, "GHz", id_str="J12", loops=[loop1], + cap=sq.Capacitor(ecj, "GHz", id_str="JC12") + )] + circuit_dict[(1, 3)] = [sq.Junction( + ej, "GHz", id_str="J13", loops=[loop1], + cap=sq.Capacitor(ecj, "GHz", id_str="JC13") + )] circuit_dict[(2, 3)] = [ - sq.Junction(ej, "GHz", id_str="J23", loops=[loop1, loop2], - cap=sq.Capacitor(ecj, "GHz", id_str="JC23")), - sq.Inductor(el, "GHz", id_str="J23", loops=[loop2])] + sq.Junction( + ej, "GHz", id_str="J23", loops=[loop1, loop2], + cap=sq.Capacitor(ecj, "GHz", id_str="JC23") + ), + sq.Inductor(el, "GHz", id_str="J23", loops=[loop2]) + ] + + circ_2 = sq.Circuit(circuit_dict) + + circ_2.set_trunc_nums([20, 15, 1]) + efreqs, _ = circ_2.diag(2) + q_freq_2 = efreqs[1] - efreqs[0] - circ = sq.Circuit(circuit_dict) + assert np.isclose(q_freq_1, 0.010243745) + assert np.isclose(q_freq_1, q_freq_2) - circ.description() From 858351e1c3de6d40f03a842b16697add36f9393f Mon Sep 17 00:00:00 2001 From: taha1373 Date: Tue, 6 Aug 2024 14:04:30 -0700 Subject: [PATCH 4/4] change loop construction in SQcircuit by adding remove_dependent_columns --- SQcircuit/circuit.py | 400 ++++++++++++++++++---------------- SQcircuit/elements.py | 30 ++- SQcircuit/functions.py | 39 +++- SQcircuit/torch_extensions.py | 48 ++-- 4 files changed, 281 insertions(+), 236 deletions(-) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index c5ebd13..2a52f5a 100755 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -93,62 +93,62 @@ def _set_matrix_rep(self) -> ndarray: This helps to construct the capacitance and susceptance matrices. """ - A = np.zeros((self.circ.n, self.circ.n)) + edge_mat = np.zeros((self.circ.n, self.circ.n)) if 0 in self.edge: i = self.edge[0] + self.edge[1] - 1 - A[i, i] = 1 + edge_mat[i, i] = 1 else: i = self.edge[0] - 1 j = self.edge[1] - 1 - A[i, i] = 1 - A[j, j] = 1 - A[i, j] = -1 - A[j, i] = -1 + edge_mat[i, i] = 1 + edge_mat[j, j] = 1 + edge_mat[i, j] = -1 + edge_mat[j, i] = -1 - return A + return edge_mat def update_circuit_loop_from_element( self, el: Union[Inductor, Junction], - B_idx: int, + b_id: int, ) -> None: """Update loop properties related to circuit from element in the edge - with its inductive index (B_idx). + with its inductive index (b_id). Parameters ---------- el: Inductive element. - B_idx: + b_id: Inductive index. """ for loop in el.loops: self.circ.add_loop(loop) - loop.add_index(B_idx) - loop.addK1(self.w) + loop.add_index(b_id) + loop.add_to_k_mat(self.w) if get_optim_mode(): self.circ.add_to_parameters(loop) def process_edge_and_update_circ( self, - B_idx: int, - W_idx: int, - K1: list, + b_id: int, + w_id: int, + k_mat: list, c_edge_mat: list, ) -> Tuple[int, list, list]: """Process the edge and update the related circuit properties. Parameters ---------- - B_idx: + b_id: Point to each row of B matrix of the circuit. - W_idx: + w_id: Point to each row of W matrix of the circuit. - K1: + k_mat: Matrix related to loop calculation c_edge_mat: edge capacitance matrix @@ -171,14 +171,14 @@ def process_edge_and_update_circ( ) self.circ.elem_keys[el.type].append( - el.get_key(self.edge, B_idx, W_idx) + el.get_key(self.edge, b_id, w_id) ) - self.update_circuit_loop_from_element(el, B_idx) + self.update_circuit_loop_from_element(el, b_id) - B_idx += 1 + b_id += 1 - K1.append(self.w) + k_mat.append(self.w) c_edge_mat.append( el.get_cap_for_flux_dist(self.circ.flux_dist) @@ -195,7 +195,7 @@ def process_edge_and_update_circ( self.processed = True - return B_idx, K1, c_edge_mat + return b_id, k_mat, c_edge_mat def _check_if_edge_is_processed(self) -> None: @@ -294,16 +294,16 @@ def __init__( # number of nodes without ground self.n: int = max(max(self.elements)) - # number of branches that contain JJ without parallel inductor. - self.countJJnoInd: int = 0 + # number of branches that contain junctions without a parallel inductor. + self.num_jun_without_ind: int = 0 self.elem_keys = { - # inductor element keys: (edge, el, B_idx) B_idx point to + # inductor element keys: (edge, el, b_id) b_id point to # each row of B matrix (external flux distribution of that element) Inductor: [], - # junction element keys: (edge, el, B_idx, W_idx) B_idx point to + # junction element keys: (edge, el, b_id, w_id) b_id point to # each row of B matrix (external flux distribution of that element) - # and W_idx point to each row of W matrix + # and w_id point to each row of W matrix Junction: [], } @@ -567,7 +567,6 @@ def add_to_parameters(self, el: Element) -> None: if (el.requires_grad or not el.is_leaf) and el not in self._parameters: self._parameters[el] = el.internal_value - def add_loop(self, loop: Loop) -> None: """Add loop to the circuit loops. @@ -653,12 +652,12 @@ def safecopy(self, save_eigs=False) -> Self: Inductor: [], Junction: [], } - for edge, el, B_idx, W_idx in self.elem_keys[Junction]: + for edge, el, b_id, w_id in self.elem_keys[Junction]: new_el = replacement_dict[el] - new_circuit.elem_keys[Junction].append((edge, new_el, B_idx, W_idx)) - for edge, el, B_idx in self.elem_keys[Inductor]: + new_circuit.elem_keys[Junction].append((edge, new_el, b_id, w_id)) + for edge, el, b_id in self.elem_keys[Inductor]: new_el = replacement_dict[el] - new_circuit.elem_keys[Inductor].append((edge, new_el, B_idx)) + new_circuit.elem_keys[Inductor].append((edge, new_el, b_id)) # Update the `.partial_mats` dictionary new_circuit.partial_mats = defaultdict(lambda: 0) @@ -678,9 +677,9 @@ def safecopy(self, save_eigs=False) -> Self: new_circuit._memory_ops[op_type] = self._memory_ops[op_type] else: new_circuit._memory_ops[op_type] = {} - for el, B_idx in self._memory_ops[op_type].keys(): - new_circuit._memory_ops[op_type][(replacement_dict[el], B_idx)] = ( - self._memory_ops[op_type][(el, B_idx)] + for el, b_id in self._memory_ops[op_type].keys(): + new_circuit._memory_ops[op_type][(replacement_dict[el], b_id)] = ( + self._memory_ops[op_type][(el, b_id)] ) # Remove old eigenvectors/values, if desired @@ -707,7 +706,7 @@ def picklecopy(self) -> Self: ---------- Copy of self with ``._toggle_fullcopy = False``. """ - # Instantiate new container + # Instantiate new container new_circuit = copy(self) # Remove large objects when saving @@ -724,93 +723,92 @@ def _get_LCWB(self) -> Tuple[ndarray, ndarray, ndarray, ndarray]: A tuple of the (capacitance, susceptance, W, B) matrices. """ - cMat = sqf.zeros((self.n, self.n), dtype=torch.float64) - lMat = sqf.zeros((self.n, self.n), dtype=torch.float64) - wMat = [] - bMat = sqf.array([]) + c_mat = sqf.zeros((self.n, self.n), dtype=torch.float64) + l_mat = sqf.zeros((self.n, self.n), dtype=torch.float64) + w_mat = [] + b_mat = sqf.array([]) self.partial_mats = defaultdict(lambda: 0) # point to each row of B matrix (external flux distribution of that # element) or count the number of inductive elements. - B_idx = 0 + b_id = 0 - # W_idx point to each row of W matrix for junctions or count the + # w_id point to each row of W matrix for junctions or count the # number of edges contain JJ - W_idx = 0 + w_id = 0 - # number of branches that contain JJ without parallel inductor. - countJJnoInd = 0 + # number of branches that contain JJ without a parallel inductor. + num_jun_without_ind = 0 - # K1 is a matrix that transfers node coordinates to the edge phase drop - # for inductive elements - K1 = [] + # k_mat is a matrix that transfers node coordinates to the edge phase + # drop for inductive elements (In the paper Appendix D we referred to it + # as W matrix) + k_mat = [] - # capacitor at each inductive elements + # capacitor at each inductive element c_edge_mat = [] for edge in self.elements.keys(): circ_edge = CircuitEdge(self, edge) - B_idx, K1, c_edge_mat = circ_edge.process_edge_and_update_circ( - B_idx, - W_idx, - K1, + b_id, k_mat, c_edge_mat = circ_edge.process_edge_and_update_circ( + b_id, + w_id, + k_mat, c_edge_mat ) if circ_edge.is_JJ_without_ind(): - countJJnoInd += 1 + num_jun_without_ind += 1 - cMat += sqf.array(circ_edge.get_eff_cap_value()) * sqf.array( + c_mat += sqf.array(circ_edge.get_eff_cap_value()) * sqf.array( circ_edge.mat_rep) - lMat += sqf.array(circ_edge.get_eff_ind_value()) * sqf.array( + l_mat += sqf.array(circ_edge.get_eff_ind_value()) * sqf.array( circ_edge.mat_rep) if circ_edge.is_JJ_in_this_edge(): - wMat.append(circ_edge.w) - W_idx += 1 + w_mat.append(circ_edge.w) + w_id += 1 - wMat = np.array(wMat) + w_mat = np.array(w_mat) try: - K1 = np.array(K1) - a = np.zeros_like(K1) - select = np.sum(K1 != a, axis=0) != 0 - # eliminate the zero columns - K1 = K1[:, select] - if K1.shape[0] == K1.shape[1]: - K1 = K1[:, 0:-1] - - X = K1.T @ np.diag(sqf.numpy(c_edge_mat)) + k_mat = np.array(k_mat) + + k_mat = sqf.remove_dependent_columns(k_mat) + + x_mat = k_mat.T @ np.diag(sqf.numpy(c_edge_mat)) for loop in self.loops: - p = np.zeros((1, B_idx)) - p[0, loop.indices] = loop.getP() - X = np.concatenate((X, p), axis=0) + g = np.zeros((1, b_id)) + g[0, loop.indices] = loop.get_g() + x_mat = np.concatenate((x_mat, g), axis=0) # number of inductive loops of the circuit n_loops = len(self.loops) if n_loops != 0: - Y = np.concatenate( - (np.zeros((B_idx - n_loops, n_loops)), np.eye(n_loops)), + i_mat = np.concatenate( + (np.zeros((b_id - n_loops, n_loops)), np.eye(n_loops)), axis=0 ) - bMat = np.linalg.inv(X) @ Y - bMat = np.around(bMat, 5) + b_mat = np.linalg.inv(x_mat) @ i_mat + b_mat = np.around(b_mat, 5) except ValueError: - print("The edge list does not specify a connected graph or " - "all inductive loops of the circuit are not specified.") + print( + "The edge list does not specify a connected graph or " + "all inductive loops of the circuit are not specified." + ) raise ValueError - self.countJJnoInd = countJJnoInd + self.num_jun_without_ind = num_jun_without_ind - return cMat, lMat, wMat, bMat + return c_mat, l_mat, w_mat, b_mat def _is_charge_mode(self, i: int) -> bool: """Check if the mode is a charge mode. @@ -856,33 +854,35 @@ def _get_and_apply_transformation_1(self) -> Tuple[ndarray, ndarray]: A tuple of the (S1, R1) matrices. """ - cMatRoot = sqrtm(sqf.numpy(self.C)) - cMatRootInv = np.linalg.inv(cMatRoot) - lMatRoot = sqrtm(sqf.numpy(self.L)) + c_mat_root = sqrtm(sqf.numpy(self.C)) + c_mat_root_inv = np.linalg.inv(c_mat_root) + l_mat_root = sqrtm(sqf.numpy(self.L)) - V, D, U = np.linalg.svd(lMatRoot @ cMatRootInv) + _, D, U = np.linalg.svd(l_mat_root @ c_mat_root_inv) # the case that there is not any inductor in the circuit if np.max(D) == 0: D = np.diag(np.eye(self.n)) - singLoc = list(range(0, self.n)) + sing_locs = list(range(0, self.n)) else: - # find the number of singularity in the circuit + # find the number of singularities in the circuit - lEig, _ = np.linalg.eig(sqf.numpy(self.L)) - numSing = len(lEig[lEig / np.max(lEig) < ACC["sing_mode_detect"]]) - singLoc = list(range(self.n - numSing, self.n)) - D[singLoc] = np.max(D) + l_mat_eigs, _ = np.linalg.eig(sqf.numpy(self.L)) + num_sings = len(l_mat_eigs[ + l_mat_eigs / np.max(l_mat_eigs) < ACC["sing_mode_detect"] + ]) + sing_locs = list(range(self.n - num_sings, self.n)) + D[sing_locs] = np.max(D) # build S1 and R1 matrix - S1 = cMatRootInv @ U.T @ np.diag(np.sqrt(D)) - R1 = np.linalg.inv(S1).T + s_1 = c_mat_root_inv @ U.T @ np.diag(np.sqrt(D)) + r_1 = np.linalg.inv(s_1).T - self._apply_transformation(S1, R1) + self._apply_transformation(s_1, r_1) - self.lTrans[singLoc, singLoc] = 0 + self.lTrans[sing_locs, sing_locs] = 0 - return S1, R1 + return s_1, r_1 @staticmethod def _independent_rows(A: ndarray) -> Tuple[List[int], List[ndarray]]: @@ -937,7 +937,7 @@ def _round_to_zero_one(self, W: ndarray) -> ndarray: rounded_W = W.copy() - if self.countJJnoInd == 0: + if self.num_jun_without_ind == 0: rounded_W[:, self.omega == 0] = 0 charge_only_W = rounded_W[:, self.omega == 0] @@ -950,7 +950,7 @@ def _round_to_zero_one(self, W: ndarray) -> ndarray: return rounded_W - def _is_Gram_Schmidt_successful(self, S) -> bool: + def _is_gram_schmidt_successful(self, S) -> bool: """Check if the Gram_Schmidt process has the sufficient accuracy for the ``S`` transformation matrix. @@ -967,13 +967,11 @@ def _is_Gram_Schmidt_successful(self, S) -> bool: is_successful = True # absolute value of the current wTrans - cur_wTrans = self.wTrans @ S - - cur_wTrans = self._round_to_zero_one(cur_wTrans) + cur_w_trans = self._round_to_zero_one(self.wTrans @ S) for j in range(self.n): if self._is_charge_mode(j): - for abs_w in np.abs(cur_wTrans[:, j]): + for abs_w in np.abs(cur_w_trans[:, j]): if abs_w != 0 and abs_w != 1: is_successful = False @@ -1000,51 +998,65 @@ def _get_and_apply_transformation_2(self) -> Tuple[ndarray, ndarray]: """ if len(self.W) != 0: + # remove independent rows + w_t = sqf.remove_dependent_columns(self.W.T) + w_reduced = w_t.T + w_reduced_transformed = w_reduced @ self.S1 + + # charge part of the reduced w_mat + w_charge = w_reduced_transformed[:, self.omega == 0].copy() + # get the charge basis part of the wTrans matrix - wQ = self.wTrans[:, self.omega == 0].copy() + # w_charge = self.wTrans[:, self.omega == 0].copy() + # number of operators represented in charge bases - nq = wQ.shape[1] + nq = w_charge.shape[1] else: nq = 0 - wQ = np.array([]) + w_charge = np.array([]) - # if we need to represent an operator in charge basis - if nq != 0 and self.countJJnoInd != 0: + # if we need to represent an operator in the charge basis + if nq != 0 and self.num_jun_without_ind != 0: # list of indices of w vectors that are independent - indList = [] + ind_lst = [] + + # random vectors to complete the basis + rnd_vecs = [] - X = [] # use Gram–Schmidt to find the linear independent rows of - # normalized wQ (wQ_norm) + # normalized w_charge (w_charge_norm) basis = [] while len(basis) != nq: if len(basis) == 0: - indList, basis = self._independent_rows(wQ) + ind_lst, basis = self._independent_rows(w_charge) else: - # to complete the basis - X = list(np.random.randn(nq - len(basis), nq)) - basisComplete = np.array(basis + X) - _, basis = self._independent_rows(basisComplete) + # random vectors to complete the basis + rnd_vecs = list( + np.linalg.norm(w_charge, 'fro') + * np.random.randn(nq - len(basis), nq) + ) + completed_basis = np.array(basis + rnd_vecs) + _, basis = self._independent_rows(completed_basis) # the second S and R matrix are: - F = np.array(list(wQ[indList, :]) + X) - S2 = block_diag(np.eye(self.n - nq), np.linalg.inv(F)) + f_mat = np.array(list(w_charge[ind_lst, :]) + rnd_vecs) + s_2 = block_diag(np.eye(self.n - nq), np.linalg.inv(f_mat)) - R2 = np.linalg.inv(S2.T) + r_2 = np.linalg.inv(s_2.T) else: - S2 = np.eye(self.n, self.n) - R2 = S2 + s_2 = np.eye(self.n, self.n) + r_2 = s_2 - if self._is_Gram_Schmidt_successful(S2): + if self._is_gram_schmidt_successful(s_2): - self._apply_transformation(S2, R2) + self._apply_transformation(s_2, r_2) self.wTrans = self._round_to_zero_one(self.wTrans) - return S2, R2 + return s_2, r_2 else: print("Gram_Schmidt process failed. Retrying...") @@ -1060,7 +1072,7 @@ def _get_and_apply_transformation_3(self) -> Tuple[ndarray, ndarray]: A tuple of the (S3, R3) transformation matrices. """ - S3 = np.eye(self.n) + s_3 = np.eye(self.n) for j in range(self.n): @@ -1081,22 +1093,22 @@ def _get_and_apply_transformation_3(self) -> Tuple[ndarray, ndarray]: # find the coefficient in wTrans for j-th mode that # has maximum alpha s = np.abs(self.wTrans[np.argmax(jth_alphas), j]) - S3[j, j] = 1 / s + s_3[j, j] = 1 / s else: # scale the uncoupled mode s = np.max(np.abs(self.S[:, j])) - S3[j, j] = 1 / s + s_3[j, j] = 1 / s else: # scale the uncoupled mode s = np.max(np.abs(self.S[:, j])) - S3[j, j] = 1 / s + s_3[j, j] = 1 / s - R3 = np.linalg.inv(S3.T) + r_3 = np.linalg.inv(s_3.T) - self._apply_transformation(S3, R3) + self._apply_transformation(s_3, r_3) - return S3, R3 + return s_3, r_3 def _transform_hamil(self) -> None: """Transform the Hamiltonian of the circuit into the charge and Fock @@ -1154,10 +1166,10 @@ def _compute_params(self) -> None: (2 * np.pi * unt.get_unit_freq())) self.descrip_vars['EL'] = [] self.descrip_vars['EL_incl'] = [] - for _, el, B_idx in self.elem_keys[Inductor]: + for _, el, b_id in self.elem_keys[Inductor]: self.descrip_vars['EL'].append(sqf.numpy(el.get_value(unt._unit_ind))) self.descrip_vars['EL_incl'].append( - np.sum(np.abs(self.descrip_vars['B'][B_idx, :])) != 0) + np.sum(np.abs(self.descrip_vars['B'][b_id, :])) != 0) self.descrip_vars['EC'] = dict() for i in range(self.descrip_vars['har_dim'], self.descrip_vars['n_modes']): @@ -1250,8 +1262,10 @@ def set_trunc_nums(self, nums: List[int]) -> None: print("set_trunc_nums called") error1 = "The input must be be a python list" assert isinstance(nums, list), error1 - error2 = ("The number of modes (length of the input) must be equal to " - "the number of nodes") + error2 = ( + "The number of modes (length of the input) must be equal to " + "the number of nodes" + ) assert len(nums) == self.n, error2 self.m = self.n*[1] @@ -1593,7 +1607,7 @@ def _get_LC_hamil(self) -> Qobj: def _get_external_flux_at_element( self, - B_idx: int, + b_id: int, torch = False ) -> Union[float, Tensor]: """ @@ -1601,7 +1615,7 @@ def _get_external_flux_at_element( Parameters ---------- - B_idx: + b_id: An integer point to each row of B matrix (external flux distribution of that element) torch: @@ -1609,11 +1623,11 @@ def _get_external_flux_at_element( Returns ---------- - The external flux at the element referenced by ``B_idx.`` + The external flux at the element referenced by ``b_id.`` """ phi_ext = sqf.zero() for i, loop in enumerate(self.loops): - phi_ext += loop.value() * self.B[B_idx, i] + phi_ext += loop.value() * self.B[b_id, i] if isinstance(phi_ext, Tensor) and not torch: return sqf.numpy(phi_ext) @@ -1629,10 +1643,10 @@ def _get_inductive_hamil(self) -> Qobj: """ H = qt.Qobj() - for edge, el, B_idx in self.elem_keys[Inductor]: + for edge, el, b_id in self.elem_keys[Inductor]: # phi = 0 - # if B_idx is not None: - phi = self._get_external_flux_at_element(B_idx) + # if b_id is not None: + phi = self._get_external_flux_at_element(b_id) # summation of the 1 over inductor values. x = np.squeeze(sqf.numpy(1 / el.get_value())) @@ -1644,24 +1658,24 @@ def _get_inductive_hamil(self) -> Qobj: H += x * phi * (unt.Phi0 / 2 / np.pi) * op / np.sqrt(unt.hbar) # save the operators for loss calculation - self._memory_ops["ind_hamil"][(el, B_idx)] = op - for _, el, B_idx, W_idx in self.elem_keys[Junction]: + self._memory_ops["ind_hamil"][(el, b_id)] = op + for _, el, b_id, w_id in self.elem_keys[Junction]: # phi = 0 - # if B_idx is not None: - phi = self._get_external_flux_at_element(B_idx) + # if b_id is not None: + phi = self._get_external_flux_at_element(b_id) EJ = np.squeeze(sqf.numpy(el.get_value())) - exp = np.exp(1j * phi) * self._memory_ops["exp"][W_idx] + exp = np.exp(1j * phi) * self._memory_ops["exp"][w_id] root_exp = np.exp(1j * phi / 2) * self._memory_ops["root_exp"][ - W_idx] + w_id] cos = (exp + exp.dag()) / 2 sin = (exp - exp.dag()) / 2j sin_half = (root_exp - root_exp.dag()) / 2j - self._memory_ops["cos"][el, B_idx] = self._squeeze_op(cos) - self._memory_ops["sin"][el, B_idx] = self._squeeze_op(sin) - self._memory_ops["sin_half"][el, B_idx] = self._squeeze_op(sin_half) + self._memory_ops["cos"][el, b_id] = self._squeeze_op(cos) + self._memory_ops["sin"][el, b_id] = self._squeeze_op(sin) + self._memory_ops["sin_half"][el, b_id] = self._squeeze_op(sin_half) H += -EJ * cos @@ -1714,7 +1728,7 @@ def charge_op(self, mode: int, basis: str = 'FC') -> Qobj: return qt.Qobj(Q_eig) - def _get_W_idx(self, el: Junction, B_idx: int) -> Optional[int]: + def _get_w_id(self, el: Junction, b_id: int) -> Optional[int]: """" Find the corresponding ``W`` matrix index given an junction and its ``B`` matrix index. @@ -1723,16 +1737,16 @@ def _get_W_idx(self, el: Junction, B_idx: int) -> Optional[int]: ---------- el: Josephson junction in circuit. - B_idx: + b_id: Index of B matrix corresponding to ``el``. Returns ---------- The corresponding ``W`` matrix index, if it exists. """ - for _, o_el, o_B_idx, W_idx in self.elem_keys[Junction]: - if o_el == el and o_B_idx == B_idx: - return W_idx + for _, o_el, o_b_id, w_id in self.elem_keys[Junction]: + if o_el == el and o_b_id == b_id: + return w_id return None @@ -1756,22 +1770,22 @@ def op(self, typ: str, keywords: Dict) -> Union[Qobj, Tensor]: The `typ` operator of the circuit specified by ``keywords``. """ if typ == 'sin_half': - B_idx = keywords['B_idx'] + b_id = keywords['b_id'] el = keywords['el'] if get_optim_mode(): - W_idx = self._get_W_idx(el, B_idx) + w_id = self._get_w_id(el, b_id) - phi = self._get_external_flux_at_element(B_idx, torch=True) + phi = self._get_external_flux_at_element(b_id, torch=True) root_exp = ( torch.exp(1j * phi / 2) - * sqf.qobj_to_tensor(self._memory_ops["root_exp"][W_idx]) + * sqf.qobj_to_tensor(self._memory_ops["root_exp"][w_id]) ) sin_half = (root_exp - sqf.dag(root_exp)) / 2j # ToDo: need to squeeze? return sin_half else: - return self._memory_ops['sin_half'][el, B_idx] + return self._memory_ops['sin_half'][el, b_id] else: raise NotImplementedError @@ -2024,7 +2038,7 @@ def _tensor_to_modes(self, tensorIndex: int) -> List[int]: Returns ---------- - indList: + ind_lst: A list of mode indices (self.n) """ @@ -2037,16 +2051,16 @@ def _tensor_to_modes(self, tensorIndex: int) -> List[int]: else: mP = [mP[0] * self.m[-1 - i]] + mP - indList = [] + ind_lst = [] indexP = tensorIndex for i in range(self.n): if i == self.n - 1: - indList.append(indexP) + ind_lst.append(indexP) continue - indList.append(int(indexP / mP[i])) + ind_lst.append(int(indexP / mP[i])) indexP = indexP % mP[i] - return indList + return ind_lst def eig_phase_coord(self, k: int, grid: Sequence[ndarray]) -> ndarray: """ @@ -2082,7 +2096,7 @@ def eig_phase_coord(self, k: int, grid: Sequence[ndarray]) -> ndarray: # decomposes the tensor product space index (i) to each mode # indices as a list - indList = self._tensor_to_modes(i) + ind_lst = self._tensor_to_modes(i) if get_optim_mode(): term = self._evecs[k][i].item() @@ -2092,7 +2106,7 @@ def eig_phase_coord(self, k: int, grid: Sequence[ndarray]) -> ndarray: for mode in range(self.n): # mode number related to that node - n = indList[mode] + n = ind_lst[mode] # For charge basis if self._is_charge_mode(mode): @@ -2334,8 +2348,10 @@ def _dec_rate_cc_np(self, states: Tuple[int, int]) -> float: Critical current dephasing rate between ``states``. """ decay = 0 - for el, B_idx in self._memory_ops['cos']: - partial_omega = self._get_partial_omega_mn(el, states=states, _B_idx=B_idx) + for el, b_id in self._memory_ops['cos']: + partial_omega = self._get_partial_omega_mn( + el, states=states, _b_id=b_id + ) A = el.A * el.get_value() decay = decay + self._dephasing(A, partial_omega) @@ -2426,8 +2442,8 @@ def dec_rate( sqf.operator_inner_product(state1, op, state2)) ** 2 if dec_type == "quasiparticle": - for el, B_idx in self._memory_ops['sin_half']: - op = self.op('sin_half', {'el': el, 'B_idx': B_idx}) + for el, b_id in self._memory_ops['sin_half']: + op = self.op('sin_half', {'el': el, 'b_id': b_id}) if np.isnan(sqf.numpy(el.Y(omega, ENV["T"]))): decay = decay + 0 else: @@ -2513,7 +2529,7 @@ def _get_quadratic_phi(self, A: ndarray) -> Qobj: def _get_partial_H( self, el: Union[Capacitor, Inductor, Junction, Loop], - _B_idx: Optional[int] = None, + _b_id: Optional[int] = None, ) -> Qobj: """ Compute the gradient of the Hamiltonian with respect to elements or @@ -2523,7 +2539,7 @@ def _get_partial_H( el: Element of a circuit that can be either ``Capacitor``, ``Inductor``, ``Junction``, or ``Loop``. - _B_idx: + _b_id: Optional integer to indicate which row of the B matrix (per-element external flux distribution) to use. This specifies which JJ of the circuit to consider specifically (ex. for @@ -2546,12 +2562,12 @@ def _get_partial_H( A = -self.S.T @ self.partial_mats[el] @ self.S partial_H += self._get_quadratic_phi(A) - for edge, el_ind, B_idx in self.elem_keys[Inductor]: + for edge, el_ind, b_id in self.elem_keys[Inductor]: if el == el_ind: - phi = self._get_external_flux_at_element(B_idx) + phi = self._get_external_flux_at_element(b_id) - partial_H += -(self._memory_ops["ind_hamil"][(el, B_idx)] + partial_H += -(self._memory_ops["ind_hamil"][(el, b_id)] / np.squeeze(sqf.numpy(el.get_value()))**2 / np.sqrt(unt.hbar) * (unt.Phi0/2/np.pi) * phi) @@ -2559,27 +2575,27 @@ def _get_partial_H( loop_idx = self.loops.index(el) - for edge, el_ind, B_idx in self.elem_keys[Inductor]: + for edge, el_ind, b_id in self.elem_keys[Inductor]: partial_H += ( - self.B[B_idx, loop_idx] - * self._memory_ops["ind_hamil"][(el_ind, B_idx)] + self.B[b_id, loop_idx] + * self._memory_ops["ind_hamil"][(el_ind, b_id)] / sqf.numpy(el_ind.get_value()) * unt.Phi0 / np.sqrt(unt.hbar) / 2 / np.pi ) - for edge, el_JJ, B_idx, W_idx in self.elem_keys[Junction]: - partial_H += (self.B[B_idx, loop_idx] * sqf.numpy(el_JJ.get_value()) - * self._memory_ops['sin'][(el_JJ, B_idx)]) + for edge, el_JJ, b_id, w_id in self.elem_keys[Junction]: + partial_H += (self.B[b_id, loop_idx] * sqf.numpy(el_JJ.get_value()) + * self._memory_ops['sin'][(el_JJ, b_id)]) elif isinstance(el, Junction): - for _, el_JJ, B_idx, W_idx in self.elem_keys[Junction]: + for _, el_JJ, b_id, w_id in self.elem_keys[Junction]: - if el == el_JJ and _B_idx is None: - partial_H += -self._memory_ops['cos'][(el, B_idx)] + if el == el_JJ and _b_id is None: + partial_H += -self._memory_ops['cos'][(el, b_id)] - elif el == el_JJ and _B_idx == B_idx: - partial_H += -self._memory_ops['cos'][(el, B_idx)] + elif el == el_JJ and _b_id == b_id: + partial_H += -self._memory_ops['cos'][(el, b_id)] return partial_H @@ -2588,7 +2604,7 @@ def get_partial_omega( el: Union[Capacitor, Inductor, Junction, Loop], m: int, subtract_ground: bool = True, - _B_idx: Optional[int] = None, + _b_id: Optional[int] = None, ) -> float: """Return the gradient of the eigen angular frequency with respect to elements or loop as ``qutip.Qobj`` format. @@ -2604,7 +2620,7 @@ def get_partial_omega( subtract_ground: If ``True``, it subtracts the ground state frequency from the desired frequency. - _B_idx: + _b_id: Optional integer to indicate which row of the B matrix (per-element external flux distribution) to use. This specifies which JJ of the circuit to consider specifically (ex. for @@ -2617,7 +2633,7 @@ def get_partial_omega( """ state_m = self._evecs[m] - partial_H = self._get_partial_H(el, _B_idx) + partial_H = self._get_partial_H(el, _b_id) partial_omega_m = sqf.operator_inner_product( state_m, partial_H, state_m ) @@ -2636,7 +2652,7 @@ def _get_partial_omega_mn( self, el: Union[Capacitor, Inductor, Junction, Loop], states: Tuple[int, int], - _B_idx: Optional[int] = None, + _b_id: Optional[int] = None, ) -> float: """Return the gradient of the eigen angular frequency with respect to elements or loop as ``qutip.Qobj`` format. Note that if @@ -2650,7 +2666,7 @@ def _get_partial_omega_mn( states: Integers indicating indices of eigenenergies to calculate the difference of. - _B_idx: + _b_id: Optional integer to indicate which row of the B matrix (external flux distribution of that element) to use. This specifies which JJ of the circuit to consider specifically (ex. for critical @@ -2665,7 +2681,7 @@ def _get_partial_omega_mn( state_m = self._evecs[states[0]] state_n = self._evecs[states[1]] - partial_H = self._get_partial_H(el, _B_idx) + partial_H = self._get_partial_H(el, _b_id) partial_omega_m = sqf.operator_inner_product( state_m, partial_H, state_m diff --git a/SQcircuit/elements.py b/SQcircuit/elements.py index eac8812..de13dbc 100755 --- a/SQcircuit/elements.py +++ b/SQcircuit/elements.py @@ -630,8 +630,8 @@ def __init__( self.A = A * 2 * np.pi # indices of inductive elements. self.indices = [] - # k1 matrix related to this specific loop - self.K1 = [] + # k matrix related to this specific loop + self.k_mat = [] if id_str is None: self.id_str = "loop" @@ -658,7 +658,7 @@ def is_leaf(self) -> bool: return self.internal_value.is_leaf def reset(self) -> None: - self.K1 = [] + self.k_mat = [] self.indices = [] def value(self, random: bool = False) -> float: @@ -693,21 +693,17 @@ def set_flux(self, value: float) -> None: def add_index(self, index): self.indices.append(index) - def addK1(self, w): - self.K1.append(w) - - def getP(self): - K1 = np.array(self.K1) - a = np.zeros_like(K1) - select = np.sum(K1 != a, axis=0) != 0 - # eliminate the zero columns - K1 = K1[:, select] - if K1.shape[0] == K1.shape[1]: - K1 = K1[:, 0:-1] - b = np.zeros((1, K1.shape[0])) + def add_to_k_mat(self, w): + self.k_mat.append(w) + + def get_g(self): + + k_mat = sqf.remove_dependent_columns(np.array(self.k_mat)) + + b = np.zeros((1, k_mat.shape[0])) b[0, 0] = 1 - p = np.linalg.inv(np.concatenate((b, K1.T), axis=0)) @ b.T - return p.T + g = np.linalg.inv(np.concatenate((b, k_mat.T), axis=0)) @ b.T + return g.T @property def internal_value(self) -> Union[float, Tensor]: diff --git a/SQcircuit/functions.py b/SQcircuit/functions.py index 1c05b5f..567e899 100755 --- a/SQcircuit/functions.py +++ b/SQcircuit/functions.py @@ -13,6 +13,36 @@ from SQcircuit.settings import get_optim_mode +def remove_zero_columns(matrix): + """Removes columns from the matrix that contain only zero values.""" + + # Identify columns that are not all zeros + non_zero_columns = np.any(matrix != 0, axis=0) + + # Filter the matrix to keep only non-zero columns + filtered_matrix = matrix[:, non_zero_columns] + + return filtered_matrix + + +def remove_dependent_columns(matrix, tol=1e-5): + """Remove dependent columns from the given matrix.""" + + # first remove the zero columns from the matrix + no_zero_col_matrix = remove_zero_columns(matrix) + + # Perform SVD + _, s, _ = np.linalg.svd(no_zero_col_matrix) + + # Identify columns with singular values above tolerance + independent_columns = np.where(np.abs(s) > tol)[0] + + # Reconstruct matrix with only independent columns + reduced_matrix = no_zero_col_matrix[:, independent_columns] + + return reduced_matrix + + def block_diag(*args: Union[ndarray, Tensor]) -> Union[ndarray, Tensor]: if get_optim_mode(): return torch.block_diag(*args) @@ -264,8 +294,9 @@ def qobj_to_tensor(qobj, dtype=torch.complex128): indices_tensor = torch.LongTensor(indices) values_tensor = torch.tensor(values, dtype=dtype) - coo_tensor = torch.sparse_coo_tensor(indices_tensor, values_tensor, torch.Size(shape), - dtype=dtype) + coo_tensor = torch.sparse_coo_tensor( + indices_tensor, values_tensor, torch.Size(shape), dtype=dtype + ) return coo_tensor @@ -289,7 +320,9 @@ def qutip(A: Union[Qobj, Tensor], dims=List[list]) -> Qobj: col_indices = indices[1, :] values = input.values().numpy() shape = tuple(input.shape) - coo_sparse = scipy.sparse.coo_matrix((values, (row_indices, col_indices)), shape=shape) + coo_sparse = scipy.sparse.coo_matrix( + (values, (row_indices, col_indices)), shape=shape + ) csr_sparse = coo_sparse.tocsr() qobj = qt.Qobj(inpt=csr_sparse) qobj.dims = dims diff --git a/SQcircuit/torch_extensions.py b/SQcircuit/torch_extensions.py index 08ec6f8..b5a988e 100755 --- a/SQcircuit/torch_extensions.py +++ b/SQcircuit/torch_extensions.py @@ -378,7 +378,7 @@ def get_B_indices( el: Union[Junction, Inductor] ) -> List[int]: """ - Return the list of ``B_idx``'s with the element ``el`` (the same element + Return the list of ``b_id``'s with the element ``el`` (the same element could be placed at multiple branches). Parameters @@ -394,13 +394,13 @@ def get_B_indices( """ B_indices = [] if isinstance(el, Junction): - for _, el_JJ, B_idx, _ in cr.elem_keys[Junction]: + for _, el_JJ, b_id, _ in cr.elem_keys[Junction]: if el_JJ is el: - B_indices.append(B_idx) + B_indices.append(b_id) elif isinstance(el, Inductor): - for _, el_ind, B_idx in cr.elem_keys[Inductor]: + for _, el_ind, b_id in cr.elem_keys[Inductor]: if el_ind is el: - B_indices.append(B_idx) + B_indices.append(b_id) return B_indices @@ -648,7 +648,7 @@ def backward(ctx, grad_output) -> Tuple[Tensor, None, None]: def partial_squared_H_EJ( cr: 'Circuit', EJ_el: Junction, - B_idx: int, + b_id: int, grad_el: SupportedGradElements ) -> qt.Qobj: """ Calculates the second derivative of the Hamiltonian of ``cr`` with @@ -660,7 +660,7 @@ def partial_squared_H_EJ( The ``Circuit`` object to differentiate the Hamiltonian of. EJ_el: The Josephson junction to differentiate with respect to. - B_idx: + b_id: The index of the ``cr.B`` matrix identifying which branch the ``EJ_el`` is on. grad_el: @@ -675,13 +675,13 @@ def partial_squared_H_EJ( return 0 loop_idx = cr.loops.index(grad_el) - return cr.B[B_idx, loop_idx] * cr._memory_ops['sin'][(EJ_el, B_idx)] + return cr.B[b_id, loop_idx] * cr._memory_ops['sin'][(EJ_el, b_id)] def partial_squared_omega_mn_EJ( cr: 'Circuit', EJ_el: Junction, - B_idx: int, + b_id: int, grad_el: SupportedGradElements, states: Tuple[int, int] ) -> float: @@ -694,7 +694,7 @@ def partial_squared_omega_mn_EJ( The ``Circuit`` object to differentiate the eigenfrequencies of. EJ_el: A Josephson junction to differentiate with respect to. - B_idx: + b_id: A number grad_el: A circuit element to differentiate with respect to. @@ -706,8 +706,8 @@ def partial_squared_omega_mn_EJ( The second derivative of the eigenfrequency difference with respect to ``EJ_el`` and ``grad_el``. """ - partial_H = cr._get_partial_H(EJ_el, _B_idx = B_idx) - partial_H_squared = partial_squared_H_EJ(cr, EJ_el, B_idx, grad_el) + partial_H = cr._get_partial_H(EJ_el, _b_id = b_id) + partial_H_squared = partial_squared_H_EJ(cr, EJ_el, b_id, grad_el) return partial_squared_omega( cr, @@ -744,16 +744,16 @@ def partial_cc_dec( dec_rate_grad = 0 # Sum over all cosine operators because each Josephson junction is # associated with exactly one. - for EJ_el, B_idx in cr._memory_ops['cos']: + for EJ_el, b_id in cr._memory_ops['cos']: partial_omega_mn = sqf.numpy(cr._get_partial_omega_mn( EJ_el, states=states, - _B_idx=B_idx + _b_id=b_id )) partial_squared_omega_mn = partial_squared_omega_mn_EJ( cr, EJ_el, - B_idx, + b_id, grad_el, states ) @@ -846,27 +846,27 @@ def partial_squared_H_phi( H_squared = 0 if isinstance(grad_el, Junction): - for B_idx in B_indices: - H_squared += cr.B[B_idx, loop_idx] * cr._memory_ops['sin'][(grad_el, B_idx)] + for b_id in B_indices: + H_squared += cr.B[b_id, loop_idx] * cr._memory_ops['sin'][(grad_el, b_id)] return H_squared elif isinstance(grad_el, Inductor): - for B_idx in B_indices: + for b_id in B_indices: H_squared += ( - cr.B[B_idx, loop_idx] + cr.B[b_id, loop_idx] / -sqf.numpy(grad_el.get_value()**2) * unt.Phi0 / np.sqrt(unt.hbar) / 2 / np.pi - * cr._memory_ops["ind_hamil"][(grad_el, B_idx)] + * cr._memory_ops["ind_hamil"][(grad_el, b_id)] ) return H_squared elif isinstance(grad_el, Loop): loop_idx_1 = cr.loops.index(loop) loop_idx_2 = cr.loops.index(grad_el) - for edge, el_JJ, B_idx, W_idx in cr.elem_keys[Junction]: + for edge, el_JJ, b_id, W_idx in cr.elem_keys[Junction]: H_squared += ( sqf.numpy(el_JJ.get_value()) - * cr.B[B_idx, loop_idx_2] - * cr.B[B_idx, loop_idx_1] - * cr._memory_ops['cos'][(el_JJ, B_idx)] + * cr.B[b_id, loop_idx_2] + * cr.B[b_id, loop_idx_1] + * cr._memory_ops['cos'][(el_JJ, b_id)] ) return H_squared else: