(*tree16是所给的基础树, cycletree[p,v](Subscript[是在tree16的第v点上挂上一个C, p]).A是取邻接矩阵,DMat是度数对角阵*)
tree16=EdgeAdd[PathGraph[Range[10]],{3\[UndirectedEdge]11,11\[UndirectedEdge]12,6\[UndirectedEdge]13,13\[UndirectedEdge]14,14\[UndirectedEdge]15,9\[UndirectedEdge]16}];
cycletree[p_Integer,v_]:=GraphUnion[tree16,Graph[Join[{v\[UndirectedEdge]17,15+p\[UndirectedEdge]v},Table[16+k\[UndirectedEdge]17+k,{k,1,p-2}]]]];
A[G_Graph]:=Normal[AdjacencyMatrix[G]];
DMat[G_Graph]:=DiagonalMatrix@VertexDegree[G];
ψ[G_Graph,t_:t,μ_:μ]:=CharacteristicPolynomial[AdjacencyMatrix[G]-μ DMat[G],t]

ClassifyDegSeq=Select[GatherBy[#,Sort[VertexDegree[#]]&],Length[#]>1&]&;
ClassifyCharPoly=Select[GatherBy[#,CharacteristicPolynomial[AdjacencyMatrix[#],t]&],Length[#]>1&]&;
ClassifyPsi=Select[GatherBy[#,ψ],Length[#]>1&]&;
(*用函数FindDegSimMat来寻找G1和G2之间是否度相似,若否则函数输出的Det为0.*)
FindDegSimMat[G1_Graph,G2_Graph]:=Module[{n=VertexCount[G1],M},If[n!=VertexCount[G2],Abort[]];M=Array[(Subscript[m, #1,#2])&,{n,n}];M=M/.Solve[{A[G1].M==M.A[G2],DMat[G1].M==M.DMat[G2]}][[1]]];
FindDegSimMatDet[G1_Graph,G2_Graph]:=Det[FindDegSimMat[G1,G2]];
FindDegSimMatDet[G1_Graph,G2_Graph]:=Det[FindDegSimMat[G1,G2]];
matrixRowSumQ=Not[MemberQ[Total/@#,0]]&;
FindDegSimMat[list_List]:=FindDegSimMat@@@Subsets[list,{2}];
FindDegSimMatDet[list_List]:=FindDegSimMatDet@@@Subsets[list,{2}];
tree16DegSimMat[p_Integer,v1_,v2_]:=FindDegSimMat[cycletree[p,v1],cycletree[p,v2]];
(*https://oeis.org/A001429*)
NumUniGraphs={1,2,5,13,33,89,240,657,1806,5026,13999,39260,110381,311465,880840,2497405,7093751,20187313,57537552,164235501,469406091,1343268050,3848223585,11035981711,31679671920,91021354454,261741776369,753265624291,2169441973139,6252511838796}[[#-2]]&;

(*例子: 一对度相似但不同构的图*)
X11=Graph[{1\[UndirectedEdge]2,3\[UndirectedEdge]4,5\[UndirectedEdge]6,6\[UndirectedEdge]7,7\[UndirectedEdge]8,8\[UndirectedEdge]9,1\[UndirectedEdge]6,2\[UndirectedEdge]7,3\[UndirectedEdge]8,4\[UndirectedEdge]9,1\[UndirectedEdge]10,3\[UndirectedEdge]10,1\[UndirectedEdge]9,2\[UndirectedEdge]5,2\[UndirectedEdge]6,3\[UndirectedEdge]7,4\[UndirectedEdge]5,4\[UndirectedEdge]8,5\[UndirectedEdge]9}];
X12=Graph[{1\[UndirectedEdge]2,1\[UndirectedEdge]5,1\[UndirectedEdge]7,1\[UndirectedEdge]8,2\[UndirectedEdge]8,2\[UndirectedEdge]9,2\[UndirectedEdge]10,3\[UndirectedEdge]4,3\[UndirectedEdge]5,3\[UndirectedEdge]6,3\[UndirectedEdge]9,4\[UndirectedEdge]6,4\[UndirectedEdge]7,4\[UndirectedEdge]10,5\[UndirectedEdge]6,5\[UndirectedEdge]9,6\[UndirectedEdge]7,7\[UndirectedEdge]8,8\[UndirectedEdge]9}];
FindDegSimMat[X11,X12]//Det//FullSimplify

rawList=Import["D:\\Projects\\sage\\DegreeSimilar\\unicyc_graph_classify_15.txt","List"];
StringTake[#,{2,-2}]&/@rawList
(*g6List=StringTake[#,{3,-3}]&/@rawList;
ImportString[#,"Graph6"]&/@g6List*)

timeUse[12]=3.9394;
timeUse[13]=10.8796;
timeUse[14]=30.8646;
timeUse[15]=82.476;
Table[((1000timeUse[k])/(NumUniGraphs[k])),{k,12,15}]

FindDegSimMat[X11,X12]

Total/@%