mathematical optimization - How can I speed up the classic GA for graph coloring?

I'm trying to compute the chromatic number of this graph (which is 28):

g = Import@"http://www.info.univ-angers.fr/pub/porumbel/graphs/dsjc250.5.col";

My genetic algorithm is getting stuck at an upper bound of 38 vertex colors:

In[] := Timing @ GAColor[g, 10, 20, 3]

Out[293]= {19.178072, {38, {28, 16, ...., 3, 22}}}

I've written the general GA implementation, but I'm using naive recombination and mutation, and my mathematica code is slow. My question is, how could I improve on this with more clever choices of Combine[] and Mutate[], as well as faster code, in the general? I'm by no means an expert here, so I'm sure there are many possible improvements both theoretically and algorithmically...

GAColor[g_Graph, PopulationSize_Integer:100, NumberOfGenerations_Integer:10, NumberOfMutants_Integer:0, mutationRadius_Integer:Automatic] :=  Module[
    {NumberOfVertices = VertexCount @ g, NumberOfBreeders, PermuteColorClasses, MutationRadius, Combine, Mutate, PopulationStep, 
    InitializePopulation, InitalPopulation, Generations, BestFitness = \[Infinity], BestColoring, GenerationsFitness, Chromatize, Adjacencies, result},

    MutationRadius = If[mutationRadius === Automatic, NumberOfVertices, mutationRadius];
    Adjacencies = Last /@ Transpose /@ GatherBy[First /@ Most[ArrayRules[AdjacencyMatrix[g]]], First];
    NumberOfBreeders = PopulationSize - NumberOfMutants;


    PermuteColorClasses[colors_, n_:1] := Module[{p},
        p = {#, Flatten @ Position[colors, #]}& /@ Union[colors];
        p[[All, 1]] = RandomSample[p[[All,1]]];
        ReplacePart[ConstantArray[0, Length[colors]], Flatten[Rule @@@ Thread[Reverse@#]& /@ p]]
    ];

    Chromatize[colorVector_] := Module[{f, h, min, co = colorVector},
        f = Function[{c,v},
            ReplacePart[c, v -> With[{ncols = c[[Adjacencies[[v]]]]},

                For[min = 1, MemberQ[ncols, min], min++]; min]
        ]];
        h = Function[{c}, Fold[f, c, RandomSample[Range[NumberOfVertices]]]];
        FixedPoint[h, co]
    ];

    Combine[colorVector1_, colorVector2_] := MapThread[RandomChoice[{#1, #2}]&, {colorVector1, colorVector2}];
    Mutate[colorVector_, mr_] := Permute[colorVector, RandomPermutation[mr]];

    PopulationStep[population_, NumberOfBreeders_] := Module[

        {fitness = Max /@ population, breeders, children, mutants},
        With[{min = Min[fitness]}, If[min < BestFitness, BestFitness = min]];
        breeders = RandomChoice[fitness -> population, NumberOfBreeders];
        children = Chromatize /@ Table[Combine @@ RandomChoice[breeders, 2], {NumberOfBreeders}
        ];
        mutants = Mutate[#, MutationRadius]& /@ RandomChoice[breeders, NumberOfMutants];
        Join[children, mutants]
    ];

    InitalPopulation = With[{color = Chromatize[RandomSample[Range @ NumberOfVertices]]},

            Table[PermuteColorClasses[color], {PopulationSize}]
    ];

    Generations = NestList[PopulationStep[#, NumberOfBreeders]&, InitalPopulation, NumberOfGenerations];
    GenerationsFitness = Map[Max, Generations, {2}];
    BestColoring = Extract[Generations, Position[GenerationsFitness, BestFitness, {2}, 1]][[1]];
    If[Or @@ (BestColoring[[First[#]]] == BestColoring[[Last[#]]]& /@ First /@ Most[ArrayRules[AdjacencyMatrix[g]]]),
        $Failed, {BestFitness, BestColoring}
    ]
]

For those with Mathematica 7 or Less

Here is code that doesn't use the version 8 Graph object, it's pretty much exactly the same:

GAColor[adjmatrix_, PopulationSize_Integer:100, NumberOfGenerations_Integer:10, NumberOfMutants_Integer:0, mutationRadius_Integer:Automatic] :=  Module[
        {NumberOfVertices = Length @ adjmatrix, NumberOfBreeders, PermuteColorClasses, MutationRadius, Combine, Mutate, PopulationStep, 
        InitializePopulation, InitalPopulation, Generations, BestFitness = \[Infinity], BestColoring, GenerationsFitness, Chromatize, Adjacencies, result},

        MutationRadius = If[mutationRadius === Automatic, NumberOfVertices, mutationRadius];
        Adjacencies = Last /@ Transpose /@ GatherBy[First /@ Most[ArrayRules[adjmatrix]], First];
        NumberOfBreeders = PopulationSize - NumberOfMutants;


        PermuteColorClasses[colors_, n_:1] := Module[{p},
            p = {#, Flatten @ Position[colors, #]}& /@ Union[colors];
            p[[All, 1]] = RandomSample[p[[All,1]]];
            ReplacePart[ConstantArray[0, Length[colors]], Flatten[Rule @@@ Thread[Reverse@#]& /@ p]]
        ];

        Chromatize[colorVector_] := Module[{f, h, min, co = colorVector},
            f = Function[{c,v},
                ReplacePart[c, v -> With[{ncols = c[[Adjacencies[[v]]]]},

                    For[min = 1, MemberQ[ncols, min], min++]; min]
            ]];
            h = Function[{c}, Fold[f, c, RandomSample[Range[NumberOfVertices]]]];
            FixedPoint[h, co]
        ];

        Combine[colorVector1_, colorVector2_] := MapThread[RandomChoice[{#1, #2}]&, {colorVector1, colorVector2}];
        Mutate[colorVector_, mr_] := Permute[colorVector, RandomPermutation[mr]];

        PopulationStep[population_, NumberOfBreeders_] := Module[

            {fitness = Max /@ population, breeders, children, mutants},
            With[{min = Min[fitness]}, If[min < BestFitness, BestFitness = min]];
            breeders = RandomChoice[fitness -> population, NumberOfBreeders];
            children = Chromatize /@ Table[Combine @@ RandomChoice[breeders, 2], {NumberOfBreeders}
            ];
            mutants = Mutate[#, MutationRadius]& /@ RandomChoice[breeders, NumberOfMutants];
            Join[children, mutants]
        ];

        InitalPopulation = With[{color = Chromatize[RandomSample[Range @ NumberOfVertices]]},

                Table[PermuteColorClasses[color], {PopulationSize}]
        ];

        Generations = NestList[PopulationStep[#, NumberOfBreeders]&, InitalPopulation, NumberOfGenerations];
        GenerationsFitness = Map[Max, Generations, {2}];
        BestColoring = Extract[Generations, Position[GenerationsFitness, BestFitness, {2}, 1]][[1]];
        If[Or @@ (BestColoring[[First[#]]] == BestColoring[[Last[#]]]& /@ First /@ Most[ArrayRules[adjmatrix]]),
            $Failed, {BestFitness, BestColoring}
        ]
    ]

Here is the sample input graph (as a compressed adjacency matrix) to test it on: http://pastebin.com/t7gnTczD

The algorithm should give a chromatic number of 28 in a few seconds. Here are the other benchmarks: http://www.info.univ-angers.fr/pub/porumbel/graphs/

Even in version 8 of Mathematica there still are no tools to compute the chromatic number or index of a graph, let alone a fast upper bound. Here is an illustration of the simulated annealing that's going on inside the algorithm:

g = Uncompress@"1:eJzt...."; (* get this string from pastebin link *)
NumberOfVertices = Length @ g;
color = RandomSample[Range[NumberOfVertices], NumberOfVertices];
n = NumberOfVertices;
A = Last /@ Transpose /@ GatherBy[First /@ Most[ArrayRules[g]], First];


NeighborComplements = Function[c,
    Module[{p, n, nc, r},
        p = Flatten @ Position[color, c];
        n = A[[p]];
        nc = Map[color[[#]]&, n, {2}];
        Thread @ {p, Complement[Range[c], #, {c}]& /@ nc}
    ]
];

Chromatic[g_, n_, col_:Range[NumberOfVertices]] := Module[{c=col, f, h, slow, fast},

    f = Function[{c, v},    
        ReplacePart[c, 
            v -> Module[{i, com = Complement[Range[n], c[[A[[v]]]]]}, 
                RandomChoice[Join[com, {c[[v]]}]]
            ]
        ]
    ];
    h = Function[{c}, Fold[f, c, RandomSample[Range[NumberOfVertices], NumberOfVertices]]];
    NestWhile[h, c, Max[#]>n&]
];


AbsoluteTiming[Monitor[color = NestWhile[Chromatic[g, n-=1, #]&, color, (color=#;n>1)&],
    ListPlot[Sort @ color, PlotRange -> All, PlotLabel -> Max[color]]]]

This is an optimization problem, and I'm sure some of you know this area intimately. When you run this code you will see a plot of the color classes which decrease slowly to around 30 different colors for the 250 vertices, however this is only a local minimum, the global minimum and chromatic number of the graph is actually 28... so my code is inefficient, if you can design a completely new function, and/or use openCL or JavaLink that is ok too...