scqubits - npm Package Compare versions

Comparing version 3.3.0 to 4.0.0

scqubits/core/circuit_input.py

+import pyparsing as pp
+import os
+
+from typing import List, Tuple
+from scqubits.utils.misc import is_string_float
+from pyparsing import Group, Opt, Or, Literal, Suppress
+import numpy as np
+import scipy as sp
+import sympy as sm
+
+# *****************************************************************
+#  OUR GRAMMAR DEFINITIONS
+# *****************************************************************
+#
+# Pattern for branch definitions in yaml file:
+#      - [ <branch_type>, <node1>, <node2>, <param>, <aux params> ]       or
+#      - [ <branch_type>, <node1>, <node2>, <param1>, <param2> ]
+#
+# where <param>:   <symbol> = <number>
+#                  <number>
+#                  <symbol>
+#
+# The last option is valid only if <symbol> has previously been
+# assigned. Optionally, <number> may be grouped with valid physical
+# unit.
+
+
+# - Ignore in parsing ********************************************
+# Mark the following as characters / character combinations not
+# to be recorded in the parsed result. I.e., the grammar may expect
+# these in various places, but they are not carried forward into
+# the parsed result
+BEG = Suppress(Literal("-") + Literal("["))
+END = Suppress("]")
+CM = Suppress(",")
+QM = Opt(Suppress('"'))  # optional quotation mark: may use JJ, or "JJ"
+
+
+# - Numbers ******************************************************
+INT = pp.common.integer  # unsigned integer
+NUM = pp.common.fnumber  # float
+
+
+# - Branch types ***************************************************
+branch_type_names = ["C", "L"]
+
+# build up dictionary of branch types
+# allow, for example, both "JJ" as well as just JJ
+BRANCH_TYPES = {name: QM + name + QM for name in branch_type_names}
+for BTYPE in BRANCH_TYPES.values():
+    BTYPE.set_results_name("branch_type")
+
+JJ_ORDER = pp.Word(pp.nums).add_condition(
+    lambda tokens: int(tokens[0]) > 0, message="Junction order must be greater than 0."
+)
+BRANCH_TYPES["JJ"] = (
+    QM + pp.Combine(pp.Word("JJ") + Opt(JJ_ORDER) + Opt(pp.Word("s"))) + QM
+)  # defining grammar to find "JJi" where i is an optional natural number
+
+
+# - Units: prefixes etc. **************************************************
+prefix_dict = {
+    "T": 1e12,  # Tera
+    "G": 1e9,  # Giga
+    "M": 1e6,  # Mega
+    "k": 1e3,  # kilo
+    "m": 1e-3,  # milli
+    "u": 1e-6,  # micro
+    "n": 1e-9,  # nano
+    "p": 1e-12,  # pico
+    "f": 1e-15,  # femto
+}
+PREFIX = pp.Char(list(prefix_dict.keys()))
+
+energy_names = ["EJ", "EC", "EL"]
+
+
+UNITS_FREQ_ENERGY = Literal("Hz") ^ Literal("J") ^ Literal("eV")
+
+UNITS = {name: Opt(PREFIX, None) for name in energy_names}
+UNITS["EJ"] += UNITS_FREQ_ENERGY ^ Literal("A")  # Ampere
+UNITS["EC"] += UNITS_FREQ_ENERGY ^ Literal("F")  # Farad
+UNITS["EL"] += UNITS_FREQ_ENERGY ^ Literal("H")  # Henry
+for name, unit in UNITS.items():
+    unit.leave_whitespace()  # allow only "kHz", not "k Hz"
+    unit.set_name(f"{name}_UNITS")
+
+# - Parameter specifications --------------------------------
+SYMBOL = pp.common.identifier
+VALUES = {name: NUM + Opt(" ") + Opt(UNITS[name], None) for name in energy_names}
+ASSIGNS = {
+    name: SYMBOL + Suppress(Literal("=")) + VALUES[name] for name in energy_names
+}
+
+PARAMS = {
+    name: Or(
+        [
+            Group(ASSIGNS[name])("ASSIGN"),
+            Group(SYMBOL)("SYMBOL"),
+            Group(VALUES[name])("VALUE"),
+        ]
+    )
+    for name in energy_names
+}  # can specify in three ways
+
+# # - Branch specifications ------------------------------------------------------
+aux_val = pp.Word(
+    pp.printables.replace("]", "").replace("[", "") + " "
+)  # allowing for numerical expressions in auxiliary params
+AUX_PARAM = Group(pp.ZeroOrMore(CM + SYMBOL + Suppress(Literal("=")) + aux_val))(
+    "AUX_PARAM"
+)
+
+order_count = pp.Empty()
+
+
+def find_jj_order(str_result: str, location: int, tokens: pp.ParseResults):
+    from scqubits.core.circuit_utils import _junction_order
+
+    JJ_TYPE = BEG + BRANCH_TYPES["JJ"]
+    JJ_TYPE.add_parse_action(lambda tokens: _junction_order(tokens[0]))
+    return JJ_TYPE.parse_string(str_result)
+
+
+order_count.set_parse_action(find_jj_order)
+
+BRANCH_JJ = (
+    BEG
+    + BRANCH_TYPES["JJ"]("BRANCH_TYPE")
+    + CM
+    + INT("node1")
+    + CM
+    + INT("node2")
+    + CM
+    + pp.counted_array(PARAMS["EJ"] + CM, int_expr=order_count)("EJ_VALS")
+    + PARAMS["EC"]("EC")
+    + AUX_PARAM
+    + END
+)
+
+BRANCH_C = (
+    BEG
+    + BRANCH_TYPES["C"]("BRANCH_TYPE")
+    + CM
+    + INT("node1")
+    + CM
+    + INT("node2")
+    + CM
+    + PARAMS["EC"]("EC")
+    + AUX_PARAM
+    + END
+)
+
+BRANCH_L = (
+    BEG
+    + BRANCH_TYPES["L"]("BRANCH_TYPE")
+    + CM
+    + INT("node1")
+    + CM
+    + INT("node2")
+    + CM
+    + PARAMS["EL"]("EL")
+    + AUX_PARAM
+    + END
+)
+
+BRANCHES = Or([BRANCH_JJ, BRANCH_C, BRANCH_L])
+
+# uncomment to create a html describing the grammar of this language
+# BRANCHES.create_diagram("branches.html")
+
+
+# - For filtering out only the code specifying the branches -----------------
+def remove_comments(code: str) -> str:
+    return pp.pythonStyleComment.suppress().transformString(code)
+
+
+def remove_branchline(code: str) -> str:
+    return Suppress(Literal("branches") + ":").transformString(code)
+
+
+def strip_empty_lines(code: str) -> str:
+    return os.linesep.join(
+        [line.lstrip() for line in code.splitlines() if line.lstrip()]
+    )
+
+
+pp.autoname_elements()
+
+
+# - Parsing and processing ParsedResults data ------------------------------
+def parse_code_line(code_line: str, _branch_count):
+    """
+
+    Args:
+        code_line (str): string describing the branch from the input file
+        _branch_count (_type_): the count of the branch in a given circuit
+
+    Returns:
+        branch_type: str
+        node_idx1: int
+        node_idx2: int
+        params: str
+        aux_params: str
+        _branch_count: int
+    """
+    pp_result = BRANCHES.parse_string(code_line)
+    branch_type = pp_result.BRANCH_TYPE[0]
+
+    branch_type, node_idx1, node_idx2, *params, aux_params = pp_result
+
+    return branch_type, node_idx1, node_idx2, params, aux_params, _branch_count
+
+
+def convert_value_to_GHz(val, units):
+    """
+    Converts a given value and units to energy in GHz. The units are given in a string in the format "pU"
+    where p is an optional multiplier prefix and U is units. For example: "pH", "nA", "fF", "eV"
+
+    Args:
+        val (float): value in given units
+        units (str): units described in a string with an optional multiplier prefix
+
+    Raises:
+        ValueError: If the unit is unknown.
+
+    Returns:
+        float: Energy in GHz
+    """
+    # all the possible units
+    # capacitance F, inductance H, current A, energy J, frequency Hz, eV
+    # returns value in GHz
+    if units is None:
+        return val  # default is GHz
+
+    prefix = 1 if units[0] is None else prefix_dict[units[0]]
+    val *= prefix
+    unit_str = units[1]
+
+    h = sp.constants.h
+    e = sp.constants.e
+    Φ0 = sp.constants.h / (2 * e)
+    if unit_str == "Hz":
+        return val * 1e-9
+    elif unit_str == "J":
+        return val / h * 1e-9
+    elif unit_str == "eV":
+        return val * 1.602176634e-19 / h * 1e-9
+    elif unit_str == "F":
+        return e**2 / (2 * val * h) * 1e-9
+    elif unit_str == "H":
+        return Φ0**2 / (val * h * (2 * np.pi) ** 2) * 1e-9
+    elif unit_str == "A":
+        return val * Φ0 / (2 * np.pi * h) * 1e-9
+    else:
+        raise ValueError(f"Unknown unit {unit_str}")
+
+
+def process_param(pattern):
+    """
+    Returns a tuple containing (symbol, value) given a pattern as detected by pyparsing.
+    Either the symbol or the value can be returned to be none, when the symbol is already assigned or no symbol is assigned to the given branch parameter.
+    """
+    name = pattern.getName()
+    if name == "ASSIGN":
+        sym = sm.symbols(pattern[0])
+        val = pattern[1]
+        units = None if pattern[-1] == None else pattern[-2:]
+        val_converted = convert_value_to_GHz(val, units)
+        return sym, val_converted
+    if name == "SYMBOL":
+        return sm.symbols(pattern[0]), None
+    if name == "VALUE":
+        units = None if pattern[-1] == None else pattern[-2:]
+        converted_val = convert_value_to_GHz(pattern[0], units)
+        return None, converted_val
+    if name == "AUX_PARAM":
+        num_of_aux_params = int(len(pattern) / 2)
+        aux_params = {}
+        for idx in range(num_of_aux_params):
+            aux_params[pattern[2 * idx]] = (
+                float(pattern[2 * idx + 1])
+                if is_string_float(pattern[2 * idx + 1])
+                else pattern[2 * idx + 1]
+            )
+        return aux_params

scqubits/core/circuit_noise.py

+# circuit_noise.py
+#
+# This file is part of scqubits: a Python package for superconducting qubits,
+# Quantum 5, 583 (2021). https://quantum-journal.org/papers/q-2021-11-17-583/
+#
+#    Copyright (c) 2019 and later, Jens Koch and Peter Groszkowski
+#    All rights reserved.
+#
+#    This source code is licensed under the BSD-style license found in the
+#    LICENSE file in the root directory of this source tree.
+############################################################################
+
+from abc import ABC
+import operator as builtin_op
+import functools
+from numpy import ndarray
+import numpy as np
+import qutip as qt
+import scipy as sp
+import re
+import copy
+
+from scqubits.core.noise import NOISE_PARAMS, NoisySystem, calc_therm_ratio
+from scqubits.core.circuit_utils import get_trailing_number, keep_terms_for_subsystem
+from scqubits.utils.misc import is_string_float, Qobj_to_scipy_csc_matrix
+import scqubits.core.units as units
+
+from types import MethodType
+from typing import Any, Callable, Dict, List, Optional, Tuple, Union
+
+import sympy as sm
+
+from scqubits.core.symbolic_circuit import Branch
+
+
+class NoisyCircuit(NoisySystem, ABC):
+    @staticmethod
+    def Q_from_branch(branch):
+        key = "Q_" + ("ind" if branch.type == "L" else "cap")
+        if key in branch.aux_params.keys():
+            Q_str = branch.aux_params[key]
+            if not is_string_float(Q_str):
+
+                def Q_func(omega, T):
+                    return eval(Q_str)
+
+                return Q_func
+            else:
+                return float(Q_str)
+        return None
+
+    def generate_methods_d_hamiltonian_d(self):
+        """
+        Generate methods which return the derivative of the Hamiltonian with respect to
+        offset charges, external fluxes and junction energies.
+        """
+
+        ext_flux_1_over_f_methods = {}
+        ng_1_over_f_methods = {}
+        cc_1_over_f_methods = {}
+
+        for param_sym in self.external_fluxes + self.offset_charges:
+
+            def param_derivative(self, param_sym=param_sym):
+                parent_instance = self.return_parent_circuit()
+                hamiltonian = parent_instance.fetch_symbolic_hamiltonian()
+                hamiltonian = parent_instance._hamiltonian_sym_for_numerics
+                hamiltonian = hamiltonian.subs("I", 1)
+                all_sym_parameters = (
+                    list(parent_instance.symbolic_params.keys())
+                    + parent_instance.external_fluxes
+                    + parent_instance.offset_charges
+                )
+                diff_sym_expr = hamiltonian.diff(param_sym)
+                # substitute all symbolic params
+                for param in all_sym_parameters:
+                    diff_sym_expr = diff_sym_expr.subs(
+                        param, getattr(parent_instance, param.name)
+                    )
+                # evaluate the expression
+                return parent_instance._evaluate_symbolic_expr(diff_sym_expr)
+
+            if param_sym in self.external_fluxes:
+                ext_flux_1_over_f_methods[
+                    f"d_hamiltonian_d_flux{get_trailing_number(param_sym.name)}"
+                ] = param_derivative
+            elif param_sym in self.offset_charges:
+                ng_1_over_f_methods[
+                    f"d_hamiltonian_d_ng{get_trailing_number(param_sym.name)}"
+                ] = param_derivative
+        ## cc noise methods
+        junction_branches = [branch for branch in self.branches if "JJ" in branch.type]
+        for idx, branch in enumerate(junction_branches):
+
+            def param_derivative(self, branch=branch):
+                return -self.junction_related_evaluation(branch, calc="dhdEJ")
+
+            cc_1_over_f_methods[f"d_hamiltonian_d_EJ{branch.id_str}"] = param_derivative
+        noise_helper_methods = {
+            **ext_flux_1_over_f_methods,
+            **ng_1_over_f_methods,
+            **cc_1_over_f_methods,
+        }
+        self.noise_helper_methods = noise_helper_methods
+        for method_name in noise_helper_methods:
+            setattr(
+                self, method_name, MethodType(noise_helper_methods[method_name], self)
+            )
+
+    def _transform_expr_to_new_variables(
+        self, expr_node_vars: sm.Expr, substitute_symbol: Optional[str] = None
+    ):
+        transformation_mat = self.transformation_matrix
+        expr_node_vars = expr_node_vars.expand()
+        num_vars = len(self.symbolic_circuit.nodes) - (1 if self.is_grounded else 0)
+        new_vars = [sm.symbols(f"θ{index}") for index in range(1, 1 + num_vars)]
+        old_vars = [sm.symbols(f"φ{index}") for index in range(1, 1 + num_vars)]
+        transformed_expr = transformation_mat.dot(new_vars)
+        for idx, var in enumerate(old_vars):
+            expr_node_vars = expr_node_vars.subs(var, transformed_expr[idx])
+
+        if substitute_symbol:
+            for var in expr_node_vars.free_symbols:
+                expr_node_vars = expr_node_vars.subs(
+                    var,
+                    sm.symbols(f"{substitute_symbol}{get_trailing_number(var.name)}"),
+                )
+        return expr_node_vars
+
+    def junction_related_evaluation(self, branch_junction: Branch, calc="dhdEJ"):
+        parent_instance = self.return_parent_circuit()
+        hamiltonian = parent_instance.fetch_symbolic_hamiltonian()
+        hamiltonian = parent_instance._hamiltonian_sym_for_numerics
+
+        for sym in parent_instance.offset_charges + list(
+            parent_instance.symbolic_params.keys()
+        ):
+            hamiltonian = hamiltonian.subs(sym, getattr(parent_instance, sym.name))
+        hamiltonian = hamiltonian.subs("I", 1)
+        branch_cos_node_expr = sm.cos(
+            sm.symbols(f"φ{branch_junction.nodes[0].index}")
+            - sm.symbols(f"φ{branch_junction.nodes[1].index}")
+        )
+        branch_cos_node_expr = branch_cos_node_expr.subs(
+            "φ0", 0
+        )  # setting ground node to zero.
+        branch_cos_expr = parent_instance._transform_expr_to_new_variables(
+            branch_cos_node_expr
+        )
+        expr_dict = hamiltonian.as_coefficients_dict()
+        for term, coefficient in expr_dict.items():
+            term_without_ext_flux = copy.copy(term)
+            for flux in parent_instance.external_fluxes:
+                term_without_ext_flux = term_without_ext_flux.subs(flux, 0)
+            if term_without_ext_flux == branch_cos_expr:
+                break
+        # substitute external flux
+        for flux in parent_instance.external_fluxes:
+            term = term.subs(flux, getattr(parent_instance, flux.name))
+        if calc == "sin_phi_qp":
+            term = term.subs(sm.cos, sm.sin)
+            term = term.subs(term.args[0], term.args[0] / 2)
+
+        # evaluate the expression
+        return parent_instance._evaluate_symbolic_expr(term)
+
+    def generate_tphi_1_over_f_methods(self):
+        """Generate methods tphi_1_over_f_{noise_type}{index} methods for
+        noise_type=['cc', 'ng', 'flux']; individual noise sources differentiated
+        using index."""
+        # calculating the rates from each of the flux sources
+        junction_branches = [branch for branch in self.branches if "JJ" in branch.type]
+        methods_noise_rates_from_flux = {}
+        methods_noise_rates_from_ng = {}
+        methods_noise_rates_from_cc = {}
+        for param_sym in self.external_fluxes + self.offset_charges + junction_branches:
+            if param_sym in self.external_fluxes:
+                diff_func_name = "d_hamiltonian_d_flux"
+                noise_type = "flux"
+            elif param_sym in self.offset_charges:
+                diff_func_name = "d_hamiltonian_d_ng"
+                noise_type = "ng"
+            if param_sym in junction_branches:
+                diff_func_name = "d_hamiltonian_d_EJ"
+                noise_type = "cc"
+            if isinstance(param_sym, sm.Expr):
+                trailing_number = get_trailing_number(param_sym.name)
+                noise_op_func = getattr(self, f"{diff_func_name}{trailing_number}")
+            elif param_sym in junction_branches:
+                trailing_number = param_sym.id_str
+                noise_op_func = getattr(self, f"{diff_func_name}{trailing_number}")
+
+            def tphi_1_over_f_func(
+                self=self,
+                A_noise: float = NOISE_PARAMS[f"A_{noise_type}"],
+                i: int = 0,
+                j: int = 1,
+                esys: Tuple[ndarray, ndarray] = None,
+                get_rate: bool = False,
+                noise_op_func=noise_op_func,
+                **kwargs,
+            ) -> float:
+                r"""
+                Calculate the 1/f dephasing time (or rate) due to critical current noise of junction associated with
+                Josephson energy :math:`EJ3`.
+
+                Parameters
+                ----------
+                A_noise:
+                    noise strength
+                i:
+                    state index that along with j defines a qubit
+                j:
+                    state index that along with i defines a qubit
+                esys:
+                    evals, evecs tuple
+                get_rate:
+                    get rate or time
+
+                Returns
+                -------
+                    decoherence time in units of :math:`2\pi ({\rm system\,\,units})`, or rate in inverse units.
+                """
+                noise_op = noise_op_func()
+                if isinstance(noise_op, qt.Qobj):
+                    noise_op = Qobj_to_scipy_csc_matrix(noise_op)
+
+                return self.tphi_1_over_f(
+                    A_noise=A_noise,
+                    i=i,
+                    j=j,
+                    noise_op=noise_op,
+                    esys=esys,
+                    get_rate=get_rate,
+                    **kwargs,
+                )
+
+            if param_sym in self.external_fluxes:
+                methods_noise_rates_from_flux[
+                    f"tphi_1_over_f_flux{trailing_number}"
+                ] = tphi_1_over_f_func
+            elif param_sym in self.offset_charges:
+                methods_noise_rates_from_ng[f"tphi_1_over_f_ng{trailing_number}"] = (
+                    tphi_1_over_f_func
+                )
+            elif param_sym in junction_branches:
+                methods_noise_rates_from_cc[f"tphi_1_over_f_cc{trailing_number}"] = (
+                    tphi_1_over_f_func
+                )
+
+        noise_methods = {
+            **methods_noise_rates_from_flux,
+            **methods_noise_rates_from_ng,
+            **methods_noise_rates_from_cc,
+        }
+        for method_name in noise_methods:
+            setattr(self, method_name, MethodType(noise_methods[method_name], self))
+
+    def generate_overall_tphi_cc(self):
+        if not any([re.match(r"tphi_1_over_f_cc\d+$", method) for method in dir(self)]):
+            return None
+
+        def tphi_1_over_f_cc(
+            self=self,
+            A_noise: float = NOISE_PARAMS["A_cc"],
+            i: int = 0,
+            j: int = 1,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+            **kwargs,
+        ) -> float:
+            tphi_times = []
+            for branch in [brnch for brnch in self.branches if "JJ" in brnch.type]:
+                tphi_times.append(
+                    getattr(self, f"tphi_1_over_f_cc{branch.id_str}")(
+                        A_noise=A_noise, i=i, j=j, esys=esys, **kwargs
+                    )
+                )
+            total_rate = sum([1 / tphi for tphi in tphi_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(self, "tphi_1_over_f_cc", MethodType(tphi_1_over_f_cc, self))
+
+    def generate_overall_tphi_flux(self):
+        if not any(
+            [re.match(r"tphi_1_over_f_flux\d+$", method) for method in dir(self)]
+        ):
+            return None
+
+        def tphi_1_over_f_flux(
+            self=self,
+            A_noise: float = NOISE_PARAMS["A_flux"],
+            i: int = 0,
+            j: int = 1,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+            **kwargs,
+        ) -> float:
+            tphi_times = []
+            for flux_sym in self.external_fluxes:
+                tphi_times.append(
+                    getattr(
+                        self, f"tphi_1_over_f_flux{get_trailing_number(flux_sym.name)}"
+                    )(
+                        A_noise=A_noise,
+                        i=i,
+                        j=j,
+                        esys=esys,
+                    )
+                )
+            total_rate = sum([1 / tphi for tphi in tphi_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(self, "tphi_1_over_f_flux", MethodType(tphi_1_over_f_flux, self))
+
+    def generate_overall_tphi_ng(self):
+        if not any([re.match(r"tphi_1_over_f_ng\d+$", method) for method in dir(self)]):
+            return None
+
+        def tphi_1_over_f_ng(
+            self=self,
+            A_noise: float = NOISE_PARAMS["A_ng"],
+            i: int = 0,
+            j: int = 1,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+            **kwargs,
+        ) -> float:
+            tphi_times = []
+            for flux_sym in self.offset_charges:
+                tphi_times.append(
+                    getattr(
+                        self, f"tphi_1_over_f_ng{get_trailing_number(flux_sym.name)}"
+                    )(
+                        A_noise=A_noise,
+                        i=i,
+                        j=j,
+                        esys=esys,
+                    )
+                )
+            total_rate = sum([1 / tphi for tphi in tphi_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(self, "tphi_1_over_f_ng", MethodType(tphi_1_over_f_ng, self))
+
+    def generate_t1_flux_bias_line_methods(self):
+        """
+        Generate methods for flux bias line t1 coherence times.
+        """
+        flux_bias_line_methods = {}
+        for flux_sym in self.external_fluxes:
+            trailing_number = get_trailing_number(flux_sym.name)
+            noise_op_method = getattr(self, f"d_hamiltonian_d_flux{trailing_number}")
+
+            def flux_bias_noise(
+                self=self,
+                i: int = 1,
+                j: int = 0,
+                M: float = NOISE_PARAMS["M"],
+                Z: Union[complex, float, Callable] = NOISE_PARAMS["R_0"],
+                T: float = NOISE_PARAMS["T"],
+                total: bool = True,
+                esys: Tuple[ndarray, ndarray] = None,
+                get_rate: bool = False,
+                noise_op_method=noise_op_method,
+            ):
+                return NoisySystem.t1_flux_bias_line(
+                    self=self,
+                    i=i,
+                    j=j,
+                    M=M,
+                    Z=Z,
+                    T=T,
+                    total=total,
+                    esys=esys,
+                    get_rate=get_rate,
+                    noise_op_method=noise_op_method,
+                )
+
+            flux_bias_line_methods[f"t1_flux_bias_line{trailing_number}"] = (
+                flux_bias_noise
+            )
+
+        for method_name in flux_bias_line_methods:
+            setattr(
+                self, method_name, MethodType(flux_bias_line_methods[method_name], self)
+            )
+
+    def generate_t1_methods(self):
+        t1_capacitive_methods = {}
+        t1_inductive_methods = {}
+        t1_charge_impedance_methods = {}
+        t1_quasiparticle_tunneling_methods = {}
+
+        for branch in self.branches:
+            if branch.type == "L":
+                t1_inductive_methods[f"t1_inductive{branch.id_str}"] = (
+                    self.wrapper_t1_inductive_capacitive(branch)
+                )
+            else:
+                t1_capacitive_methods[f"t1_capacitive{branch.id_str}"] = (
+                    self.wrapper_t1_inductive_capacitive(branch)
+                )
+            # # quasiparticle noise
+            # if "JJ" in branch.type:
+            #     t1_quasiparticle_tunneling_methods[
+            #         f"t1_quasiparticle_tunneling{branch.id_str}"
+            #     ] = self.wrapper_t1_quasiparticle_tunneling(branch)
+            # quasiparticle noise methods are not included yet
+        noise_methods = {
+            **t1_capacitive_methods,
+            **t1_inductive_methods,
+            **t1_charge_impedance_methods,
+        }
+        for method_name in noise_methods:
+            setattr(self, method_name, MethodType(noise_methods[method_name], self))
+        # self._data.update(t1_quasiparticle_tunneling_methods)
+
+    def wrapper_t1_quasiparticle_tunneling(self, branch: Branch):
+        def t1_quasiparticle_tunneling(
+            self=self,
+            i: int = 1,
+            j: int = 0,
+            Y_qp: Union[float, Callable] = None,
+            x_qp: float = NOISE_PARAMS["x_qp"],
+            T: float = NOISE_PARAMS["T"],
+            Delta: float = NOISE_PARAMS["Delta"],
+            total: bool = True,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+        ) -> float:
+            return NoisySystem.t1_quasiparticle_tunneling(
+                self=self,
+                i=i,
+                j=j,
+                Y_qp=Y_qp,
+                x_qp=x_qp,
+                T=T,
+                Delta=Delta,
+                total=total,
+                esys=esys,
+                get_rate=get_rate,
+                noise_op=self.junction_related_evaluation(branch, calc="sin_phi_qp"),
+            )
+
+        return t1_quasiparticle_tunneling
+
+    def wrapper_t1_charge_impedance(self, branch: Branch):
+        def t1_charge_impedance(
+            self,
+            i: int = 1,
+            j: int = 0,
+            Z: Union[float, Callable] = NOISE_PARAMS["R_0"],
+            T: float = NOISE_PARAMS["T"],
+            total: bool = True,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+            branch=branch,
+        ) -> float:
+            parent_circuit = self.return_parent_circuit()
+            branch_var_expr = parent_circuit.symbolic_circuit._branch_sym_expr(
+                branch, return_charge=False if branch.type == "L" else True
+            )
+
+            if branch.type != "L":
+                branch_param = (
+                    branch.parameters["EC"]
+                    if branch.type == "C"
+                    else branch.parameters["ECJ"]
+                )
+            else:
+                branch_param = branch.parameters["EL"]
+            if isinstance(branch_param, sm.Expr):
+                branch_param = getattr(parent_circuit, branch_param.name)
+
+            return NoisySystem.t1_charge_impedance(
+                self=self,
+                i=i,
+                j=j,
+                Z=Z,
+                T=T,
+                total=total,
+                esys=esys,
+                get_rate=get_rate,
+                noise_op=parent_circuit._evaluate_symbolic_expr(branch_var_expr),
+            )
+
+        return t1_charge_impedance
+
+    def wrapper_t1_inductive_capacitive(
+        self,
+        branch: Branch,
+    ):
+        if branch.type != "L":
+
+            def t1_method(
+                self,
+                i: int = 1,
+                j: int = 0,
+                Q_cap: Union[float, Callable] = None,
+                T: float = NOISE_PARAMS["T"],
+                total: bool = True,
+                esys: Tuple[ndarray, ndarray] = None,
+                get_rate: bool = False,
+                branch: Branch = branch,
+            ) -> float:
+                parent_circuit = self.return_parent_circuit()
+                branch_charge_expr = parent_circuit.symbolic_circuit._branch_sym_expr(
+                    branch, return_charge=True
+                )
+
+                branch_param = (
+                    branch.parameters["EC"]
+                    if branch.type == "C"
+                    else branch.parameters["ECJ"]
+                )
+                if isinstance(branch_param, sm.Expr):
+                    branch_param = getattr(parent_circuit, branch_param.name)
+
+                return NoisySystem.t1_capacitive(
+                    self=self,
+                    i=i,
+                    j=j,
+                    Q_cap=Q_cap or self.Q_from_branch(branch),
+                    T=T,
+                    total=total,
+                    esys=esys,
+                    get_rate=get_rate,
+                    noise_op=parent_circuit._evaluate_symbolic_expr(branch_charge_expr),
+                    branch_params=branch_param,
+                )
+
+        else:
+
+            def t1_method(
+                self,
+                i: int = 1,
+                j: int = 0,
+                Q_ind: Union[float, Callable] = None,
+                T: float = NOISE_PARAMS["T"],
+                total: bool = True,
+                esys: Tuple[ndarray, ndarray] = None,
+                get_rate: bool = False,
+                branch: Branch = branch,
+            ) -> float:
+                parent_circuit = self.return_parent_circuit()
+                branch_var_expr = parent_circuit.symbolic_circuit._branch_sym_expr(
+                    branch
+                )
+
+                branch_param = branch.parameters["EL"]
+
+                if isinstance(branch_param, sm.Expr):
+                    branch_param = getattr(parent_circuit, branch_param.name)
+
+                return NoisySystem.t1_inductive(
+                    self=self,
+                    i=i,
+                    j=j,
+                    Q_ind=Q_ind or self.Q_from_branch(branch),
+                    T=T,
+                    total=total,
+                    esys=esys,
+                    get_rate=get_rate,
+                    noise_op=parent_circuit._evaluate_symbolic_expr(branch_var_expr),
+                    branch_params=branch_param,
+                )
+
+        return t1_method
+
+    def generate_overall_t1_quasiparticle_tunneling(self):
+        if not any(
+            [
+                re.match(r"t1_quasiparticle_tunneling\d+$", method)
+                for method in dir(self)
+            ]
+        ):
+            return None
+        if self.is_purely_harmonic:
+            return None
+
+        def t1_quasiparticle_tunneling(
+            self=self,
+            i: int = 1,
+            j: int = 0,
+            Y_qp: Union[float, Callable] = None,
+            x_qp: float = NOISE_PARAMS["x_qp"],
+            T: float = NOISE_PARAMS["T"],
+            Delta: float = NOISE_PARAMS["Delta"],
+            total: bool = True,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+        ) -> float:
+            t1_times = []
+            for branch in [b for b in self.branches if "JJ" in b.type]:
+                t1_times.append(
+                    getattr(self, f"t1_quasiparticle_tunneling{branch.id_str}")(
+                        i=i,
+                        j=j,
+                        Y_qp=Y_qp,
+                        x_qp=x_qp,
+                        T=T,
+                        Delta=Delta,
+                        total=total,
+                        esys=esys,
+                    )
+                )
+            total_rate = sum([1 / t1 for t1 in t1_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(
+            self,
+            "t1_quasiparticle_tunneling",
+            MethodType(t1_quasiparticle_tunneling, self),
+        )
+
+    def generate_overall_t1_inductive(self):
+        if not any([re.match(r"t1_inductive\d+$", method) for method in dir(self)]):
+            return None
+
+        def t1_method(
+            self,
+            i: int = 1,
+            j: int = 0,
+            Q_ind: Union[float, Callable] = None,
+            T: float = NOISE_PARAMS["T"],
+            total: bool = True,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+        ) -> float:
+            t1_times = []
+            parent_circuit = self.return_parent_circuit()
+            for branch in [b for b in parent_circuit.branches if b.type == "L"]:
+                t1_times.append(
+                    getattr(parent_circuit, f"t1_inductive{branch.id_str}")(
+                        i=i,
+                        j=j,
+                        Q_ind=Q_ind or self.Q_from_branch(branch),
+                        T=T,
+                        total=total,
+                        esys=esys,
+                    )
+                )
+            total_rate = sum([1 / t1 for t1 in t1_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(self, "t1_inductive", MethodType(t1_method, self))
+
+    def generate_overall_t1_capacitive(self):
+        if not any([re.match(r"t1_capacitive\d+$", method) for method in dir(self)]):
+            return None
+
+        def t1_method(
+            self,
+            i: int = 1,
+            j: int = 0,
+            Q_cap: Union[float, Callable] = None,
+            T: float = NOISE_PARAMS["T"],
+            total: bool = True,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+        ) -> float:
+            t1_times = []
+            parent_circuit = self.return_parent_circuit()
+            for branch in [b for b in parent_circuit.branches if b.type != "L"]:
+                t1_times.append(
+                    getattr(parent_circuit, f"t1_capacitive{branch.id_str}")(
+                        i=i,
+                        j=j,
+                        Q_cap=Q_cap or self.Q_from_branch(branch),
+                        T=T,
+                        total=total,
+                        esys=esys,
+                    )
+                )
+            total_rate = sum([1 / t1 for t1 in t1_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(self, "t1_capacitive", MethodType(t1_method, self))
+
+    def generate_overall_t1_charge_impedance(self):
+        if not any(
+            [re.match(r"t1_charge_impedance\d+$", method) for method in dir(self)]
+        ):
+            return None
+
+        def t1_method(
+            self=self,
+            i: int = 1,
+            j: int = 0,
+            Z: Union[float, Callable] = NOISE_PARAMS["R_0"],
+            T: float = NOISE_PARAMS["T"],
+            total: bool = True,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+        ) -> float:
+            t1_times = []
+            parent_circuit = self.return_parent_circuit()
+            for branch in [b for b in parent_circuit.branches if b.type != "L"]:
+                t1_times.append(
+                    getattr(parent_circuit, f"t1_charge_impedance{branch.id_str}")(
+                        i=i,
+                        j=j,
+                        Z=Z,
+                        T=T,
+                        total=total,
+                        esys=esys,
+                    )
+                )
+            total_rate = sum([1 / t1 for t1 in t1_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(self, "t1_charge_impedance", MethodType(t1_method, self))
+
+    def generate_overall_t1_flux_bias_line(self):
+        if not any(
+            [re.match(r"t1_flux_bias_line\d+$", method) for method in dir(self)]
+        ):
+            return None
+
+        def t1_flux_bias_line(
+            self=self,
+            i: int = 1,
+            j: int = 0,
+            M: float = NOISE_PARAMS["M"],
+            Z: Union[complex, float, Callable] = NOISE_PARAMS["R_0"],
+            T: float = NOISE_PARAMS["T"],
+            total: bool = True,
+            esys: Tuple[ndarray, ndarray] = None,
+            get_rate: bool = False,
+        ) -> float:
+            t1_times = []
+            for external_flux_sym in self.external_fluxes:
+                t1_times.append(
+                    getattr(
+                        self,
+                        f"t1_flux_bias_line{get_trailing_number(external_flux_sym.name)}",
+                    )(
+                        i=i,
+                        j=j,
+                        M=M,
+                        Z=Z,
+                        T=T,
+                        total=total,
+                        esys=esys,
+                    )
+                )
+            total_rate = sum([1 / t1 for t1 in t1_times])
+            if get_rate:
+                return total_rate
+            return 1 / total_rate if total_rate != 0 else np.inf
+
+        setattr(self, "t1_flux_bias_line", MethodType(t1_flux_bias_line, self))
+
+    def generate_noise_methods(self):
+        self._frozen = False
+        self.generate_methods_d_hamiltonian_d()
+        self.generate_tphi_1_over_f_methods()
+        self.generate_t1_flux_bias_line_methods()
+        self.generate_t1_methods()
+
+        self.generate_overall_tphi_cc()
+        self.generate_overall_tphi_flux()
+        self.generate_overall_tphi_ng()
+        self.generate_overall_t1_capacitive()
+        self.generate_overall_t1_charge_impedance()
+        self.generate_overall_t1_inductive()
+        self.generate_overall_t1_flux_bias_line()
+        self.generate_overall_t1_quasiparticle_tunneling()
+        self._noise_methods_generated = True
+        self._frozen = True

optional-requirements.txt

@@ -6,4 +6,3 @@ h5py>=2.10
 pathos>=0.3.0
-dill
 traitlets
 typing_extensions

PKG-INFO

 Metadata-Version: 2.1
 Name: scqubits
-Version: 3.3.0
+Version: 4.0.0
 Summary: scqubits: superconducting qubits in Python
@@ -5,0 +5,0 @@ Home-page: https://scqubits.readthedocs.io

requirements.txt

@@ -5,7 +5,7 @@ cycler
 numpy>=1.14.2
-pyyaml
 qutip>=4.3.1
-scipy>=1.1.0
+scipy>=1.5
+dill
 sympy
 tqdm
 typing_extensions

scqubits.egg-info/PKG-INFO

 Metadata-Version: 2.1
 Name: scqubits
-Version: 3.3.0
+Version: 4.0.0
 Summary: scqubits: superconducting qubits in Python
@@ -5,0 +5,0 @@ Home-page: https://scqubits.readthedocs.io

scqubits.egg-info/requires.txt

@@ -5,5 +5,5 @@ cycler
 numpy>=1.14.2
-pyyaml
 qutip>=4.3.1
-scipy>=1.1.0
+scipy>=1.5
+dill
 sympy
@@ -10,0 +10,0 @@ tqdm

scqubits.egg-info/SOURCES.txt

@@ -20,2 +20,5 @@ LICENSE
 scqubits/core/circuit.py
+scqubits/core/circuit_input.py
+scqubits/core/circuit_noise.py
+scqubits/core/circuit_routines.py
 scqubits/core/circuit_utils.py
@@ -22,0 +25,0 @@ scqubits/core/constants.py

scqubits/__init__.py

@@ -43,2 +43,3 @@ # scqubits: superconducting qubits in Python
 from scqubits.core.storage import DataStore, SpectrumData
+from scqubits.core.symbolic_circuit import SymbolicCircuit
 from scqubits.core.transmon import Transmon, TunableTransmon
@@ -95,2 +96,11 @@ from scqubits.core.units import (
 
+# Import of custom-circuit modules needs to take place after other imports to
+# avoid circular import issues
+from scqubits.core.circuit import Circuit
+from scqubits.core.circuit_utils import (
+    truncation_template,
+    assemble_circuit,
+    assemble_transformation_matrix,
+)
+
 # version
@@ -97,0 +107,0 @@ try:

scqubits/core/central_dispatch.py

@@ -38,2 +38,3 @@ # central_dispatch.py
     "PARAMETERSWEEP_UPDATE",
+    "CIRCUIT_UPDATE",
 ]
@@ -40,0 +41,0 @@

scqubits/core/circuit_utils.py

@@ -14,6 +14,8 @@ # circuit-utils.py
 import re
-from typing import TYPE_CHECKING, Any, Callable, List, Union
+from typing import TYPE_CHECKING, Any, Callable, List, Union, Optional, Tuple, Dict
 
 import numpy as np
+import scipy as sp
 import sympy as sm
+import qutip as qt
 from numpy import ndarray
@@ -24,3 +26,9 @@ from scipy import sparse
 from scqubits.core import discretization as discretization
-from scqubits.utils.misc import flatten_list_recursive
+from scqubits.utils.misc import (
+    flatten_list_recursive,
+    is_string_float,
+    unique_elements_in_list,
+    Qobj_to_scipy_csc_matrix,
+)
+from scqubits.core import circuit_input
 
@@ -31,2 +39,75 @@ if TYPE_CHECKING:
 
+def _junction_order(branch_type: str) -> int:
+    """
+    Returns the order of the branch if it is a JJ branch,
+    if the order is n its energy is given by cos(phi) + cos(2*phi) + ... + cos(n*phi)
+
+    Args:
+        branch (_type_): Branch
+
+    Raises:
+        ValueError: when the branch is not a josephson junction
+
+    Returns:
+        _type_: int, order of the josephson junction
+    """
+    if "JJ" not in branch_type:
+        raise ValueError("The branch is not a JJ branch")
+    if len(branch_type) > 2:
+        if (
+            branch_type[2] == "s"
+        ):  # adding "JJs" which is a junction with sawtooth current phase relationship
+            return 1
+        return int(branch_type[2:])
+    else:
+        return 1
+
+
+def sawtooth_operator(x: Union[ndarray, csc_matrix]):
+    """
+    Returns the operator evaluated using applying the sawtooth_potential function on the
+    diagonal elements of the operator x
+
+    Args:
+        x (Union[ndarray, csc_matrix]): argument of the sawtooth operator in the Hamiltonian
+    """
+    diagonal_elements = sawtooth_potential(x.diagonal())
+
+    operator = sp.sparse.dia_matrix(
+        (diagonal_elements, 0), shape=(len(diagonal_elements), len(diagonal_elements))
+    )
+    return operator.tocsc()
+
+
+# def sawtooth_potential(x: float) -> float:
+#     """
+#     Is the function which returns the potential of a sawtooth junction,
+#     i.e. a junction with a sawtooth current phase relationship, only in the discretized phi basis.
+#     """
+#     x_rel = (x - np.pi) % (2*np.pi) - np.pi
+#     return (x_rel)**2/(np.pi)**2 # normalized to have a maximum of 1
+
+
+def sawtooth_potential(phi_pts):
+    # definition from Andras
+    skewness = 0.99
+    N = 1000
+    V = np.zeros_like(phi_pts)
+    for idx in range(1, N + 1):
+        V += (skewness + 1) * (-skewness) ** (idx - 1) * np.cos(idx * phi_pts) / idx**2
+    return -V
+
+
+def _capactiance_variable_for_branch(branch_type: str):
+    """
+    Returns the parameter name that stores the capacitance of the branch
+    """
+    if "C" in branch_type:
+        return "EC"
+    elif "JJ" in branch_type:
+        return "ECJ"
+    else:
+        raise ValueError("Branch type is not a capacitor or a JJ")
+
+
 def truncation_template(
@@ -214,3 +295,3 @@ system_hierarchy: list, individual_trunc_dim: int = 6, combined_trunc_dim: int = 30
     sin_op = sparse.dia_matrix((pt_count, pt_count))
-    diag_elements = np.cos(grid.make_linspace())
+    diag_elements = np.sin(grid.make_linspace())
     sin_op.setdiag(diag_elements)
@@ -240,9 +321,11 @@ return sin_op.tocsc()
 
-def _exp_i_theta_operator(ncut) -> csc_matrix:
+def _exp_i_theta_operator(ncut, prefactor=1) -> csc_matrix:
     r"""
     Operator :math:`\cos(\theta)`, acting only on the `\theta` Hilbert subspace.
     """
+    # if type(prefactor) != int:
+    #     raise ValueError("Prefactor must be an integer")
     dim_theta = 2 * ncut + 1
     matrix = sparse.dia_matrix(
-        (np.ones(dim_theta), [-1]),
+        (np.ones(dim_theta), [-prefactor]),
         shape=(dim_theta, dim_theta),
@@ -348,10 +431,5 @@ ).tocsc()
     def operator_func(self: "Subsystem"):
-        if not self.hierarchical_diagonalization:
-            return self._kron_operator(
-                inner_op(self.grids_dict_for_discretized_extended_vars()[index]), index
-            )
-        else:
-            return self.identity_wrap_for_hd(
-                inner_op(self.grids_dict_for_discretized_extended_vars()[index]), index
-            ).data.tocsc()
+        return self._kron_operator(
+            inner_op(self.discretized_grids_dict_for_vars()[index]), index
+        )
 
@@ -361,10 +439,40 @@ return operator_func
 
-def operator_func_factory(inner_op: Callable, index: int) -> Callable:
-    def operator_func(self):
-        if not self.hierarchical_diagonalization:
+def hierarchical_diagonalization_func_factory(symbol_name: str) -> Callable:
+    def operator_func(self: "Subsystem"):
+        return Qobj_to_scipy_csc_matrix(self.get_operator_by_name(symbol_name))
+
+    return operator_func
+
+
+def keep_terms_for_subsystem(sym_expr, subsys, substitute_zero=False):
+    if substitute_zero:
+        for var_sym in sym_expr.free_symbols:
+            sym_expr = sym_expr.subs(var_sym, 0)
+        return sym_expr
+    terms = sym_expr.as_ordered_terms()
+    for term in terms:
+        var_indices = [
+            get_trailing_number(sym_var.name) for sym_var in list(term.free_symbols)
+        ]
+        if len(set(var_indices) & set(subsys.dynamic_var_indices)) == 0:
+            sym_expr = sym_expr - term
+    return sym_expr
+
+
+def operator_func_factory(
+    inner_op: Callable, index: int, op_type: Optional[str] = None
+) -> Callable:
+    def operator_func(self, op_type=op_type):
+        prefactor = None
+        if self.ext_basis == "harmonic":
+            if op_type in ["position", "sin", "cos"]:
+                prefactor = self.osc_lengths[index] / (2**0.5)
+            elif op_type == "momentum":
+                prefactor = 1 / (self.osc_lengths[index] * 2**0.5)
+        if prefactor:
+            return self._kron_operator(
+                inner_op(self.cutoffs_dict()[index], prefactor=prefactor), index
+            )
+        else:
             return self._kron_operator(inner_op(self.cutoffs_dict()[index]), index)
-        else:
-            return self.identity_wrap_for_hd(
-                inner_op(self.cutoffs_dict()[index]), index
-            ).data.tocsc()
 
@@ -416,1 +524,327 @@ return operator_func
     return sparse_mat**n
+
+
+def round_symbolic_expr(expr: sm.Expr, number_of_digits: int) -> sm.Expr:
+    rounded_expr = expr.expand()
+    for term in sm.preorder_traversal(expr.expand()):
+        if isinstance(term, sm.Float):
+            rounded_expr = rounded_expr.subs(term, round(term, number_of_digits))
+    return rounded_expr
+
+
+def yaml_like_out_with_pp(circuit_yaml):
+    import pyparsing as pp
+
+    code = circuit_input.remove_comments(circuit_yaml)
+    code = circuit_input.remove_branchline(code)
+    code = circuit_input.strip_empty_lines(code)
+
+    bp = [
+        circuit_input.BRANCHES.parse_string(branch)
+        for branch in [branch for branch in code.splitlines()]
+    ]
+
+    yaml_like_out = []
+    for branch in bp:
+        yaml_branch = []
+        for idx, param in enumerate(branch):
+            if isinstance(param, pp.ParseResults):
+                parse_type = param.getName()
+                if parse_type == "ASSIGN" or parse_type == "AUX_PARAM":
+                    yaml_branch.append(
+                        param[0] + "=" + "".join([str(x) for x in param[1:] if x])
+                    )
+                elif parse_type == "SYMBOL":
+                    yaml_branch.append(str(param[0]))
+                else:
+                    yaml_branch.append("".join([str(x) for x in param if x]))
+            else:
+                yaml_branch.append(param)
+        yaml_like_out.append(yaml_branch)
+    return yaml_like_out
+
+
+def assemble_circuit(
+    circuit_list: List[str],
+    couplers: str,
+    rename_parameters=False,
+) -> Tuple[str, List[Dict[int, int]]]:
+    """
+    Assemble a yaml string for a large circuit that are made of smaller
+    sub-circuits and coupling elements. This method takes a list of Sub-circuit yaml
+    strings as the first argument, and a yaml string that characterize the coupler
+    branches as the second argument. For example, if one wish to make a yaml string for
+    a circuit consist of a grounded fluxonium capacitively coupled to another fluxonium,
+    then one need to define:
+
+    circuit_1 = '''
+    branches:
+    - [C, 0, 1, EC = 1]
+    - [JJ, 0, 1, EJ = 20, ECJ = 3]
+    - [L, 0, 1, EL = 10]
+    '''
+    circuit_2 = '''
+    branches:
+    - [C, 0, 1, EC = 3]
+    - [JJ, 0, 1, EJ = 1, ECJ = 2]
+    - [L, 0, 1, EL = 0.5]
+    '''
+    circuit_list = [circuit_1, circuit_2]
+    couplers = '''
+    branches:
+    - [C, 1: 1, 2: 1, E_coup = 1]
+    '''
+
+    If rename_parameters argument set to False, the resulting yaml string for the assembled
+    composite circuit is:
+    branches:
+    - [C, 0, 1, EC = 1]
+    - [JJ, 0, 1, EJ = 20, ECJ = 3]
+    - [L, 0, 1, EL = 10]
+    - [C, 0, 2, EC]
+    - [JJ, 0, 2, EJ, ECJ]
+    - [L, 0, 2, EL]
+    - [C, 1, 2, E_coup = 1]
+
+    If rename_parameters argument set to True, the resulting yaml string for the assembled
+    composite circuit is:
+    branches:
+    - [C, 0, 1, EC_1 = 1]
+    - [JJ, 0, 1, EJ_1 = 20, ECJ_1 = 3]
+    - [L, 0, 1, EL_1 = 10]
+    - [C, 0, 2, EC_2 = 3]
+    - [JJ, 0, 2, EJ_2 = 1, ECJ_2 = 2]
+    - [L, 0, 2, EL_2 = 0.5]
+    - [C, 1, 2, E_coup_12 = 1]
+
+    The yaml strings for each coupler circuit follow the syntax of input strings used in the
+    custom circuit module, whereas the syntax for coupler branches is different. Each coupler
+    branch is represented by:
+    <branch-type>, <subcircuit-index>:<node-index-of-the-circuit>, <subcircuit-index>:<node-index-of-the-circuit>, <param-1> [, <param-2>]]
+
+    All the grounded sub-circuit share the same ground node in the composite circuit.
+    The parameter symbols are global, i.e. the same parameter symbol that appears in different
+    subcircuit yaml strings will be treated as one parameter in the composite circuit. The
+    symbolic parameters are only initialized once, with the value specified at the first
+    instance of appearance (notice the initial value for EC in the above example).
+
+    Parameters
+    ----------
+    circuit_list:
+        A list of yaml strings encoding branches of sub-circuits.
+    couplers:
+        A yaml string that encodes information of coupler branches
+    rename_parameters:
+        If set to True, parameters in the sub-circuits will be renamed as <original-parameter-name>_<sub-circuit-index>
+        parameters in the couplers will be renamed as <original-parameter-name>_<sub-circuit-1-index><sub-circuit-2-index>
+
+    Returns
+    -------
+        A yaml string for the composite circuit, which can be used as the input for the custom
+        circuit module, and a list of dictionaries which provides the mapping between the
+        node indices of the sub-circuits (keys) and those of the composite circuit (values).
+    """
+    # identify numbers of subcircuits
+    subcircuit_number = len(circuit_list)
+    subcircuit_branches_list = []
+    subcircuit_nodes_list = []
+    subcircuit_is_grounded_list = []
+    subcircuit_node_number_list = []
+    subcircuit_node_index_dict_list = []
+    for circuit_yaml in circuit_list:
+        # load subcircuit yaml strings
+        subcircuit_branches = yaml_like_out_with_pp(circuit_yaml)
+        # append the dictionary for each subcircuit
+        subcircuit_branches_list.append(subcircuit_branches)
+        # for each subcircuit, extract their node indices
+        subcircuit_nodes = [
+            [subcircuit_branch[1], subcircuit_branch[2]]
+            for subcircuit_branch in subcircuit_branches
+        ]
+        subcircuit_nodes = [
+            *unique_elements_in_list(flatten_list_recursive(subcircuit_nodes))
+        ]
+        # add node indices of each subcircuit to a single list
+        subcircuit_nodes_list.append(subcircuit_nodes)
+        # add total node number of each subcircuit to a single list
+        subcircuit_node_number_list.append(len(subcircuit_nodes))
+        subcircuit_is_grounded_list.append(True if (0 in subcircuit_nodes) else False)
+    # generate a dictionary for each subcircuit which has subcircuit node indices as keys and
+    # the assembled circuit node indices as values
+    node_index_offset = 0
+    for subcircuit_index in range(subcircuit_number):
+        subcircuit_node_index_dict = {}
+        for subcircuit_node_index in subcircuit_nodes_list[subcircuit_index]:
+            subcircuit_node_index_dict[subcircuit_node_index] = (
+                subcircuit_node_index + node_index_offset
+            )
+        if subcircuit_is_grounded_list[subcircuit_index]:
+            subcircuit_node_index_dict[0] = 0
+        node_index_offset += (
+            subcircuit_node_number_list[subcircuit_index]
+            - subcircuit_is_grounded_list[subcircuit_index]
+        )
+        subcircuit_node_index_dict_list.append(subcircuit_node_index_dict)
+    # create new yaml string for the composite circuit
+    composite_circuit_yaml = "\nbranches:\n"
+    # initialize parameter dictionary
+    param_dict = {}
+    # write all the subcircuit branch info into the composite circuit yaml,
+    # converting their node indices
+    for subcircuit_index in range(subcircuit_number):
+        for subcircuit_branch in subcircuit_branches_list[subcircuit_index]:
+            composite_circuit_yaml += " - ["
+            # identify branch type
+            branch_type = subcircuit_branch[0]
+            composite_circuit_yaml += branch_type + " ,"
+            # include the converted first node index
+            composite_circuit_yaml += (
+                str(
+                    subcircuit_node_index_dict_list[subcircuit_index][
+                        subcircuit_branch[1]
+                    ]
+                )
+                + " ,"
+            )
+            # include the converted second node index
+            composite_circuit_yaml += (
+                str(
+                    subcircuit_node_index_dict_list[subcircuit_index][
+                        subcircuit_branch[2]
+                    ]
+                )
+                + " ,"
+            )
+            # identify parameter numbers
+            num_params = 2 if "JJ" in branch_type else 1
+            # include parameters
+            for word in subcircuit_branch[3 : 3 + num_params]:
+                if not is_string_float(word):
+                    if not rename_parameters:
+                        if len(word.split("=")) == 2:
+                            param_str, init_val = word.split("=")
+                            param_str, init_val = param_str.strip(), float(
+                                init_val.strip()
+                            )
+                            # if the parameter is already initialized, the subsequent initialization
+                            # is neglected
+                            if param_str in param_dict:
+                                composite_circuit_yaml += str(param_str) + ", "
+                            else:
+                                composite_circuit_yaml += str(word) + ", "
+                                param_dict[param_str] = init_val
+                        elif len(word.split("=")) == 1:
+                            composite_circuit_yaml += str(word) + ", "
+                    else:
+                        if len(word.split("=")) == 2:
+                            param_str, init_val = word.split("=")
+                            param_str, init_val = param_str.strip(), float(
+                                init_val.strip()
+                            )
+                            composite_circuit_yaml += (
+                                param_str
+                                + "_"
+                                + str(subcircuit_index + 1)
+                                + " = "
+                                + str(init_val)
+                                + ", "
+                            )
+                        elif len(word.split("=")) == 1:
+                            composite_circuit_yaml += (
+                                str(word.strip())
+                                + "_"
+                                + str(subcircuit_index + 1)
+                                + ", "
+                            )
+                else:
+                    composite_circuit_yaml += str(word) + ", "
+            composite_circuit_yaml += "]\n"
+    # add coupling branches to the composite circuit yaml string
+    # load coupler yaml strings
+    coupler_branches = yaml_like_out_with_pp(couplers)
+    for coupler_branch in coupler_branches:
+        composite_circuit_yaml += " - ["
+        branch_type = coupler_branch[0]
+        composite_circuit_yaml += branch_type + ", "
+        # include the converted first node index
+        composite_circuit_yaml += (
+            str(
+                subcircuit_node_index_dict_list[list(coupler_branch[1].keys())[0] - 1][
+                    list(coupler_branch[1].values())[0]
+                ]
+            )
+            + " ,"
+        )
+        # include the converted second node index
+        composite_circuit_yaml += (
+            str(
+                subcircuit_node_index_dict_list[list(coupler_branch[2].keys())[0] - 1][
+                    list(coupler_branch[2].values())[0]
+                ]
+            )
+            + " ,"
+        )
+        # identify parameter numbers
+        num_params = 2 if "JJ" in branch_type else 1
+        # include parameters
+        for word in coupler_branch[3 : 3 + num_params]:
+            if not is_string_float(word):
+                if not rename_parameters:
+                    if len(word.split("=")) == 2:
+                        param_str, init_val = word.split("=")
+                        param_str, init_val = param_str.strip(), float(init_val.strip())
+                        # if the parameter is already initialized, the subsequent initialization
+                        # is neglected
+                        if param_str in param_dict:
+                            composite_circuit_yaml += str(param_str) + ", "
+                        else:
+                            composite_circuit_yaml += str(word) + ", "
+                            param_dict[param_str] = init_val
+                    elif len(word.split("=")) == 1:
+                        composite_circuit_yaml += str(word) + ", "
+                else:
+                    if len(word.split("=")) == 2:
+                        param_str, init_val = word.split("=")
+                        param_str, init_val = param_str.strip(), float(init_val.strip())
+                        composite_circuit_yaml += (
+                            param_str
+                            + "_"
+                            + str(list(coupler_branch[1].keys())[0])
+                            + str(list(coupler_branch[2].keys())[0])
+                            + " = "
+                            + str(init_val)
+                            + ", "
+                        )
+                    elif len(word.split("=")) == 1:
+                        composite_circuit_yaml += (
+                            str(word.strip())
+                            + "_"
+                            + str(list(coupler_branch[1].keys())[0])
+                            + str(list(coupler_branch[2].keys())[0])
+                            + ", "
+                        )
+            else:
+                composite_circuit_yaml += str(word) + ", "
+        composite_circuit_yaml += "]\n"
+    return composite_circuit_yaml, subcircuit_node_index_dict_list
+
+
+def assemble_transformation_matrix(
+    transformation_matrix_list: List[ndarray],
+) -> ndarray:
+    """
+    Assemble a transformation matrix for a large circuit that are made of smaller sub-circuits
+    and coupling elements. This method takes a list of sub-circuit transformation matrices as
+    the argument.
+
+    Parameters
+    ----------
+    transformation_matrix_list:
+        A list of transformation matrices as numpy ndarray.
+
+    Returns
+    -------
+        A numpy ndarray for the composite circuit.
+    """
+
+    return sp.linalg.block_diag(*transformation_matrix_list)

scqubits/core/cos2phi_qubit.py

@@ -144,3 +144,3 @@ # cos2phi_qubit.py
 
-        def spectral_density1(omega):
+        def spectral_density1(omega, T):
             therm_ratio = calc_therm_ratio(omega, T)
@@ -162,3 +162,3 @@ s = (
 
-        def spectral_density2(omega):
+        def spectral_density2(omega, T):
             therm_ratio = calc_therm_ratio(omega, T)
@@ -265,3 +265,3 @@ s = (
 
-        def spectral_density1(omega):
+        def spectral_density1(omega, T):
             therm_ratio = calc_therm_ratio(omega, T)
@@ -282,3 +282,3 @@ s1 = (
 
-        def spectral_density2(omega):
+        def spectral_density2(omega, T):
             therm_ratio = calc_therm_ratio(omega, T)
@@ -385,3 +385,3 @@ s2 = (
 
-        def spectral_density(omega):
+        def spectral_density(omega, T):
             therm_ratio = calc_therm_ratio(omega, T)
@@ -471,2 +471,3 @@ s = (
     """
+
     EJ = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
@@ -473,0 +474,0 @@ ECJ = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")

scqubits/core/descriptors.py

@@ -55,3 +55,3 @@ # descriptors.py
     inner_object_name:
-        Used, e.g., in FulLZeroPi where an inner-object property is to be set.
+        Used, e.g., in FullZeroPi where an inner-object property is to be set.
     attr_name:
@@ -68,2 +68,4 @@ custom attribute name to be used (default: name from defining property in
         attr_name: Optional[str] = None,
+        fget=None,
+        fset=None,
     ) -> None:
@@ -73,2 +75,4 @@ self.event = event
         self.attr_name = attr_name
+        self.setter = fset
+        self.getter = fget
 
@@ -87,3 +91,8 @@ def __set_name__(self, owner, name: str) -> None:
             return getattr(inner_instance, self.attr_name)
-        return instance.__dict__[self.attr_name]
+        if self.getter is None:
+            return instance.__dict__[
+                self.attr_name
+            ]  # cannot use getattr, otherwise recursion
+        else:
+            return self.getter(instance)
 
@@ -97,7 +106,16 @@ def __set__(self, instance: DispatchClient, value: TargetType) -> None:
             assert self.attr_name
-            if self.attr_name not in instance.__dict__:
-                instance.__dict__[self.attr_name] = value
+            if (
+                self.attr_name not in instance.__dict__
+                and f"_{self.attr_name}" not in instance.__dict__
+            ):
+                if self.setter is None:
+                    instance.__dict__[self.attr_name] = value
+                else:
+                    self.setter(instance, value, name=self.attr_name)
                 # Rely on inner_instance.attr_name to do the broadcasting.
             else:
-                instance.__dict__[self.attr_name] = value
+                if self.setter is None:
+                    instance.__dict__[self.attr_name] = value
+                else:
+                    self.setter(instance, value, name=self.attr_name)
                 instance.broadcast(self.event)

scqubits/core/fluxonium.py

@@ -72,2 +72,3 @@ # fluxonium.py
     """
+
     EJ = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
@@ -74,0 +75,0 @@ EC = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")

scqubits/core/hilbert_space.py

@@ -37,3 +37,2 @@ # hilbert_space.py
 from numpy import ndarray
-from qutip import Qobj
 from scipy.sparse import csc_matrix, dia_matrix
@@ -69,2 +68,3 @@
 from scqubits.utils.typedefs import OscillatorList, QuantumSys, QubitList
+from scqubits.core.qubit_base import QubitBaseClass
 
@@ -139,3 +139,3 @@
         bare_esys: Optional[Dict[int, ndarray]] = None,
-    ) -> Qobj:
+    ) -> qt.Qobj:
         """
@@ -158,3 +158,3 @@ Returns the full Hamiltonian of the interacting quantum system described by the
         """
-        hamiltonian = cast(Qobj, self.g_strength)
+        hamiltonian = cast(qt.Qobj, self.g_strength)
         id_wrapped_ops = self.id_wrap_all_ops(
@@ -174,3 +174,3 @@ self.operator_list, subsystem_list, bare_esys=bare_esys
         bare_esys: Optional[Dict[int, ndarray]] = None,
-    ) -> List[Qobj]:
+    ) -> List[qt.Qobj]:
         """
@@ -260,3 +260,4 @@ Returns a list of identity-wrapped operators, one for each operator in
         operator_list: List[Tuple[int, str, Union[ndarray, csc_matrix, dia_matrix]]],
-        const: Optional[Dict[str, Union[float, complex]]] = None,
+        id_wrapped_operator_list: Optional[List[Tuple[str, callable]]] = None,
+        const: Optional[Dict[str, Union[float, complex, QubitBaseClass]]] = None,
         add_hc: bool = False,
@@ -276,2 +277,3 @@ ) -> None:
         self.operator_list = operator_list
+        self.id_wrapped_operator_list = id_wrapped_operator_list or []
         self.const = const or {}
@@ -308,6 +310,6 @@ self.add_hc = add_hc
     def run_string_code(
-        self, expression: str, idwrapped_ops_by_name: Dict[str, Qobj]
-    ) -> Qobj:
+        self, expression: str, idwrapped_ops_by_name: Dict[str, qt.Qobj]
+    ) -> qt.Qobj:
         expression = self.parse_qutip_functions(expression)
-        idwrapped_ops_by_name["Qobj"] = Qobj
+        idwrapped_ops_by_name["Qobj"] = qt.Qobj
 
@@ -324,3 +326,3 @@ main = importlib.import_module("__main__")
         bare_esys: Optional[Dict[int, ndarray]] = None,
-    ) -> Dict[str, Qobj]:
+    ) -> Dict[str, qt.Qobj]:
         idwrapped_ops_by_name = {}
@@ -341,3 +343,3 @@ for subsys_index, name, op in self.operator_list:
         bare_esys: Optional[Dict[int, ndarray]] = None,
-    ) -> Qobj:
+    ) -> qt.Qobj:
         """
@@ -356,2 +358,8 @@ Parameters
         )
+        idwrapped_ops_by_name.update(
+            {
+                item[0]: item[1](bare_esys=bare_esys)
+                for item in self.id_wrapped_operator_list
+            }
+        )
         hamiltonian = self.run_string_code(self.expr, idwrapped_ops_by_name)
@@ -428,6 +436,9 @@ if not self.add_hc:
         self._osc_subsys_list = [
-            subsys for subsys in self if isinstance(subsys, osc.Oscillator)
+            subsys
+            for subsys in self
+            if isinstance(subsys, osc.Oscillator)
+            or (hasattr(subsys, "is_purely_harmonic") and subsys.is_purely_harmonic)
         ]
         self._qbt_subsys_list = [
-            subsys for subsys in self if not isinstance(subsys, osc.Oscillator)
+            subsys for subsys in self if subsys not in self._osc_subsys_list
         ]
@@ -454,4 +465,3 @@
     @overload
-    def __getitem__(self, key: int) -> QuantumSys:
-        ...
+    def __getitem__(self, key: int) -> QuantumSys: ...
 
@@ -461,4 +471,3 @@ @overload
         self, key: str
-    ) -> Union[QuantumSys, InteractionTerm, InteractionTermStr]:
-        ...
+    ) -> Union[QuantumSys, InteractionTerm, InteractionTermStr]: ...
 
@@ -549,3 +558,3 @@ def __getitem__(
         """
-        init_parameters = self._init_params
+        init_parameters = self._init_params.copy()
         init_parameters.remove("ignore_low_overlap")
@@ -659,2 +668,10 @@ init_parameters.append("_ignore_low_overlap")
         for subsys_index, subsys in enumerate(self):
+            # generate bare_esys for the subsystem as well if necessary
+            if (
+                hasattr(subsys, "hierarchical_diagonalization")
+                and subsys.hierarchical_diagonalization
+            ):
+                subsys.hilbert_space.generate_bare_esys(
+                    update_subsystem_indices=subsys.affected_subsystem_indices
+                )
             # diagonalizing only those subsystems present in update_subsystem_indices
@@ -683,3 +700,2 @@ if subsys_index in update_subsystem_indices:
         )
-
         return bare_esys_dict
@@ -801,3 +817,3 @@
         bare_esys: Optional[Dict[int, ndarray]] = None,
-    ) -> Qobj:
+    ) -> qt.Qobj:
         """
@@ -819,3 +835,5 @@ Parameters
 
-    def bare_hamiltonian(self, bare_esys: Optional[Dict[int, ndarray]] = None) -> Qobj:
+    def bare_hamiltonian(
+        self, bare_esys: Optional[Dict[int, ndarray]] = None
+    ) -> qt.Qobj:
         """
@@ -833,3 +851,7 @@ Parameters
         """
-        bare_hamiltonian = Qobj(0)
+        bare_hamiltonian = (
+            qt.Qobj(0, dims=[self.subsystem_dims] * 2)
+            if qt.__version__ >= "5.0.0"
+            else qt.Qobj(0)
+        )
         for subsys_index, subsys in enumerate(self):
@@ -845,3 +867,3 @@ if bare_esys is not None and subsys_index in bare_esys:
         self, bare_esys: Optional[Dict[int, ndarray]] = None
-    ) -> Qobj:
+    ) -> qt.Qobj:
         """
@@ -862,7 +884,11 @@ Returns the interaction Hamiltonian, based on the interaction terms specified
         if not self.interaction_list:
-            return Qobj(0)
+            return (
+                qt.Qobj(0, dims=[self.subsystem_dims] * 2)
+                if qt.__version__ >= "5.0.0"
+                else qt.Qobj(0)
+            )
 
         operator_list = []
         for term in self.interaction_list:
-            if isinstance(term, Qobj):
+            if isinstance(term, qt.Qobj):
                 operator_list.append(term)
@@ -881,3 +907,3 @@ elif isinstance(term, (InteractionTerm, InteractionTermStr)):
 
-    def diag_hamiltonian(self, subsystem: QuantumSys, evals: ndarray = None) -> Qobj:
+    def diag_hamiltonian(self, subsystem: QuantumSys, evals: ndarray = None) -> qt.Qobj:
         """Returns a `qutip.Qobj` which has the eigenenergies of the object `subsystem`
@@ -897,3 +923,3 @@ on the diagonal.
             evals = subsystem.eigenvals(evals_count=evals_count)
-        diag_qt_op = qt.Qobj(inpt=np.diagflat(evals[0:evals_count]))  # type:ignore
+        diag_qt_op = qt.Qobj(np.diagflat(evals[0:evals_count]))  # type:ignore
         return spec_utils.identity_wrap(diag_qt_op, subsystem, self.subsystem_list)
@@ -905,3 +931,3 @@
 
-    def diag_operator(self, diag_elements: ndarray, subsystem: QuantumSys) -> Qobj:
+    def diag_operator(self, diag_elements: ndarray, subsystem: QuantumSys) -> qt.Qobj:
         """For given diagonal elements of a diagonal operator in `subsystem`, return
@@ -924,3 +950,3 @@ the `Qobj` operator for the full Hilbert space (perform wrapping in
 
-    def hubbard_operator(self, j: int, k: int, subsystem: QuantumSys) -> Qobj:
+    def hubbard_operator(self, j: int, k: int, subsystem: QuantumSys) -> qt.Qobj:
         """Hubbard operator :math:`|j\\rangle\\langle k|` for system `subsystem`
@@ -939,3 +965,3 @@
 
-    def annihilate(self, subsystem: QuantumSys) -> Qobj:
+    def annihilate(self, subsystem: QuantumSys) -> qt.Qobj:
         """Annihilation operator a for `subsystem`
@@ -1054,3 +1080,4 @@
         for idx, evec in enumerate(self._data["evecs"][0]):
-            phase = spec_utils.extract_phase(evec.data.toarray())
+            array = utils.Qobj_to_scipy_csc_matrix(evec)
+            phase = spec_utils.extract_phase(array)
             self._data["evecs"][0][idx] = evec * np.exp(-1j * phase)
@@ -1065,3 +1092,3 @@
         **kwargs,
-    ) -> Qobj:
+    ) -> qt.Qobj:
         """
@@ -1113,4 +1140,6 @@ Express a subsystem operator in the dressed eigenbasis of the full system
         dressed_evecs = self._data["evecs"][0]
-        dressed_op_data = id_wrapped_op.transform(dressed_evecs).data.toarray()
-        dressed_op_truncated = Qobj(
+        dressed_op_data = utils.Qobj_to_scipy_csc_matrix(
+            id_wrapped_op.transform(dressed_evecs)
+        )
+        dressed_op_truncated = qt.Qobj(
             dressed_op_data[0:truncated_dim, 0:truncated_dim],
@@ -1179,5 +1208,5 @@ dims=[[truncated_dim], [truncated_dim]],
         if "expr" in kwargs:
-            interaction: Union[
-                InteractionTerm, InteractionTermStr
-            ] = self._parse_interactiontermstr(**kwargs)
+            interaction: Union[InteractionTerm, InteractionTermStr] = (
+                self._parse_interactiontermstr(**kwargs)
+            )
         elif "qobj" in kwargs:
@@ -1214,8 +1243,20 @@ interaction = self._parse_qobj(**kwargs)
         operator_list = []
+        id_wrapped_operator_list = []
         for key in kwargs.keys():
+            if callable(kwargs[key][1]) and not hasattr(kwargs[key][1], "__self__"):
+                id_wrapped_operator_list.append(kwargs[key])
+                continue
             if re.match(r"op\d+$", key) is None:
                 raise TypeError("Unexpected keyword argument {}.".format(key))
             operator_list.append(self._parse_str_based_op(kwargs[key]))
+        if id_wrapped_operator_list == []:
+            id_wrapped_operator_list = None
 
-        return InteractionTermStr(expr, operator_list, const=const, add_hc=add_hc)
+        return InteractionTermStr(
+            expr,
+            operator_list,
+            id_wrapped_operator_list=id_wrapped_operator_list,
+            const=const,
+            add_hc=add_hc,
+        )
 
@@ -1238,5 +1279,5 @@ def _parse_interactionterm(self, **kwargs) -> InteractionTerm:
     @staticmethod
-    def _parse_qobj(**kwargs) -> Qobj:
+    def _parse_qobj(**kwargs) -> qt.Qobj:
         op = kwargs["qobj"]
-        if len(kwargs) > 1 or not isinstance(op, Qobj):
+        if len(kwargs) > 1 or not isinstance(op, qt.Qobj):
             raise TypeError("Cannot interpret specified operator {}".format(op))
@@ -1243,0 +1284,0 @@ return kwargs["qobj"]

scqubits/core/noise.py

@@ -23,2 +23,3 @@ # noise.py
 import scipy as sp
+import qutip as qt
 
@@ -30,2 +31,3 @@ from matplotlib.axes import Axes
 from scipy.sparse import csc_matrix
+from sympy import csc
 
@@ -35,2 +37,3 @@ import scqubits.core.units as units
 import scqubits.utils.plotting as plotting
+from scqubits.utils.misc import Qobj_to_scipy_csc_matrix
 
@@ -1146,2 +1149,3 @@ from scqubits.core.storage import SpectrumData
         spectral_density: Callable,
+        T: float = NOISE_PARAMS["T"],
         total: bool = True,
@@ -1175,4 +1179,6 @@ esys: Tuple[ndarray, ndarray] = None,
             noise operator
+        T:
+            Temperature defined in Kelvin
         spectral_density:
-            defines a spectral density, must take one argument: `omega`
+            defines a spectral density, must take two arguments: `omega` and `T`
             (assumed to be in units of `2 \pi * <system units>`)
@@ -1225,5 +1231,5 @@ total:
         s = (
-            spectral_density(omega) + spectral_density(-omega)
+            spectral_density(omega, T) + spectral_density(-omega, T)
             if total
-            else spectral_density(omega)
+            else spectral_density(omega, T)
         )
@@ -1253,2 +1259,4 @@
         get_rate: bool = False,
+        noise_op: Optional[Union[ndarray, csc_matrix, qt.Qobj]] = None,
+        branch_params: Optional[dict] = None,
     ) -> float:
@@ -1293,3 +1301,3 @@ r"""
             # See Smith et al (2020)
-            def q_cap_fun(omega):
+            def q_cap_fun(omega, T):
                 return (
@@ -1304,6 +1312,6 @@ 1e6
 
-            def q_cap_fun(omega):
+            def q_cap_fun(omega, T):
                 return Q_cap
 
-        def spectral_density(omega):
+        def spectral_density(omega, T):
             therm_ratio = calc_therm_ratio(omega, T)
@@ -1313,4 +1321,4 @@ s = (
                 * 8
-                * self.EC
-                / q_cap_fun(omega)
+                * (branch_params if branch_params else self.EC)
+                / q_cap_fun(omega, T)
                 * (1 / np.tanh(0.5 * np.abs(therm_ratio)))
@@ -1324,3 +1332,9 @@ / (1 + np.exp(-therm_ratio))
 
-        noise_op = self.n_operator()  # type: ignore
+        noise_op = noise_op or self.n_operator()  # type: ignore
+        if not isinstance(noise_op, (ndarray, csc_matrix, qt.Qobj)):
+            raise AttributeError(
+                "The type of the matrix noise_op is invalid. It should be an instance of ndarray, csc_matrix or qutip Qobj."
+            )
+        if isinstance(noise_op, (qt.Qobj)):
+            noise_op = Qobj_to_scipy_csc_matrix(noise_op)
 
@@ -1331,2 +1345,3 @@ return self.t1(
             noise_op=noise_op,
+            T=T,
             spectral_density=spectral_density,
@@ -1347,2 +1362,3 @@ total=total,
         get_rate: bool = False,
+        noise_op: Optional[Union[ndarray, csc_matrix, qt.Qobj]] = None,
     ) -> float:
@@ -1383,3 +1399,3 @@ r"""Noise due to charge coupling to an impedance (such as a transmission line).
 
-        def spectral_density(omega):
+        def spectral_density(omega, T):
             # Note, our definition of Q_c is different from Zhang et al (2020) by a
@@ -1398,3 +1414,9 @@ # factor of 2
 
-        noise_op = self.n_operator()  # type: ignore
+        noise_op = noise_op or self.n_operator()  # type: ignore
+        if not isinstance(noise_op, (ndarray, csc_matrix, qt.Qobj)):
+            raise AttributeError(
+                "The type of the matrix noise_op is invalid. It should be an instance of ndarray, csc_matrix or qutip Qobj."
+            )
+        if isinstance(noise_op, (qt.Qobj)):
+            noise_op = Qobj_to_scipy_csc_matrix(noise_op)
 
@@ -1405,2 +1427,3 @@ return self.t1(
             noise_op=noise_op,
+            T=T,
             spectral_density=spectral_density,
@@ -1422,2 +1445,3 @@ total=total,
         get_rate: bool = False,
+        noise_op_method: Optional[Callable] = None,
     ) -> float:
@@ -1462,3 +1486,3 @@ r"""Noise due to a bias flux line.
 
-        def spectral_density(omega, Z=Z):
+        def spectral_density(omega, T, Z=Z):
             """
@@ -1484,4 +1508,5 @@ Our definitions assume that the noise_op is dH/dflux.
 
-        noise_op = self.d_hamiltonian_d_flux()  # type: ignore
-
+        noise_op = (noise_op_method or self.d_hamiltonian_d_flux)()  # type: ignore
+        if isinstance(noise_op, qt.Qobj):
+            noise_op = Qobj_to_scipy_csc_matrix(noise_op)
         return self.t1(
@@ -1491,2 +1516,3 @@ i=i,
             noise_op=noise_op,
+            T=T,
             spectral_density=spectral_density,
@@ -1507,2 +1533,4 @@ total=total,
         get_rate: bool = False,
+        noise_op: Optional[Union[ndarray, csc_matrix, qt.Qobj]] = None,
+        branch_params: Optional[dict] = None,
     ) -> float:
@@ -1547,3 +1575,3 @@ r"""
             # See Smith et al (2020)
-            def q_ind_fun(omega):
+            def q_ind_fun(omega, T):
                 therm_ratio = abs(calc_therm_ratio(omega, T))
@@ -1570,11 +1598,11 @@ therm_ratio_500MHz = calc_therm_ratio(
 
-            def q_ind_fun(omega):
+            def q_ind_fun(omega, T):
                 return Q_ind
 
-        def spectral_density(omega):
+        def spectral_density(omega, T):
             therm_ratio = calc_therm_ratio(omega, T)
             s = (
                 2
-                * self.EL
-                / q_ind_fun(omega)
+                * (branch_params if branch_params else self.EL)
+                / q_ind_fun(omega, T)
                 * (1 / np.tanh(0.5 * np.abs(therm_ratio)))
@@ -1588,3 +1616,9 @@ / (1 + np.exp(-therm_ratio))
 
-        noise_op = self.phi_operator()  # type: ignore
+        noise_op = noise_op or self.phi_operator()  # type: ignore
+        if not isinstance(noise_op, (ndarray, csc_matrix, qt.Qobj)):
+            raise AttributeError(
+                "The type of the matrix noise_op is invalid. It should be an instance of ndarray, csc_matrix or qutip Qobj."
+            )
+        if isinstance(noise_op, (qt.Qobj)):
+            noise_op = Qobj_to_scipy_csc_matrix(noise_op)
 
@@ -1595,2 +1629,3 @@ return self.t1(
             noise_op=noise_op,
+            T=T,
             spectral_density=spectral_density,
@@ -1613,2 +1648,3 @@ total=total,
         get_rate: bool = False,
+        noise_op: Optional[Union[ndarray, csc_matrix, qt.Qobj]] = None,
     ) -> float:
@@ -1655,3 +1691,3 @@ r"""Noise due to quasiparticle tunneling across a Josephson junction.
 
-            def y_qp_fun(omega):
+            def y_qp_fun(omega, T):
                 """
@@ -1692,6 +1728,6 @@ Based on Eq. S23 in the appendix of Smith et al (2020).
 
-            def y_qp_fun(omega):
+            def y_qp_fun(omega, T):
                 return Y_qp
 
-        def spectral_density(omega):
+        def spectral_density(omega, T):
             """Based on Eq. 19 in Smith et al (2020)."""
@@ -1703,3 +1739,3 @@ therm_ratio = calc_therm_ratio(omega, T)
                 * omega
-                * complex(y_qp_fun(omega)).real
+                * complex(y_qp_fun(omega, T)).real
                 * (1 / np.tanh(0.5 * therm_ratio))
@@ -1713,4 +1749,11 @@ / (1 + np.exp(-therm_ratio))
         # transform the operator using phi -> phi + 2*pi*flux
-        noise_op = self.sin_phi_operator(alpha=0.5, beta=0.5 * (2 * np.pi * self.flux))  # type: ignore
+        noise_op = noise_op or self.sin_phi_operator(alpha=0.5, beta=0.5 * (2 * np.pi * self.flux))  # type: ignore
 
+        if not isinstance(noise_op, (ndarray, csc_matrix, qt.Qobj)):
+            raise AttributeError(
+                "The type of the matrix noise_op is invalid. It should be an instance of ndarray, csc_matrix or qutip Qobj."
+            )
+        if isinstance(noise_op, (qt.Qobj)):
+            noise_op = Qobj_to_scipy_csc_matrix(noise_op)
+
         return self.t1(
@@ -1720,2 +1763,3 @@ i=i,
             noise_op=noise_op,
+            T=T,
             spectral_density=spectral_density,
@@ -1722,0 +1766,0 @@ total=total,

scqubits/core/operators.py

@@ -169,4 +169,46 @@ # operators.py
 
-def ia_minus_iadag_sparse(
+def cos_theta_harmonic(
     dimension: int, prefactor: Union[float, complex, None] = None
+) -> ndarray:
+    """Operator matrix for cos(prefactor(a+a^dag)) of size dimension x dimension in
+    sparse matrix representation.
+
+    Parameters
+    ----------
+    dimension:
+        matrix size
+    prefactor:
+        prefactor multiplying the number operator matrix
+        (if not given, this defaults to 1)
+
+    Returns
+    -------
+        prefactor (a + a^dag) as ndarray, size dimension x dimension
+    """
+    return sp.linalg.cosm(a_plus_adag_sparse(dimension, prefactor=prefactor).toarray())
+
+
+def sin_theta_harmonic(
+    dimension: int, prefactor: Union[float, complex, None] = None
+) -> ndarray:
+    """Operator matrix for sin(prefactor(a+a^dag)) of size dimension x dimension in
+    sparse matrix representation.
+
+    Parameters
+    ----------
+    dimension:
+        matrix size
+    prefactor:
+        prefactor multiplying the number operator matrix
+        (if not given, this defaults to 1)
+
+    Returns
+    -------
+        prefactor (a + a^dag) as ndarray, size dimension x dimension
+    """
+    return sp.linalg.sinm(a_plus_adag_sparse(dimension, prefactor=prefactor).toarray())
+
+
+def iadag_minus_ia_sparse(
+    dimension: int, prefactor: Union[float, complex, None] = None
 ) -> csc_matrix:
@@ -190,7 +232,7 @@ """Operator matrix for prefactor(ia-ia^dag) of size dimension x dimension as
     return prefactor * (
-        1j * annihilation_sparse(dimension) - 1j * creation_sparse(dimension)
+        1j * creation_sparse(dimension) - 1j * annihilation_sparse(dimension)
     )
 
 
-def ia_minus_iadag(
+def iadag_minus_ia(
     dimension: int, prefactor: Union[float, complex, None] = None
@@ -213,3 +255,3 @@ ) -> ndarray:
     """
-    return ia_minus_iadag_sparse(dimension, prefactor=prefactor).toarray()
+    return iadag_minus_ia_sparse(dimension, prefactor=prefactor).toarray()
 
@@ -216,0 +258,0 @@

scqubits/core/oscillator.py

@@ -92,2 +92,3 @@ # oscillator.py
     """
+
     E_osc = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
@@ -211,2 +212,3 @@ l_osc = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
     """
+
     K = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
@@ -252,5 +254,3 @@
         """
-        evals = [
-            (self.E_osc + self.K) * n - self.K * n**2 for n in range(evals_count)
-        ]
+        evals = [(self.E_osc + self.K) * n - self.K * n**2 for n in range(evals_count)]
         return np.asarray(evals)

scqubits/core/param_sweep.py

@@ -14,2 +14,3 @@ # param_sweep.py
 import copy
+import dill
 import functools
@@ -42,2 +43,3 @@ import inspect
 
+import scqubits as scq
 import scqubits.core.central_dispatch as dispatch
@@ -494,4 +496,3 @@ import scqubits.core.descriptors as descriptors
         param_indices: Optional[NpIndices] = None,
-    ) -> SpectrumData:
-        ...
+    ) -> SpectrumData: ...
 
@@ -509,4 +510,3 @@ @overload
         param_indices: Optional[NpIndices] = None,
-    ) -> Tuple[List[Tuple[StateLabel, StateLabel]], List[NamedSlotsNdarray]]:
-        ...
+    ) -> Tuple[List[Tuple[StateLabel, StateLabel]], List[NamedSlotsNdarray]]: ...
 
@@ -994,2 +994,4 @@ def transitions(
         for interactionterm in self._hilbertspace.interaction_list:
+            if isinstance(interactionterm, (ndarray, Qobj, csc_matrix)):
+                return True
             for idx_operator in interactionterm.operator_list:
@@ -1035,4 +1037,5 @@ if isinstance(idx_operator[1], (ndarray, Qobj, csc_matrix)):
         if self._deepcopy:
-            stored_hilbertspace = copy.deepcopy(self.hilbertspace)
-            self._hilbertspace = copy.deepcopy(self.hilbertspace)
+            instance_str = dill.dumps(self.hilbertspace)
+            stored_hilbertspace = dill.loads(instance_str)
+            self._hilbertspace = dill.loads(instance_str)
         else:
@@ -1042,3 +1045,7 @@ self.cause_dispatch()
 
-        self._data["bare_evals"], self._data["bare_evecs"] = self._bare_spectrum_sweep()
+        (
+            self._data["bare_evals"],
+            self._data["bare_evecs"],
+            self._data["circuit_esys"],
+        ) = self._bare_spectrum_sweep()
         if not self._bare_only:
@@ -1073,13 +1080,30 @@ self._data["evals"], self._data["evecs"] = self._dressed_spectrum_sweep()
         bare_evecs = np.empty((self.subsystem_count,), dtype=object)
+        # creating data arrays for subsystems, to store the esys for all subsystems when HD is used
+        circuit_esys = np.empty((self.subsystem_count,), dtype=object)
 
         for subsys_index, subsystem in enumerate(self.hilbertspace):
             bare_esys = self._subsys_bare_spectrum_sweep(subsystem)
+            if (
+                hasattr(subsystem, "hierarchical_diagonalization")
+                and subsystem.hierarchical_diagonalization
+            ):
+                evals = np.empty_like(bare_esys[..., 0])
+                evecs = np.empty_like(bare_esys[..., 0])
+                for array_index, esys in np.ndenumerate(bare_esys[..., 0]):
+                    evals[array_index] = esys[0]
+                    evecs[array_index] = esys[1]
+            else:
+                evals = bare_esys[..., 0]
+                evecs = bare_esys[..., 1]
             bare_evals[subsys_index] = NamedSlotsNdarray(
-                np.asarray(bare_esys[..., 0].tolist()),
+                np.asarray(evals.tolist()),
                 self._parameters.paramvals_by_name,
             )
             bare_evecs[subsys_index] = NamedSlotsNdarray(
-                np.asarray(bare_esys[..., 1].tolist()),
+                np.asarray(evecs.tolist()),
                 self._parameters.paramvals_by_name,
             )
+            circuit_esys[subsys_index] = (
+                bare_esys  # when param =(p0, p1, p2, ...), subsys i esys is circuit_esys[i][p0, p1, p3, ...]
+            )
 
@@ -1089,2 +1113,5 @@ return (
             NamedSlotsNdarray(bare_evecs, {"subsys": np.arange(self.subsystem_count)}),
+            NamedSlotsNdarray(
+                circuit_esys, {"subsys": np.arange(self.subsystem_count)}
+            ),
         )
@@ -1099,2 +1126,5 @@
         update_func(self, *paramval_tuple)
+        # use the Circuit method to return esys for all the subsystems when HD is used
+        if isinstance(subsystem, (scq.Circuit, scq.core.circuit.Subsystem)):
+            return subsystem.generate_bare_eigensys()
         evals, evecs = subsystem.eigensys(evals_count=subsystem.truncated_dim)
@@ -1178,3 +1208,2 @@ esys_array = np.empty(shape=(2,), dtype=object)
         update_func(self, *paramval_tuple)
-
         assert self._data is not None
@@ -1188,2 +1217,15 @@ bare_esys: Dict[int, List[ndarray]] = {
         }
+        # update the lookuptables for subsystems using hierarchical diagonalization
+        for subsys_index, subsys in enumerate(hilbertspace.subsystem_list):
+            if (
+                hasattr(subsys, "hierarchical_diagonalization")
+                and subsys.hierarchical_diagonalization
+            ):
+                subsys.set_bare_eigensys(
+                    self._data["circuit_esys"][subsys_index][paramindex_tuple]
+                )
+        if hasattr(
+            hilbertspace.subsystem_list[0], "parent"
+        ):  # update necessary interactions and attributes
+            hilbertspace.subsystem_list[0].parent.update(calculate_bare_esys=False)
 
@@ -1190,0 +1232,0 @@ evals, evecs = hilbertspace.eigensys(

scqubits/core/qubit_base.py

@@ -324,4 +324,3 @@ # qubit_base.py
         return_spectrumdata: "Literal[False]" = False,
-    ) -> ndarray:
-        ...
+    ) -> ndarray: ...
 
@@ -334,4 +333,3 @@ @overload
         return_spectrumdata: "Literal[True]",
-    ) -> SpectrumData:
-        ...
+    ) -> SpectrumData: ...
 
@@ -389,4 +387,3 @@ def eigenvals(
         return_spectrumdata: "Literal[False]" = False,
-    ) -> Tuple[ndarray, ndarray]:
-        ...
+    ) -> Tuple[ndarray, ndarray]: ...
 
@@ -399,4 +396,3 @@ @overload
         return_spectrumdata: "Literal[True]",
-    ) -> SpectrumData:
-        ...
+    ) -> SpectrumData: ...
 
@@ -533,4 +529,3 @@ def eigensys(
         return_datastore: "Literal[False]" = False,
-    ) -> ndarray:
-        ...
+    ) -> ndarray: ...
 
@@ -545,4 +540,3 @@ @overload
         return_datastore: "Literal[True]",
-    ) -> DataStore:
-        ...
+    ) -> DataStore: ...
 
@@ -596,2 +590,4 @@ def matrixelement_table(
         setattr(self, param_name, paramval)
+        if hasattr(self, "hierarchical_diagonalization"):
+            self.update()
         return self.eigensys(evals_count=evals_count)
@@ -603,2 +599,4 @@
         setattr(self, param_name, paramval)
+        if hasattr(self, "hierarchical_diagonalization"):
+            self.update()
         return self.eigenvals(evals_count)
@@ -605,0 +603,0 @@

scqubits/core/spec_lookup.py

@@ -50,8 +50,6 @@ # spec_lookup.py
 
-    def __getitem__(self, key: Any) -> Any:
-        ...
+    def __getitem__(self, key: Any) -> Any: ...
 
     @property
-    def hilbertspace(self) -> "HilbertSpace":
-        ...
+    def hilbertspace(self) -> "HilbertSpace": ...
 
@@ -58,0 +56,0 @@

scqubits/core/transmon.py

@@ -77,2 +77,3 @@ # transmon.py
     """
+
     EJ = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
@@ -643,2 +644,3 @@ EC = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
     """
+
     EJmax = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
@@ -686,4 +688,3 @@ d = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
         return self.EJmax * np.sqrt(
-            np.cos(np.pi * self.flux) ** 2
-            + self.d**2 * np.sin(np.pi * self.flux) ** 2
+            np.cos(np.pi * self.flux) ** 2 + self.d**2 * np.sin(np.pi * self.flux) ** 2
         )
@@ -690,0 +691,0 @@

scqubits/core/zeropi_full.py

@@ -113,2 +113,3 @@ # zeropi_full.py
     """
+
     EJ = descriptors.WatchedProperty(
@@ -319,3 +320,3 @@ float, "QUANTUMSYSTEM_UPDATE", inner_object_name="_zeropi"
     def esys_method_options(self, value: Union[dict, None] = None) -> None:
-        self._zeropi.esys_merthod_options = value
+        self._zeropi.esys_method_options = value
 
@@ -336,3 +337,3 @@ @property
     def evals_method_options(self, value: Union[dict, None] = None) -> None:
-        self._zeropi.evals_merthod_options = value
+        self._zeropi.evals_method_options = value
 
@@ -339,0 +340,0 @@ def hamiltonian(

scqubits/core/zeropi.py

@@ -110,2 +110,3 @@ # zeropi.py
    """
+
     EJ = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")
@@ -112,0 +113,0 @@ EL = descriptors.WatchedProperty(float, "QUANTUMSYSTEM_UPDATE")

scqubits/io_utils/fileio_qutip.py

@@ -36,3 +36,7 @@ # fileio_qutip.py
             [
-                qt.Qobj(inpt=evec, dims=qobj_dims, shape=qobj_shape, type="ket")
+                (
+                    qt.Qobj(evec, type="ket")
+                    if qt.__version__ >= "5.0.0"
+                    else qt.Qobj(evec, dims=qobj_dims, shape=qobj_shape, type="ket")
+                )
                 for evec in evec_array
@@ -39,0 +43,0 @@ ],

scqubits/io_utils/fileio_serializers.py

@@ -360,3 +360,8 @@ # fileio_serializers.py
 def ndarray_deserialize(iodata: "IOData") -> ndarray:
-    return np.asarray(list_deserialize(iodata), dtype=object)
+    # changing instead of np.asarray(, dtype=object), try a = [np.zeros(13), np.zeros((13, 13))]; np.asarray(a, dtype=object)
+    data_list = list_deserialize(iodata)
+    data_array = np.empty(len(data_list), dtype=object)
+    for idx, arr in enumerate(data_list):
+        data_array[idx] = arr
+    return data_array
 
@@ -363,0 +368,0 @@

scqubits/settings.py

@@ -16,2 +16,3 @@ """
 """
+
 # settings.py
@@ -18,0 +19,0 @@ #

scqubits/tests/test_circuit_plot.py

@@ -10,5 +10,3 @@ import os
 
-circ = scq.Circuit(
-    DATADIR + "circuit_zeropi.yaml", ext_basis="harmonic", initiate_sym_calc=True
-)
+circ = scq.Circuit(DATADIR + "circuit_zeropi.yaml", ext_basis="harmonic")
 system_hierarchy = [[1, 3], [2]]
@@ -29,13 +27,16 @@
 )
+circ.update()
 
 
-def test_plot_wf():
-    circ.plot_wavefunction(which=0, var_indices=(2, 3))
+class TestCircuitPlot:
+    @staticmethod
+    def test_plot_wf():
+        circ.plot_wavefunction(which=0, var_indices=(2, 3))
 
-
-def test_plot_potential():
-    circ.plot_potential(
-        θ1=circ._default_grid_phi.make_linspace(),
-        θ2=circ._default_grid_phi.make_linspace(),
-        θ3=0,
-    )
+    @staticmethod
+    def test_plot_potential():
+        circ.plot_potential(
+            θ1=circ._default_grid_phi.make_linspace(),
+            θ2=circ._default_grid_phi.make_linspace(),
+            θ3=0,
+        )

scqubits/tests/test_circuit.py

@@ -39,8 +39,3 @@ # test_circuit.py
         """
-        REFERENCE = (
-            "<bound method Printable.__str__ of "
-            "EJ*cos(θ1 - 1.0*θ3) + EJ*cos(-(2πΦ_{1}) + θ1 + θ3) + "
-            "6.25625*\\dot{θ_1}**2 + 25.0*\\dot{θ_2}**2 + 0.00625*\\dot{θ_3}**2 - "
-            "0.036*θ2**2 - 0.004*θ2*θ3 - 0.009*θ3**2>"
-        )
+        REFERENCE = "<bound method Printable.__str__ of EJ*cos(θ1 + θ3) + EJ*cos((2πΦ_{1}) + θ1 - 1.0*θ3) + 6.25625*\\dot{θ_1}**2 + 25.0*\\dot{θ_2}**2 + 0.00625*\\dot{θ_3}**2 - 0.036*θ2**2 - 0.004*θ2*θ3 - 0.009*θ3**2>"
 
@@ -92,35 +87,34 @@ zero_pi = scq.Circuit(zp_yaml, from_file=False, ext_basis="discretized")
     @staticmethod
-    def test_zero_pi_harmonic():
+    def test_circuit_with_symbolic_hamiltonian():
         """
-        Test for symmetric zero-pi in harmonic oscillator basis.
+        Test for initiating Circuit module with symbolic Hamiltonian.
         """
-        zp_yaml = """# zero-pi circuit
-        branches:
-        - ["JJ", 1, 2, 10, 20]
-        - ["JJ", 3, 4, 10, 20]
-        - ["L", 2, 3, 0.008]
-        - ["L", 4, 1, 0.008]
-        - ["C", 1, 3, 0.02]
-        - ["C", 2, 4, 0.02]
-        """
-        circ = scq.Circuit(zp_yaml, from_file=False, ext_basis="harmonic")
-        circ.cutoff_n_1 = 30
-        circ.cutoff_ext_2 = 30
-        circ.cutoff_ext_3 = 80
-        circ.configure(system_hierarchy=[[1, 3], [2]], subsystem_trunc_dims=[30, 20])
-        circ.cutoff_ext_3 = 200
-        sym_zp = circ.subsystems[0]
-        eigensys = sym_zp.eigensys()
-        eigs = eigensys[0]
-        eig_ref = np.array(
-            [
-                -3.69858244,
-                -3.69261899,
-                -2.90463196,
-                -2.89989473,
-                -2.81204032,
-                -2.81003324,
-            ]
+        import sympy as sm
+
+        sym_hamiltonian = sm.parse_expr(
+            "0.25*θ3**2 + 2.0*Q3**2 + 0.790697674419*Q2**2 + 0.45*θ2**2 + 7.674418604651*n1**2 + 7.674418604651*ng1**2 - 1.0*cos(θ1) + 0.5*θ2*θ3 + 1.395348837209*Q2*n1 + 1.395348837209*Q2*ng1 + 15.348837209302*n1*ng1"
         )
-        assert np.allclose(eigs, eig_ref)
+        circ = scq.Circuit(
+            input_string=None,
+            symbolic_hamiltonian=sym_hamiltonian,
+            symbolic_param_dict={"ng1": 0},
+            ext_basis="harmonic",
+        )
+        circ.configure(
+            transformation_matrix=np.array([[1, 0, 0], [0, 1, 0], [0, 1, 1]])
+        )
+        circ.cutoff_n_1 = 20
+        circ.cutoff_ext_2 = 20
+        circ.cutoff_ext_3 = 20
+        circ.configure(
+            system_hierarchy=[[1], [[2], [3]]],
+            subsystem_trunc_dims=[20, [50, [10, 10]]],
+        )
+        # new_circ.configure(system_hierarchy=[[1], [2, 3]], subsystem_trunc_dims=[20, 30])
+        circ.ng1 = 0.5
+        eigs = circ.eigenvals()
+        eigs_ref = np.array(
+            [2.51547879, 3.00329327, 3.5556228, 3.57568727, 4.13233136, 4.29671029]
+        )
+        assert np.allclose(eigs, eigs_ref)
 
@@ -132,31 +126,26 @@ @staticmethod
         )
-        DFC = scq.Circuit(
-            DATADIR + "circuit_DFC.yaml",
-            ext_basis="harmonic",
-            initiate_sym_calc=False,
-            basis_completion="canonical",
+        inp_yaml = """
+        branches:
+        - [JJ, 0, 1, 1, 15]
+        - [C, 1, 2, 2]
+        - [L, 2, 0, 0.4]
+        - [C, 2, 0, 0.2]
+        - [C, 2, 3, 0.5]
+        - [L, 3, 0, 0.5]
+        # - [JJ, 3, 0, EJ=0, 1e5]
+        """
+        circ = scq.Circuit(inp_yaml, from_file=False, ext_basis="discretized")
+        circ.configure(
+            transformation_matrix=np.array([[1, 0, 0], [0, 1, 0], [0, 1, 1]])
         )
-
-        closure_branches = [DFC.branches[0], DFC.branches[4], DFC.branches[-1]]
-        system_hierarchy = [[[1], [3]], [2], [4]]
-        subsystem_trunc_dims = [[34, [6, 6]], 6, 6]
-
-        DFC.configure(
-            closure_branches=closure_branches,
-            system_hierarchy=system_hierarchy,
-            subsystem_trunc_dims=subsystem_trunc_dims,
+        circ.cutoff_n_1 = 20
+        circ.cutoff_ext_2 = 10
+        circ.cutoff_ext_3 = 10
+        circ.configure(system_hierarchy=[[1], [2, 3]], subsystem_trunc_dims=[20, 30])
+        circ.ng1 = 0.5
+        eigs = circ.eigenvals()
+        generated_eigs = eigs - eigs[0]
+        ref_eigs = np.array(
+            [0.0, 0.48790869, 1.04058606, 1.06037218, 1.61763356, 1.78158506]
         )
-
-        DFC._Φ1 = 0.5 + 0.01768
-        DFC._Φ2 = -0.2662
-        DFC._Φ3 = -0.5 + 0.01768
-
-        DFC._cutoff_ext_1 = 110
-        DFC._cutoff_ext_2 = 110
-        DFC._cutoff_ext_3 = 110
-        DFC._cutoff_ext_4 = 110
-
-        DFC.EJ = 4.6
-        eigs = DFC.eigenvals()
-        generated_eigs = eigs - eigs[0]
         assert np.allclose(generated_eigs, ref_eigs)
@@ -186,12 +175,12 @@
 
-        DFC._Φ1 = 0.5 + 0.01768
-        DFC._Φ2 = -0.2662
-        DFC._Φ3 = -0.5 + 0.01768
+        DFC.Φ1 = 0.5 + 0.01768
+        DFC.Φ2 = -0.2662
+        DFC.Φ3 = -0.5 + 0.01768
 
-        DFC._cutoff_ext_1 = 110
-        DFC._cutoff_ext_2 = 110
-        DFC._cutoff_ext_3 = 110
-        DFC._cutoff_ext_4 = 110
+        DFC.cutoff_ext_1 = 110
+        DFC.cutoff_ext_2 = 110
+        DFC.cutoff_ext_3 = 110
+        DFC.cutoff_ext_4 = 110
+        DFC.update()
 
-        DFC.EJ = 4.6
         eigs = DFC.eigenvals()
@@ -217,2 +206,3 @@ generated_eigs = eigs - eigs[0]
         circ.EJ = 0.01
+        circ.update()
         eigs_ref = np.array(
@@ -232,3 +222,22 @@ [
     @staticmethod
-    def test_param_sweep(num_cpus):
+    def test_noisy_circuit(num_cpus):
+        yaml_inp = f"""branches:
+        - [JJ, 1, 2, EJ=6.8, 1]
+        - [L, 1, 2, 0.2]
+        """
+        circ = scq.Circuit(
+            yaml_inp,
+            from_file=False,
+            ext_basis="harmonic",
+            use_dynamic_flux_grouping=True,
+        )
+        circ.cutoff_ext_1 = 300
+        circ.Φ1 = 0.5
+        circ.configure(generate_noise_methods=True)
+        coherence_times = np.array([circ.t1_effective(), circ.t2_effective()])
+        ref = np.array([3319890.8160632304, 5385675.324726781])
+        assert np.allclose(coherence_times, ref)
+
+    @staticmethod
+    def test_get_spectrum_vs_paramvals(num_cpus):
         DFC = scq.Circuit(
@@ -251,13 +260,30 @@ DATADIR + "circuit_DFC.yaml",
 
-        DFC._Φ1 = 0.5 + 0.01768
-        DFC._Φ2 = -0.2662
-        DFC._Φ3 = -0.5 + 0.01768
+        DFC.Φ1 = 0.5 + 0.01768
+        DFC.Φ2 = -0.2662
+        DFC.Φ3 = -0.5 + 0.01768
 
-        DFC._cutoff_ext_1 = 110
-        DFC._cutoff_ext_2 = 110
-        DFC._cutoff_ext_3 = 110
-        DFC._cutoff_ext_4 = 110
+        DFC.cutoff_ext_1 = 110
+        DFC.cutoff_ext_2 = 110
+        DFC.cutoff_ext_3 = 110
+        DFC.cutoff_ext_4 = 110
+        DFC.update()
+        DFC.get_spectrum_vs_paramvals("Φ1", np.linspace(0, 1, 11), num_cpus=num_cpus)
 
-        DFC.EJ = 4.6
+        paramvals_by_name = {
+            "Φ1": np.linspace(0.4, 0.5, 6),
+            "Φ2": np.linspace(0.4, 0.5, 3),
+        }
 
-        DFC.get_spectrum_vs_paramvals("Φ1", np.linspace(0, 1, 11), num_cpus=num_cpus)
+        # update Hilbert space function
+        def update_hilbertspace(Φ1, Φ2):
+            DFC.Φ1 = Φ1
+            DFC.Φ2 = Φ2
+            DFC.update()
+
+        ps = scq.ParameterSweep(
+            hilbertspace=DFC.hilbert_space,
+            paramvals_by_name=paramvals_by_name,
+            update_hilbertspace=update_hilbertspace,
+            evals_count=6,
+            num_cpus=num_cpus,
+        )

scqubits/tests/test_units.py

@@ -29,12 +29,12 @@ # test_units.py
 
-    def test_units_warning(self):
-        scq.set_units("GHz")
-        qubit = Transmon(EJ=0.5, EC=12.0, ng=0.3, ncut=150)
-        # Expect a warning when changing units since a QuantumSystem is present
-        with pytest.warns(UserWarning):
-            scq.set_units("MHz")
+    # def test_units_warning(self):
+    #     scq.set_units("GHz")
+    #     qubit = Transmon(EJ=0.5, EC=12.0, ng=0.3, ncut=150)
+    #     # Expect a warning when changing units since a QuantumSystem is present
+    #     with pytest.warns(UserWarning):
+    #         scq.set_units("MHz")
 
-        # Do not expect warning after deleting the only QuantumSystem
-        del qubit
-        scq.set_units("kHz")
+    #     # Do not expect warning after deleting the only QuantumSystem
+    #     del qubit
+    #     scq.set_units("kHz")
 
@@ -41,0 +41,0 @@ def test_units_auxiliary(self):

scqubits/ui/gui_custom_widgets.py

@@ -371,8 +371,8 @@ # gui_custom_widgets.py
             self.plot_settings_dialog = plot_settings_dialog
-            self.panel_by_btn: OrderedDict[
-                v.Btn, ClosablePlotPanel
-            ] = collections.OrderedDict()
-            self.panel_by_id: OrderedDict[
-                str, ClosablePlotPanel
-            ] = collections.OrderedDict()
+            self.panel_by_btn: OrderedDict[v.Btn, ClosablePlotPanel] = (
+                collections.OrderedDict()
+            )
+            self.panel_by_id: OrderedDict[str, ClosablePlotPanel] = (
+                collections.OrderedDict()
+            )
             self.toggle_switches_by_plot_id = toggle_switches_by_plot_id
@@ -379,0 +379,0 @@

scqubits/ui/gui_defaults.py

@@ -155,3 +155,3 @@ # gui_defaults.py
     "EL": EL_range,
-    "ECJ": EC_range,
+    "ECJ": {"v_min": STEP, "v_max": 25.0},
     "dEJ": {"v_min": 0.0, "v_max": 1.0},
@@ -283,3 +283,2 @@ "dCJ": {"v_min": 0.0, "v_max": 1.0},
                 "EJ": 6.0,
-                "EJ": 6.0,
                 "ECJ": 2.28,
@@ -286,0 +285,0 @@ "EC": 0.184,

scqubits/ui/hspace_widget.py

@@ -302,11 +302,11 @@ # hspace_widget.py
             value = self.op1subsys_widget[self.current_interaction_idx].v_model
-            self.op1_ddown_widget[
-                self.current_interaction_idx
-            ].items = self.possible_operators(value)
+            self.op1_ddown_widget[self.current_interaction_idx].items = (
+                self.possible_operators(value)
+            )
 
         def on_op_subsys2_selected(*args):
             value = self.op2subsys_widget[self.current_interaction_idx].v_model
-            self.op2_ddown_widget[
-                self.current_interaction_idx
-            ].items = self.possible_operators(value)
+            self.op2_ddown_widget[self.current_interaction_idx].items = (
+                self.possible_operators(value)
+            )
 
@@ -313,0 +313,0 @@ self.op1subsys_widget[-1].observe(on_op_subsys1_selected, names="v_model")

scqubits/utils/misc.py

@@ -187,2 +187,16 @@ # misc.py
 
+def check_sync_status_circuit(func: Callable) -> Callable:
+    @functools.wraps(func)
+    def wrapper(self, *args, **kwargs):
+        # update the circuit if necessary
+        if (self._user_changed_parameter) or (
+            self.hierarchical_diagonalization
+            and (self._out_of_sync or len(self.affected_subsystem_indices) > 0)
+        ):
+            self.update()
+        return func(self, *args, **kwargs)
+
+    return wrapper
+
+
 def check_lookup_exists(func: Callable) -> Callable:
@@ -238,5 +252,17 @@ @functools.wraps(func)
     # is a Qobj ket.
-    return np.asarray(qobj_ket.data.todense())
+    return (
+        qobj_ket.data.as_ndarray()
+        if qt.__version__ >= "5.0.0"
+        else qobj_ket.data.toarray()
+    )
 
 
+def Qobj_to_scipy_csc_matrix(qobj_array: qt.Qobj) -> sp.sparse.csc_matrix:
+    return (
+        qobj_array.to("csr").data.as_scipy().tocsc()
+        if qt.__version__ >= "5.0.0"
+        else qobj_array.data.tocsc()
+    )
+
+
 def get_shape(lst, shape=()):
@@ -386,3 +412,3 @@ """
 def list_intersection(list1: list, list2: list) -> list:
-    return list(set(list1) & set(list2))
+    return [item for item in list1 if item in list2]
 
@@ -429,2 +455,20 @@
 
+def unique_elements_in_list(list_object: list) -> list:
+    """
+    Returns a list of all the unique elements in the list
+
+    Parameters
+    ----------
+    list_object :
+        A list of any objects
+    """
+    unique_list = []
+    [
+        unique_list.append(element)
+        for element in list_object
+        if element not in unique_list
+    ]
+    return unique_list
+
+
 def number_of_lists_in_list(list_object: list) -> int:
@@ -431,0 +475,0 @@ """

scqubits/utils/spectrum_utils.py

@@ -23,3 +23,3 @@ # spectrum_utils.py
 from qutip import Qobj
-from scipy.sparse import csc_matrix, dia_matrix
+from scipy.sparse import csc_matrix, dia_matrix, csr_matrix
 
@@ -34,2 +34,3 @@ import scqubits.settings as settings
 from scqubits.utils.typedefs import QuantumSys
+from scqubits.utils.misc import Qobj_to_scipy_csc_matrix
 
@@ -98,3 +99,3 @@
     if position is None:
-        halfway_position = len(complex_array) // 2
+        halfway_position = (complex_array.shape[0]) // 2
         flattened_position = np.argmax(
@@ -163,3 +164,3 @@ np.abs(complex_array[:halfway_position])
     if isinstance(state1, qt.Qobj):
-        vec1 = state1.data.toarray()
+        vec1 = Qobj_to_scipy_csc_matrix(state1)
     else:
@@ -169,3 +170,3 @@ vec1 = state1
     if isinstance(state2, qt.Qobj):
-        vec2 = state2.data.toarray()
+        vec2 = Qobj_to_scipy_csc_matrix(state2)
     else:
@@ -321,4 +322,4 @@ vec2 = state2
         operator_matrixelements = get_matrixelement_table(operator, evecs)
-        return qt.Qobj(inpt=operator_matrixelements)
-    return qt.Qobj(inpt=operator[:dim, :dim])
+        return qt.Qobj(operator_matrixelements)
+    return qt.Qobj(operator[:dim, :dim])
 
@@ -336,3 +337,3 @@
     operator_matrixelements = subsystem.matrixelement_table(operator, evecs=evecs)
-    return qt.Qobj(inpt=operator_matrixelements)
+    return qt.Qobj(operator_matrixelements)
 
@@ -348,3 +349,3 @@
         return operator
-    if isinstance(operator, (np.ndarray, csc_matrix, dia_matrix)):
+    if isinstance(operator, (np.ndarray, csc_matrix, csr_matrix, dia_matrix)):
         return convert_matrix_to_qobj(operator, subsystem, op_in_eigenbasis, evecs)
@@ -351,0 +352,0 @@ if isinstance(operator, str):

scqubits/version.py

 # THIS FILE IS GENERATED FROM scqubits SETUP.PY
-short_version = '3.3.0'
-version = '3.3.0'
+short_version = '4.0.0'
+version = '4.0.0'
 release = True

setup.py

@@ -20,2 +20,3 @@ """scqubits: superconducting qubits in Python
 """
+
 #
@@ -55,4 +56,4 @@ # This file is part of scqubits.
 # version information about scqubits goes here
-MAJOR = 3
-MINOR = 3
+MAJOR = 4
+MINOR = 0
 MICRO = 0
@@ -59,0 +60,0 @@ ISRELEASED = True

Alert delta unavailable

Currently unable to show alert delta for PyPI packages.

Improved metrics

Total package byte prevSize

8638136

1.75%

Number of package files

141

2.17%

Lines of code

30815

12.2%