class BuilderBase(object):
    """
    Class for building the system for stationary transport calculations.

    For descriptions of all attributes use help(Builder).
    """

    def __init__(self,
                 nsingle=0, hsingle={}, coulomb={},
                 nleads=0, tleads={}, mulst={}, tlst={}, dband={},
                 indexing=None, kpnt=None,
                 kerntype='Pauli', symq=True, norm_row=0, solmethod=None,
                 itype=0, dqawc_limit=10000, mfreeq=False, phi0_init=None,
                 mtype_qd=complex, mtype_leads=complex,
                 symmetry=None, herm_hs=True, herm_c=False, m_less_n=True):

        self._init_copy_data(locals())
        self._init_validate_data()
        self._init_set_globals()
        self._init_set_approach_class()
        self._init_create_setup()
        self._init_create_appr()

    def _init_copy_data(self, data):
        self.data = ModelParameters(data)

    def _init_validate_data(self):
        data = self.data
        data.itype = validate_itype(data.itype, data.kerntype)
        data.kerntype = validate_kerntype(data.kerntype)
        data.indexing, data.symmetry = validate_indexing(data.indexing,
                                          data.symmetry,
                                          data.kerntype)

    def _init_set_globals(self):
        self.globals = globals()

    def _init_set_approach_class(self):
        kerntype = self.data.kerntype
        if isinstance(kerntype, str):
            approach_string = kerntype[0].capitalize() + kerntype[1:]
            self.Approach = self.globals['Approach'+approach_string]
        elif issubclass(kerntype, Approach):
            self.Approach = kerntype
            self.kerntype = self.Approach.kerntype

    def _init_create_setup(self):
        data = self.data
        self.funcp = FunctionProperties(symq=data.symq, norm_row=data.norm_row, solmethod=data.solmethod,
                                        itype=data.itype, dqawc_limit=data.dqawc_limit,
                                        mfreeq=data.mfreeq, phi0_init=data.phi0_init,
                                        mtype_qd=data.mtype_qd, mtype_leads=data.mtype_leads,
                                        kpnt=data.kpnt, dband=data.dband)

        icn = self.Approach.indexing_class_name
        self.si = self.globals[icn](data.nsingle, data.indexing, data.symmetry)
        self.qd = QuantumDot(data.hsingle, data.coulomb, self.si,
                             data.herm_hs, data.herm_c, data.m_less_n, data.mtype_qd)
        self.leads = LeadsTunneling(data.nleads, data.tleads, self.si,
                                    data.mulst, data.tlst, data.dband, data.mtype_leads)

    def _init_create_appr(self):
        self.appr = self.Approach(self)

    def __getattr__(self, item):
        sub_class_str = attribute_map.get(item)
        if sub_class_str is None:
            return super(BuilderBase, self).__getattribute__(item)
        else:
            sub_class = getattr(self, attribute_map[item])
            return getattr(sub_class, item)

    def __setattr__(self, item, value):
        sub_class_str = attribute_map.get(item)
        if sub_class_str is None:
            super(BuilderBase, self).__setattr__(item, value)
        else:
            sub_class = getattr(self, attribute_map[item])
            setattr(sub_class, item, value)

    def change_si(self):
        si = self.si
        icn = self.Approach.indexing_class_name
        # noinspection PyTypeHints
        if not isinstance(si, self.globals[icn]):
            self.si = self.globals[icn](si.nsingle, si.indexing, si.symmetry, si.nleads)
            self.qd.si = self.si
            self.leads.si = self.si
        self.si.states_changed = True

    # kerntype
    def get_kerntype(self):
        return self.appr.kerntype

    def set_kerntype(self, value):
        if isinstance(value, str):
            if self.appr.kerntype != value:
                approach_string = value[0].capitalize() + value[1:]
                self.Approach = self.globals['Approach'+approach_string]
                self.change_si()
                self.appr = self.Approach(self)
        else:
            if issubclass(value, Approach):
                self.Approach = value
                self.change_si()
                self.appr = self.Approach(self)
    kerntype = property(get_kerntype, set_kerntype)

    # indexing
    def get_indexing(self):
        return self.si.indexing

    def set_indexing(self, value):
        if self.si.indexing != value:
            print('WARNING: Cannot change indexing from \''+self.si.indexing+'\' to \''+value
                  + '\' Consider constructing a new system using \'indexing='+value+'\'')
    indexing = property(get_indexing, set_indexing)

    # dband
    def get_dband(self):
        return self.funcp.dband

    def set_dband(self, value):
        self.funcp.dband = value
        self.leads.change(dlst=value)
    dband = property(get_dband, set_dband)

    def add(self, hsingle=None, coulomb=None, tleads=None, mulst=None, tlst=None, dlst=None):
        """
        Adds the values to the specified dictionaries and correspondingly redefines
        relevant many-body properties of the system.

        Parameters
        ----------
        hsingle,coulomb,tleads,mulst,tlst,dlst : dict
            Dictionaries describing what values to add.
            For example, tleads[(lead, state)] = value to add.
        """
        if not (hsingle is None and coulomb is None):
            self.qd.add(hsingle, coulomb)
        if not (tleads is None and mulst is None and tlst is None and dlst is None):
            self.leads.add(tleads, mulst, tlst, dlst)

    def change(self, hsingle=None, coulomb=None, tleads=None, mulst=None, tlst=None, dlst=None):
        """
        Changes the values of the specified dictionaries and correspondingly redefines
        relevant many-body properties of the system.

        Parameters
        ----------
        hsingle,coulomb,tleads,mulst,tlst,dlst : dict
            Dictionaries describing what values to change.
            For example, tleads[(lead, state)] = value to change.
        """
        if not (hsingle is None and coulomb is None):
            self.qd.change(hsingle, coulomb)
        if not (tleads is None and mulst is None and tlst is None and dlst is None):
            self.leads.change(tleads, mulst, tlst, dlst)

    def get_phi0(self, b, bp):
        """
        Get the reduced density matrix element corresponding to
        many-body states b and bp.
        """
        return get_phi0(self, b, bp)

    def get_phi1(self, l, c, b):
        """
        Get the energy integrated current amplitudes corresponding to
        lead l and many-body states c and b.
        """
        return get_phi1(self, l, c, b)

    def print_state(self, b, eps=0.0, prntq=True, filename=None, separator=''):
        """
        Prints properties of given many-body eigenstate of the quantum dot Hamiltonian
        """
        print_state(self, b, eps, prntq, filename, separator)

    def print_all_states(self, filename, eps=0.0, separator='', mode='w'):
        """
        Prints properties of all many-body eigenstates to a file.
        """
        print_all_states(self, filename, eps, separator, mode)

    def sort_eigenstates(self, srt=None):
        """
        Sort many-body states of the system by given order of properties.
        """
        sort_eigenstates(self, srt=srt)

    def remove_coherences(self, dE):
        """
        Remove the coherences with energy difference larger than dE.
        """
        remove_coherences(self, dE)

    def remove_states(self, dE):
        """
        Remove the states with energy dE larger than the ground state
        for the transport calculations.
        """
        remove_states(self, dE)

    def remove_fock_states(self, lin_state_indices):
        """Remove Fock states."""
        self.si.remove_fock_states(lin_state_indices)
        self.qd._init_hamiltonian()
        self.leads._init_coupling()

    def use_all_states(self):
        """
        Use all states for the transport calculations.
        """
        use_all_states(self) class ModelParameters(object):
class BuilderBase(object):
class BuilderManyBody(BuilderBase):
    def __init__(self, params):
    def __init__(self,
                 nsingle=0, hsingle={}, coulomb={},
                 nleads=0, tleads={}, mulst={}, tlst={}, dband={},
                 indexing=None, kpnt=None,
                 kerntype='Pauli', symq=True, norm_row=0, solmethod=None,
                 itype=0, dqawc_limit=10000, mfreeq=False, phi0_init=None,
                 mtype_qd=complex, mtype_leads=complex,
                 symmetry=None, herm_hs=True, herm_c=False, m_less_n=True):
    def _init_copy_data(self, data):
    def _init_validate_data(self):
    def _init_set_globals(self):
    def _init_set_approach_class(self):
    def _init_create_setup(self):
    def _init_create_appr(self):
    def __getattr__(self, item):
    def __setattr__(self, item, value):
    def change_si(self):
    def get_kerntype(self):
    def set_kerntype(self, value):
    def get_indexing(self):
    def set_indexing(self, value):
    def get_dband(self):
    def set_dband(self, value):
    def add(self, hsingle=None, coulomb=None, tleads=None, mulst=None, tlst=None, dlst=None):
    def change(self, hsingle=None, coulomb=None, tleads=None, mulst=None, tlst=None, dlst=None):
    def get_phi0(self, b, bp):
    def get_phi1(self, l, c, b):
    def print_state(self, b, eps=0.0, prntq=True, filename=None, separator=''):
    def print_all_states(self, filename, eps=0.0, separator='', mode='w'):
    def sort_eigenstates(self, srt=None):
    def remove_coherences(self, dE):
    def remove_states(self, dE):
    def remove_fock_states(self, lin_state_indices):
    def use_all_states(self):
    def __init__(self,
                 Ea=None, Na=[0], Tba=None,
                 mulst={}, tlst={}, dband={}, kpnt=None,
                 kerntype='Pauli', symq=True, norm_row=0, solmethod=None,
                 itype=0, dqawc_limit=10000, mfreeq=False, phi0_init=None,
                 mtype_qd=complex, mtype_leads=complex,
                 symmetry=None, herm_hs=True, herm_c=False, m_less_n=True):
    def _init_state_indexing(self, Na, Ea): class BuilderManyBody(BuilderBase):
    """
    Class for building the system for stationary transport calculations,
    using many-body states as an input.

    For missing descriptions of attributes use help(Builder).

    Attributes
    ----------
    Ea : array
        nmany by 1 array containing many-body Hamiltonian eigenvalues.
    Na : array
        nmany by 1 array containing particle numbers of many-body states.
    Tba : array
        nleads by nmany by nmany array, which contains many-body tunneling amplitude matrix.
    """

    def __init__(self,
                 Ea=None, Na=[0], Tba=None,
                 mulst={}, tlst={}, dband={}, kpnt=None,
                 kerntype='Pauli', symq=True, norm_row=0, solmethod=None,
                 itype=0, dqawc_limit=10000, mfreeq=False, phi0_init=None,
                 mtype_qd=complex, mtype_leads=complex,
                 symmetry=None, herm_hs=True, herm_c=False, m_less_n=True):

        nleads = Tba.shape[0] if Tba is not None else 0

        # noinspection PyPep8
        BuilderBase.__init__(self,
            nleads=nleads, mulst=mulst, tlst=tlst, dband=dband, kpnt=kpnt,
            kerntype=kerntype, symq=symq, norm_row=norm_row, solmethod=solmethod,
            itype=itype, dqawc_limit=dqawc_limit, mfreeq=mfreeq, phi0_init=phi0_init,
            mtype_qd=mtype_qd, mtype_leads=mtype_leads,
            symmetry=symmetry, herm_hs=herm_hs, herm_c=herm_c, m_less_n=m_less_n,
            indexing='charge')

        self._init_state_indexing(Na, Ea)

        self.qd.Ea = Ea
        self.leads.Tba = Tba

    def _init_state_indexing(self, Na, Ea):
        Na = np.array(Na, dtype=int)
        nmin, nmax = Na.min(), Na.max()
        ncharge = nmax-nmin+1
        nmany = len(Ea) if Ea is not None else 0

        statesdm = [[] for _ in range(ncharge)]
        for i in range(nmany):
            statesdm[Na[i]].append(i)

        self.Na = Na
        self.si.nmany = nmany
        self.si.ncharge = ncharge
        self.si.shiftlst0 = np.zeros(ncharge+1, dtype=longnp)
        self.si.shiftlst1 = np.zeros(ncharge, dtype=longnp)
        self.si.lenlst = np.zeros(ncharge, dtype=longnp)
        self.si.dictdm = np.zeros(nmany, dtype=longnp)
        self.si.set_statesdm(statesdm) class StateIndexingDM(StateIndexing):
    """
    Class for indexing density matrix elements.
    Separates the off-diagonal density matrix elements into real and imaginary parts.
    Derived from StateIndexing.

    Attributes
    ----------
    ndm0 : int
        Number of density matrix elements to zeroth order corresponding to statesdm,
        when symmetries are used.
    ndm0r : int
        Number of unique density matrix elements to zeroth order corresponding to statesdm,
        when symmetries are used.
        The real and imaginary parts of the off-diagonal matrix element are considered
        as separate entities.
    ndm0_ : int
        Number of density matrix elements to zeroth order corresponding to statesdm,
        when no symmetries are used.
    ndm0_tot : int
        Total number of density matrix elements to zeroth order for given nsingle,
        when no symmetries are used.
    ndm1, ndm1_ : int, int
        Number of density matrix elements to first order corresponding to statesdm,
        when no symmetries are used.
    ndm1_tot : int
        Total number of density matrix elements to first order for given nsingle,
        when no symmetries are used.
    npauli_ : int
        Number of diagonal density matrix elements to zeroth order corresponding to statesdm,
        when no symmetries are used.
    npauli : int
        Number of diagonal density matrix elements to zeroth order corresponding to statesdm,
        when symmetries are used.
    statesdm : list
        List of states considered in calculations involving density matrix.
    shiftlst0 : numpy array
        Array necessary for labeling zeroth order density matrix elements. See method get_ind_dm0().
    shiftlst1 : numpy array
        Array necessary for labeling first order density matrix elements. See method get_ind_dm1().
    lenlst : numpy array
        Array containing number of many-body states for each charge. lenlst[charge]=number of states.
    dictdm : numpy array
        Array giving enumeration of considered many-body states in statesdm.
    inddm0 : dict
        Dictionary between index of density matrix element and pair of many-body states (a,b).
    mapdm0 : list
        List showing which density matrix elements are mapped to each other due to symmetries.
        For example mapdm0[1]=1, mapdm0[2]=1, shows that 1 and 2 are equivalent.
        Also mapdm0[ind]=-1 shows that the state is excluded.
    booldm0 : list
        List giving density matrix elements which will be used.
        Other elements will be mapped to elements, which have booldm0[ind]=True.
        From example for mapdm0 we have booldm0[1]=True, booldm0[2]=False.
    conjdm0 : list
        List showing, which density matrix elements are complex conjugate and are not unique.
    """

    def __init__(self, nsingle, indexing='Lin', symmetry=None, nleads=0):
        """
        Initialization of the StateIndexingDM class

        Parameters
        ----------
        nsingle : int
            Number of single-particle states.
        indexing : str
            String determining type of the indexing. Possible values are 'Lin', 'charge', 'sz', 'ssq'.
        symmetry : str
            String determining that the states will be augmented by the symmetry.
        """
        StateIndexing.__init__(self, nsingle, indexing, symmetry, nleads)
        self.ndm0_tot = int(factorial(2*self.nsingle)/factorial(self.nsingle)**2)
        self.ndm1_tot = int(self.nsingle/(self.nsingle+1)*self.ndm0_tot)
        #
        self.shiftlst0 = np.zeros(self.ncharge+1, dtype=longnp)
        self.shiftlst1 = np.zeros(self.ncharge, dtype=longnp)
        self.lenlst = np.zeros(self.ncharge, dtype=longnp)
        self.dictdm = np.zeros(self.nmany, dtype=longnp)
        #
        self.statesdm, self.mapdm0, self.mapdm1 = None, None, None
        self.set_statesdm(self.chargelst)

    def set_statesdm(self, statesdm):
        """
        Sets the considered many-body states for density matrix calculations.

        Parameters
        ----------
        statesdm: list
            List containing indices of many-body state under consideration.
        """
        self.statesdm = list(statesdm)
        self.statesdm.append([])
        self.ndm0_, self.ndm1_, self.npauli_ = 0, 0, 0
        for j1 in range(self.ncharge):
            self.npauli_ += len(statesdm[j1])
            self.ndm0_ += len(statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.ndm1_ += len(statesdm[j1])*len(statesdm[j1+1])
        self.set_dictdm()

        self.states_changed = True

    def set_dictdm(self):
        """
        Makes dictdm, shiftlst0, and shiftlst1 necessary for density-matrix element indexing.
        """
        for j1 in range(self.ncharge):
            self.lenlst[j1] = len(self.statesdm[j1])
            self.shiftlst0[j1+1] = self.shiftlst0[j1] + len(self.statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.shiftlst1[j1+1] = (self.shiftlst1[j1]
                                        + len(self.statesdm[j1])*len(self.statesdm[j1+1]))
            counter = 0
            for j2 in self.statesdm[j1]:
                self.dictdm[j2] = counter
                counter += 1
        self.set_mapdm()

    def set_mapdm(self):
        """
        Reduce the number of diagonal matrix elements by using symmetries.
        """
        # noinspection PyShadowingNames
        def add_elem(counter, b, bp, charge, dictq=True, conjq=True):
            bbp = self.get_ind_dm0(b, bp, charge, maptype=0)
            self.mapdm0[bbp] = counter
            self.booldm0[bbp] = dictq
            self.conjdm0[bbp] = conjq
            if dictq:
                self.inddm0.update({counter: (b, bp)})
            return counter+1
        #
        self.mapdm0 = np.ones(self.ndm0_, dtype=longnp)*(-1)
        self.booldm0 = np.zeros(self.ndm0_, dtype=boolnp)
        self.conjdm0 = np.zeros(self.ndm0_, dtype=boolnp)
        self.inddm0 = {}
        counter = 0
        # Diagonal density matrix elements
        for charge in range(self.ncharge):
            for b in self.statesdm[charge]:
                if self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    if b_ssq == -b_sz:
                        add_elem(counter, b, b, charge)
                        for sz in range(b_sz+2, b_ssq+1, 2):
                            b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                            add_elem(counter, b1, b1, charge, dictq=False)
                        counter = counter+1
                else:
                    counter = add_elem(counter, b, b, charge)
        self.npauli = counter
        # Off-diagonal density matrix elements
        for charge in range(self.ncharge):
            for b, bp in itertools.combinations(self.statesdm[charge], 2):
                if self.indexing == 'sz':
                    if self.ind_qn[b][1] == self.ind_qn[bp][1]:
                        add_elem(counter, bp, b, charge, dictq=False, conjq=False)
                        counter = add_elem(counter, b, bp, charge)
                elif self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    (bp_ch, bp_sz, bp_ssq, bp_alpha) = self.ind_qn[bp]
                    if (b_sz, b_ssq) == (bp_sz, bp_ssq):
                        if b_ssq == -b_sz:
                            add_elem(counter, bp, b, charge, dictq=False, conjq=False)
                            add_elem(counter, b, bp, charge)
                            for sz in range(b_sz+2, b_ssq+1, 2):
                                b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                                b1p = self.qn_ind[(bp_ch, sz, bp_ssq, bp_alpha)]
                                add_elem(counter, b1p, b1, charge, dictq=False, conjq=False)
                                add_elem(counter, b1, b1p, charge, dictq=False)
                            counter = counter+1
                else:
                    add_elem(counter, bp, b, charge, dictq=False, conjq=False)
                    counter = add_elem(counter, b, bp, charge)
        self.ndm0 = counter
        self.ndm0r = self.npauli+2*(self.ndm0-self.npauli)
        self.ndm1 = self.ndm1_

    def get_ind_dm0(self, b, bp, charge, maptype=1):
        """
        Get the index of zeroth order density matrix element.

        Parameters
        ----------
        b,bp : int
            Indices of many-body states within the same charge state.
        charge : int
            Charge of b and bp.
        maptype : int
            Determines what kind of mapping to use when returning the index.

        Returns
        -------
        int
            Index of the zeroth order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[b]
        # j = self.dictdm[bp]
        # shiftas = self.shiftlst0[charge]
        # print('index is', l*i + j + shiftas)
        if maptype == 0:
            return (self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge])
        elif maptype == 1:
            return (self.mapdm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])
        elif maptype == 2:
            return (self.booldm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])
        elif maptype == 3:
            return (self.conjdm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])

    def get_ind_dm1(self, c, b, bcharge):
        """
        Get the index of first order density matrix element.

        Parameters
        ----------
        c,b : int
            Indices of many-body states differing by single charge, charge(c)>charge(b).
        bcharge : int
            Charge of many-body state b.

        Returns
        -------
        int
            Index of the first order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[c]
        # j = self.dictdm[b]
        # shiftas = self.shiftlst1[charge]
        # print('index is', l*i + j + shiftas)
        #       b
        #     -----
        #   c |   |
        #     -----
        return self.lenlst[bcharge]*self.dictdm[c] + self.dictdm[b] + self.shiftlst1[bcharge]

    def remove_fock_states(self, lin_state_indices):
        StateIndexing.remove_fock_states(self, lin_state_indices)
        self.set_statesdm(self.chargelst) class StateIndexingDMc(StateIndexing):
    """
    Class for indexing density matrix elements.
    Does not separates the off-diagonal density matrix elements into real and imaginary parts,
    and treats Phi[0]_{b,bp} and Phi[0]_{bp,b} as separate entities.
    Derived from StateIndexing.

    Attributes
    ----------
    Same as in StateIndexingDM.
    """

    def __init__(self, nsingle, indexing='Lin', symmetry=None, nleads=0):
        """
        Initialization of the StateIndexingDMc class

        Parameters
        ----------
        nsingle : int
            Number of single-particle states.
        indexing : str
            String determining type of the indexing. Possible values are 'Lin', 'charge', 'sz', 'ssq'.
        symmetry : str
            String determining that the states will be augmented by the symmetry.
        """
        StateIndexing.__init__(self, nsingle, indexing, symmetry, nleads)
        self.ndm0_tot = int(factorial(2*self.nsingle)/factorial(self.nsingle)**2)
        self.ndm1_tot = int(self.nsingle/(self.nsingle+1)*self.ndm0_tot)
        #
        self.shiftlst0 = np.zeros(self.ncharge+1, dtype=longnp)
        self.shiftlst1 = np.zeros(self.ncharge, dtype=longnp)
        self.lenlst = np.zeros(self.ncharge, dtype=longnp)
        self.dictdm = np.zeros(self.nmany, dtype=longnp)
        #
        self.statesdm, self.mapdm0, self.mapdm1 = None, None, None
        self.set_statesdm(self.chargelst)

    def set_statesdm(self, statesdm):
        """
        Sets the considered many-body states for density matrix calculations.

        Parameters
        ----------
        statesdm: list
            List containing indices of many-body state under consideration.
        """
        self.statesdm = statesdm
        self.statesdm.append([])
        self.ndm0_, self.ndm1_, self.npauli_ = 0, 0, 0
        for j1 in range(self.ncharge):
            self.npauli_ += len(statesdm[j1])
            self.ndm0_ += len(statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.ndm1_ += len(statesdm[j1])*len(statesdm[j1+1])
        self.set_dictdm()

        self.states_changed = True

    def set_dictdm(self):
        """
        Makes dictdm, shiftlst0, and shiftlst1 necessary for density-matrix element indexing.
        """
        for j1 in range(self.ncharge):
            self.lenlst[j1] = len(self.statesdm[j1])
            self.shiftlst0[j1+1] = self.shiftlst0[j1] + len(self.statesdm[j1])**2
            if j1 < self.ncharge-1:
                self.shiftlst1[j1+1] = (self.shiftlst1[j1]
                                        + len(self.statesdm[j1])*len(self.statesdm[j1+1]))
            counter = 0
            for j2 in self.statesdm[j1]:
                self.dictdm[j2] = counter
                counter += 1
        self.set_mapdm()

    def set_mapdm(self):
        """
        Reduce the number of diagonal matrix elements by using symmetries.
        """
        # noinspection PyShadowingNames
        def add_elem(counter, b, bp, charge, dictq=True):
            bbp = self.get_ind_dm0(b, bp, charge, maptype=0)
            self.mapdm0[bbp] = counter
            self.booldm0[bbp] = dictq
            if dictq:
                self.inddm0.update({counter: (b, bp)})
            return counter+1
        #
        self.mapdm0 = np.ones(self.ndm0_, dtype=longnp)*(-1)
        self.booldm0 = np.zeros(self.ndm0_, dtype=boolnp)
        self.conjdm0 = None
        self.inddm0 = {}
        counter = 0
        # Diagonal density matrix elements
        for charge in range(self.ncharge):
            for b in self.statesdm[charge]:
                if self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    if b_ssq == -b_sz:
                        add_elem(counter, b, b, charge)
                        for sz in range(b_sz+2, b_ssq+1, 2):
                            b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                            add_elem(counter, b1, b1, charge, dictq=False)
                        counter = counter+1
                else:
                    counter = add_elem(counter, b, b, charge)
        self.npauli = counter
        # self.npauli = 0
        # Off-diagonal density matrix elements
        for charge in range(self.ncharge):
            # for b, bp in itertools.product(self.statesdm[charge], self.statesdm[charge]):
            for b, bp in itertools.permutations(self.statesdm[charge], 2):
                if self.indexing == 'sz':
                    if self.ind_qn[b][1] == self.ind_qn[bp][1]:
                        counter = add_elem(counter, b, bp, charge)
                elif self.indexing == 'ssq':
                    (b_ch,  b_sz,  b_ssq,  b_alpha) = self.ind_qn[b]
                    (bp_ch, bp_sz, bp_ssq, bp_alpha) = self.ind_qn[bp]
                    if (b_sz, b_ssq) == (bp_sz, bp_ssq):
                        if b_ssq == -b_sz:
                            add_elem(counter, b, bp, charge)
                            for sz in range(b_sz+2, b_ssq+1, 2):
                                b1 = self.qn_ind[(b_ch,  sz, b_ssq,  b_alpha)]
                                b1p = self.qn_ind[(bp_ch, sz, bp_ssq, bp_alpha)]
                                add_elem(counter, b1, b1p, charge, dictq=False)
                            counter = counter+1
                else:
                    counter = add_elem(counter, b, bp, charge)
        self.ndm0 = counter
        self.ndm0r = self.npauli+2*(self.ndm0-self.npauli)
        self.ndm1 = self.ndm1_

    def get_ind_dm0(self, b, bp, charge, maptype=1):
        """
        Get the index of zeroth order density matrix element.

        Parameters
        ----------
        b,bp : int
            Indices of many-body states within the same charge state.
        charge : int
            Charge of b and bp.
        maptype : int
            Determines what kind of mapping to use when returning the index.

        Returns
        -------
        int
            Index of the zeroth order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[b]
        # j = self.dictdm[bp]
        # shiftas = self.shiftlst0[charge]
        # print('index is', l*i + j + shiftas)
        if maptype == 0:
            return self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp] + self.shiftlst0[charge]
        elif maptype == 1:
            return (self.mapdm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])
        elif maptype == 2:
            return (self.booldm0[self.lenlst[charge]*self.dictdm[b] + self.dictdm[bp]
                    + self.shiftlst0[charge]])

    def get_ind_dm1(self, c, b, bcharge):
        """
        Get the index of first order density matrix element.

        Parameters
        ----------
        c,b : int
            Indices of many-body states differing by single charge, charge(c)>charge(b).
        bcharge : int
            Charge of many-body state b.

        Returns
        -------
        int
            Index of the first order density matrix element.
        """
        # l = len(self.statesdm[charge])
        # i = self.dictdm[c]
        # j = self.dictdm[b]
        # shiftas = self.shiftlst1[charge]
        # print('index is', l*i + j + shiftas)
        #       b
        #     -----
        #   c |   |
        #     -----
        return self.lenlst[bcharge]*self.dictdm[c] + self.dictdm[b] + self.shiftlst1[bcharge]

    def remove_fock_states(self, lin_state_indices):
        StateIndexing.remove_fock_states(self, lin_state_indices)
        self.set_statesdm(self.chargelst) class PhononBaths(object):
    """
    Class for defining electron-phonon couplings from baths to the quantum dot.

    Attributes
    ----------
    nbaths : int
        Number of the phonon baths.
    velph : list, dict, or ndarray
        Dictionary, list or numpy array defining single-particle electron-phonon couplings.
        numpy array has to be nbaths by nsingle by nsingle.
    si : StateIndexing
        StateIndexing or StateIndexingDM object.
    tlst_ph : dict, list or ndarray
        Dictionary, list or numpy array containing temperatures of the phonon baths.
    dlst_ph : dict, list or ndarray
        Dictionary, list or numpy array containing bandwidths of the phonon baths.
    bath_func : list
        List of length nbaths containing density of states functions for the phonon baths.
    mtype : type
        Defines type of Vbbp0 and Vbbp matrices. For example, float, complex, etc.
    Vbbp0 : ndarray
        nmany by nmany array, which contains many-body electron-phonon coupling matrix in Fock basis.
    Vbbp : list
        nmany by nmany array, which contains many-body electron-phonon coupling matrix,
        which is used in calculations.
    """

    def __init__(self, nbaths, velph, si, tlst_ph, dlst_ph, bath_func=None, mtype=complex):
        """Initialization of the LeadsTunneling class."""
        si.nbaths = nbaths
        #
        self.si = si
        self.velph = make_velph_dict(velph, si)
        self.si.nbaths = nbaths
        self.tlst_ph = make_array(None, tlst_ph, si, nbaths, False)
        self.dlst_ph = make_array_dlst(None, dlst_ph, si, nbaths, False)
        self.bath_func = bath_func
        self.mtype = mtype
        self._init_coupling()

    def _init_coupling(self):
        self.Vbbp0 = elph_construct_Vbbp(self, self.velph)
        self.Vbbp = np.array(self.Vbbp0)

    def add(self, velph=None, tlst_ph=None, dlst_ph=None, updateq=True, lstq=True):
        """
        Adds a value to single particle electron-phonon couplings
        and correspondingly redefines many-body matrix Vbbp.

        Parameters
        ----------
        velph : dict
            Dictionary describing what values to add.
            For example, velph[(bath, state1, state2)] = value to add.
        tlst_ph : dict, list or ndarray
            Dictionary, list or numpy array containing temperatures of the phonon baths.
        dlst_ph : dict, list or ndarray
            Dictionary, list or numpy array containing bandwidths of the phonon baths.
        updateq : bool
            Specifies if the values of the single particle couplings will be updated.
            The many-body couplings Vbbp will be updates in either case.
        lstq : bool
            Determines if the values will be added to tlst_ph, dlst_ph.
        """
        if lstq:
            if tlst_ph is not None:
                self.tlst_ph += make_array(None, tlst_ph, self.si, self.si.nbaths, False)
            if dlst_ph is not None:
                self.dlst_ph += make_array_dlst(None, dlst_ph, self.si, self.si.nbaths, False)
        if velph is not None:
            if updateq:
                velph = make_velph_dict(velph, self.si)
                for j0 in velph:
                    if j0 in self.velph:
                        self.velph[j0] += velph[j0]
                    else:
                        self.velph.update({j0: velph[j0]})
            self.Vbbp0 = elph_construct_Vbbp(self, velph, self.Vbbp0)

    def change(self, velph=None, tlst_ph=None, dlst_ph=None):
        """
        Changes the values of the single particle electron-phonon couplings
        and correspondingly redefines many-body coupling matrix Vbbp.

        Parameters
        ----------
        velph : dict
            Dictionary describing which electron-phonon couplings to change.
            For example, velph[(bath, state1, state2)] = the new value.
        tlst_ph : dict, list or ndarray
            Dictionary, list or numpy array containing temperatures of the phonon baths.
        dlst_ph : dict, list or ndarray
            Dictionary, list or numpy array containing bandwidths of the phonon baths.
        """
        if tlst_ph is not None:
            self.tlst_ph[:] = make_array(self.tlst_ph, tlst_ph, self.si, self.si.nbaths, False)
        if dlst_ph is not None:
            self.dlst_ph[:] = make_array_dlst(self.dlst_ph, dlst_ph, self.si, self.si.nbaths, False)
        #
        if velph is not None:
            velph = make_velph_dict(velph, self.si, True)
            # Find the differences from the previous electron-phonon coupling
            velph_add = {}
            for j0 in velph:
                if j0 in self.velph:
                    velph_diff = velph[j0]-self.velph[j0]
                    if velph_diff != 0:
                        velph_add.update({j0: velph_diff})
                        self.velph[j0] += velph_diff
                else:
                    velph_diff = velph[j0]
                    if velph_diff != 0:
                        velph_add.update({j0: velph_diff})
                        self.velph.update({j0: velph_diff})
            # Add the differences
            self.add(velph_add, updateq=False, lstq=False)

    def rotate(self, vecslst, indexing=None):
        """
        Rotates electron-phonon coupling matrix Vbbp0 in Fock basis to Vbbp,
        which is in eigenstate basis of the quantum dot.

        Parameters
        ----------
        vecslst : list of ndarrays
            List of size ncharge containing arrays defining eigenvector matrices for given charge.
        indexing : string
            Specifies what kind of rotation procedure to use. Default is si.indexing.
        """
        self.Vbbp[:] = elph_rotate_Vbbp(self.Vbbp0, vecslst, self.si, indexing, self.mtype)

    def use_Vbbp0(self):
        """
        Sets the Vbbp matrix for calculation to Vbbp0 in the Fock basis.
        """
        self.Vbbp[:] = self.Vbbp0

    def update_Vbbp0(self, nbaths, velph, mtype=complex):
        """
        Updates the Vbbp0 in the Fock basis using new single-particle coupling amplitudes.

        Parameters
        ----------
        nbaths : int
            Number of the phonon baths.
        velph : dict
            The new single-particle electron-phonon couplings. See attribute velph.
        mtype : type
            Defines type of Vbbp0 and Vbbp matrices. For example, float, complex, etc.
        """
        si = self.si
        si.nbaths = nbaths
        self.velph = make_velph_dict(velph, si)
        self.mtype = mtype
        self.Vbbp0 = elph_construct_Vbbp(self, velph) def validate_itype_ph(itype_ph):
    if itype_ph not in {0, 2}:
        print("WARNING: itype_ph needs to be 0, or 2. Using default itype=0.")
        itype_ph = 0
    return itype_ph class ApproachPauli(ApproachElPh):

    kerntype = 'pyPauli'

    def get_kern_size(self):
        return self.si.npauli

    def prepare_arrays(self):
        ApproachPauliBase.prepare_arrays(self)
        nbaths, ndm0 = self.si_elph.nbaths, self.si_elph.ndm0
        self.paulifct_elph = np.zeros((nbaths, ndm0), dtype=doublenp)

    def clean_arrays(self):
        ApproachPauliBase.clean_arrays(self)
        self.paulifct_elph.fill(0.0)

    def generate_fct(self):
        ApproachPauliBase.generate_fct(self)

        E, Vbbp = self.qd.Ea, self.baths.Vbbp
        si, kh = self.si_elph, self.kernel_handler
        ncharge, nbaths, statesdm = si.ncharge, si.nbaths, si.statesdm

        func_pauli = FuncPauliElPh(self.baths.tlst_ph, self.baths.dlst_ph,
                                   self.baths.bath_func, self.funcp.eps_elph)

        paulifct = self.paulifct_elph
        for charge in range(ncharge):
            # The diagonal elements b=bp are excluded, because they do not contribute
            for b, bp in itertools.permutations(statesdm[charge], 2):
                bbp_bool = si.get_ind_dm0(b, bp, charge, maptype=2)
                if not bbp_bool:
                    continue
                bbp = si.get_ind_dm0(b, bp, charge)
                Ebbp = E[b]-E[bp]
                for l in range(nbaths):
                    xbbp = 0.5*(Vbbp[l, b, bp]*Vbbp[l, b, bp].conjugate() +
                                Vbbp[l, bp, b].conjugate()*Vbbp[l, bp, b]).real
                    func_pauli.eval(Ebbp, l)
                    paulifct[l, bbp] = xbbp*func_pauli.val

    def generate_kern(self):
        ApproachPauliBase.generate_kern(self)

    def generate_coupling_terms(self, b, bp, bcharge):
        ApproachPauliBase.generate_coupling_terms(self, b, bp, bcharge)

        paulifct = self.paulifct_elph
        si, si_elph, kh = self.si, self.si_elph, self.kernel_handler
        nbaths, statesdm = si.nbaths, si.statesdm

        acharge = bcharge

        bb = si.get_ind_dm0(b, b, bcharge)
        for a in statesdm[acharge]:
            aa = si.get_ind_dm0(a, a, acharge)
            ab = si_elph.get_ind_dm0(a, b, bcharge)
            ba = si_elph.get_ind_dm0(b, a, acharge)
            if aa == -1 or ba == -1:
                continue

            fctm, fctp = 0, 0
            for l in range(nbaths):
                fctm -= paulifct[l, ab]
                fctp += paulifct[l, ba]
            kh.set_matrix_element_pauli(fctm, fctp, bb, aa)

    def generate_current(self):
        ApproachPauliBase.generate_current(self) class ApproachLindblad(ApproachElPh):

    kerntype = 'pyLindblad'

    def prepare_arrays(self):
        ApproachLindbladBase.prepare_arrays(self)
        Vbbp, mtype = self.baths.Vbbp, self.baths.mtype
        tLbbp_shape = Vbbp.shape + (2,)
        self.tLbbp = np.zeros(tLbbp_shape, dtype=mtype)

    def clean_arrays(self):
        ApproachLindbladBase.clean_arrays(self)
        self.tLbbp.fill(0.0)

    def generate_fct(self):
        ApproachLindbladBase.generate_fct(self)

        Vbbp, E = self.baths.Vbbp, self.qd.Ea,
        si = self.si
        ncharge, nbaths, statesdm = si.ncharge, si.nbaths, si.statesdm

        func_pauli = FuncPauliElPh(self.baths.tlst_ph, self.baths.dlst_ph,
                                   self.baths.bath_func, self.funcp.eps_elph)

        tLbbp = self.tLbbp
        # Diagonal elements
        for l in range(nbaths):
            func_pauli.eval(0., l)
            for charge in range(ncharge):
                for b in si.statesdm[charge]:
                    tLbbp[l, b, b, 0] = np.sqrt(0.5*func_pauli.val)*Vbbp[l, b, b]
                    tLbbp[l, b, b, 1] = tLbbp[l, b, b, 0].conjugate()
        # Off-diagonal elements
        for charge in range(ncharge):
            for b, bp in itertools.permutations(statesdm[charge], 2):
                Ebbp = E[b]-E[bp]
                for l in range(nbaths):
                    func_pauli.eval(Ebbp, l)
                    tLbbp[l, b, bp, 0] = np.sqrt(0.5*func_pauli.val)*Vbbp[l, b, bp]
                    tLbbp[l, b, bp, 1] = np.sqrt(0.5*func_pauli.val)*Vbbp[l, bp, b].conjugate()
        self.tLbbp = tLbbp

    def generate_coupling_terms(self, b, bp, bcharge):
        ApproachLindbladBase.generate_coupling_terms(self, b, bp, bcharge)

        tLbbp = self.tLbbp
        si, kh = self.si, self.kernel_handler
        nbaths, statesdm = si.nbaths, si.statesdm

        acharge = bcharge

        # --------------------------------------------------
        # Here letter convention is not used
        # For example, the label `a' has the same charge as the label `b'
        for a, ap in itertools.product(statesdm[acharge], statesdm[acharge]):
            if kh.is_included(a, ap, acharge):
                fct_aap = 0
                for (l, q) in itertools.product(range(si.nbaths), range(2)):
                    fct_aap += tLbbp[l, b, a, q]*tLbbp[l, bp, ap, q].conjugate()
                kh.set_matrix_element(1j*fct_aap, b, bp, bcharge, a, ap, acharge)
        # --------------------------------------------------
        for bpp in statesdm[bcharge]:
            if kh.is_included(bpp, bp, bcharge):
                fct_bppbp = 0
                for a in statesdm[acharge]:
                    for (l, q) in itertools.product(range(si.nbaths), range(2)):
                        fct_bppbp += -0.5*tLbbp[l, a, b, q].conjugate()*tLbbp[l, a, bpp, q]
                kh.set_matrix_element(1j*fct_bppbp, b, bp, bcharge, bpp, bp, bcharge)
            # --------------------------------------------------
            if kh.is_included(b, bpp, bcharge):
                fct_bbpp = 0
                for a in statesdm[acharge]:
                    for (l, q) in itertools.product(range(si.nbaths), range(2)):
                        fct_bbpp += -0.5*tLbbp[l, a, bpp, q].conjugate()*tLbbp[l, a, bp, q]
                kh.set_matrix_element(1j*fct_bbpp, b, bp, bcharge, b, bpp, bcharge)
        # --------------------------------------------------

    def generate_current(self):
        ApproachLindbladBase.generate_current(self) class Approach1vN(ApproachElPh):

    kerntype = 'py1vN'

    def prepare_arrays(self):
        Approach1vNBase.prepare_arrays(self)
        nbaths, ndm0 = self.si_elph.nbaths, self.si_elph.ndm0
        self.w1fct = np.zeros((nbaths, ndm0, 2), dtype=complexnp)

    def clean_arrays(self):
        Approach1vNBase.clean_arrays(self)
        self.w1fct.fill(0.0)

    def generate_fct(self):
        Approach1vNBase.generate_fct(self)

        E, si = self.qd.Ea, self.si_elph
        ncharge, nbaths, statesdm = si.ncharge, si.nbaths, si.statesdm

        func_1vN_elph = Func1vNElPh(self.baths.tlst_ph, self.baths.dlst_ph,
                                    self.funcp.itype_ph, self.funcp.dqawc_limit,
                                    self.baths.bath_func,
                                    self.funcp.eps_elph)

        w1fct = self.w1fct
        # Diagonal elements
        for l in range(nbaths):
            func_1vN_elph.eval(0., l)
            for charge in range(ncharge):
                for b in statesdm[charge]:
                    bb = si.get_ind_dm0(b, b, charge)
                    bb_bool = si.get_ind_dm0(b, b, charge, maptype=2)
                    if bb != -1 and bb_bool:
                        w1fct[l, bb, 0] = func_1vN_elph.val0 - 0.5j*func_1vN_elph.val0.imag
                        w1fct[l, bb, 1] = func_1vN_elph.val1 - 0.5j*func_1vN_elph.val1.imag
        # Off-diagonal elements
        for charge in range(ncharge):
            for b, bp in itertools.permutations(statesdm[charge], 2):
                bbp = si.get_ind_dm0(b, bp, charge)
                bbp_bool = si.get_ind_dm0(b, bp, charge, maptype=2)
                if bbp != -1 and bbp_bool:
                    Ebbp = E[b]-E[bp]
                    for l in range(nbaths):
                        func_1vN_elph.eval(Ebbp, l)
                        w1fct[l, bbp, 0] = func_1vN_elph.val0
                        w1fct[l, bbp, 1] = func_1vN_elph.val1

    def generate_coupling_terms(self, b, bp, bcharge):
        Approach1vNBase.generate_coupling_terms(self, b, bp, bcharge)

        Vbbp, w1fct = self.baths.Vbbp, self.w1fct
        si, si_elph, kh = self.si, self.si_elph, self.kernel_handler
        nbaths, statesdm = si.nbaths, si.statesdm

        acharge = bcharge
        ccharge = bcharge

        # --------------------------------------------------
        for a, ap in itertools.product(statesdm[acharge], statesdm[acharge]):
            if kh.is_included(a, ap, acharge):
                bpa = si_elph.get_ind_dm0(bp, a, acharge)
                bap = si_elph.get_ind_dm0(b, ap, acharge)
                if bpa == -1 or bap == -1:
                    continue
                fct_aap = 0
                for l in range(nbaths):
                    gamma_ba_bpap = 0.5*(Vbbp[l, b, a]*Vbbp[l, bp, ap].conjugate() +
                                         Vbbp[l, a, b].conjugate()*Vbbp[l, ap, bp])
                    fct_aap += gamma_ba_bpap*(w1fct[l, bpa, 0].conjugate() - w1fct[l, bap, 0])
                kh.set_matrix_element(fct_aap, b, bp, bcharge, a, ap, acharge)
        # --------------------------------------------------
        for bpp in statesdm[bcharge]:
            if kh.is_included(bpp, bp, bcharge):
                fct_bppbp = 0
                for a in statesdm[acharge]:
                    bpa = si_elph.get_ind_dm0(bp, a, acharge)
                    if bpa == -1:
                        continue
                    for l in range(nbaths):
                        gamma_ba_bppa = 0.5*(Vbbp[l, b, a]*Vbbp[l, bpp, a].conjugate() +
                                             Vbbp[l, a, b].conjugate()*Vbbp[l, a, bpp])
                        fct_bppbp += gamma_ba_bppa*w1fct[l, bpa, 1].conjugate()
                for c in statesdm[ccharge]:
                    cbp = si_elph.get_ind_dm0(c, bp, bcharge)
                    if cbp == -1:
                        continue
                    for l in range(nbaths):
                        gamma_bc_bppc = 0.5*(Vbbp[l, b, c]*Vbbp[l, bpp, c].conjugate() +
                                             Vbbp[l, c, b].conjugate()*Vbbp[l, c, bpp])
                        fct_bppbp += gamma_bc_bppc*w1fct[l, cbp, 0]
                kh.set_matrix_element(fct_bppbp, b, bp, bcharge, bpp, bp, bcharge)
            # --------------------------------------------------
            if kh.is_included(b, bpp, bcharge):
                fct_bbpp = 0
                for a in statesdm[acharge]:
                    ba = si_elph.get_ind_dm0(b, a, acharge)
                    if ba == -1:
                        continue
                    for l in range(nbaths):
                        gamma_abpp_abp = 0.5*(Vbbp[l, a, bpp].conjugate()*Vbbp[l, a, bp] +
                                              Vbbp[l, bpp, a]*Vbbp[l, bp, a].conjugate())
                        fct_bbpp += -gamma_abpp_abp*w1fct[l, ba, 1]
                for c in statesdm[ccharge]:
                    cb = si_elph.get_ind_dm0(c, b, bcharge)
                    if cb == -1:
                        continue
                    for l in range(nbaths):
                        gamma_cbpp_cbp = 0.5*(Vbbp[l, c, bpp].conjugate()*Vbbp[l, c, bp] +
                                              Vbbp[l, bpp, c]*Vbbp[l, bp, c].conjugate())
                        fct_bbpp += -gamma_cbpp_cbp*w1fct[l, cb, 0].conjugate()
                kh.set_matrix_element(fct_bbpp, b, bp, bcharge, b, bpp, bcharge)
        # --------------------------------------------------
        for c, cp in itertools.product(statesdm[ccharge], statesdm[ccharge]):
            if kh.is_included(c, cp, ccharge):
                cbp = si_elph.get_ind_dm0(c, bp, bcharge)
                cpb = si_elph.get_ind_dm0(cp, b, bcharge)
                if cbp == -1 or cpb == -1:
                    continue
                fct_ccp = 0
                for l in range(nbaths):
                    gamma_bc_bpcp = 0.5*(Vbbp[l, b, c]*Vbbp[l, bp, cp].conjugate() +
                                         Vbbp[l, c, b].conjugate()*Vbbp[l, cp, bp])
                    fct_ccp += gamma_bc_bpcp*(w1fct[l, cbp, 1] - w1fct[l, cpb, 1].conjugate())
                kh.set_matrix_element(fct_ccp, b, bp, bcharge, c, cp, ccharge)
        # --------------------------------------------------

    def generate_current(self):
        Approach1vNBase.generate_current(self) class ApproachRedfield(ApproachElPh):

    kerntype = 'pyRedfield'

    def prepare_arrays(self):
        Approach1vN.prepare_arrays(self)

    def clean_arrays(self):
        Approach1vN.clean_arrays(self)

    def generate_fct(self):
        Approach1vN.generate_fct(self)

    def generate_coupling_terms(self, b, bp, bcharge):
        ApproachRedfieldBase.generate_coupling_terms(self, b, bp, bcharge)

        Vbbp, w1fct = self.baths.Vbbp, self.w1fct
        si, si_elph, kh = self.si, self.si_elph, self.kernel_handler
        nbaths, statesdm = si.nbaths, si.statesdm

        acharge = bcharge
        ccharge = bcharge

        # --------------------------------------------------
        for a, ap in itertools.product(statesdm[acharge], statesdm[acharge]):
            if kh.is_included(a, ap, acharge):
                bpap = si_elph.get_ind_dm0(bp, ap, acharge)
                ba = si_elph.get_ind_dm0(b, a, acharge)
                if bpap == -1 or ba == -1:
                    continue
                fct_aap = 0
                for l in range(nbaths):
                    gamma_ba_bpap = 0.5*(Vbbp[l, b, a]*Vbbp[l, bp, ap].conjugate() +
                                         Vbbp[l, a, b].conjugate()*Vbbp[l, ap, bp])
                    fct_aap += gamma_ba_bpap*(w1fct[l, bpap, 0].conjugate() - w1fct[l, ba, 0])
                kh.set_matrix_element(fct_aap, b, bp, bcharge, a, ap, acharge)
        # --------------------------------------------------
        for bpp in statesdm[bcharge]:
            if kh.is_included(bpp, bp, bcharge):
                fct_bppbp = 0
                for a in statesdm[acharge]:
                    bppa = si_elph.get_ind_dm0(bpp, a, acharge)
                    if bppa == -1:
                        continue
                    for l in range(nbaths):
                        gamma_ba_bppa = 0.5*(Vbbp[l, b, a]*Vbbp[l, bpp, a].conjugate() +
                                             Vbbp[l, a, b].conjugate()*Vbbp[l, a, bpp])
                        fct_bppbp += +gamma_ba_bppa*w1fct[l, bppa, 1].conjugate()
                for c in statesdm[ccharge]:
                    cbpp = si_elph.get_ind_dm0(c, bpp, bcharge)
                    if cbpp == -1:
                        continue
                    for l in range(nbaths):
                        gamma_bc_bppc = 0.5*(Vbbp[l, b, c]*Vbbp[l, bpp, c].conjugate() +
                                             Vbbp[l, c, b].conjugate()*Vbbp[l, c, bpp])
                        fct_bppbp += +gamma_bc_bppc*w1fct[l, cbpp, 0]
                kh.set_matrix_element(fct_bppbp, b, bp, bcharge, bpp, bp, bcharge)
            # --------------------------------------------------
            if kh.is_included(b, bpp, bcharge):
                fct_bbpp = 0
                for a in statesdm[acharge]:
                    bppa = si_elph.get_ind_dm0(bpp, a, acharge)
                    if bppa == -1:
                        continue
                    for l in range(nbaths):
                        gamma_abpp_abp = 0.5*(Vbbp[l, a, bpp].conjugate()*Vbbp[l, a, bp] +
                                              Vbbp[l, bpp, a]*Vbbp[l, bp, a].conjugate())
                        fct_bbpp += -gamma_abpp_abp*w1fct[l, bppa, 1]
                for c in statesdm[ccharge]:
                    cbpp = si_elph.get_ind_dm0(c, bpp, bcharge)
                    if cbpp == -1:
                        continue
                    for l in range(nbaths):
                        gamma_cbpp_cbp = 0.5*(Vbbp[l, c, bpp].conjugate()*Vbbp[l, c, bp] +
                                              Vbbp[l, bpp, c]*Vbbp[l, bp, c].conjugate())
                        fct_bbpp += -gamma_cbpp_cbp*w1fct[l, cbpp, 0].conjugate()
                kh.set_matrix_element(fct_bbpp, b, bp, bcharge, b, bpp, bcharge)
        # --------------------------------------------------
        for c, cp in itertools.product(statesdm[ccharge], statesdm[ccharge]):
            if kh.is_included(c, cp, ccharge):
                cpbp = si_elph.get_ind_dm0(cp, bp, bcharge)
                cb = si_elph.get_ind_dm0(c, b, bcharge)
                if cpbp == -1 or cb == -1:
                    continue
                fct_ccp = 0
                for l in range(nbaths):
                    gamma_bc_bpcp = 0.5*(Vbbp[l, b, c]*Vbbp[l, bp, cp].conjugate() +
                                         Vbbp[l, c, b].conjugate()*Vbbp[l, cp, bp])
                    fct_ccp += gamma_bc_bpcp*(w1fct[l, cpbp, 1] - w1fct[l, cb, 1].conjugate())
                kh.set_matrix_element(fct_ccp, b, bp, bcharge, c, cp, ccharge)
        # --------------------------------------------------

    def generate_current(self):
        ApproachRedfieldBase.generate_current(self) class Approach2vN(ApproachBase2vN):

    kerntype = 'py2vN'

    #region Preparation

    def prepare_kern(self):
        if self.is_zeroth_iteration:
            self.prepare_energy_grid()
            self.prepare_distribution()

        ApproachBase2vN.prepare_kern(self)

    def prepare_energy_grid(self):
        """"Prepares the energy grid."""
        self.make_Ek_grid()

        # Define the extended grid Ek_grid_ext for calculations outside the bandwidth
        # Here funcp.emin, funcp.emax are defined
        self.determine_emin_emax()

        # Here Ek_grid_ext, funcp.kpnt_left, funcp.kpnt_right are defined
        self.determine_grid_ext()

    def make_Ek_grid(self):
        """Make an energy grid on which 2vN equations are solved. """
        if self.funcp.kpnt is None:
            raise ValueError('kpnt needs to be specified.')
        if self.si.nleads > 0:
            dmin = np.min(self.leads.dlst)
            dmax = np.max(self.leads.dlst)
            Ek_grid, kpnt = self.Ek_grid, self.funcp.kpnt
            if Ek_grid[0] != dmin or Ek_grid[-1] != dmax or Ek_grid.shape[0] != kpnt:
                self.funcp.dmin = dmin
                self.funcp.dmax = dmax
                self.Ek_grid = np.linspace(dmin, dmax, kpnt)
                #
                if self.niter != -1:
                    print("WARNING: Ek_grid has changed. Restarting the calculation.")
                    self.restart()
                #
                if ((dmin * np.ones(self.si.nleads)).tolist() != self.leads.dlst.T[0].tolist() or
                    (dmax * np.ones(self.si.nleads)).tolist() != self.leads.dlst.T[1].tolist()):
                    print("WARNING: The bandwidth and Ek_grid for all leads will be the same: from " +
                          "dmin=" + str(dmin) + " to dmax=" + str(dmax) + ".")

    def determine_emin_emax(self):
        """
        Finds how much the Ek_grid has to be expanded
        above the bandwidth D of the leads.

        Parameters
        ----------
        funcp.emin : float
            (Modifies) Minimal energy in the updated Ek_grid.
        funcp.emax : float
            (Modifies) Maximal energy in the updated Ek_grid.
        """
        E, si = self.qd.Ea, self.si
        dmin, dmax = self.funcp.dmin, self.funcp.dmax

        lst = [dmin, dmax]
        for charge in range(si.ncharge):
            for b, bp in itertools.product(si.statesdm[charge], si.statesdm[charge]):
                lst.append(dmax-E[b]+E[bp])
                lst.append(dmin-E[b]+E[bp])
        for charge in range(si.ncharge-2):
            for d, b in itertools.product(si.statesdm[charge+2], si.statesdm[charge]):
                lst.append(dmin+E[d]-E[b])
                lst.append(dmax+E[d]-E[b])

        self.funcp.emax = max(lst)
        self.funcp.emin = min(lst)

    def determine_grid_ext(self):
        """
        Expands the Ek_grid from emin to emax.

        Parameters
        ----------
        Ek_grid_ext : ndarray
            (Modifies) Extended Ek_grid from emin to emax.
        funcp.kpnt_left : int
            (Modifies) Number of points Ek_grid is extended to the left.
        funcp.kpnt_right : int
            (Modifies) Number of points Ek_grid is extended to the right.
        """
        emin_, emax_ = self.funcp.emin, self.funcp.emax
        dmin, dmax = self.funcp.dmin, self.funcp.dmax
        ext_fct = self.funcp.ext_fct
        Ek_grid = self.Ek_grid

        step = Ek_grid[1]-Ek_grid[0]
        emin = ext_fct*(emin_-dmin)+dmin
        emax = ext_fct*(emax_-dmax)+dmax
        ext_left = np.sort(-np.arange(-dmin+step, -emin+step, step))
        ext_right = np.arange(dmax+step, emax+step, step)

        self.Ek_grid_ext = np.concatenate((ext_left, Ek_grid, ext_right))
        self.funcp.kpnt_left, self.funcp.kpnt_right = len(ext_left), len(ext_right)

    def prepare_distribution(self):
        """"
        Prepares the Fermi distribution functions
        and their Hilbert transform on the energy grid.

        Parameters
        ----------
        fkp : ndarray
            (Modifies) nleads by len(Ek_grid_ext) numpy array containing Fermi function.
        fkm : ndarray
            (Modifies) nleads by len(Ek_grid_ext) numpy array containing 1-Fermi function.
        hfkp : ndarray
            (Modifies) Hilbert transform of fkp.
        hfkm : ndarray
            (Modifies) Hilbert transform of fkm.
        """
        # Generate the Fermi functions on the grid
        # This is necessary to generate only if Ek_grid, mulst, or tlst are changed
        nleads = self.si.nleads
        mulst, tlst = self.leads.mulst, self.leads.tlst
        Ek_grid = self.Ek_grid
        Eklen = len(Ek_grid)

        self.fkp = np.zeros((nleads, Eklen), dtype=doublenp)
        for l in range(nleads):
            self.fkp[l] = 1/(np.exp((Ek_grid - mulst[l])/tlst[l]) + 1)
        self.fkm = 1-self.fkp
        self.fkp, self.hfkp = get_htransf_fk(self.fkp, self.funcp)
        self.fkm, self.hfkm = get_htransf_fk(self.fkm, self.funcp)

    #endregion Preparation

    #region Generation

    def iterate(self):
        terms_calculator = TermsCalculator2vN(self)
        terms_calculator.iterate()

    def determine_phi1_phi0(self):
        """
        Integrates phi1k over energy.

        Parameters
        ----------
        phi1_phi0 : ndarray
            (Modifies) Numpy array with dimensions (nleads, ndm1, ndm0),
            corresponding to energy integrated Phi[1](k).
        e_phi1_phi0 : ndarray
            (Modifies) Numpy array with dimensions (nleads, ndm1, ndm0),
            corresponding to energy integrated Ek*Phi[1](k).
        """
        phi1k, Ek_grid = self.phi1k, self.Ek_grid
        # Get integrated Phi[1]_{cb} in terms of Phi[0]_{bb'}
        phi1_phi0 = self.phi1_phi0
        e_phi1_phi0 = self.e_phi1_phi0

        Eklen = len(Ek_grid)
        for i in range(Eklen):
            # Trapezoidal rule
            if i == 0:
                dx = Ek_grid[i+1] - Ek_grid[i]
            elif i == Eklen-1:
                dx = Ek_grid[i] - Ek_grid[i-1]
            else:
                dx = Ek_grid[i+1] - Ek_grid[i-1]
            phi1_phi0 += 1/2*dx*phi1k[i]
            e_phi1_phi0 += 1/2*dx*Ek_grid[i]*phi1k[i]

    def generate_kern(self):
        """
        From Phi[1](k) generate equations containing just Phi[0].

        Parameters
        ----------
        kern : ndarray
            (Modifies) Kernel for master equation involving just Phi[0].
        bvec : ndarray
            (Modifies) Right hand side column vector for master equation.
        """
        phi1_phi0, E, Tba, si = self.phi1_phi0, self.qd.Ea, self.leads.Tba, self.si
        ncharge, nleads, ndm0, statesdm = si.ncharge, si.nleads, si.ndm0, si.statesdm

        # Integrated Phi[1]_{bc} in terms of phi1_phi0
        shuffle = np.zeros((ndm0, ndm0), dtype=intnp)
        for charge in range(ncharge):
            for b, bp in itertools.combinations_with_replacement(statesdm[charge], 2):
                bbp = si.get_ind_dm0(b, bp, charge)
                bpb = si.get_ind_dm0(bp, b, charge)
                shuffle[bbp, bpb] = 1
                shuffle[bpb, bbp] = 1
        phi1_phi0_conj = np.dot(np.conjugate(phi1_phi0), shuffle)

        kern = self.kern
        for charge in range(ncharge):
            acharge = charge-1
            bcharge = charge
            ccharge = charge+1
            for b, bp in itertools.product(statesdm[bcharge], statesdm[bcharge]):
                bbp = si.get_ind_dm0(b, bp, bcharge)
                kern[bbp, bbp] = E[b]-E[bp]

                for a1 in statesdm[acharge]:
                    bpa1 = si.get_ind_dm1(bp, a1, acharge)
                    ba1 = si.get_ind_dm1(b, a1, acharge)
                    for l in range(nleads):
                        kern[bbp] += +Tba[l, b, a1]*phi1_phi0_conj[l, bpa1]
                        kern[bbp] += -phi1_phi0[l, ba1]*Tba[l, a1, bp]

                for c1 in statesdm[ccharge]:
                    c1bp = si.get_ind_dm1(c1, bp, bcharge)
                    c1b = si.get_ind_dm1(c1, b, bcharge)
                    for l in range(nleads):
                        kern[bbp] += +Tba[l, b, c1]*phi1_phi0[l, c1bp]
                        kern[bbp] += -phi1_phi0_conj[l, c1b]*Tba[l, c1, bp]

    def generate_current(self):
        """
        Finds the currents for 2vN approach after an approximate
        solution for integral 2vN equations was found.

        Parameters
        ----------
        phi1 : ndarray
            (Modifies) nleads by ndm1 numpy array, which gives energy integrated
            Phi[1](k) after the values for Phi[0] were inserted.
        current : ndarray
            (Modifies) Values of the current having nleads entries.
        energy_current : ndarray
            (Modifies) Values of the energy current having nleads entries.
        heat_current : ndarray
            (Modifies) Values of the heat current having nleads entries.
        """
        phi1_phi0, e_phi1_phi0 = self.phi1_phi0, self.e_phi1_phi0
        phi0, Tba, si = self.phi0, self.leads.Tba, self.si
        ncharge, nleads, ndm1, statesdm = si.ncharge, si.nleads, si.ndm1, si.statesdm

        phi1 = self.phi1
        h_phi1 = self.h_phi1
        current = self.current
        energy_current = self.energy_current

        for l in range(nleads):
            phi1[l] = np.dot(phi1_phi0[l], phi0)
            h_phi1[l] = np.dot(e_phi1_phi0[l], phi0)

        for charge in range(ncharge-1):
            ccharge = charge+1
            bcharge = charge
            for c, b in itertools.product(statesdm[ccharge], statesdm[bcharge]):
                cb = si.get_ind_dm1(c, b, bcharge)

                for l in range(nleads):
                    current_l = phi1[l, cb]*Tba[l, b, c]
                    energy_current_l = h_phi1[l, cb]*Tba[l, b, c]

                    current[l] += -2*current_l.imag
                    energy_current[l] += -2*energy_current_l.imag

        self.heat_current[:] = energy_current - current*self.leads.mulst

    #endregion Generation
from .builder_base import BuilderBase
from .builder_base import attribute_map
from .builder_base import BuilderManyBody
from ..indexing import StateIndexingDM
from ..indexing import StateIndexingDMc
from ..baths import PhononBaths
from .validation import validate_itype_ph
from ..approach.elph.pauli import ApproachPauli as ApproachPyPauli
from ..approach.elph.lindblad import ApproachLindblad as ApproachPyLindblad
from ..approach.elph.neumann1 import Approach1vN as ApproachPy1vN
from ..approach.elph.redfield import ApproachRedfield as ApproachPyRedfield
from ..approach.base.neumann2 import Approach2vN as ApproachPy2vN
    from ..approach.elph.c_pauli import ApproachPauli
    from ..approach.elph.c_lindblad import ApproachLindblad
    from ..approach.elph.c_redfield import ApproachRedfield
    from ..approach.elph.c_neumann1 import Approach1vN
    from ..approach.base.c_neumann2 import Approach2vN
"""Module containing Builder_elph and Builder_many_body_elph classes."""



# -----------------------------------------------------------
# Python modules


# Cython compiled modules

try:
except ImportError:
    print("WARNING: Cannot import Cython compiled modules for the approaches (builder_elph.py).")
    ApproachPauli = ApproachPyPauli

    ApproachLindblad = ApproachPyLindblad