From e3db7da27891d7b417e8b9b7042e0b5a06143423 Mon Sep 17 00:00:00 2001 From: sambonkov Date: Tue, 28 May 2024 17:17:07 -0700 Subject: [PATCH 1/7] Allow loop values to be tensors --- SQcircuit/circuit.py | 4 ++-- SQcircuit/elements.py | 23 ++++++++++++++++++++++- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index f5b1882..ddfe5e2 100644 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -1054,7 +1054,7 @@ def description( loopTxt = txt.loops() + txt.tab() for i in range(len(self.loops)): - phiExt = self.loops[i].value() / 2 / np.pi + phiExt = sqf.numpy(self.loops[i].value()) / 2 / np.pi loopTxt += txt.phiExt(i + 1) + txt.tPi() + txt.eq() + str( phiExt) + txt.tab() @@ -1480,7 +1480,7 @@ def _get_external_flux_at_element(self, B_idx: int) -> float: """ phi_ext = 0.0 for i, loop in enumerate(self.loops): - phi_ext += loop.value() * self.B[B_idx, i] + phi_ext += sqf.numpy(loop.value()) * self.B[B_idx, i] return phi_ext diff --git a/SQcircuit/elements.py b/SQcircuit/elements.py index 7622151..2d1ce4a 100644 --- a/SQcircuit/elements.py +++ b/SQcircuit/elements.py @@ -626,10 +626,15 @@ def __init__( self, value: float = 0, A: float = 1e-6, + requires_grad: bool = False, id_str: Optional[str] = None ) -> None: - self.lpValue = value * 2 * np.pi + self.set_flux(value) + + if requires_grad: + self.requires_grad = requires_grad + self.A = A * 2 * np.pi # indices of inductive elements. self.indices = [] @@ -669,6 +674,9 @@ def set_flux(self, value: float) -> None: value: The external flux value """ + if get_optim_mode(): + value = torch.as_tensor(value) + self.lpValue = value * 2 * np.pi def add_index(self, index): @@ -690,6 +698,19 @@ def getP(self): p = np.linalg.inv(np.concatenate((b, K1.T), axis=0)) @ b.T return p.T + @property + def requires_grad(self) -> bool: + raise_optim_error_if_needed() + + return self.lpValue.requires_grad + + @requires_grad.setter + def requires_grad(self, f: bool) -> None: + + raise_optim_error_if_needed() + + self.lpValue.requires_grad = f + class Charge: """Class that contains the charge island properties. From 7615858b77d317e9cdb376debc590eb131fc9b8f Mon Sep 17 00:00:00 2001 From: sambonkov Date: Sun, 2 Jun 2024 17:03:41 -0700 Subject: [PATCH 2/7] Support flux eigenfrequency gradient. --- SQcircuit/circuit.py | 24 ++++++++--- SQcircuit/elements.py | 8 ++++ SQcircuit/tests/conftest.py | 41 ++++++++++++++++++ SQcircuit/tests/test_grad.py | 79 ++++++++++++++++++++++++++++++++--- SQcircuit/torch_extensions.py | 9 ++-- 5 files changed, 146 insertions(+), 15 deletions(-) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index ddfe5e2..f3fa4ef 100644 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -136,6 +136,9 @@ def update_circuit_loop_from_element( loop.add_index(B_idx) loop.addK1(self.w) + if get_optim_mode(): + self.circ.add_to_parameters(loop) + def process_edge_and_update_circ( self, B_idx: int, @@ -412,18 +415,29 @@ def safecopy(self, save_eigs=False): # Instantiate new container new_circuit = copy(self) - # Explicitly copy any non-leaf tensors - # (these don't implement a __deepcopy__ method) + # Explicitly copy all elements to new circuit, in particular those with + # non-leaf `.internal_value`s ((these don't implement a `._deepcopy__` + # method) if get_optim_mode(): new_circuit.C = new_circuit.C.detach() new_circuit.L = new_circuit.L.detach() + new_loops: List[Loop] = [] + replacement_dict: Dict[Union[Loop, Element], Union[Loop, Element]] = {} + for loop in self.loops: + new_loop = copy(loop) + new_loop.internal_value = loop.internal_value.detach().clone() + new_loops.append(new_loop) + replacement_dict[loop] = new_loop + new_circuit.loops = new_loops + new_elements = defaultdict(list) - replacement_dict = dict() for edge in self.elements: for el in self.elements[edge]: new_el = copy(el) new_el.internal_value = el.internal_value.detach().clone() + if hasattr(el, 'loop'): + new_el.loop = replacement_dict[el.loop] new_elements[edge].append(new_el) replacement_dict[el] = new_el @@ -544,7 +558,7 @@ def add_to_parameters(self, el: Element) -> None: """Add elements with ``requires_grad=True`` to parameters. """ - if el.requires_grad: + if el.requires_grad and el not in self._parameters: self._parameters[el] = el.internal_value def add_loop(self, loop: Loop) -> None: @@ -2497,4 +2511,4 @@ def get_all_circuit_elements(self): def flatten(l): return [item for sublist in l for item in sublist] elements = flatten(list(self.elements.values())) - return elements \ No newline at end of file + return elements diff --git a/SQcircuit/elements.py b/SQcircuit/elements.py index 2d1ce4a..aa1ae60 100644 --- a/SQcircuit/elements.py +++ b/SQcircuit/elements.py @@ -711,6 +711,14 @@ def requires_grad(self, f: bool) -> None: self.lpValue.requires_grad = f + @property + def internal_value(self) -> Union[float, Tensor]: + return self.lpValue + + @internal_value.setter + def internal_value(self, v: Union[float, Tensor]) -> None: + self.lpValue = v + class Charge: """Class that contains the charge island properties. diff --git a/SQcircuit/tests/conftest.py b/SQcircuit/tests/conftest.py index 5b737eb..d2eb0d2 100644 --- a/SQcircuit/tests/conftest.py +++ b/SQcircuit/tests/conftest.py @@ -103,6 +103,36 @@ def create_transmon_torch(trunc_num): return circuit_torch +def create_flux_transmon_numpy(trunc_num, phi_ext): + cap_value, ind_value, Q = 7.746, 5, 1e6 + cap_unit, ind_unit = 'fF', 'GHz' + disorder = 1.1 + + set_optim_mode(False) + loop = Loop(phi_ext) + C_numpy = Capacitor(cap_value, cap_unit, Q=Q) + J1_numpy = Junction(ind_value, ind_unit, loops=[loop]) + J2_numpy = Junction(ind_value * disorder, ind_unit, loops=[loop]) + circuit_numpy = Circuit({(0, 1): [C_numpy, J1_numpy, J2_numpy], }) + circuit_numpy.set_trunc_nums([trunc_num, ]) + return circuit_numpy + + +def create_flux_transmon_torch(trunc_num, phi_ext): + cap_value, ind_value, Q = 7.746, 5, 1e6 + cap_unit, ind_unit = 'fF', 'GHz' + disorder = 1.1 + + set_optim_mode(True) + loop_torch = Loop(phi_ext, requires_grad=True) + C_torch = Capacitor(cap_value, cap_unit, Q=Q, requires_grad=True) + J1_torch = Junction(ind_value, ind_unit, loops=[loop_torch], requires_grad=True) + J2_torch = Junction(ind_value * disorder, ind_unit, loops=[loop_torch], requires_grad=True) + circuit_torch = Circuit({(0, 1): [C_torch, J1_torch, J2_torch], }) + circuit_torch.set_trunc_nums([trunc_num, ]) + return circuit_torch + + def create_fluxonium_numpy(trunc_num, phi_ext=0): set_optim_mode(False) loop = Loop(phi_ext) @@ -123,3 +153,14 @@ def create_fluxonium_torch(trunc_num, phi_ext=0): circuit_torch = Circuit({(0, 1): [C_torch, L_torch, JJ_torch], }, flux_dist='junctions') circuit_torch.set_trunc_nums([trunc_num, ]) return circuit_torch + + +def create_fluxonium_torch_flux(trunc_num, phi_ext=0): + set_optim_mode(True) + loop = Loop(phi_ext, requires_grad=True) + C_torch = Capacitor(3.6, 'GHz', Q=1e6, requires_grad=True) + L_torch = Inductor(0.46, 'GHz', Q=500e6, loops=[loop], requires_grad=True) + JJ_torch = Junction(10.2, 'GHz', A=1e-7, x=3e-06, loops=[loop], requires_grad=True) + circuit_torch = Circuit({(0, 1): [C_torch, L_torch, JJ_torch], }, flux_dist='junctions') + circuit_torch.set_trunc_nums([trunc_num, ]) + return circuit_torch diff --git a/SQcircuit/tests/test_grad.py b/SQcircuit/tests/test_grad.py index 880723b..8253413 100644 --- a/SQcircuit/tests/test_grad.py +++ b/SQcircuit/tests/test_grad.py @@ -16,7 +16,10 @@ create_fluxonium_numpy, create_fluxonium_torch, create_transmon_numpy, - create_transmon_torch + create_transmon_torch, + create_flux_transmon_numpy, + create_flux_transmon_torch, + create_fluxonium_torch_flux, ) @@ -30,9 +33,9 @@ def max_ratio(a, b): return np.max([np.abs(b / a), np.abs(a / b)]) -def function_grad_test(circuit_numpy, +def function_grad_test(circuit_numpy: Circuit, function_numpy, - circuit_torch, + circuit_torch: Circuit, function_torch, num_eigenvalues=20, delta=1e-4): @@ -60,7 +63,6 @@ def function_grad_test(circuit_numpy, set_optim_mode(True) circuit_torch.diag(num_eigenvalues) tensor_val = function_torch(circuit_torch) - optimizer = torch.optim.SGD(circuit_torch.parameters, lr=1) tensor_val.backward() for edge_idx, elements_by_edge in enumerate(circuit_numpy.elements.values()): @@ -101,7 +103,7 @@ def function_grad_test(circuit_numpy, print(f"edge element: {edge_element}") print(f"value: {edge_element._value}") print(f"value grad: {edge_element._value.grad}") - grad_torch = edge_elements_torch[element_idx]._value.grad.detach().numpy() + grad_torch = edge_elements_torch[element_idx].internal_value.grad.detach().numpy() # TODO: Modify Element class so that following gradient scaling is not necessary if isinstance(element_numpy, Capacitor) and element_numpy.value_unit in unt.freq_list: @@ -116,7 +118,38 @@ def function_grad_test(circuit_numpy, assert np.sign(grad_torch) == np.sign(grad_numpy) assert max_ratio(grad_torch, grad_numpy) <= 1 + tolerance - optimizer.zero_grad() + for loop_idx, loop in enumerate(circuit_numpy.loops): + set_optim_mode(False) + + loop_flux = loop.internal_value / 2 / np.pi + # Calculate f(x+delta) + loop.set_flux(loop_flux + delta) + circuit_numpy.diag(num_eigenvalues) + val_plus = function_numpy(circuit_numpy) + + # Calculate f(x-delta) + loop.set_flux(loop_flux - delta) + circuit_numpy.diag(num_eigenvalues) + val_minus = function_numpy(circuit_numpy) + + # Calculate gradient, being careful for flux insensitivity + if val_plus == val_minus: + grad_numpy = 0 + else: + grad_numpy = (val_plus - val_minus) / (2 * delta * 2 * np.pi) + + # Reset circuit + loop.set_flux(loop_flux) + + grad_torch = circuit_torch.loops[loop_idx].internal_value.grad + if grad_torch is not None: + print(f"grad torch: {grad_torch}, grad numpy: {grad_numpy}") + grad_torch = grad_torch.detach().numpy() + assert np.sign(grad_torch) == np.sign(grad_numpy) + assert max_ratio(grad_torch, grad_numpy) <= 1 + tolerance + + set_optim_mode(True) + circuit_torch.zero_parameters_grad() def first_eigendifference_numpy(circuit): return circuit._efreqs[1] - circuit._efreqs[0] @@ -141,6 +174,40 @@ def test_omega_transmon(): first_eigendifference_torch) set_optim_mode(False) +def test_omega_flux_transmon(): + """Verify gradient of first eigendifference omega_1-omega_0 + in transmon circuit with linearized value.""" + + # Create circuits + + for flux_point in [1+1e-2, 0.3, 0.5 - 1e-3, 0.5+1e-3, 0.7, 1-1e-2]: + print('flux point:', flux_point) + circuit_numpy = create_flux_transmon_numpy(trunc_num, flux_point) + circuit_torch = create_flux_transmon_torch(trunc_num, flux_point) + + function_grad_test(circuit_numpy, + first_eigendifference_numpy, + circuit_torch, + first_eigendifference_torch) + set_optim_mode(False) + +def test_omega_flux_fluxonium(): + """Verify gradient of first eigendifference omega_1-omega_0 + in transmon circuit with linearized value.""" + + # Create circuits + + for flux_point in [1+1e-2, 0.3, 0.5 - 1e-3, 0.5+1e-3, 0.7, 1-1e-2]: + print('flux point:', flux_point) + circuit_numpy = create_fluxonium_numpy(trunc_num, flux_point) + circuit_torch = create_fluxonium_torch_flux(trunc_num, flux_point) + + function_grad_test(circuit_numpy, + first_eigendifference_numpy, + circuit_torch, + first_eigendifference_torch) + set_optim_mode(False) + def test_T1_transmon(): """Compare gradient of T1 decoherence due to capacitive, inductive, and diff --git a/SQcircuit/torch_extensions.py b/SQcircuit/torch_extensions.py index d472979..7223bdd 100644 --- a/SQcircuit/torch_extensions.py +++ b/SQcircuit/torch_extensions.py @@ -66,7 +66,7 @@ def dec_rate_flux_torch(circuit: 'Circuit', states: Tuple[int, int]): class EigenSolver(Function): @staticmethod - def forward(ctx, + def forward(ctx, element_tensors: Tensor, circuit: 'Circuit', n_eig: int) -> Tensor: @@ -103,14 +103,14 @@ def forward(ctx, @staticmethod @once_differentiable def backward(ctx, grad_output) -> Tuple[Tensor]: - # Break grad_output into eigenvalue sub-tensor and eigenvector + # Break grad_output into eigenvalue sub-tensor and eigenvector # sub-tensor elements = list(ctx.circuit._parameters.keys()) m, n, l = ( len(ctx.circuit.parameters), # number of parameters ctx.n_eig, # number of eigenvalues (grad_output.shape[1] - 1), # length of eigenvectors - ) + ) grad_output_eigenvalue = grad_output[:, 0] grad_output_eigenvector = grad_output[:, 1:] @@ -121,7 +121,8 @@ def backward(ctx, grad_output) -> Tuple[Tensor]: # Compute backward pass for eigenvalues partial_omega[el_idx, eigen_idx] = ctx.circuit.get_partial_omega( el=elements[el_idx], - m=eigen_idx, subtract_ground=False + m=eigen_idx, + subtract_ground=False ) partial_tensor = torch.squeeze(torch.as_tensor( ctx.circuit.get_partial_vec( From 39e32b8f97560423e679a55e3aab66301412e292 Mon Sep 17 00:00:00 2001 From: sambonkov Date: Sun, 2 Jun 2024 18:03:59 -0700 Subject: [PATCH 3/7] Support elements being placed on multiple branches in T2 grad --- SQcircuit/torch_extensions.py | 36 ++++++++++++++++++++++------------- 1 file changed, 23 insertions(+), 13 deletions(-) diff --git a/SQcircuit/torch_extensions.py b/SQcircuit/torch_extensions.py index 7223bdd..84ed069 100644 --- a/SQcircuit/torch_extensions.py +++ b/SQcircuit/torch_extensions.py @@ -227,20 +227,25 @@ def partial_dephasing_rate( ) -def get_B_idx( +def get_B_indices( cr: 'Circuit', el: Union[Junction, Inductor] -): +) -> list[int]: + """ + Return the list of `B_idx`'s with the element `el` (the same element + could be placed at multiple branches) + """ + B_indices = [] if isinstance(el, Junction): for _, el_JJ, B_idx, _ in cr.elem_keys[Junction]: if el_JJ is el: - return B_idx + B_indices.append(B_idx) elif isinstance(el, Inductor): for _, el_ind, B_idx in cr.elem_keys[Inductor]: if el_ind is el: - return B_idx + B_indices.append(B_idx) - return None + return B_indices ############################################################################### @@ -584,17 +589,22 @@ def partial_squared_H_phi( return 0 loop_idx = cr.loops.index(loop) - B_idx = get_B_idx(cr, grad_el) + B_indices = get_B_indices(cr, grad_el) + H_squared = 0 if isinstance(grad_el, Junction): - return cr.B[B_idx, loop_idx] * cr._memory_ops['sin'][(grad_el, B_idx)] + for B_idx in B_indices: + H_squared += cr.B[B_idx, loop_idx] * cr._memory_ops['sin'][(grad_el, B_idx)] + return H_squared elif isinstance(grad_el, Inductor): - return ( - cr.B[B_idx, 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)] - ) + for B_idx in B_indices: + H_squared += ( + cr.B[B_idx, 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)] + ) + return H_squared else: raise NotImplementedError From a43f14482d80d84de6747da7fbb66fcdac940ff2 Mon Sep 17 00:00:00 2001 From: sambonkov Date: Mon, 3 Jun 2024 19:05:01 -0700 Subject: [PATCH 4/7] Add T1 support to phi_ext backprop --- SQcircuit/circuit.py | 51 +++++++++++++++++++++++++++++++----- SQcircuit/functions.py | 5 ++++ SQcircuit/tests/test_grad.py | 34 +++++++++++++++++++----- 3 files changed, 77 insertions(+), 13 deletions(-) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index f3fa4ef..b97fc06 100644 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -1482,7 +1482,11 @@ def _get_LC_hamil(self) -> Qobj: return LC_hamil - def _get_external_flux_at_element(self, B_idx: int) -> float: + def _get_external_flux_at_element( + self, + B_idx: int, + torch = False + ) -> Union[float, Tensor]: """ Return the external flux at an inductive element. @@ -1492,11 +1496,15 @@ def _get_external_flux_at_element(self, B_idx: int) -> float: An integer point to each row of B matrix (external flux distribution of that element) """ - phi_ext = 0.0 + phi_ext = sqf.zero() for i, loop in enumerate(self.loops): - phi_ext += sqf.numpy(loop.value()) * self.B[B_idx, i] + print(loop.value()) + phi_ext += loop.value() * self.B[B_idx, i] - return phi_ext + if isinstance(phi_ext, Tensor) and not torch: + return sqf.numpy(phi_ext) + else: + return phi_ext def _get_inductive_hamil(self) -> Qobj: @@ -1581,6 +1589,37 @@ def charge_op(self, mode: int, basis: str = 'FC') -> Qobj: return qt.Qobj(Q_eig) + def _get_W_idx(self, my_el: Junction, my_B_idx: int) -> Optional[int]: + for _, el, B_idx, W_idx in self.elem_keys[Junction]: + if el == my_el and B_idx == my_B_idx: + return W_idx + + return None + + def op(self, name: str, keywords: Dict): + """ + Returns the `name` operator of the circuit, specified by `keywords`, + as a sparse torch matrix. + """ + if name == 'sin_half': + B_idx = keywords['B_idx'] + el = keywords['el'] + if get_optim_mode(): + W_idx = self._get_W_idx(el, B_idx) + + phi = self._get_external_flux_at_element(B_idx, torch=True) + root_exp = ( + torch.exp(1j * phi / 2) + * sqf.qobj_to_tensor(self._memory_ops["root_exp"][W_idx]) + ) + + sin_half = (root_exp - sqf.dag(root_exp)) / 2j + return sin_half ## TO CHECK: need to squeeze? + else: + return self._memory_ops['sin_half'][el, B_idx] + else: + raise NotImplementedError + def diag_np( self, n_eig: int @@ -2127,8 +2166,8 @@ def dec_rate( sqf.operator_inner_product(state1, op, state2)) ** 2 if dec_type == "quasiparticle": - for el, _ in self._memory_ops['sin_half']: - op = self._memory_ops['sin_half'][(el, _)] + for el, B_idx in self._memory_ops['sin_half']: + op = self.op('sin_half', {'el': el, 'B_idx': B_idx}) if np.isnan(sqf.numpy(el.Y(omega, ENV["T"]))): decay = decay + 0 else: diff --git a/SQcircuit/functions.py b/SQcircuit/functions.py index aefeab8..1c05b5f 100644 --- a/SQcircuit/functions.py +++ b/SQcircuit/functions.py @@ -105,6 +105,11 @@ def init_op(size): # return torch.sparse_coo_tensor(size = size, dtype=torch.complex128) return qt.Qobj() +def zero(dtype=torch.complex128): + if get_optim_mode(): + return torch.tensor(0, dtype=dtype) + return 0 + def zeros(shape, dtype=torch.complex128): if get_optim_mode(): diff --git a/SQcircuit/tests/test_grad.py b/SQcircuit/tests/test_grad.py index 8253413..26b25f9 100644 --- a/SQcircuit/tests/test_grad.py +++ b/SQcircuit/tests/test_grad.py @@ -122,6 +122,8 @@ def function_grad_test(circuit_numpy: Circuit, set_optim_mode(False) loop_flux = loop.internal_value / 2 / np.pi + # assert 0 + delta < loop_flux < 1 - delta + # Calculate f(x+delta) loop.set_flux(loop_flux + delta) circuit_numpy.diag(num_eigenvalues) @@ -143,6 +145,7 @@ def function_grad_test(circuit_numpy: Circuit, grad_torch = circuit_torch.loops[loop_idx].internal_value.grad if grad_torch is not None: + print(f"loop #: {loop_idx}") print(f"grad torch: {grad_torch}, grad numpy: {grad_numpy}") grad_torch = grad_torch.detach().numpy() assert np.sign(grad_torch) == np.sign(grad_numpy) @@ -209,6 +212,13 @@ def test_omega_flux_fluxonium(): set_optim_mode(False) +def T1_inv(circuit): + return ( + circuit.dec_rate('capacitive', (0, 1)) + + circuit.dec_rate('inductive', (0, 1)) + + circuit.dec_rate('quasiparticle', (0, 1)) + ) + def test_T1_transmon(): """Compare gradient of T1 decoherence due to capacitive, inductive, and quasiparticle noise in transmon circuit with linearized value.""" @@ -217,13 +227,6 @@ def test_T1_transmon(): circuit_numpy = create_transmon_numpy(trunc_num) circuit_torch = create_transmon_torch(trunc_num) - def T1_inv(circuit): - return ( - circuit.dec_rate('capacitive', (0, 1)) - + circuit.dec_rate('inductive', (0, 1)) - + circuit.dec_rate('quasiparticle', (0, 1)) - ) - function_grad_test(circuit_numpy, T1_inv, circuit_torch, @@ -231,6 +234,23 @@ def T1_inv(circuit): delta=1e-6) set_optim_mode(False) +def test_T1_transmon_flux(): + """Compare gradient of T1 decoherence due to capacitive, inductive, and + quasiparticle noise in transmon circuit with linearized value.""" + + # Create circuits + for flux_point in [1+1e-2, 0.3, 0.5 - 1e-3, 0.5+1e-3, 0.7, 1-1e-2]: + print('flux point:', flux_point) + circuit_numpy = create_flux_transmon_numpy(trunc_num, flux_point) + circuit_torch = create_flux_transmon_torch(trunc_num, flux_point) + + function_grad_test(circuit_numpy, + T1_inv, + circuit_torch, + T1_inv, + delta=1e-6) + set_optim_mode(False) + def test_grad_multiple_steps(): """Sample ability of PyTorch to successfully update circuit parameters and From acb25629d0d7aaf57ea62fcd9efe3f7fc3f790fd Mon Sep 17 00:00:00 2001 From: sambonkov Date: Mon, 3 Jun 2024 20:23:58 -0700 Subject: [PATCH 5/7] Add phi_ext to T2 --- SQcircuit/tests/test_grad.py | 55 +++++++++++++++++++++++++++ SQcircuit/torch_extensions.py | 71 ++++++++++++++++++++++++++--------- 2 files changed, 109 insertions(+), 17 deletions(-) diff --git a/SQcircuit/tests/test_grad.py b/SQcircuit/tests/test_grad.py index 26b25f9..2f3e764 100644 --- a/SQcircuit/tests/test_grad.py +++ b/SQcircuit/tests/test_grad.py @@ -518,3 +518,58 @@ def test_T2_flux(): num_eigenvalues=50, delta=1e-6 ) + +def test_T2_cc_phi_ext(): + flux_points = [1e-2, 0.25, 0.5 - 1e-2, 0.5 + 1e-2, 0.75] + + for phi_ext in flux_points: + set_optim_mode(False) + circuit_numpy = create_fluxonium_numpy(trunc_num, phi_ext) + + # Create torch circuit + set_optim_mode(True) + circuit_torch = create_fluxonium_torch_flux(trunc_num, phi_ext) + + function_grad_test( + circuit_numpy, + lambda cr: cr.dec_rate('cc', states=(0, 1)), + circuit_torch, + lambda cr: cr.dec_rate('cc', states=(0, 1)), + num_eigenvalues=50, + delta=1e-6 + ) + +def test_T2_charge_phi_ext(): + charge_offsets = [1e-2, 0.2, 0.4, 0.5 - 1e-2, 0.5 + 1e-2, 0.6, 0.8, 1-1e-2] + phi_ext = 0.5-1e-3 + + for ng in charge_offsets: + circuit_numpy = create_flux_transmon_numpy(trunc_num, phi_ext) + circuit_numpy.set_charge_offset(1, ng) + + circuit_torch = create_flux_transmon_torch(trunc_num, phi_ext) + circuit_torch.set_charge_offset(1, ng) + + function_grad_test( + circuit_numpy, + lambda cr: cr.dec_rate('charge', states=(0, 1)), + circuit_torch, + lambda cr: cr.dec_rate('charge', states=(0, 1)), + num_eigenvalues=50, + ) + +def test_T2_flux_phi_ext(): + flux_points = [1e-2, 0.25, 0.5 - 1e-2, 0.5 + 1e-2, 0.75] + + for phi_ext in flux_points: + circuit_numpy = create_fluxonium_numpy(trunc_num, phi_ext) + circuit_torch = create_fluxonium_torch_flux(trunc_num, phi_ext) + + function_grad_test( + circuit_numpy, + lambda cr: cr.dec_rate('flux', states=(0, 1)), + circuit_torch, + lambda cr: cr.dec_rate('flux', states=(0, 1)), + num_eigenvalues=50, + delta=1e-6 + ) diff --git a/SQcircuit/torch_extensions.py b/SQcircuit/torch_extensions.py index 84ed069..c7ce68f 100644 --- a/SQcircuit/torch_extensions.py +++ b/SQcircuit/torch_extensions.py @@ -12,14 +12,14 @@ from SQcircuit import functions as sqf from SQcircuit import units as unt - +SupportedGradElements = Union[Capacitor, Inductor, Junction, Loop] ############################################################################### # Interface to custom torch functions ############################################################################### def eigencircuit(circuit: 'Circuit', n_eig: int): - """Given a circuit, computes the concatenated the tensor including both + """Given a circuit, computes the concatenated tensor including both eigenvalues and eigenvectors of a circuit via torch `Function`. Parameters @@ -145,7 +145,7 @@ def backward(ctx, grad_output) -> Tuple[Tensor]: def partial_squared_omega( cr: 'Circuit', - grad_el: Element, + grad_el: SupportedGradElements, partial_H: qt.Qobj, partial_H_squared: qt.Qobj, states: Tuple[int, int] @@ -171,7 +171,7 @@ def partial_squared_omega( """ m, n = states - + state_m = cr.evecs[m] partial_state_m = cr.get_partial_vec(grad_el, m) state_n = cr.evecs[n] @@ -198,10 +198,10 @@ def partial_squared_omega( def partial_dephasing_rate( - A, - partial_A, - partial_omega_mn, - partial_squared_omega_mn + A: float, + partial_A: float, + partial_omega_mn: float, + partial_squared_omega_mn: float ): """ Calculate the derivative of the dephasing rate with noise amplitude `A` @@ -281,7 +281,7 @@ def partial_H_ng( def partial_squared_H_ng( cr: 'Circuit', charge_idx: int, - grad_el: Union[Capacitor, Inductor, Junction] + grad_el: SupportedGradElements ): """ Calculates the second derivative of the Hamiltonian of `cr` with respect to the gate charge on the `charge_idx` charge mode and `grad_el`. @@ -344,7 +344,7 @@ def partial_omega_ng( def partial_squared_omega_mn_ng( cr: 'Circuit', charge_idx: int, - grad_el: Union[Capacitor, Inductor, Junction], + grad_el: SupportedGradElements, states: Tuple[int, int] ): """ Calculates the second derivative of the difference between the `m`th @@ -378,7 +378,7 @@ def partial_squared_omega_mn_ng( def partial_charge_dec( cr: 'Circuit', - grad_el: Element, + grad_el: SupportedGradElements, states: Tuple[int, int] ): """ @@ -452,12 +452,38 @@ def backward(ctx, grad_output) -> Tuple[Tensor, None, None]: # Critical current noise ############################################################################### +def partial_squared_H_EJ( + cr: 'Circuit', + EJ_el: Junction, + B_idx: int, + grad_el: SupportedGradElements +): + """ Calculates the second derivative of the Hamiltonian of `cr` with + respect to the Josephson energy of`EJ_el` charge mode and `grad_el`. + + Parameters + ---------- + cr: + The `Circuit` object to differentiate the Hamiltonian of. + EJ_el: + The Josephson junction to differentiate with respect to. + B_idx: + The index of the `cr.B` matrix identifying which branch the + `EJ_el` is on. + grad_el: + The circuit element to differentiate with respect to. + """ + if not isinstance(grad_el, Loop): + return 0 + + loop_idx = cr.loops.index(grad_el) + return cr.B[B_idx, loop_idx] * cr._memory_ops['sin'][(EJ_el, B_idx)] def partial_squared_omega_mn_EJ( cr: 'Circuit', EJ_el: Junction, B_idx: int, - grad_el: Union[Capacitor, Inductor, Junction], + grad_el: SupportedGradElements, states: Tuple[int, int] ): """ Calculates the second derivative of the difference between the `m`th @@ -477,7 +503,7 @@ def partial_squared_omega_mn_EJ( The numbers `(m, n)` of the eigenfrequencies to differentiate. """ partial_H = cr._get_partial_H(EJ_el, _B_idx = B_idx) - partial_H_squared = 0 + partial_H_squared = partial_squared_H_EJ(cr, EJ_el, B_idx, grad_el) return partial_squared_omega( cr, @@ -490,7 +516,7 @@ def partial_squared_omega_mn_EJ( def partial_cc_dec( cr: 'Circuit', - grad_el: Element, + grad_el: SupportedGradElements, states: Tuple[int, int] ): """ Calculate the derivative of the critical current decoherence between the @@ -571,7 +597,7 @@ def backward(ctx, grad_output) -> Tuple[Tensor, None, None]: def partial_squared_H_phi( cr: 'Circuit', loop: Loop, - grad_el: Union[Capacitor, Inductor, Junction] + grad_el: SupportedGradElements ): """ Calculates the second derivative of the Hamiltonian of `cr` with respect to the external flux through `loop` and `grad_el` @@ -605,6 +631,17 @@ def partial_squared_H_phi( * cr._memory_ops["ind_hamil"][(grad_el, B_idx)] ) 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]: + 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)] + ) + return H_squared else: raise NotImplementedError @@ -612,7 +649,7 @@ def partial_squared_H_phi( def partial_squared_omega_mn_phi( cr: 'Circuit', loop: Loop, - grad_el: Union[Capacitor, Inductor, Junction], + grad_el: SupportedGradElements, states: Tuple[int, int] ): """ Calculate the second derivative of the difference between the `m`th @@ -643,7 +680,7 @@ def partial_squared_omega_mn_phi( def partial_flux_dec( cr: 'Circuit', - grad_el: Element, + grad_el: SupportedGradElements, states: Tuple[int, int] ): """ Calculate the derivative of the flux decoherence between the states From 693cc3baf34abc19ad2876832e7188b6587b37ef Mon Sep 17 00:00:00 2001 From: sambonkov Date: Mon, 3 Jun 2024 20:40:53 -0700 Subject: [PATCH 6/7] Bugfix: load circuit tests in correct optim_mode --- SQcircuit/tests/conftest.py | 1 + SQcircuit/tests/test_qubit_grad.py | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/SQcircuit/tests/conftest.py b/SQcircuit/tests/conftest.py index d2eb0d2..50d9416 100644 --- a/SQcircuit/tests/conftest.py +++ b/SQcircuit/tests/conftest.py @@ -25,6 +25,7 @@ def setup_class(cls): def test_transform_process(self): # load the data + set_optim_mode(False) data = SQdata.load(DATADIR + "/" + self.fileName) # build the new circuit based on data circuit parameters diff --git a/SQcircuit/tests/test_qubit_grad.py b/SQcircuit/tests/test_qubit_grad.py index 86992c6..e69949e 100644 --- a/SQcircuit/tests/test_qubit_grad.py +++ b/SQcircuit/tests/test_qubit_grad.py @@ -81,7 +81,7 @@ def _get_fluxonium_efreqs( def test_fluxonium_grad(): - + sq.set_optim_mode(False) c_val = 3 l_val = 0.353 j_val = 10.2 @@ -193,7 +193,7 @@ def _get_zeropi_efreqs( def test_zeropi_grad(): - + sq.set_optim_mode(False) c_val = 129 cj_val = 1.93 l_val = 1.257 From 0516bf7c12242dc9ff99644560a461ba781237d3 Mon Sep 17 00:00:00 2001 From: sambonkov Date: Sat, 8 Jun 2024 16:46:12 -0700 Subject: [PATCH 7/7] Try to add test for JJL --- SQcircuit/circuit.py | 2 +- SQcircuit/tests/conftest.py | 46 ++++++++++++++++++++++++++++++++++++ SQcircuit/tests/test_grad.py | 19 +++++++++++++++ 3 files changed, 66 insertions(+), 1 deletion(-) diff --git a/SQcircuit/circuit.py b/SQcircuit/circuit.py index b97fc06..bbfc7ea 100644 --- a/SQcircuit/circuit.py +++ b/SQcircuit/circuit.py @@ -2382,7 +2382,7 @@ def _get_partial_omega_mn( partial_omega_n = sqf.operator_inner_product(state_n, partial_H, state_n) partial_omega_mn = partial_omega_m - partial_omega_n - assert sqf.imag(partial_omega_mn)/sqf.real(partial_omega_mn) < 1e-6 + # assert sqf.imag(partial_omega_mn)/sqf.real(partial_omega_mn) < 1e-6 return sqf.real(partial_omega_mn) diff --git a/SQcircuit/tests/conftest.py b/SQcircuit/tests/conftest.py index 50d9416..7456036 100644 --- a/SQcircuit/tests/conftest.py +++ b/SQcircuit/tests/conftest.py @@ -165,3 +165,49 @@ def create_fluxonium_torch_flux(trunc_num, phi_ext=0): circuit_torch = Circuit({(0, 1): [C_torch, L_torch, JJ_torch], }, flux_dist='junctions') circuit_torch.set_trunc_nums([trunc_num, ]) return circuit_torch + +def create_JJL_numpy(trunc_num, phi_ext): + set_optim_mode(False) + Cs = [3.6, 1.2, 2.0] + Js = [10.2, 6.8] + Ls = [0.46] + loop = Loop(phi_ext) + C1 = Capacitor(Cs[0], 'GHz', requires_grad=False) + C2 = Capacitor(Cs[1], 'GHz', requires_grad=False) + C3 = Capacitor(Cs[2], 'GHz', requires_grad=False) + L = Inductor(Ls[0], 'GHz', loops=[loop], requires_grad=False) + J1 = Junction(Js[0], 'GHz',loops=[loop], requires_grad=False) + J2 = Junction(Js[1], 'GHz',loops=[loop], requires_grad=False) + + circuit_numpy = Circuit( + {(0, 1): [C1, J1], + (1, 2): [C2, J2], + (2, 0): [L, C3] + }, + flux_dist='junctions' + ) + circuit_numpy.set_trunc_nums([trunc_num, trunc_num]) + return circuit_numpy + +def create_JJL_torch(trunc_num, phi_ext): + set_optim_mode(True) + Cs = [3.6, 1.2, 2.0] + Js = [10.2, 6.8] + Ls = [0.46] + loop = Loop(phi_ext, requires_grad=False) + C1 = Capacitor(Cs[0], 'GHz', requires_grad=True) + C2 = Capacitor(Cs[1], 'GHz', requires_grad=True) + C3 = Capacitor(Cs[2], 'GHz', requires_grad=True) + L = Inductor(Ls[0], 'GHz', loops=[loop], requires_grad=True) + J1 = Junction(Js[0], 'GHz',loops=[loop], requires_grad=True) + J2 = Junction(Js[1], 'GHz',loops=[loop], requires_grad=True) + + circuit_torch = Circuit( + {(0, 1): [C1, J1], + (1, 2): [C2, J2], + (2, 0): [L, C3] + }, + flux_dist='junctions' + ) + circuit_torch.set_trunc_nums([trunc_num, trunc_num]) + return circuit_torch diff --git a/SQcircuit/tests/test_grad.py b/SQcircuit/tests/test_grad.py index 2f3e764..400b6a2 100644 --- a/SQcircuit/tests/test_grad.py +++ b/SQcircuit/tests/test_grad.py @@ -20,6 +20,8 @@ create_flux_transmon_numpy, create_flux_transmon_torch, create_fluxonium_torch_flux, + create_JJL_numpy, + create_JJL_torch, ) @@ -573,3 +575,20 @@ def test_T2_flux_phi_ext(): num_eigenvalues=50, delta=1e-6 ) + +# def test_T2_flux_JJL(): +# flux_points = [0.5] #, 0.25, 0.5 - 1e-2, 0.5 + 1e-2, 0.75] + +# for phi_ext in flux_points: +# print('phi_ext', phi_ext) +# circuit_numpy = create_JJL_numpy(45, phi_ext) +# circuit_torch = create_JJL_torch(45, phi_ext) + +# function_grad_test( +# circuit_numpy, +# first_eigendifference_numpy, #lambda cr: cr.dec_rate('flux', states=(0, 1)), +# circuit_torch, +# first_eigendifference_torch, #lambda cr: cr.dec_rate('flux', states=(0, 1)), +# num_eigenvalues=50, +# delta=1e-4 +# )