[Bug report] sym_external_fluxes() doesn't work for the system of coupled fluxonia

[AlesyaSokol]

Documentation check

If the bug manifests in unexpected behavior (as opposed to a crash), confirm that you have consulted the
API documentation

• I have checked the API documentation.
• I could not locate relevant information in the documentation or information is missing.

Describe the bug

I make a system of capacitively coupled fluxoniums and try to see how it assigned external fluxes:

EC1 = 1.*2*np.pi
EJ1 = 4*2*np.pi
EL1 = 0.55*2*np.pi
EC2 = 1.5*2*np.pi
EJ2 = 5*2*np.pi
EL2 = 0.5*2*np.pi
ECg = 10000*2*np.pi

coupled_fl_yaml = """
branches:
- ["JJ", 1, 2, """+str(EJ1)+""", """+str(EC1)+"""]
- ["L", 1, 2, """+str(EL1)+"""]
- ["C", 2, 3, """+str(ECg)+"""]
- ["JJ", 3, 4, """ +str(EJ2)+ """, """+str(EC2)+"""]
- ["L", 3, 4, """+str(EL2)+"""]
"""

coupled_fl = scq.Circuit.from_yaml(coupled_fl_yaml, from_file=False, ext_basis="harmonic")

coupled_fl.Φ1 = 0.495
coupled_fl.Φ2 = 0.495

coupled_fl.sym_external_fluxes()

However, it gives an error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Input In [12], in <cell line: 16>()
     13 coupled_fl.Φ2 = 0.495
     15 # print(coupled_fl.sym_hamiltonian())
---> 16 coupled_fl.sym_external_fluxes()

File ~/miniconda3/envs/quantum/lib/python3.9/site-packages/scqubits/core/circuit.py:2900, in Circuit.sym_external_fluxes(self)
   2890 def sym_external_fluxes(self) -> Dict[sm.Expr, Tuple["Branch", List["Branch"]]]:
   2891     """
   2892     Method returns a dictionary of Human readable external fluxes with associated
   2893     branches and loops (represented as lists of branches) for the current instance
   (...)
   2898         branches and loops
   2899     """
-> 2900     return {
   2901         self._make_expr_human_readable(self.external_fluxes[ibranch]): (
   2902             self.closure_branches[ibranch],
   2903             self.symbolic_circuit._find_loop(self.closure_branches[ibranch]),
   2904         )
   2905         for ibranch in range(len(self.external_fluxes))
   2906     }

File ~/miniconda3/envs/quantum/lib/python3.9/site-packages/scqubits/core/circuit.py:2903, in <dictcomp>(.0)
   2890 def sym_external_fluxes(self) -> Dict[sm.Expr, Tuple["Branch", List["Branch"]]]:
   2891     """
   2892     Method returns a dictionary of Human readable external fluxes with associated
   2893     branches and loops (represented as lists of branches) for the current instance
   (...)
   2898         branches and loops
   2899     """
   2900     return {
   2901         self._make_expr_human_readable(self.external_fluxes[ibranch]): (
   2902             self.closure_branches[ibranch],
-> 2903             self.symbolic_circuit._find_loop(self.closure_branches[ibranch]),
   2904         )
   2905         for ibranch in range(len(self.external_fluxes))
   2906     }

File ~/miniconda3/envs/quantum/lib/python3.9/site-packages/scqubits/core/symbolic_circuit.py:1470, in SymbolicCircuit._find_loop(self, closure_branch)
   1468 gen_last_same_ancestor = 0
   1469 for igen in range(min(gen_1, gen_2)):
-> 1470     if ancestors_1[igen].id == ancestors_2[igen].id:
   1471         gen_last_same_ancestor = igen
   1472     elif ancestors_1[igen].id != ancestors_2[igen].id:

IndexError: list index out of range

Expected behavior

No errors

To Reproduce

See above

OS and version used (please complete the following information):

• OS: Linux
• scqubits version 3.0.0
• Python version 3.9.12

Additional context

If I try to output symbolic Hamiltonian, it randomly gives different results from just repeated execution of the same code. In particular, the command

print(coupled_fl.sym_hamiltonian())

randomly gives either

15.708*Q1**2 - 6.283*Q1*Q2 + 15.708*Q2**2 + 1.571*Φ1**2 - 3.142*Φ1*θ1 + 3.142*Φ1*θ2 + 1.728*Φ2**2 - 3.456*Φ2*θ1 - 3.456*Φ2*θ2 + 3.299*θ1**2 + 0.314*θ1*θ2 + 3.299*θ2**2 - 31.416*cos(θ1 - 1.0*θ2) - 25.133*cos(θ1 + θ2)

or

15.708*Q1**2 - 6.283*Q1*Q2 + 15.708*Q2**2 + 1.728*Φ1**2 - 3.456*Φ1*θ1 - 3.456*Φ1*θ2 + 1.571*Φ2**2 - 3.142*Φ2*θ1 + 3.142*Φ2*θ2 + 3.299*θ1**2 + 0.314*θ1*θ2 + 3.299*θ2**2 - 31.416*cos(θ1 - 1.0*θ2) - 25.133*cos(θ1 + θ2)

One can see that it rearranges external fluxes. I assume that by itself it should not be a problem if it would be possible to know in advance which flux is assigned to which loop. That's why I tried to use sym_external_fluxes(), which doesn't work

[ZhaoTianPu]

Many thanks for reporting this bug, and sorry for replying late. We confirm we can reproduce it, and the fix will be released with the updated code release soon.

[saipavanc]

Thanks again for reporting the bug. About your additional comments, it seemed like the external fluxes were getting swapped with every run. I was able to track the issue to how we were finding the closure branches, where I used set objects. But, sets in Python do not generally preserve the order in an iterable, which makes the ordering random at times. I have replaced the code with an equivalent list comprehension which does preserve ordering.

Also, just a suggestion: as Sympy uses its own printing internally, its preferable to call coupled_fl.sym_hamiltonian() instead of print(coupled_fl.sym_hamiltonian()). This also ensures that the resultant expression is more readable.

[jkochNU]

This has been addressed with 71fa720 and released with v3.0.1. Issue #145 with external fluxes is still open, but is separate, I believe.