How to perform below operation in Qiskit?

Question

I want to implement the below equation in Qiskit.

$(A \otimes B).\rho.(B^\dagger \otimes A^\dagger)$ where $\rho$ is a density matrix and $A$ and $B$ are CNOT gates. $$ A=\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \end{bmatrix}$$ $$B=\begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & 1 & 0 \end{bmatrix}. $$ The final $pf$ I got using the codes below is an operator, not a density matrix, as $dot()$ returns an operator. Hence, I can not use $\texttt{state_fidelity()}$ on it. The dimension of $\rho$ is 16. I can not use $evolve()$ since it is not any normal operator usage. Please let me know if there is any other method I can implement the above operation. I will be grateful.

p1=[[1,0,0,0],[0,0,0,1],[0,0,1,0],[0,1,0,0]]        
p2=[[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]]        
pa1=qi.Operator(p1)     
pa2=qi.Operator(p2)         
pa3=(pa2.transpose()).tensor(pa1.transpose())
pf=((pa1.tensor(pa2)).dot(rho)).dot(pa3)    

where qi is :

import qiskit.quantum_info as qi

Answer 1

The constructor of Qiskit's DensityMatrix class accepts any object with a to_matrix() method. And since Operator class has a to_matrix(), you can simply convert pf to a DensityMatrix that can be used with state_fidelity():

F = state_fidelity(qi.DensityMatrix(pf), rho2)