diff --git a/pkg/obigraph/clique.go b/pkg/obigraph/clique.go index 3297f18..4631368 100644 --- a/pkg/obigraph/clique.go +++ b/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 -// current clique plus the remaining candidates cannot reach target. 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. +// The search is a depth first search with backtracking, which comes in two +// flavours that the restarts alternate: +// +// - the focused search prunes every branch that cannot reach target, which +// is the most efficient way to reach it but gives poor cliques when target +// 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, - filter: filter, - target: target, + ctx: ctx, + filter: filter, + 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 - filter NeighborFilter - target int - publish func([]int32) + ctx context.Context + filter NeighborFilter + target int + 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 { - if depth+len(candidates)-i < s.target { + var colors []int32 + 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 { - s.bestLen = len(s.current) + if len(s.current) > int(s.bestSize.Load()) { 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:] - if depth+1+len(rest) < s.target { + rest := candidates[:i] + 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) diff --git a/pkg/obigraph/clique_test.go b/pkg/obigraph/clique_test.go index 0d9f0ad..c71e2b0 100644 --- a/pkg/obigraph/clique_test.go +++ b/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) + } +} diff --git a/pkg/obikmer/kmer_dissimilarity_test.go b/pkg/obikmer/kmer_dissimilarity_test.go index 94a92e6..1c52306 100644 --- a/pkg/obikmer/kmer_dissimilarity_test.go +++ b/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) + } + } +} diff --git a/pkg/obitools/obioligotag/obioligotag.go b/pkg/obitools/obioligotag/obioligotag.go index bd6dbd7..f23326e 100644 --- a/pkg/obitools/obioligotag/obioligotag.go +++ b/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(), - Seed: 1, - Progress: progress.Update, + Workers: obidefault.ParallelWorkers(), + Seed: 1, + 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))