1# This file is part of QuTiP: Quantum Toolbox in Python.
2#
3#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
4#    All rights reserved.
5#
6#    Redistribution and use in source and binary forms, with or without
7#    modification, are permitted provided that the following conditions are
8#    met:
9#
10#    1. Redistributions of source code must retain the above copyright notice,
11#       this list of conditions and the following disclaimer.
12#
13#    2. Redistributions in binary form must reproduce the above copyright
14#       notice, this list of conditions and the following disclaimer in the
15#       documentation and/or other materials provided with the distribution.
16#
17#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
18#       of its contributors may be used to endorse or promote products derived
19#       from this software without specific prior written permission.
20#
21#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
24#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32###############################################################################
33from collections.abc import Iterable
34import warnings
35from copy import deepcopy
36
37import numpy as np
38from scipy.interpolate import CubicSpline
39
40from qutip.qobj import Qobj
41from qutip.qobjevo import QobjEvo
42from qutip.operators import identity
43from qutip.qip.operations.gates import expand_operator, globalphase
44from qutip.tensor import tensor
45from qutip.mesolve import mesolve
46from qutip.mcsolve import mcsolve
47from qutip.qip.circuit import QubitCircuit
48from qutip.qip.noise import (
49    Noise, RelaxationNoise, DecoherenceNoise,
50    ControlAmpNoise, RandomNoise, process_noise)
51from qutip.qip.pulse import Pulse, Drift, _merge_qobjevo, _fill_coeff
52
53
54__all__ = ['Processor']
55
56
57class Processor(object):
58    """
59    A simulator of a quantum device based on the QuTiP solver
60    :func:`qutip.mesolve`.  It is defined by the available driving Hamiltonian
61    and the decoherence time for each component systems.  The processor can
62    simulate the evolution under the given control pulses. Noisy evolution is
63    supported by :class:`.Noise` and can be added to the processor.
64
65    Parameters
66    ----------
67    N: int
68        The number of component systems.
69
70    t1: list or float, optional
71        Characterize the decoherence of amplitude damping for
72        each qubit. A list of size `N` or a float for all qubits.
73
74    t2: list of float, optional
75        Characterize the decoherence of dephasing for
76        each qubit. A list of size `N` or a float for all qubits.
77
78    dims: list, optional
79        The dimension of each component system.
80        Default value is a
81        qubit system of ``dim=[2,2,2,...,2]``
82
83    spline_kind: str, optional
84        Type of the coefficient interpolation. Default is "step_func"
85        Note that they have different requirement for the length of ``coeff``.
86
87        - "step_func":
88          The coefficient will be treated as a step function.  E.g.
89          ``tlist=[0,1,2]`` and ``coeff=[3,2]``, means that the coefficient is
90          3 in t=[0,1) and 2 in t=[2,3). It requires
91          ``len(coeff)=len(tlist)-1`` or ``len(coeff)=len(tlist)``, but in the
92          second case the last element of `coeff` has no effect.
93
94        - "cubic": Use cubic interpolation for the coefficient. It requires
95          ``len(coeff)=len(tlist)``
96
97    Attributes
98    ----------
99    N: int
100        The number of component systems.
101
102    pulses: list of :class:`.Pulse`
103        A list of control pulses of this device
104
105    t1: float or list
106        Characterize the decoherence of amplitude damping of
107        each qubit.
108
109    t2: float or list
110        Characterize the decoherence of dephasing for
111        each qubit.
112
113    noise: :class:`.Noise`, optional
114        A list of noise objects. They will be processed when creating the
115        noisy :class:`qutip.QobjEvo` from the processor or run the simulation.
116
117    drift: :class:`qutip.qip.pulse.Drift`
118        A `Drift` object representing the drift Hamiltonians.
119
120    dims: list
121        The dimension of each component system.
122        Default value is a
123        qubit system of ``dim=[2,2,2,...,2]``
124
125    spline_kind: str
126        Type of the coefficient interpolation.
127        See parameters of :class:`.Processor` for details.
128    """
129    def __init__(self, N, t1=None, t2=None,
130                 dims=None, spline_kind="step_func"):
131        self.N = N
132        self.pulses = []
133        self.t1 = t1
134        self.t2 = t2
135        self.noise = []
136        self.drift = Drift()
137        if dims is None:
138            self.dims = [2] * N
139        else:
140            self.dims = dims
141        self.pulse_mode = "discrete"
142        self.spline_kind = spline_kind
143
144    @property
145    def num_qubits(self):
146        return self.N
147
148    @num_qubits.setter
149    def num_qubits(self, value):
150        self.N = value
151
152    def add_drift(self, qobj, targets, cyclic_permutation=False):
153        """
154        Add a drift Hamiltonians. The drift Hamiltonians are intrinsic
155        of the quantum system and cannot be controlled by external field.
156
157        Parameters
158        ----------
159        qobj: :class:`qutip.Qobj`
160            The drift Hamiltonian.
161        targets: list
162            The indices of the target qubits
163            (or subquantum system of other dimensions).
164        """
165        if not isinstance(qobj, Qobj):
166            raise TypeError("The drift Hamiltonian must be a qutip.Qobj.")
167        if not qobj.isherm:
168            raise ValueError("The drift Hamiltonian must be Hermitian.")
169
170        num_qubits = len(qobj.dims[0])
171        if targets is None:
172            targets = list(range(num_qubits))
173        if not isinstance(targets, list):
174            targets = [targets]
175        if cyclic_permutation:
176            for i in range(self.N):
177                temp_targets = [(t + i) % self.N for t in targets]
178                self.drift.add_drift(qobj, temp_targets)
179        else:
180            self.drift.add_drift(qobj, targets)
181
182    def add_control(self, qobj, targets=None, cyclic_permutation=False,
183                    label=None):
184        """
185        Add a control Hamiltonian to the processor. It creates a new
186        :class:`.Pulse`
187        object for the device that is turned off
188        (``tlist = None``, ``coeff = None``). To activate the pulse, one
189        can set its `tlist` and `coeff`.
190
191        Parameters
192        ----------
193        qobj: :class:`qutip.Qobj`
194            The Hamiltonian for the control pulse..
195
196        targets: list, optional
197            The indices of the target qubits
198            (or subquantum system of other dimensions).
199
200        cyclic_permutation: bool, optional
201            If true, the Hamiltonian will be expanded for
202            all cyclic permutation of the target qubits.
203
204        label: str, optional
205            The label (name) of the pulse
206        """
207        # Check validity of ctrl
208        if not isinstance(qobj, Qobj):
209            raise TypeError("The control Hamiltonian must be a qutip.Qobj.")
210        if not qobj.isherm:
211            raise ValueError("The control Hamiltonian must be Hermitian.")
212
213        num_qubits = len(qobj.dims[0])
214        if targets is None:
215            targets = list(range(num_qubits))
216        if not isinstance(targets, list):
217            targets = [targets]
218        if cyclic_permutation:
219            for i in range(self.N):
220                temp_targets = [(t + i) % self.N for t in targets]
221                if label is not None:
222                    temp_label = label + "_" + str(temp_targets)
223                temp_label = label
224                self.pulses.append(
225                    Pulse(qobj, temp_targets, spline_kind=self.spline_kind,
226                          label=temp_label))
227        else:
228            self.pulses.append(
229                Pulse(qobj, targets, spline_kind=self.spline_kind, label=label)
230                )
231
232    def find_pulse(self, pulse_name):
233        if isinstance(pulse_name, str):
234            try:
235                return self.pulses[self.pulse_dict[pulse_name]]
236            except (KeyError):
237                raise KeyError(
238                    "Pulse name {} undefined. "
239                    "Please define it in the attribute "
240                    "`pulse_dict`.".format(pulse_name))
241        elif isinstance(pulse_name, int):
242            return self.pulses[pulse_name]
243        else:
244            raise TypeError(
245                "pulse_name is either a string or an integer, not "
246                "{}".format(type(pulse_name))
247                )
248
249    @property
250    def ctrls(self):
251        """
252        A list of Hamiltonians of all pulses.
253        """
254        result = []
255        for pulse in self.pulses:
256            result.append(pulse.get_ideal_qobj(self.dims))
257        return result
258
259    @property
260    def coeffs(self):
261        """
262        A list of the coefficients for all control pulses.
263        """
264        if not self.pulses:
265            return None
266        coeffs_list = [pulse.coeff for pulse in self.pulses]
267        return coeffs_list
268
269    @coeffs.setter
270    def coeffs(self, coeffs_list):
271        for i, coeff in enumerate(coeffs_list):
272            self.pulses[i].coeff = coeff
273
274    @property
275    def pulse_mode(self):
276        if self.spline_kind == "step_func":
277            return "discrete"
278        elif self.spline_kind == "cubic":
279            return "continuous"
280        else:
281            raise ValueError(
282                "Saved spline_kind not understood.")
283
284    @pulse_mode.setter
285    def pulse_mode(self, mode):
286        if mode == "discrete":
287            spline_kind = "step_func"
288        elif mode == "continuous":
289            spline_kind = "cubic"
290        else:
291            raise ValueError(
292                "Pulse mode must be either discrete or continuous.")
293
294        self.spline_kind = spline_kind
295        for pulse in self.pulses:
296            pulse.spline_kind = spline_kind
297
298    def get_full_tlist(self, tol=1.0e-10):
299        """
300        Return the full tlist of the ideal pulses.
301        If different pulses have different time steps,
302        it will collect all the time steps in a sorted array.
303
304        Returns
305        -------
306        full_tlist: array-like 1d
307            The full time sequence for the ideal evolution.
308        """
309        full_tlist = [pulse.tlist
310                      for pulse in self.pulses if pulse.tlist is not None]
311        if not full_tlist:
312            return None
313        full_tlist = np.unique(np.sort(np.hstack(full_tlist)))
314        # account for inaccuracy in float-point number
315        full_tlist = np.concatenate(
316            (full_tlist[:1], full_tlist[1:][np.diff(full_tlist) > tol]))
317        return full_tlist
318
319    def get_full_coeffs(self, full_tlist=None):
320        """
321        Return the full coefficients in a 2d matrix form.
322        Each row corresponds to one pulse. If the `tlist` are
323        different for different pulses, the length of each row
324        will be same as the `full_tlist` (see method
325        `get_full_tlist`). Interpolation is used for
326        adding the missing coefficient according to `spline_kind`.
327
328        Returns
329        -------
330        coeffs: array-like 2d
331            The coefficients for all ideal pulses.
332        """
333        # TODO add tests
334        self._is_pulses_valid()
335        if not self.pulses:
336            return np.array((0, 0), dtype=float)
337        if full_tlist is None:
338            full_tlist = self.get_full_tlist()
339        coeffs_list = []
340        for pulse in self.pulses:
341            if pulse.tlist is None and pulse.coeff is None:
342                coeffs_list.append(np.zeros(len(full_tlist)))
343                continue
344            if not isinstance(pulse.coeff, (bool, np.ndarray)):
345                raise ValueError(
346                    "get_full_coeffs only works for "
347                    "NumPy array or bool coeff.")
348            if isinstance(pulse.coeff, bool):
349                if pulse.coeff:
350                    coeffs_list.append(np.ones(len(full_tlist)))
351                else:
352                    coeffs_list.append(np.zeros(len(full_tlist)))
353                continue
354            if self.spline_kind == "step_func":
355                arg = {"_step_func_coeff": True}
356                coeffs_list.append(
357                    _fill_coeff(pulse.coeff, pulse.tlist, full_tlist, arg))
358            elif self.spline_kind == "cubic":
359                coeffs_list.append(
360                    _fill_coeff(pulse.coeff, pulse.tlist, full_tlist, {}))
361            else:
362                raise ValueError("Unknown spline kind.")
363        return np.array(coeffs_list)
364
365    def set_all_tlist(self, tlist):
366        """
367        Set the same `tlist` for all the pulses.
368
369        Parameters
370        ----------
371        tlist: array-like, optional
372            A list of time at which the time-dependent coefficients are
373            applied. See :class:`.Pulse` for detailed information`
374        """
375        if isinstance(tlist, list) and len(tlist) == len(self.pulses):
376            for i, pulse in enumerate(self.pulses):
377                pulse.tlist = tlist[i]
378        else:
379            for pulse in self.pulses:
380                pulse.tlist = tlist
381
382    def add_pulse(self, pulse):
383        """
384        Add a new pulse to the device.
385
386        Parameters
387        ----------
388        pulse: :class:`.Pulse`
389            `Pulse` object to be added.
390        """
391        if isinstance(pulse, Pulse):
392            if pulse.spline_kind is None:
393                pulse.spline_kind = self.spline_kind
394            self.pulses.append(pulse)
395        else:
396            raise ValueError("Invalid input, pulse must be a Pulse object")
397
398    def remove_pulse(self, indices=None, label=None):
399        """
400        Remove the control pulse with given indices.
401
402        Parameters
403        ----------
404        indices: int or list of int
405            The indices of the control Hamiltonians to be removed.
406        label: str
407            The label of the pulse
408        """
409        if indices is not None:
410            if not isinstance(indices, Iterable):
411                indices = [indices]
412            indices.sort(reverse=True)
413            for ind in indices:
414                del self.pulses[ind]
415        else:
416            for ind, pulse in enumerate(self.pulses):
417                if pulse.label == label:
418                    del self.pulses[ind]
419
420    def _is_pulses_valid(self):
421        """
422        Check if the pulses are in the correct shape.
423
424        Returns: bool
425            If they are valid or not
426        """
427        for i, pulse in enumerate(self.pulses):
428            if pulse.coeff is None or isinstance(pulse.coeff, bool):
429                # constant pulse
430                continue
431            if pulse.tlist is None:
432                raise ValueError(
433                    "Pulse id={} is invalid. "
434                    "Please define a tlist for the pulse.".format(i))
435            if pulse.tlist is not None and pulse.coeff is None:
436                raise ValueError(
437                    "Pulse id={} is invalid. "
438                    "Please define a coeff for the pulse.".format(i))
439            coeff_len = len(pulse.coeff)
440            tlist_len = len(pulse.tlist)
441            if pulse.spline_kind == "step_func":
442                if coeff_len == tlist_len-1 or coeff_len == tlist_len:
443                    pass
444                else:
445                    raise ValueError(
446                        "The length of tlist and coeff of the pulse "
447                        "labelled {} is invalid. "
448                        "It's either len(tlist)=len(coeff) or "
449                        "len(tlist)-1=len(coeff) for coefficients "
450                        "as step function".format(i))
451            else:
452                if coeff_len == tlist_len:
453                    pass
454                else:
455                    raise ValueError(
456                        "The length of tlist and coeff of the pulse "
457                        "labelled {} is invalid. "
458                        "It should be either len(tlist)=len(coeff)".format(i))
459        return True
460
461    def add_noise(self, noise):
462        """
463        Add a noise object to the processor
464
465        Parameters
466        ----------
467        noise: :class:`.Noise`
468            The noise object defined outside the processor
469        """
470        if isinstance(noise, Noise):
471            self.noise.append(noise)
472        else:
473            raise TypeError("Input is not a Noise object.")
474
475    def save_coeff(self, file_name, inctime=True):
476        """
477        Save a file with the control amplitudes in each timeslot.
478
479        Parameters
480        ----------
481        file_name: string
482            Name of the file.
483
484        inctime: bool, optional
485            True if the time list should be included in the first column.
486        """
487        self._is_pulses_valid()
488        coeffs = np.array(self.get_full_coeffs())
489        if inctime:
490            shp = coeffs.T.shape
491            data = np.empty((shp[0], shp[1] + 1), dtype=np.float64)
492            data[:, 0] = self.get_full_tlist()
493            data[:, 1:] = coeffs.T
494        else:
495            data = coeffs.T
496
497        np.savetxt(file_name, data, delimiter='\t', fmt='%1.16f')
498
499    def read_coeff(self, file_name, inctime=True):
500        """
501        Read the control amplitudes matrix and time list
502        saved in the file by `save_amp`.
503
504        Parameters
505        ----------
506        file_name: string
507            Name of the file.
508
509        inctime: bool, optional
510            True if the time list in included in the first column.
511
512        Returns
513        -------
514        tlist: array_like
515            The time list read from the file.
516
517        coeffs: array_like
518            The pulse matrix read from the file.
519        """
520        data = np.loadtxt(file_name, delimiter='\t')
521        if not inctime:
522            self.coeffs = data.T
523            return self.coeffs
524        else:
525            tlist = data[:, 0]
526            self.set_all_tlist(tlist)
527            self.coeffs = data[:, 1:].T
528            return self.get_full_tlist, self.coeffs
529
530    def get_noisy_pulses(self, device_noise=False, drift=False):
531        """
532        It takes the pulses defined in the `Processor` and
533        adds noise according to `Processor.noise`. It does not modify the
534        pulses saved in `Processor.pulses` but returns a new list.
535        The length of the new list of noisy pulses might be longer
536        because of drift Hamiltonian and device noise. They will be
537        added to the end of the pulses list.
538
539        Parameters
540        ----------
541        device_noise: bool, optional
542            If true, include pulse independent noise such as single qubit
543            Relaxation. Default is False.
544        drift: bool, optional
545            If true, include drift Hamiltonians. Default is False.
546
547        Returns
548        -------
549        noisy_pulses: list of :class:`.Pulse`
550            A list of noisy pulses.
551        """
552        pulses = deepcopy(self.pulses)
553        noisy_pulses = process_noise(
554            pulses, self.noise, self.dims, t1=self.t1, t2=self.t2,
555            device_noise=device_noise)
556        if drift:
557            noisy_pulses += [self.drift]
558        return noisy_pulses
559
560    def get_qobjevo(self, args=None, noisy=False):
561        """
562        Create a :class:`qutip.QobjEvo` representation of the evolution.
563        It calls the method `get_noisy_pulses` and create the `QobjEvo`
564        from it.
565
566        Parameters
567        ----------
568        args: dict, optional
569            Arguments for :class:`qutip.QobjEvo`
570        noisy: bool, optional
571            If noise are included. Default is False.
572
573        Returns
574        -------
575        qobjevo: :class:`qutip.QobjEvo`
576            The :class:`qutip.QobjEvo` representation of the unitary evolution.
577        c_ops: list of :class:`qutip.QobjEvo`
578            A list of lindblad operators is also returned. if ``noisy==Flase``,
579            it is always an empty list.
580        """
581        # TODO test it for non array-like coeff
582        # check validity
583        self._is_pulses_valid()
584
585        if args is None:
586            args = {}
587        else:
588            args = args
589        # set step function
590
591        if not noisy:
592            dynamics = self.pulses
593        else:
594            dynamics = self.get_noisy_pulses(
595                device_noise=True, drift=True)
596
597        qu_list = []
598        c_ops = []
599        for pulse in dynamics:
600            if noisy:
601                qu, new_c_ops = pulse.get_noisy_qobjevo(dims=self.dims)
602                c_ops += new_c_ops
603            else:
604                qu = pulse.get_ideal_qobjevo(dims=self.dims)
605            qu_list.append(qu)
606
607        final_qu = _merge_qobjevo(qu_list)
608        final_qu.args.update(args)
609
610        # bring all c_ops to the same tlist, won't need it in QuTiP 5
611        full_tlist = self.get_full_tlist()
612        temp = []
613        for c_op in c_ops:
614            temp.append(_merge_qobjevo([c_op], full_tlist))
615        c_ops = temp
616
617        if noisy:
618            return final_qu, c_ops
619        else:
620            return final_qu, []
621
622    def run_analytically(self, init_state=None, qc=None):
623        """
624        Simulate the state evolution under the given `qutip.QubitCircuit`
625        with matrice exponentiation. It will calculate the propagator
626        with matrix exponentiation and return a list of :class:`qutip.Qobj`.
627        This method won't include noise or collpase.
628
629        Parameters
630        ----------
631        qc: :class:`.QubitCircuit`, optional
632            Takes the quantum circuit to be implemented. If not given, use
633            the quantum circuit saved in the processor by ``load_circuit``.
634
635        init_state: :class:`qutip.Qobj`, optional
636            The initial state of the qubits in the register.
637
638        Returns
639        -------
640        evo_result: :class:`qutip.Result`
641            An instance of the class
642            :class:`qutip.Result` will be returned.
643        """
644        if init_state is not None:
645            U_list = [init_state]
646        else:
647            U_list = []
648        tlist = self.get_full_tlist()
649        # TODO replace this by get_complete_coeff
650        coeffs = self.get_full_coeffs()
651        for n in range(len(tlist)-1):
652            H = sum([coeffs[m, n] * self.ctrls[m]
653                    for m in range(len(self.ctrls))])
654            dt = tlist[n + 1] - tlist[n]
655            U = (-1j * H * dt).expm()
656            U = self.eliminate_auxillary_modes(U)
657            U_list.append(U)
658
659        try:  # correct_global_phase are defined for ModelProcessor
660            if self.correct_global_phase and self.global_phase != 0:
661                U_list.append(globalphase(self.global_phase, N=self.N))
662        except AttributeError:
663            pass
664
665        return U_list
666
667    def run(self, qc=None):
668        """
669        Calculate the propagator of the evolution by matrix exponentiation.
670        This method won't include noise or collpase.
671
672        Parameters
673        ----------
674        qc: :class:`.QubitCircuit`, optional
675            Takes the quantum circuit to be implemented. If not given, use
676            the quantum circuit saved in the processor by `load_circuit`.
677
678        Returns
679        -------
680        U_list: list
681            The propagator matrix obtained from the physical implementation.
682        """
683        if qc:
684            self.load_circuit(qc)
685        return self.run_analytically(qc=qc, init_state=None)
686
687    def run_state(self, init_state=None, analytical=False, states=None,
688                  noisy=True, solver="mesolve", **kwargs):
689        """
690        If `analytical` is False, use :func:`qutip.mesolve` to
691        calculate the time of the state evolution
692        and return the result. Other arguments of mesolve can be
693        given as keyword arguments.
694
695        If `analytical` is True, calculate the propagator
696        with matrix exponentiation and return a list of matrices.
697        Noise will be neglected in this option.
698
699        Parameters
700        ----------
701        init_state: Qobj
702            Initial density matrix or state vector (ket).
703
704        analytical: bool
705            If True, calculate the evolution with matrices exponentiation.
706
707        states: :class:`qutip.Qobj`, optional
708            Old API, same as init_state.
709
710        solver: str
711            "mesolve" or "mcsolve"
712
713        **kwargs
714            Keyword arguments for the qutip solver.
715
716        Returns
717        -------
718        evo_result: :class:`qutip.Result`
719            If ``analytical`` is False,  an instance of the class
720            :class:`qutip.Result` will be returned.
721
722            If ``analytical`` is True, a list of matrices representation
723            is returned.
724        """
725        if states is not None:
726            warnings.warn(
727                "states will be deprecated and replaced by init_state",
728                DeprecationWarning)
729        if init_state is None and states is None:
730            raise ValueError("Qubit state not defined.")
731        elif init_state is None:
732            # just to keep the old parameters `states`,
733            # it is replaced by init_state
734            init_state = states
735        if analytical:
736            if kwargs or self.noise:
737                raise warnings.warn(
738                    "Analytical matrices exponentiation"
739                    "does not process noise or"
740                    "any keyword arguments.")
741            return self.run_analytically(init_state=init_state)
742
743        # kwargs can not contain H or tlist
744        if "H" in kwargs or "tlist" in kwargs:
745            raise ValueError(
746                "`H` and `tlist` are already specified by the processor "
747                "and can not be given as a keyword argument")
748
749        # construct qobjevo for unitary evolution
750        if "args" in kwargs:
751            noisy_qobjevo, sys_c_ops = self.get_qobjevo(
752                    args=kwargs["args"], noisy=noisy)
753        else:
754            noisy_qobjevo, sys_c_ops = self.get_qobjevo(noisy=noisy)
755
756        # add collpase operators into kwargs
757        if "c_ops" in kwargs:
758            if isinstance(kwargs["c_ops"], (Qobj, QobjEvo)):
759                kwargs["c_ops"] += [kwargs["c_ops"]] + sys_c_ops
760            else:
761                kwargs["c_ops"] += sys_c_ops
762        else:
763            kwargs["c_ops"] = sys_c_ops
764
765        # choose solver:
766        if solver == "mesolve":
767            evo_result = mesolve(
768                H=noisy_qobjevo, rho0=init_state,
769                tlist=noisy_qobjevo.tlist, **kwargs)
770        elif solver == "mcsolve":
771            evo_result = mcsolve(
772                H=noisy_qobjevo, psi0=init_state,
773                tlist=noisy_qobjevo.tlist, **kwargs)
774
775        return evo_result
776
777    def load_circuit(self, qc):
778        """
779        Translate an :class:`.QubitCircuit` to its
780        corresponding Hamiltonians. (Defined in subclasses)
781        """
782        raise NotImplementedError("Use the function in the sub-class")
783
784    def eliminate_auxillary_modes(self, U):
785        """
786        Eliminate the auxillary modes like the cavity modes in cqed.
787        (Defined in subclasses)
788        """
789        return U
790
791    def get_operators_labels(self):
792        """
793        Get the labels for each Hamiltonian.
794        It is used in the method``plot_pulses``.
795        It is a 2-d nested list, in the plot,
796        a different color will be used for each sublist.
797        """
798        label_list = []
799        for pulse in self.pulses:
800            label_list.append(pulse.label)
801        return [label_list]
802
803    def plot_pulses(
804            self, title=None, figsize=(12, 6), dpi=None,
805            show_axis=False, rescale_pulse_coeffs=True,
806            num_steps=1000):
807        """
808        Plot the ideal pulse coefficients.
809
810        Parameters
811        ----------
812        title: str, optional
813            Title for the plot.
814
815        figsize: tuple, optional
816            The size of the figure.
817
818        dpi: int, optional
819            The dpi of the figure.
820
821        show_axis: bool, optional
822            If the axis are shown.
823
824        rescale_pulse_coeffs: bool, optional
825            Rescale the hight of each pulses.
826
827        num_steps: int, optional
828            Number of time steps in the plot.
829
830        Returns
831        -------
832        fig: matplotlib.figure.Figure
833            The `Figure` object for the plot.
834
835        ax: matplotlib.axes._subplots.AxesSubplot
836            The axes for the plot.
837
838        Notes
839        -----
840        ``plot_pulses`` only works for array_like coefficients
841        """
842        import matplotlib.pyplot as plt
843        import matplotlib.gridspec as gridspec
844        color_list = plt.rcParams['axes.prop_cycle'].by_key()['color']
845
846        # create a axis for each pulse
847        fig = plt.figure(figsize=figsize, dpi=dpi)
848        grids = gridspec.GridSpec(len(self.pulses), 1)
849        grids.update(wspace=0., hspace=0.)
850
851        tlist = np.linspace(0., self.get_full_tlist()[-1], num_steps)
852        dt = tlist[1] - tlist[0]
853
854        # make sure coeffs start and end with zero, for ax.fill
855        tlist = np.hstack(([-dt*1.e-20], tlist, [tlist[-1] + dt*1.e-20]))
856        coeffs = []
857        for pulse in self.pulses:
858            coeffs.append(_pulse_interpolate(pulse, tlist))
859
860        pulse_ind = 0
861        axis = []
862        for i, label_group in enumerate(self.get_operators_labels()):
863            for j, label in enumerate(label_group):
864                grid = grids[pulse_ind]
865                ax = plt.subplot(grid)
866                axis.append(ax)
867                ax.fill(tlist, coeffs[pulse_ind], color_list[i], alpha=0.7)
868                ax.plot(tlist, coeffs[pulse_ind], color_list[i])
869                if rescale_pulse_coeffs:
870                    ymax = np.max(np.abs(coeffs[pulse_ind])) * 1.1
871                else:
872                    ymax = np.max(np.abs(coeffs)) * 1.1
873                if ymax != 0.:
874                    ax.set_ylim((-ymax, ymax))
875
876                # disable frame and ticks
877                if not show_axis:
878                    ax.set_xticks([])
879                    ax.spines['bottom'].set_visible(False)
880                ax.spines['top'].set_visible(False)
881                ax.spines['right'].set_visible(False)
882                ax.spines['left'].set_visible(False)
883                ax.set_yticks([])
884                ax.set_ylabel(label, rotation=0)
885                pulse_ind += 1
886                if i == 0 and j == 0 and title is not None:
887                    ax.set_title(title)
888        fig.tight_layout()
889        return fig, axis
890
891
892def _pulse_interpolate(pulse, tlist):
893    """
894    A function that calls Scipy interpolation routine. Used for plotting.
895    """
896    if pulse.tlist is None and pulse.coeff is None:
897        coeff = np.zeros(len(tlist))
898        return coeff
899    if isinstance(pulse.coeff, bool):
900        if pulse.coeff:
901            coeff = np.ones(len(tlist))
902        else:
903            coeff = np.zeros(len(tlist))
904        return coeff
905    coeff = pulse.coeff
906    if len(coeff) == len(pulse.tlist)-1:  # for discrete pulse
907        coeff = np.concatenate([coeff, [0]])
908
909    from scipy import interpolate
910    if pulse.spline_kind == "step_func":
911        kind = "previous"
912    else:
913        kind = "cubic"
914    inter = interpolate.interp1d(
915        pulse.tlist, coeff, kind=kind,
916        bounds_error=False, fill_value=0.0)
917    return inter(tlist)
918