can you find the total chromatic number (edge and vertices) of a graph?

edit

You can solve this problem using mixed linear integer prrogramming, as follows:

• loop over the number n of colors
• for each such n, add n binary variables to each vertex and to each edge: bv[v,c] and be[e,c], where v is a vertex, e is an edge, and 0<=c<=n-1 is an integer. The meaning is: bv[v,c] is 1 iff the vertex v has color c, be[e,c] is 1 iff the edge e has color c.
• add the following constraints:
  • for each vertex v, the sum of bv[v,c] for c<=n-1 is equal to 1 (meaning each vertex get exactely one color)
  • for each edge, the sum of be[e,c] for c<=n-1 is equal to 1 (meaning each edge get exactly one color)
  • for each vertex v, for each c<=n-1, bv[v,c] plus the sum of be[e,c] where e is an edge adjacent to v is at most 1 (coloring condition around vertices)
  • for each edge e = (v,w), for each c<=n-1, bv[v,c]+bv[w,c]+be[e,c] is at most 1 (coloring condition around edges)
• try to solve, if you succeed, you get your coloring, if you get an exception, increase n by 1 and iterate

EDIT Here is a possible implementation;

sage: from sage.numerical.mip import MIPSolverException
sage: def total_chromatic_number(G, certificate=False):
....:     nmax = len(G) + len(G.edges()) # trivial upper bound on the number of colors.
....:     for n in range(1,nmax+1):
....:         p = MixedIntegerLinearProgram()
....:         bv = p.new_variable(binary=True)
....:         be = p.new_variable(binary=True)
....:         for v in G.vertices():
....:             p.add_constraint(sum(bv[v,c] for c in range(n)) == 1)
....:         for e in G.edges(labels=False):
....:             p.add_constraint(sum(be[e,c] for c in range(n)) == 1)
....:         for v in G.vertices():
....:             for c in range(n):
....:                 p.add_constraint(bv[v,c] + sum(be[e,c] for e in G.edges_incident(v, labels=False)) <= 1)
....:         for v,w in G.edges(labels=False):
....:             for c in range(n):
....:                 p.add_constraint(bv[v,c] + bv[w,c] + be[(v,w),c] <= 1)
....:         try:    
....:             p.solve()
....:             if certificate:
....:                 bv_sol = p.get_values(bv)
....:                 be_sol = p.get_values(be)
....:                 coloration = {}
....:                 for v in G.vertices():
....:                     for c in range(n):
....:                         if bv_sol[v,c] == 1:
....:                             coloration[v] = c
....:                 for e in G.edges(labels=False):
....:                     for c in range(n):
....:                         if be_sol[e,c] == 1:
....:                             coloration[e] = c
....:                 return coloration
....:             else:
....:                 return n
....:         except MIPSolverException:
....:             pass

We can test in on the Petersen graph, and see that its total chromatic number is 4:

sage: G = graphs.PetersenGraph()
sage: total_chromatic_number(G)
4

And get an admissible total coloration:

sage: total_chromatic_number(G, certificate=True)
{0: 0,
 1: 1,
 2: 3,
 3: 0,
 4: 1,
 5: 3,
 6: 0,
 7: 0,
 8: 2,
 9: 2,
 (0, 1): 3,
 (0, 4): 2,
 (0, 5): 1,
 (1, 2): 0,
 (1, 6): 2,
 (2, 3): 2,
 (2, 7): 1,
 (3, 4): 3,
 (3, 8): 1,
 (4, 9): 0,
 (5, 7): 2,
 (5, 8): 0,
 (6, 8): 3,
 (6, 9): 1,
 (7, 9): 3}