Reformat to standard Python coding practices

JJArray/junction_array.py

@@ -8,21 +8,20 @@
 
 
 class junction_array(scq.Circuit):
-    def __init__(
-        self, EJs: list, ECs: list, converge: bool = False, Ncut: int = 5
-    ) -> None:
+    def __init__(self, ejs: list, ecs: list, converge: bool = False, ncut: int = 5,
+                 monitor=True) -> None:
         """
 
 
         Parameters
         ----------
-        EJs :
+        ejs :
             list of Josephson energies
-        ECs :
+        ecs :
             list of charging energies
         converge :
-            boolean option to converge Ncut for a sample of the flux and ngs parameter space
-        Ncut :
+            boolean option to converge ncut for a sample of the flux and ngs parameter space
+        ncut :
             cutoff value, the same for all junctions
 
         """
@@ -31,34 +30,39 @@ def __init__(
         self.ECs = ECs
         self.N = len(EJs) - 1
 
-        self.Ncut = None
+        self.ncut = None
         self.symmetric_basis_transformation = None
 
-        self.circuit_setup(Ncut, converge)
+        self.circuit_setup(ncut, converge, monitor)
 
     def circuit_setup(self, Ncut: int, converge: bool) -> None:
 
-        JJcirc_yaml = "branches:\n"
-
-        N_junctions = self.N + 1
-        for n in range(N_junctions):
-            JJcirc_yaml = (
-                JJcirc_yaml
-                + '- ["JJ", '
-                + str(n)
-                + ", "
-                + str((n + 1) % N_junctions)
-                + ", EJ"
-                + str(n % N_junctions)
-                + str((n + 1) % N_junctions)
-                + " = "
-                + str(self.EJs[n % N_junctions])
-                + ", EC"
-                + str(n % N_junctions)
-                + str((n + 1) % N_junctions)
-                + " = "
-                + str(self.ECs[n % N_junctions])
-                + "]\n"
+        :rtype: object
+        :param ncut: 
+        :param converge: 
+        :param monitor: 
+        """
+        jj_circ_yaml = "branches:\n"
+
+        n_junctions = self.N + 1
+        for n in range(n_junctions):
+            jj_circ_yaml = (
+                    jj_circ_yaml
+                    + '- ["JJ", '
+                    + str(n)
+                    + ", "
+                    + str((n + 1) % n_junctions)
+                    + ", EJ"
+                    + str(n % n_junctions)
+                    + str((n + 1) % n_junctions)
+                    + " = "
+                    + str(self.EJs[n % n_junctions])
+                    + ", EC"
+                    + str(n % n_junctions)
+                    + str((n + 1) % n_junctions)
+                    + " = "
+                    + str(self.ECs[n % n_junctions])
+                    + "]\n"
             )
 
         super().__init__(JJcirc_yaml, from_file=False)
@@ -71,34 +75,40 @@ def circuit_setup(self, Ncut: int, converge: bool) -> None:
         )
 
         if converge:
-            self.Ncut = self.converge_cutoff()
+            self.ncut = self.converge_cutoff(monitor)
         else:
-            self.Ncut = Ncut
-            self.set_cutoff(Ncut)
+            self.ncut = ncut
+            self.set_cutoff(ncut)
 
-    def set_cutoff(self, Ncut: int) -> None:
+    def set_cutoff(self, ncut: int) -> None:
         """
         Parameters
         ----------
-        Ncut :
+        ncut :
             cutoff value, the same for all junctions
         """
         for nc in self.cutoff_names:
-            self.__dict__["_" + nc] = Ncut
+            self.__dict__["_" + nc] = ncut
 
-    def converge_cutoff(self):
+    def converge_cutoff(self, monitor: object = True) -> object:
+        """
 
-        Ncut = 1
+        :rtype: object
+        :return: 
+        :param monitor: 
+        :return: 
+        """
+        ncut = 1
         for nc in self.cutoff_names:
-            self.__dict__["_" + nc] = Ncut
+            self.__dict__["_" + nc] = ncut
 
-        nvals = np.min([5, (2 * Ncut + 1) ** self.N])
-        ntries = 10
+        nvals = np.min([5, (2 * ncut + 1) ** self.N])
+        ntries = 2
         grid = np.linspace(0, 1, 5)
-        grid = grid[0 : len(grid) - 1]
+        grid = grid[0: len(grid) - 1]
         sample = np.concatenate(
             (
-                self.cartesian(tuple([grid for r in range(self.N + 1)])),
+                self.cartesian(tuple([grid for _ in range(self.N + 1)])),
                 np.random.rand(ntries, self.N + 1),
             )
         )
@@ -134,9 +144,7 @@ def converge_cutoff(self):
                     Ncut = Ncut + 1
                     print("Changed at sample:", (Ncut, sample[smp]))
 
-        print("Final Ncut value:", Ncut)
-
-        return Ncut
+        return ncut
 
     def cartesian(self, arrays, out=None):
         arrays = [np.asarray(x) for x in arrays]
@@ -151,77 +159,77 @@ def cartesian(self, arrays, out=None):
         if arrays[1:]:
             self.cartesian(arrays[1:], out=out[0:m, 1:])
             for j in range(1, arrays[0].size):
-                out[j * m : (j + 1) * m, 1:] = out[0:m, 1:]
+                out[j * m: (j + 1) * m, 1:] = out[0:m, 1:]
         return out
 
-    def GenerateMomentumBasis(self):
-        dms = (2 * self.Ncut + 1) ** self.N
+    def generate_momentum_basis(self):
+        """
+
+        :return:
+        :rtype: object
+        """
+        dms = (2 * self.ncut + 1) ** self.N
         lst = list(range(dms))
         vbs = []
 
         while len(lst) > 0:
 
             bas = [lst[0]]
-            tmp = self.TransInd(bas[-1], self.N, self.Ncut)
+            tmp = self.trans_ind(bas[-1], self.N, self.ncut)
 
             while (tmp - bas[0]) != 0:
                 bas.append(tmp)
-                tmp = self.TransInd(bas[-1], self.N, self.Ncut)
+                tmp = self.trans_ind(bas[-1], self.N, self.ncut)
 
             vbs.append(bas)
 
             lst = list(set(lst) - set(bas))
 
         return vbs
 
-    def amps(self, cycle, k):
-        N0 = len(cycle)
-        seq = np.arange(N0)
-        return (1 / np.sqrt(N0)) * np.exp(-1j * 2 * np.pi * (k) * seq)
-
-    def ChargeToTranslation(self):
+    def charge_to_translation(self):
 
         # Check if vbs file has been generated, if not then generate it
         flnms = os.listdir()
         flag = 0
         for fname in flnms:
-            if "vbs_" + str((self.N, self.Ncut)) + ".npz" in fname:
+            if "vbs_" + str((self.N, self.ncut)) + ".npz" in fname:
                 flag = 1
         if flag == 0:
-            vbs = np.array(self.GenerateMomentumBasis(), dtype=object)
-            np.savez_compressed("vbs_" + str((self.N, self.Ncut)), vbs=vbs)
+            vbs = np.array(self.generate_momentum_basis(), dtype=object)
+            np.savez_compressed("vbs_" + str((self.N, self.ncut)), vbs=vbs)
         else:
             loaded = np.load(
-                "vbs_" + str((self.N, self.Ncut)) + ".npz", allow_pickle=True
+                "vbs_" + str((self.N, self.ncut)) + ".npz", allow_pickle=True
             )
             vbs = loaded["vbs"].tolist()
 
-        DIMS = [
-            [2 * self.Ncut + 1 for r in range(self.N)],
-            [2 * self.Ncut + 1 for r in range(self.N)],
+        dims = [
+            [2 * self.ncut + 1 for _ in range(self.N)],
+            [2 * self.ncut + 1 for _ in range(self.N)],
         ]
 
         mms = np.array([n / self.N for n in range(self.N)])
         lns = np.unique(list(map(len, vbs)))
 
-        m_seqs = [[] for r in range(self.N)]
+        m_seqs = [[] for _ in range(self.N)]
         for N0 in lns:
             tmp = np.array([n / N0 for n in range(N0)])
             for k in range(self.N):
                 if mms[k] in tmp:
                     m_seqs[k].append(N0)
 
-        V = []
+        v_op = []
         sm = 0
-        print("Full Hilbert space dimension = ", (2 * self.Ncut + 1) ** self.N)
+        print("Full Hilbert space dimension = ", (2 * self.ncut + 1) ** self.N)
         for n in tqdm(range(self.N)):
             row = []
             col = []
             data = []
             ind = 0
             for state in tqdm(range(len(vbs))):
                 if len(vbs[state]) in m_seqs[n]:
-                    vec = self.amps(vbs[state], n / self.N)
+                    vec = _amps(vbs[state], n / self.N)
                     row.append(vbs[state])
                     col.append((ind * np.ones(len(vec), dtype=int)).tolist())
                     data.append(vec.tolist())
@@ -231,35 +239,35 @@ def ChargeToTranslation(self):
             data = [item for sublist in data for item in sublist]
             print("Sector " + str(n) + "=", ind)
             sm = sm + ind
-            V.append(
+            v_op.append(
                 qt.Qobj(
                     csr_matrix(
-                        (data, (row, col)), shape=((2 * self.Ncut + 1) ** self.N, ind)
+                        (data, (row, col)), shape=((2 * self.ncut + 1) ** self.N, ind)
                     ),
-                    dims=DIMS,
+                    dims=dims,
                 )
             )
             # print(V[-1].shape)
         print("Sum of sector dimensions = ", sm)
-        if sm != (2 * self.Ncut + 1) ** self.N:
-            np.savetxt("WARNING_" + str([self.N, self.Ncut]) + ".txt", [1])
-        self.symmetric_basis_transformation = V
-
-    def SectorDiagonalization(self, Nvals):
-        self.ChargeToTranslation()
-        DIMS = ((2 * self.Ncut + 1) * np.ones((2, self.N), dtype=int)).tolist()
-        H = qt.Qobj(self.hamiltonian(), dims=DIMS)
+        if sm != (2 * self.ncut + 1) ** self.N:
+            np.savetxt("WARNING_" + str([self.N, self.ncut]) + ".txt", [1])
+        self.symmetric_basis_transformation = v_op
+
+    def sector_diagonalization(self, nvals):
+        self.charge_to_translation()
+        dims = ((2 * self.ncut + 1) * np.ones((2, self.N), dtype=int)).tolist()
+        h = qt.Qobj(self.hamiltonian(), dims=dims)
         # print(H.shape)
 
         data = []
         for k in range(len(self.symmetric_basis_transformation)):
             # print(self.symmetric_basis_transformation[k].dims)
-            HF0 = (
-                self.symmetric_basis_transformation[k].dag()
-                * H
-                * self.symmetric_basis_transformation[k]
+            hf0 = (
+                    self.symmetric_basis_transformation[k].dag()
+                    * h
+                    * self.symmetric_basis_transformation[k]
             )
-            evals, evecs = HF0.eigenstates(eigvals=Nvals, sparse=True)
+            evals, evecs = hf0.eigenstates(eigvals=nvals, sparse=True)
             data.append(
                 [
                     evals,
@@ -272,9 +280,9 @@ def SectorDiagonalization(self, Nvals):
 
         return data
 
-    def SortedDiagonalization(self, Nvals):
+    def sorted_diagonalization(self, nvals):
 
-        tst = self.SectorDiagonalization(Nvals)
+        tst = self.sector_diagonalization(nvals)
 
         vecs = []
         eigs = []
@@ -291,16 +299,13 @@ def SortedDiagonalization(self, Nvals):
 
         return eigs, vecs, tst
 
-    def TransInd(self, s, N, Ncut):
+    def trans_ind(self, s, number_junctions, ncut):
         return int(
-            self.ind2occ(s, N - 1, N, Ncut) * ((2 * Ncut + 1) ** (N - 1))
-            - (self.ind2occ(s, N - 1, N, Ncut) / (2 * Ncut + 1))
-            + s / (2 * Ncut + 1)
+            _ind2occ(s, number_junctions - 1, number_junctions, ncut) * ((2 * ncut + 1) ** (number_junctions - 1))
+            - (_ind2occ(s, number_junctions - 1, number_junctions, ncut) / (2 * ncut + 1))
+            + s / (2 * ncut + 1)
         )
 
-    def ind2occ(self, s, r, N, Ncut):
-        return int((s // ((2 * Ncut + 1) ** (N - r - 1))) % (2 * Ncut + 1))
-
     def hamiltonian_one_mode_model(self, ph_points=200):
 
         ph_list = np.linspace(-self.N * np.pi, self.N * np.pi, ph_points)
@@ -310,16 +315,16 @@ def hamiltonian_one_mode_model(self, ph_points=200):
         d2d2phi[ph_points - 1, 0] = 1
         d2d2phi = qt.Qobj(d2d2phi) / (dphi * dphi)
 
-        EC = self.ECs[0] / (1 + self.ECs[0] / (np.mean(self.ECs[1:]) * self.N))
+        ec = self.ECs[0] / (1 + self.ECs[0] / (np.mean(self.ECs[1:]) * self.N))
         h = qt.Qobj(
-            -4 * EC * d2d2phi
-            - (self.N) * np.mean(self.EJs[1:]) * np.diag(np.cos(ph_list / (self.N)))
+            -4 * ec * d2d2phi
+            - self.N * np.mean(self.EJs[1:]) * np.diag(np.cos(ph_list / self.N))
             - self.EJs[0] * np.diag(np.cos(ph_list + 2 * np.pi * self.Φ1))
         )
 
         return h
 
-    def paramSweep(self, param):
+    def param_sweep(self, param):
         ng_list = np.linspace(-2, 2, 220)
         self.plot_evals_vs_paramvals(
             param, ng_list, evals_count=7, num_cpus=1, subtract_ground=True