restart;

with(GraphTheory):

claw:=CompleteGraph(1,3):
g:=CompleteGraph(2,2,2,2):
#DrawGraph(claw,size=[200,200]);
#DrawGraph(g,size=[200,200]);
IsSubgraphInducedIsomorphic(claw,g);

G1:=PathGraph(3):
G2:=Graph([1,2,3]):AddEdge(G2,{{3,1},{1,2}}):
DrawGraph(G1,size=[200,200]);
DrawGraph(G2,size=[200,200]);

IsSubgraphInducedIsomorphic(G1,G2,matching=true);

G1:=PathGraph(14); # longest for HyperGraph(5);
G2:=SpecialGraphs:-HypercubeGraph(5);
#DrawGraph(G1);
#DrawGraph(G2);

t,m:=IsSubgraphInducedIsomorphic(G1,G2,matching=true);
IsIsomorphic(G1,InducedSubgraph(G2,map(rhs,m)));

G1:=CycleGraph(14); # longest for HyperGraph(5);
G2:=SpecialGraphs:-HypercubeGraph(5);
#DrawGraph(G1);
#DrawGraph(G2);

t,m:=IsSubgraphInducedIsomorphic(G1,G2,matching=true);
IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

CodeTools:-Usage(IsSubgraphInducedIsomorphic(CycleGraph(15),G2,matching=true));

randomize(740472561204):
G2 := RandomGraphs:-RandomGraph(10000, 0.1);
NumberOfEdges(G2);
p:=combinat:-randcomb(Vertices(G2),30):
G1 := RelabelVertices(InducedSubgraph(G2,p),[$1..nops(p)]);

t,m:=CodeTools:-Usage(IsSubgraphInducedIsomorphic(G1,G2,matching=true));
IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

randomize(22807444796519):
G2 := RandomGraphs:-RandomGraph(50, 0.5);
p:=combinat:-randcomb(Vertices(G2),11):
G1 := RelabelVertices(InducedSubgraph(G2,p),[$1..nops(p)]);

t,m:=CodeTools:-Usage(IsSubgraphInducedIsomorphic(G1,G2,matching=true));
IsIsomorphic(G1,InducedSubgraph(G2,rhs~(m)));

IsSubgraphInducedIsomorphic := proc(G1::Graph, G2::Graph, {matching::truefalse:=false})
        description "Determines if undirected loopless G1 is isomorphic to an induced subgraph of G2",
                        "and if so optionally returns the matching function for the vertices";
        local n1,n2,D1,W1,V1,A1orig,A1,D2,W2,V2,A2,perm,MM,Morder,r,Vd,V1remaining,
                depth,allV1,allV2,Mout,Mstate,expandmatching,Premaining,
                M,M1,M2,T1,T1sorted,T2,T1t,T2t,P,pair,consistent,cut,i,j;
        D1,W1,V1,A1orig := op(1..4,G1); # directedness,weightedness,vertices,neighbors
        D2,W2,V2,A2     := op(1..4,G2);
        if D1 <> 'undirected' or D2 <> 'undirected' then error "both graphs must be undirected" end if;
        n1 := nops(V1):                # number of vertices in G1
        n2 := nops(V2):
        if orseq(i in A1orig[i], i = 1..n1) or orseq(i in A2[i], i=1..n2) then error "graph(s) have self-loops" end if;
        allV1 := {$1..n1};                # all vertices in G1
        allV2 := {$1..n2};
        # simple rejections based on number of vertices and max degrees  
        if n1 > n2 then return false end if;
        if max(map(nops,A1orig)) > max(map(nops,A2)) then return false end if;
        # First presort vertices of G1 in special BFS order
        # We are not using labels so can presort in decreasing degree order outside outermost while loop
        V1remaining := sort([$1..n1], (a,b) -> nops(A1orig[a]) >= nops(A1orig[b]));
        Morder:=Array([]);
        MM:={};
        while numelems(V1remaining) <> 0 do # loop for each connected component
                r := V1remaining[1]; # root for BFS tree for this component
                Vd := [r]; # depth = 0
                while nops(Vd) > 0 do # for each depth
                        # sort according to decreasing connections to M and
                        # otherwise decreasing degree
                        if nops(Vd) > 1 then
                                Vd:=sort(Vd, proc(a,b)
                                                        local na := nops(A1orig[a] intersect MM),
                                                                 nb := nops(A1orig[b] intersect MM);
                                                        na > nb  or na = nb and nops(A1orig[a]) >= nops(A1orig[b])
                                                   end proc)
                        end if;
                        # add these to the order and update
                        Morder ,= Vd[];
                        MM := MM union {Vd[]};
                        V1remaining := remove(has, V1remaining, Vd);
                        # go one level deeper by finding new neighbors of Vd;
                        Vd:=convert(`union`(seq(A1orig[Vd])) minus MM, list);
                end do;
        end do;        
        perm:=[seq(Morder)]; # order of vertices we want.
        # We now reorder and relabel the graph so the chosen order becomes the natural order
        A1 := A1orig[perm];
        A1 := eval(A1, perm =~ [$1..n1]);
        # Start of main part of algorithm
        # current matching M is represented as a set of equations
        # that map G1 vertices (lhs) to G2 vertices (rhs), e.g., (2=4,3=5}
        #
        # expandmatching works out M1,M2,T1,T1sorted,T2,T1t,T2t for a given M
        expandmatching := proc(M) option remember;
                local
                        M1 := map(lhs,M), # set of verts in G1 in matching
                        M2 := map(rhs,M), # set of verts in G2 in matching
                        T1 := {seq(op(A1[i]), i in M1)} minus M1, # set of neighbors of M1 not in M1
                        # sort T1 in order of decreasing numbers of neighbors in M1
                        T1sorted := sort([T1[]],(a,b)->nops(A1[a] intersect M1) >= nops(A1[b] intersect M1)),
                        T2 := {seq(op(A2[i]), i in M2)} minus M2, # set of neighbors of M2 not in M2
                        # T1t = T1tilde are further away vertices
                        T1t := allV1 minus M1 minus T1,
                        T2t := allV2 minus M2 minus T2;
                M1,M2,T1,T1sorted,T2,T1t,T2t
        end proc;
        # Premaining is table indexed by a matching that gives
        # the remaining addition possibilities (P) not yet tried for this matching
        Premaining := table('sparse');
        M := {};
        Mout := 0;
        depth := 0;
        # main loop - each iteration assesses a candidate state M
        # in a depth-first search
        while depth >= 0 do
                M1,M2,T1,T1sorted,T2,T1t,T2t := expandmatching(M);
                if M1 = allV1 then # success
                        Mout := M;
                        break
                end if;
                # find remaining candidate pairs for addition to this state
                if Premaining[M] = 0 then # new M
                        if nops(T1) <> 0 then
                                P := {seq([T1sorted[1],j],j in T2)};
                        elif         nops(allV1 minus M1) > 0 then        
                                P := {seq([(allV1 minus M1)[1],j], j in (allV2 minus M2))};
                        else
                                P:={}
                        end if;
                        Premaining[M] := P;
                else
                        P := Premaining[M];
                end if;
                # check out possibilities and go deeper if one is found
                while numelems(P) > 0 do
                        pair := P[1];
                        P := P minus {pair};
                        Premaining[M] := P;
                        i,j := pair[];
                        # consistency check for [i,j]
                        # neighbors of i in M1 must correspond to neighbors of j in M2
                        consistent:=evalb(eval(A1[i] intersect M1, M) = (A2[j] intersect M2));
                        # cut rule
                        cut := evalb(nops(A2[j] intersect T2) < nops(A1[i] intersect T1)
                                or nops(A2[j] intersect T2t) < nops(A1[i] intersect T1t));
                        if consistent and not cut then
                                Mstate[depth] := M; # save state
                                depth := depth + 1;  
                                M := M union {i = j};
                                next 2
                        end if
                end do;
                # no possibilities, backtrack
                depth := depth - 1;
                M := Mstate[depth];
        end do;
        forget(expandmatching);
        # return values
        if Mout <> 0 then   
                if matching then
                        true, map(eq->V1[perm[lhs(eq)]]=V2[rhs(eq)], Mout)
                else
                        true
                    end if;
          else
                false
        end if;
end proc: