# spectrum_utils.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.
############################################################################
import cmath
from typing import TYPE_CHECKING, Callable, List, Optional, Tuple, Union
import numpy as np
import qutip as qt
import scipy as sp
from numpy import ndarray
from qutip import Qobj
from scipy.sparse import csc_matrix, dia_matrix, csr_matrix
import scqubits.settings as settings
if TYPE_CHECKING:
from scqubits import Oscillator, ParameterSweep, SpectrumData
from scqubits.core.qubit_base import QubitBaseClass
from scqubits.io_utils.fileio_qutip import QutipEigenstates
from scqubits.utils.typedefs import QuantumSys
from scqubits.utils.misc import Qobj_to_scipy_csc_matrix
[docs]def eigsh_safe(*args, **kwargs):
"""Wrapper method for `scipy.sparse.linalg.eigsh` which ensures the following.
1. Always use the same "random" starting vector v0. Otherwise, results show
random behavior (small deviations between different runs, problem for pytests)
2. Test for degenerate eigenvalues. If there are any, need to orthogonalize the
eigenvectors properly."""
mat_size = args[0].shape[0]
kwargs["v0"] = settings.RANDOM_ARRAY[:mat_size]
if kwargs.get("return_eigenvectors"):
evals, evecs = sp.sparse.linalg.eigsh(*args, **kwargs)
if has_degeneracy(evals):
evecs, _ = sp.linalg.qr(evecs, mode="economic")
return evals, evecs
return sp.sparse.linalg.eigsh(*args, **kwargs)
[docs]def has_degeneracy(evals: ndarray) -> bool:
evals_rightpad = np.pad(evals, (0, 1))
evals_leftpad = np.pad(evals, (1, 0))
evals_neighbor_diffs = evals_leftpad - evals_rightpad
return np.isclose(np.min(np.abs(evals_neighbor_diffs)), 0)
[docs]def order_eigensystem(
evals: np.ndarray, evecs: np.ndarray
) -> Tuple[np.ndarray, np.ndarray]:
"""Takes eigenvalues and corresponding eigenvectors and orders them (in place)
according to the eigenvalues (from smallest to largest; real valued eigenvalues
are assumed). Compare http://stackoverflow.com/questions/22806398.
Parameters
----------
evals:
array of eigenvalues
evecs:
array containing eigenvectors; evecs[:, 0] is the first eigenvector etc.
"""
ordered_evals_indices = evals.argsort() # sort manually
evals[:] = evals[ordered_evals_indices]
evecs[:] = evecs[:, ordered_evals_indices]
return evals, evecs
[docs]def standardize_phases(complex_array: np.ndarray) -> np.ndarray:
"""Uses `extract_phase` to obtain global phase from `array` and returns
standardized array with global phase factor standardized.
Parameters
----------
complex_array:
complex
"""
phase = extract_phase(complex_array)
std_array = complex_array * np.exp(-1j * phase)
return std_array
[docs]def standardize_sign(real_array: np.ndarray) -> np.ndarray:
"""Standardizes the sign of a real-valued wavefunction by calculating the sign of
the sum of all amplitudes up to the wavefunctions mid-position and making it
positive.
Summing up to the midpoint only is to address the danger that the sum is
actually zero, which may be the case for odd wavefunctions taken over an interval
centered at zero.
"""
halfway_position = len(real_array) // 2
return np.sign(np.sum(real_array[:halfway_position])) * real_array
# -Matrix elements and operators (outside qutip) --------------------------------------
[docs]def matrix_element(
state1: Union[np.ndarray, qt.Qobj],
operator: Union[np.ndarray, csc_matrix, qt.Qobj],
state2: Union[np.ndarray, qt.Qobj],
) -> Union[float, complex]:
"""Calculate the matrix element `<state1|operator|state2>`.
Parameters
----------
state1:
state vector/ket
state2:
state vector/ket
operator:
representation of an operator
Returns
-------
matrix element
"""
if isinstance(operator, qt.Qobj):
op_matrix = operator.data
else:
op_matrix = operator
if isinstance(state1, qt.Qobj):
vec1 = Qobj_to_scipy_csc_matrix(state1)
else:
vec1 = state1
if isinstance(state2, qt.Qobj):
vec2 = Qobj_to_scipy_csc_matrix(state2)
else:
vec2 = state2
return vec1.conj().T @ op_matrix @ vec2
[docs]def get_matrixelement_table(
operator: Union[np.ndarray, csc_matrix, dia_matrix, qt.Qobj],
state_table: Union[np.ndarray, qt.Qobj],
) -> np.ndarray:
"""Calculates a table of matrix elements.
Parameters
----------
operator:
operator with respect to which matrix elements are to be calculated
state_table:
list or array of numpy arrays representing the states `|v0>, |v1>, ...`
Note: `state_table` is expected to be in scipy's `eigsh` transposed form.
Returns
-------
table of matrix elements
"""
if isinstance(operator, qt.Qobj):
states_in_columns = state_table.T
else:
states_in_columns = state_table
mtable = states_in_columns.conj().T @ operator @ states_in_columns
return mtable
[docs]def closest_dressed_energy(
bare_energy: float, dressed_energy_vals: np.ndarray
) -> float:
"""For a given bare energy value, this returns the closest lying dressed energy
value from an array.
Parameters
----------
bare_energy:
bare energy value
dressed_energy_vals:
array of dressed-energy values
Returns
-------
element from `dressed_energy_vals` closest to `bare_energy`
"""
index = (np.abs(dressed_energy_vals - bare_energy)).argmin()
return dressed_energy_vals[index]
[docs]def get_eigenstate_index_maxoverlap(
eigenstates_qobj: "QutipEigenstates",
reference_state_qobj: qt.Qobj,
return_overlap: bool = False,
) -> Union[int, Tuple[int, float], None]:
"""For given list of qutip states, find index of the state that has largest
overlap with the qutip ket `reference_state_qobj`. If `|overlap|` is smaller than
0.5, return None.
Parameters
----------
eigenstates_qobj:
as obtained from qutip `.eigenstates()`
reference_state_qobj:
specific reference state
return_overlap:
set to true if the value of largest overlap should be also returned
(default value = False)
Returns
-------
index of eigenstate from `eigenstates_Qobj` with the largest overlap with the
`reference_state_qobj`, None if `|overlap|<0.5`
"""
overlaps = np.asarray(
[
eigenstates_qobj[j].overlap(reference_state_qobj)
for j in range(len(eigenstates_qobj))
]
)
max_overlap = np.max(np.abs(overlaps))
if max_overlap < 0.5:
return None
index = (np.abs(overlaps)).argmax()
if return_overlap:
return index, np.abs(overlaps[index])
return index
[docs]def absorption_spectrum(spectrum_data: "SpectrumData") -> "SpectrumData":
"""Takes spectral data of energy eigenvalues and returns the absorption spectrum
relative to a state of given index. Calculated by subtracting from eigenenergies
the energy of the select state. Resulting negative frequencies, if the reference
state is not the ground state, are omitted.
"""
assert isinstance(spectrum_data.energy_table, ndarray)
spectrum_data.energy_table = spectrum_data.energy_table.clip(min=0.0) # type:ignore
return spectrum_data
[docs]def emission_spectrum(spectrum_data: "SpectrumData") -> "SpectrumData":
"""Takes spectral data of energy eigenvalues and returns the emission spectrum
relative to a state of given index. The resulting "upwards" transition
frequencies are calculated by subtracting from eigenenergies the energy of the
select state, and multiplying the result by -1. Resulting negative frequencies,
corresponding to absorption instead, are omitted.
"""
assert isinstance(spectrum_data.energy_table, ndarray)
spectrum_data.energy_table *= -1.0
spectrum_data.energy_table = spectrum_data.energy_table.clip(min=0.0) # type:ignore
return spectrum_data
[docs]def convert_evecs_to_ndarray(evecs_qutip: ndarray) -> np.ndarray:
"""Takes a qutip eigenstates array, as obtained with .eigenstates(), and converts
it into a pure numpy array.
Parameters
----------
evecs_qutip:
ndarray of eigenstates in qt.Qobj format
Returns
-------
converted eigenstate data
"""
evals_count = len(evecs_qutip)
dimension = evecs_qutip[0].shape[0]
evecs_ndarray = np.empty((evals_count, dimension), dtype=np.complex_)
for index, eigenstate in enumerate(evecs_qutip):
evecs_ndarray[index] = eigenstate.full()[:, 0]
return evecs_ndarray
[docs]def convert_matrix_to_qobj(
operator: Union[np.ndarray, csc_matrix, dia_matrix],
subsystem: Union["QubitBaseClass", "Oscillator"],
op_in_eigenbasis: bool,
evecs: Optional[np.ndarray],
) -> qt.Qobj:
dim = subsystem.truncated_dim
if op_in_eigenbasis is False:
if evecs is None:
_, evecs = subsystem.eigensys(evals_count=dim)
operator_matrixelements = get_matrixelement_table(operator, evecs)
return qt.Qobj(operator_matrixelements)
return qt.Qobj(operator[:dim, :dim])
[docs]def convert_opstring_to_qobj(
operator: str,
subsystem: Union["QubitBaseClass", "Oscillator"],
evecs: Optional[np.ndarray],
) -> qt.Qobj:
dim = subsystem.truncated_dim
if evecs is None:
_, evecs = subsystem.eigensys(evals_count=dim)
operator_matrixelements = subsystem.matrixelement_table(operator, evecs=evecs)
return qt.Qobj(operator_matrixelements)
[docs]def convert_operator_to_qobj(
operator: Union[np.ndarray, csc_matrix, dia_matrix, qt.Qobj, str],
subsystem: Union["QubitBaseClass", "Oscillator"],
op_in_eigenbasis: bool,
evecs: Optional[np.ndarray],
) -> qt.Qobj:
if isinstance(operator, qt.Qobj):
return operator
if isinstance(operator, (np.ndarray, csc_matrix, csr_matrix, dia_matrix)):
return convert_matrix_to_qobj(operator, subsystem, op_in_eigenbasis, evecs)
if isinstance(operator, str):
return convert_opstring_to_qobj(operator, subsystem, evecs)
raise TypeError("Unsupported operator type: ", type(operator))
[docs]def generate_target_states_list(
sweep: "ParameterSweep", initial_state_labels: Tuple[int, ...]
) -> List[Tuple[int, ...]]:
"""Based on a bare state label (i1, i2, ...) with i1 being the excitation level
of subsystem 1, i2 the excitation level of subsystem 2 etc., generate a list of
new bare state labels. These bare state labels correspond to target states
reached from the given initial one by single-photon qubit transitions. These are
transitions where one of the qubit excitation levels increases at a time. There
are no changes in oscillator photon numbers.
Parameters
----------
sweep:
initial_state_labels:
bare-state labels of the initial state whose energy is supposed to be subtracted
from the spectral data
"""
target_states_list = []
for qbt_subsys in sweep.qbt_subsys_list: # iterate through qubit subsys_list
assert qbt_subsys.truncated_dim is not None
subsys_index = sweep._hilbertspace.get_subsys_index(qbt_subsys)
initial_qbt_state = initial_state_labels[subsys_index]
for state_label in range(initial_qbt_state + 1, qbt_subsys.truncated_dim):
# for given qubit subsystem, generate target labels by increasing that qubit
# excitation level
target_labels = list(initial_state_labels)
target_labels[subsys_index] = state_label
target_states_list.append(tuple(target_labels))
return target_states_list
[docs]def recast_esys_mapdata(
esys_mapdata: List[Tuple[np.ndarray, np.ndarray]]
) -> Tuple[np.ndarray, List[np.ndarray]]:
"""
Takes data generated by a map of eigensystem calls and returns the eigenvalue and
eigenstate tables
Returns
-------
eigenvalues and eigenvectors
"""
paramvals_count = len(esys_mapdata)
eigenenergy_table = np.asarray(
[esys_mapdata[index][0] for index in range(paramvals_count)]
)
eigenstate_table = [esys_mapdata[index][1] for index in range(paramvals_count)]
return eigenenergy_table, eigenstate_table
[docs]def identity_wrap(
operator: Union[str, ndarray, Qobj, Callable],
subsystem: "QuantumSys",
subsys_list: List["QuantumSys"],
op_in_eigenbasis: bool = False,
evecs: Optional[ndarray] = None,
) -> Qobj:
"""Takes the `operator` belonging to `subsystem` and "wraps" it in identities.
The full Hilbert space is taken to consist of all subsystems given as
`subsys_list`. `subsystem` must be one element in that list. For each of the
other subsystems in the list, an identity operator of the correct dimension is
generated and inserted into the appropriate Kronecker product "sandwiching" the
operator.
Parameters
----------
operator:
operator acting in Hilbert space of `subsystem`; if str, then this should be an
operator name in the subsystem, typically not in eigenbasis
subsystem:
subsystem where diagonal operator is defined
subsys_list:
list of all subsystems relevant to the Hilbert space.
op_in_eigenbasis:
whether `operator` is given in the `subsystem` eigenbasis; otherwise,
`operator` is assumed to be in the internal QuantumSystem basis. This
argument is ignored if `operator` is given as a Qobj.
evecs:
internal `QuantumSystem` eigenstates, used to convert `operator` into eigenbasis
Returns
-------
operator in the full Hilbert space (as specified by `subsystem_list`). This
operator is expressed in the bare product basis consisting of the energy
eigenstates of each subsystem (unless `operator` is provided as a `Qobj`,
in which case no conversion takes place).
"""
if not isinstance(operator, qt.Qobj) and callable(operator):
try:
operator = operator(energy_esys=(None, evecs))
except TypeError:
operator = operator()
op_in_eigenbasis = True
subsys_operator = convert_operator_to_qobj(
operator, subsystem, op_in_eigenbasis, evecs # type:ignore
)
operator_identitywrap_list = [
qt.operators.qeye(the_subsys.truncated_dim) for the_subsys in subsys_list
]
subsystem_index = subsys_list.index(subsystem)
operator_identitywrap_list[subsystem_index] = subsys_operator
return qt.tensor(operator_identitywrap_list)