Implement robust clique search configuration and testing

The changes introduce new configuration options for clique searching, state tracking for search outcomes, and logic for handling unreachable targets. Tests are added to verify search modes, failure conditions, and the correctness of coloring bounds.

pkg/obigraph/clique.go

@@ -4,6 +4,7 @@ import (
 	"context"
 	"math/rand"
 	"sync"
+	"sync/atomic"
 )
 
 // NeighborFilter defines an implicit graph over the vertices 0..n-1.
@@ -25,8 +26,25 @@ type CliqueOptions struct {
 	// larger one is found. Calls are serialized and made while a lock is held:
 	// it must return quickly and must not call FindClique.
 	Progress func(best int)
+
+	// ColorLimit enables the coloring bound for the candidate lists of at most
+	// ColorLimit vertices (0 disables it). The vertices of a list are split into
+	// classes of pairwise non-adjacent vertices: a clique takes at most one
+	// vertex per class, so the number of classes bounds its size. Computing the
+	// classes costs up to the square of the list length in calls to the filter,
+	// which is why only short lists are colored.
+	ColorLimit int
+
+	// mode forces the kind of search (tests only): 0 alternates, modeFocused
+	// and modeImproving restrict the search to one of them.
+	mode int
 }
 
+const (
+	modeFocused   = 1
+	modeImproving = 2
+)
+
 // CliqueResult is the outcome of FindClique.
 type CliqueResult struct {
 	// Clique is the largest clique found (a clique of size >= target if Reached).
@@ -36,6 +54,9 @@ type CliqueResult struct {
 	// Exhausted is true if the whole search space has been explored. If Reached
 	// is false, it proves that no clique of the target size exists.
 	Exhausted bool
+	// Maximum is true if Clique is proven to be a maximum clique of the graph:
+	// no larger clique exists.
+	Maximum bool
 }
 
 // FindClique looks for a clique of at least target vertices in the graph of n
@@ -43,12 +64,28 @@ type CliqueResult struct {
 // search space is exhausted, or when ctx is done; in the last case, the largest
 // clique met so far is returned.
 //
-// The search is a depth first search with backtracking, pruned when the
+// The search is a depth first search with backtracking, which comes in two
-// current clique plus the remaining candidates cannot reach target. The first
+// flavours that the restarts alternate:
-// pass of the first worker follows the natural order of the vertices (greedy
+//
-// "lexicographic" construction); every other pass restarts from a random
+//   - the focused search prunes every branch that cannot reach target, which
-// order, with a work budget that doubles at each restart, so that the search
+//     is the most efficient way to reach it but gives poor cliques when target
-// is complete in the long run.
+//     is out of reach;
+//   - the improving search only prunes the branches that cannot beat the best
+//     clique found so far (by any worker), so it keeps improving that clique
+//     whatever target is. When it ends, the best clique is a maximum one.
+//
+// When a focused search proves that no clique of target vertices exists, every
+// search turns to the improving flavour, so that the largest clique is returned
+// (and proved to be so) instead of an arbitrary one.
+//
+// On short candidate lists (see CliqueOptions.ColorLimit), a coloring bound
+// prunes the branches that cannot be large enough, and orders the vertices so
+// that the most constrained ones are tried first.
+//
+// The first pass of the first worker follows the natural order of the vertices
+// (greedy "lexicographic" construction); every other pass restarts from a
+// random order, with a work budget that doubles at each restart, so that the
+// search is complete in the long run.
 func FindClique(ctx context.Context, n int, filter NeighborFilter, target int, opts CliqueOptions) CliqueResult {
 	if target <= 0 {
 		return CliqueResult{Reached: true, Exhausted: true}
@@ -64,14 +101,20 @@ func FindClique(ctx context.Context, n int, filter NeighborFilter, target int, o
 	var (
 		mutex     sync.Mutex
 		best      []int32
+		bestSize  atomic.Int32
 		exhausted bool
+		maximum   bool
 		wg        sync.WaitGroup
+		// set when a focused search has proved that target cannot be reached:
+		// every search then turns to finding the largest clique
+		unreachable atomic.Bool
 	)
 
 	publish := func(clique []int32) {
 		mutex.Lock()
 		if len(clique) > len(best) {
 			best = append(best[:0], clique...)
+			bestSize.Store(int32(len(best)))
 			if opts.Progress != nil {
 				opts.Progress(len(best))
 			}
@@ -82,9 +125,12 @@ func FindClique(ctx context.Context, n int, filter NeighborFilter, target int, o
 	for w := range workers {
 		wg.Go(func() {
 			s := &cliqueSearch{
-				ctx:    ctx,
+				ctx:         ctx,
-				filter: filter,
+				filter:      filter,
-				target: target,
+				target:      target,
+				bestSize:    &bestSize,
+				colorLimit:  opts.ColorLimit,
+				unreachable: &unreachable,
 				publish: func(c []int32) {
 					publish(c)
 					if len(c) >= target {
@@ -93,15 +139,24 @@ func FindClique(ctx context.Context, n int, filter NeighborFilter, target int, o
 				},
 			}
 			rng := rand.New(rand.NewSource(opts.Seed + int64(w)))
+			// the search starts from the end of the list: lowest vertex first
 			candidates := make([]int32, n)
 			for i := range candidates {
-				candidates[i] = int32(i)
+				candidates[i] = int32(n - 1 - i)
 			}
 
 			for restart := 0; ctx.Err() == nil; restart++ {
 				if w > 0 || restart > 0 {
 					rng.Shuffle(n, func(i, j int) { candidates[i], candidates[j] = candidates[j], candidates[i] })
 				}
+				switch opts.mode {
+				case modeFocused:
+					s.improving = false
+				case modeImproving:
+					s.improving = true
+				default:
+					s.improving = unreachable.Load() || (w+restart)%2 == 1
+				}
 				s.budget = int64(n) * int64(target) << min(restart, 20)
 				s.work = 0
 				s.aborted = false
@@ -110,8 +165,17 @@ func FindClique(ctx context.Context, n int, filter NeighborFilter, target int, o
 				if !s.aborted && ctx.Err() == nil {
 					mutex.Lock()
 					exhausted = true
+					if s.improving {
+						maximum = true
+					}
 					mutex.Unlock()
-					cancel()
+
+					if s.improving || opts.mode == modeFocused {
+						cancel()
+					} else {
+						// target is out of reach: look for the largest clique
+						unreachable.Store(true)
+					}
 				}
 			}
 		})
@@ -122,47 +186,143 @@ func FindClique(ctx context.Context, n int, filter NeighborFilter, target int, o
 		Clique:    best,
 		Reached:   len(best) >= target,
 		Exhausted: exhausted,
+		Maximum:   maximum,
 	}
 }
 
 type cliqueSearch struct {
-	ctx     context.Context
+	ctx        context.Context
-	filter  NeighborFilter
+	filter     NeighborFilter
-	target  int
+	target     int
-	publish func([]int32)
+	bestSize   *atomic.Int32 // size of the best clique found by any worker
+	publish    func([]int32)
+	colorLimit int
+
+	unreachable *atomic.Bool // target is proved out of reach
+
+	improving bool // see FindClique: the kind of the current pass
 
 	current []int32   // the clique under construction
 	buffers [][]int32 // candidate list of each depth
-	bestLen int
+	sorted  [][]int32 // candidate list of each depth, sorted by color
+	colors  [][]int32 // color (from 1) of each vertex of sorted
+	classes [][]int32 // scratch space of the coloring
+	probe   []int32   // scratch space of the coloring
 
 	budget  int64 // maximum work (number of candidates tested) for this restart
 	work    int64
 	aborted bool // budget or context exhausted before the end of the exploration
 }
 
+// threshold is the size a clique must reach to be worth exploring: the target
+// for the focused search, one more than the best clique for the improving one.
+func (s *cliqueSearch) threshold() int {
+	if s.improving {
+		return min(s.target, int(s.bestSize.Load())+1)
+	}
+	return s.target
+}
+
+// mustStop tells whether the current pass has to be interrupted: its work
+// budget is spent, or the pass must stop whatever its budget (see mustAbort).
+func (s *cliqueSearch) mustStop() bool {
+	return s.work > s.budget || s.mustAbort()
+}
+
+// mustAbort tells whether the pass has to be interrupted whatever its budget:
+// the context is done, or it is a focused pass and its target has been proved
+// out of reach.
+func (s *cliqueSearch) mustAbort() bool {
+	return s.ctx.Err() != nil || (!s.improving && s.unreachable.Load())
+}
+
+// colorSort splits candidates into classes of pairwise non-adjacent vertices
+// (greedy sequential coloring) and returns them class after class, with the
+// color of each vertex. Any clique takes at most one vertex of each class.
+//
+// The coloring can be long: it stops, and sets s.aborted, as soon as the search
+// has to be aborted. Its cost is not charged to the work budget of the pass,
+// which limits the exploration of the tree, not the bound computed at each node.
+func (s *cliqueSearch) colorSort(depth int, candidates []int32) (sorted, colors []int32) {
+	k := 0
+	for _, v := range candidates {
+		if s.mustAbort() {
+			s.aborted = true
+			return nil, nil
+		}
+		placed := false
+		for c := 0; c < k; c++ {
+			// v joins the class if it is adjacent to none of its members
+			s.probe = s.filter(v, s.classes[c], s.probe[:0])
+			if len(s.probe) == 0 {
+				s.classes[c] = append(s.classes[c], v)
+				placed = true
+				break
+			}
+		}
+		if !placed {
+			if k == len(s.classes) {
+				s.classes = append(s.classes, nil)
+			}
+			s.classes[k] = append(s.classes[k][:0], v)
+			k++
+		}
+	}
+
+	for len(s.sorted) <= depth {
+		s.sorted = append(s.sorted, nil)
+		s.colors = append(s.colors, nil)
+	}
+	sorted, colors = s.sorted[depth][:0], s.colors[depth][:0]
+	for c := 0; c < k; c++ {
+		for _, v := range s.classes[c] {
+			sorted = append(sorted, v)
+			colors = append(colors, int32(c+1))
+		}
+	}
+	s.sorted[depth], s.colors[depth] = sorted, colors
+	return sorted, colors
+}
+
+// search extends the current clique of depth vertices with vertices taken in
+// candidates, all adjacent to the vertices of the current clique. The list is
+// read from its end: the vertex at position i is combined with the vertices
+// before it, so that every clique is met once.
 func (s *cliqueSearch) search(depth int, candidates []int32) {
-	for i := range candidates {
+	var colors []int32
-		if depth+len(candidates)-i < s.target {
+	if s.colorLimit > 0 && len(candidates) > 1 && len(candidates) <= s.colorLimit {
+		candidates, colors = s.colorSort(depth, candidates)
+		if s.aborted {
 			return
 		}
-		if s.work > s.budget || s.ctx.Err() != nil {
+	}
+
+	for i := len(candidates) - 1; i >= 0; i-- {
+		// largest clique the vertices 0..i can make
+		bound := i + 1
+		if colors != nil {
+			bound = int(colors[i])
+		}
+		if depth+bound < s.threshold() {
+			return
+		}
+		if s.mustStop() {
 			s.aborted = true
 			return
 		}
 
 		v := candidates[i]
 		s.current = append(s.current[:depth], v)
-		if len(s.current) > s.bestLen {
+		if len(s.current) > int(s.bestSize.Load()) {
-			s.bestLen = len(s.current)
 			s.publish(s.current)
-			if s.bestLen >= s.target {
+			if len(s.current) >= s.target {
 				s.aborted = true // the search is over: do not report an exhaustion
 				return
 			}
 		}
 
-		rest := candidates[i+1:]
+		rest := candidates[:i]
-		if depth+1+len(rest) < s.target {
+		if depth+1+len(rest) < s.threshold() {
 			continue
 		}
 		for len(s.buffers) <= depth+1 {
@@ -172,7 +332,7 @@ func (s *cliqueSearch) search(depth int, candidates []int32) {
 		s.buffers[depth+1] = next
 		s.work += int64(len(rest))
 
-		if depth+1+len(next) < s.target {
+		if depth+1+len(next) < s.threshold() {
 			continue
 		}
 		s.search(depth+1, next)

pkg/obigraph/clique_test.go

@@ -161,3 +161,89 @@ func TestFindCliqueProgress(t *testing.T) {
 		}
 	}
 }
+
+// The improving search must end with a proven maximum clique, even when the
+// target is out of reach; the focused one only proves the absence of a clique
+// of the target size.
+func TestFindCliqueModes(t *testing.T) {
+	rng := rand.New(rand.NewSource(5))
+	for trial := 0; trial < 40; trial++ {
+		n := 8 + rng.Intn(10)
+		adj := randomGraph(rng, n, 0.3+0.5*rng.Float64())
+		maxSize := bruteMaxClique(adj)
+		filter := matrixFilter(adj)
+
+		for _, workers := range []int{1, 3} {
+			for _, limit := range []int{0, 5, 1 << 20} {
+				r := FindClique(context.Background(), n, filter, maxSize+1,
+					CliqueOptions{Workers: workers, ColorLimit: limit, mode: modeImproving})
+				if r.Reached || !r.Exhausted || !r.Maximum {
+					t.Fatalf("trial %d limit %d: improving search should prove the maximum: %+v", trial, limit, r)
+				}
+				if len(r.Clique) != maxSize || !isClique(adj, r.Clique) {
+					t.Fatalf("trial %d limit %d: expected a maximum clique of %d, got %v", trial, limit, maxSize, r.Clique)
+				}
+
+				r = FindClique(context.Background(), n, filter, maxSize+1,
+					CliqueOptions{Workers: workers, ColorLimit: limit, mode: modeFocused})
+				if r.Reached || !r.Exhausted || r.Maximum {
+					t.Fatalf("trial %d limit %d: focused search proves absence only: %+v", trial, limit, r)
+				}
+
+				r = FindClique(context.Background(), n, filter, maxSize,
+					CliqueOptions{Workers: workers, ColorLimit: limit, mode: modeImproving})
+				if !r.Reached || len(r.Clique) < maxSize || !isClique(adj, r.Clique) {
+					t.Fatalf("trial %d limit %d: target %d not reached by the improving search: %+v", trial, limit, maxSize, r)
+				}
+			}
+		}
+	}
+}
+
+// When the target is out of reach, the default search must end with a proven
+// maximum clique, even if the coloring bound proves the absence of the target
+// before any clique has been found.
+func TestFindCliqueUnreachableTargetGivesMaximum(t *testing.T) {
+	rng := rand.New(rand.NewSource(6))
+	for trial := 0; trial < 40; trial++ {
+		n := 8 + rng.Intn(10)
+		adj := randomGraph(rng, n, 0.3+0.5*rng.Float64())
+		maxSize := bruteMaxClique(adj)
+
+		for _, workers := range []int{1, 3} {
+			for _, limit := range []int{0, 5, 1 << 20} {
+				r := FindClique(context.Background(), n, matrixFilter(adj), maxSize+1,
+					CliqueOptions{Workers: workers, ColorLimit: limit})
+				if r.Reached || !r.Exhausted || !r.Maximum {
+					t.Fatalf("trial %d workers %d limit %d: expected a proven maximum: %+v", trial, workers, limit, r)
+				}
+				if len(r.Clique) != maxSize || !isClique(adj, r.Clique) {
+					t.Fatalf("trial %d workers %d limit %d: expected a maximum clique of %d, got %v",
+						trial, workers, limit, maxSize, r.Clique)
+				}
+			}
+		}
+	}
+}
+
+// A long coloring must not delay the deadline.
+func TestFindCliqueColoringHonoursDeadline(t *testing.T) {
+	const n = 3000
+	adj := randomGraph(rand.New(rand.NewSource(7)), n, 0.5)
+	base := matrixFilter(adj)
+	slow := func(v int32, candidates, out []int32) []int32 {
+		time.Sleep(200 * time.Microsecond)
+		return base(v, candidates, out)
+	}
+
+	ctx, cancel := context.WithTimeout(context.Background(), 300*time.Millisecond)
+	defer cancel()
+	start := time.Now()
+	r := FindClique(ctx, n, slow, 100, CliqueOptions{Workers: 2, ColorLimit: 1 << 20})
+	if elapsed := time.Since(start); elapsed > 2*time.Second {
+		t.Fatalf("deadline not honoured by the coloring: %v", elapsed)
+	}
+	if r.Reached || !isClique(adj, r.Clique) {
+		t.Fatalf("unexpected result %+v", r)
+	}
+}

pkg/obikmer/kmer_dissimilarity_test.go

@@ -111,3 +111,29 @@ func TestParseKmerMetric(t *testing.T) {
 		t.Fatal("unknown metric must be refused")
 	}
 }
+
+// With a target out of reach (17 4-mers at Hamming distance >= 3 do not exist,
+// 16 do), the search must return the best family it can build, not a poor one
+// stuck by the pruning on the target, and the coloring bound must prove that
+// the target is out of reach.
+func TestFindClique_BestWhenTargetOutOfReach(t *testing.T) {
+	const k, dmin = 4, 3
+	kmers := allKmers(k)
+	filter := obikmer.NewKmerDissimilarityFilter(kmers, k, obikmer.HammingMetric, dmin)
+
+	for _, limit := range []int{0, 1024} {
+		ctx, cancel := context.WithTimeout(context.Background(), time.Second)
+		r := obigraph.FindClique(ctx, len(kmers), filter, 17, obigraph.CliqueOptions{Workers: 2, ColorLimit: limit})
+		cancel()
+		if r.Reached {
+			t.Fatal("17 4-mers at Hamming distance >= 3 cannot exist (sphere-packing bound 19, true maximum 16)")
+		}
+		if len(r.Clique) < 16 {
+			t.Fatalf("limit %d: expected the family of 16 to be found, got %d (exhausted %v, maximum %v)",
+				limit, len(r.Clique), r.Exhausted, r.Maximum)
+		}
+		if limit > 0 && !r.Exhausted {
+			t.Fatalf("limit %d: the coloring bound should prove that 17 is out of reach", limit)
+		}
+	}
+}

pkg/obitools/obioligotag/obioligotag.go

@@ -16,6 +16,12 @@ import (
 	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq"
 )
 
+// cliqueColorLimit is the length of the longest list of candidates for which the
+// search computes its coloring bound. Measured on 7-mers and 8-mers, the larger
+// it is the better the families found in a given time, with little gain beyond
+// 1024 (while the cost of a coloring grows with the square of the length).
+const cliqueColorLimit = 1024
+
 // buildConstraints translates the command line options into constraints on
 // oligonucleotides of the given size.
 func buildConstraints(size int) (obikmer.KmerConstraints, error) {
@@ -166,21 +172,25 @@ func CLIOligoFamily() ([]string, error) {
 	progress := startCliqueProgress(target)
 	result := obigraph.FindClique(ctx, len(kmers), filter, target,
 		obigraph.CliqueOptions{
-			Workers:  obidefault.ParallelWorkers(),
+			Workers:    obidefault.ParallelWorkers(),
-			Seed:     1,
+			Seed:       1,
-			Progress: progress.Update,
+			ColorLimit: cliqueColorLimit,
+			Progress:   progress.Update,
 		})
 	progress.Stop(result.Reached)
 
 	switch {
 	case result.Reached:
 		log.Infof("family of %d oligonucleotides found", len(result.Clique))
+	case result.Maximum:
+		log.Warnf("no family of %d oligonucleotides exists with these constraints: "+
+			"the largest one has %d oligonucleotides (proven maximum)", target, len(result.Clique))
 	case result.Exhausted:
 		log.Warnf("no family of %d oligonucleotides exists with these constraints: "+
 			"the best family found has %d oligonucleotides", target, len(result.Clique))
 	default:
 		log.Warnf("time limit reached before finding a family of %d oligonucleotides: "+
-			"the best family found has %d oligonucleotides", target, len(result.Clique))
+			"the best family found has %d oligonucleotides (not proven maximum)", target, len(result.Clique))
 	}
 
 	family := make([]uint64, len(result.Clique))