diff --git a/pkg/obialign/fastlcs_bench_test.go b/pkg/obialign/fastlcs_bench_test.go new file mode 100644 index 0000000..1c4c554 --- /dev/null +++ b/pkg/obialign/fastlcs_bench_test.go @@ -0,0 +1,111 @@ +package obialign + +import ( + "fmt" + "math/rand" + "testing" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" +) + +// randomDNA returns a random ACGT sequence of length n. +func randomDNA(rng *rand.Rand, n int) []byte { + const nuc = "ACGT" + s := make([]byte, n) + for i := range s { + s[i] = nuc[rng.Intn(4)] + } + return s +} + +// mutate returns a copy of s in which each position is, with probability +// rate, substituted, deleted or followed by an inserted base (equiprobable). +func mutate(rng *rand.Rand, s []byte, rate float64) []byte { + const nuc = "ACGT" + out := make([]byte, 0, len(s)+len(s)/10) + for _, c := range s { + if rng.Float64() >= rate { + out = append(out, c) + continue + } + switch rng.Intn(3) { + case 0: + out = append(out, nuc[rng.Intn(4)]) + case 1: // deletion + case 2: + out = append(out, c, nuc[rng.Intn(4)]) + } + } + return out +} + +// benchPairs builds a fixed set of pairs. "similar" pairs differ by about +// 2% of mutated positions, "random" pairs are unrelated sequences. +func benchPairs(n, length int, similar bool) [][2]*obiseq.BioSequence { + rng := rand.New(rand.NewSource(42)) + pairs := make([][2]*obiseq.BioSequence, n) + for i := range pairs { + a := randomDNA(rng, length) + var b []byte + if similar { + b = mutate(rng, a, 0.02) + } else { + b = randomDNA(rng, length) + } + pairs[i] = [2]*obiseq.BioSequence{ + obiseq.NewBioSequence("a", a, ""), + obiseq.NewBioSequence("b", b, ""), + } + } + return pairs +} + +// BenchmarkLegacyLCS measures the historical implementation (the one +// FastLCSScore used before BandedLCS), kept as the "before" reference. +func BenchmarkLegacyLCS(b *testing.B) { + for _, similar := range []bool{true, false} { + kind := "random" + if similar { + kind = "similar" + } + for _, length := range []int{100, 500, 2000} { + pairs := benchPairs(64, length, similar) + for _, maxErr := range []int{-1, 5, length / 20, length / 5} { + name := fmt.Sprintf("%s/len=%d/maxErr=%d", kind, length, maxErr) + b.Run(name, func(b *testing.B) { + var buffer []uint64 + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + p := pairs[i%len(pairs)] + FastLCSEGFScoreByte(p[0].Sequence(), p[1].Sequence(), maxErr, false, &buffer) + } + }) + } + } + } +} + +// BenchmarkBandedLCS runs on the same pairs as BenchmarkLegacyLCS. +func BenchmarkBandedLCS(b *testing.B) { + for _, similar := range []bool{true, false} { + kind := "random" + if similar { + kind = "similar" + } + for _, length := range []int{100, 500, 2000} { + pairs := benchPairs(64, length, similar) + for _, maxErr := range []int{-1, 5, length / 20, length / 5} { + name := fmt.Sprintf("%s/len=%d/maxErr=%d", kind, length, maxErr) + b.Run(name, func(b *testing.B) { + b.ReportAllocs() + b.ResetTimer() + for i := 0; i < b.N; i++ { + p := pairs[i%len(pairs)] + BandedLCS(p[0].Sequence(), p[1].Sequence(), maxErr) + } + }) + } + } + } +} diff --git a/pkg/obialign/fastlcsegf.go b/pkg/obialign/fastlcsegf.go index bce5d06..691f79b 100644 --- a/pkg/obialign/fastlcsegf.go +++ b/pkg/obialign/fastlcsegf.go @@ -2,6 +2,7 @@ package obialign import ( "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" + log "github.com/sirupsen/logrus" ) // FastLCSEGFScoreByte calculates the score of the Longest Common Subsequence (LCS) between two byte slices. @@ -355,27 +356,35 @@ func FastLCSEGFScore(seqA, seqB *obiseq.BioSequence, maxError int, buffer *[]uin return FastLCSEGFScoreByte(seqA.Sequence(), seqB.Sequence(), maxError, true, buffer) } -// FastLCSScore calculates the score of the longest common subsequence between two bio sequences. +// FastLCSScore calculates the score of the longest common subsequence between two bio sequences, +// and the length of the shortest alignment achieving it (the distance is alilen - lcs). // -// if maxError > 0, the maximum allowed error between the sequences is maxError. -// Otherwise, no error checking is done. -// If the actual number of errors is larger than maxError, -1 is returned for both values. +// If maxError >= 0, the result is exact when its distance is at most maxError, +// and (-1, -1) is returned otherwise. If maxError < 0, no limit is applied. +// See BandedLCS, which does the work. // // The score matrix is: -// - Matching: 1 +// - Matching: 1 (bases are compared with the IUPAC rules) // - Mismatch or gap: 0 // +// The program stops with a fatal error if a sequence is longer than +// MaxBandedLCSLength. +// // Parameters: // - seqA: The first bio sequence. // - seqB: The second bio sequence. // - maxError: The maximum allowed error between the sequences. If set to -1, no limit is applied. -// - buffer: A pointer to a uint64 slice to store intermediate results. If nil, a new slice is created. +// - buffer: Unused, kept for compatibility; the working memory is pooled. // // Returns: // - The score of the longest common subsequence. // - The length of the shortest alignment corresponding to the LCS. func FastLCSScore(seqA, seqB *obiseq.BioSequence, maxError int, buffer *[]uint64) (int, int) { - score, alilen, _ := FastLCSEGFScoreByte(seqA.Sequence(), seqB.Sequence(), maxError, false, buffer) + score, alilen, err := BandedLCS(seqA.Sequence(), seqB.Sequence(), maxError) + if err != nil { + log.Fatalf("FastLCSScore on sequences %s (%d bp) and %s (%d bp): %v (maximum %d bp)", + seqA.Id(), seqA.Len(), seqB.Id(), seqB.Len(), err, MaxBandedLCSLength) + } return score, alilen } diff --git a/pkg/obialign/lcsband.go b/pkg/obialign/lcsband.go new file mode 100644 index 0000000..3c0c095 --- /dev/null +++ b/pkg/obialign/lcsband.go @@ -0,0 +1,175 @@ +package obialign + +import ( + "errors" + "sync" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" +) + +// MaxBandedLCSLength is the longest sequence BandedLCS accepts. The score and +// the alignment length are packed on 16 bits each, and the alignment length +// can reach the sum of both sequence lengths. +const MaxBandedLCSLength = 1<<15 - 1 + +// ErrSequenceTooLong is returned by BandedLCS when a sequence is longer than +// MaxBandedLCSLength. +var ErrSequenceTooLong = errors.New("obialign: sequence too long for BandedLCS") + +// A DP cell is packed on 32 bits as +// +// score<<16 | (lcsLenBase - alignment length) +// +// so that the maximum of two cells ranks them by score first and, for equal +// scores, by shortest alignment. The three moves then reduce to +// +// cell(i,j) = max(diag + match<<16, up, left) - 1 +// +// since every move lengthens the alignment by one and only a diagonal match +// adds one to the score. The value 0 stands for an unavailable cell (outside +// the band): it loses against any real cell. Every cell of the band has at +// least one real neighbour (the band is connected and contains (0,0)), so +// the subtraction never wraps for a cell that is used; the scalar reference +// still clamps the maximum to 1 as a safety net. +const ( + lcsScoreUnit = uint32(1) << 16 + lcsLenBase = uint32(0xFFFF) +) + +// bandedLCSWorkspace holds the scratch memory of one alignment. +type bandedLCSWorkspace struct { + aRev, bCode []uint32 // encoded bases (aRev is reversed), SIMD kernel only + rows [3][]uint32 // antidiagonal (SIMD) or row (scalar) buffers + widened int // number of alignments that needed a second pass +} + +var bandedLCSPool = sync.Pool{New: func() any { return new(bandedLCSWorkspace) }} + +func growUint32(s []uint32, n int) []uint32 { + if cap(s) < n { + return make([]uint32, n) + } + return s[:n] +} + +// iupacMask gives for each byte its 4-bit nucleotide set (A=1, C=2, G=4, +// T=8), or 0 when the byte is not an IUPAC letter. For two bytes of non-zero +// mask, SameIUPACNuc(a, b) is exactly (mask[a] & mask[b]) != 0. +var iupacMask = buildIUPACMask() + +func buildIUPACMask() (m [256]uint8) { + probes := [4]byte{'a', 'c', 'g', 't'} + for c := 0; c < 256; c++ { + if !(c >= 'a' && c <= 'z' || c >= 'A' && c <= 'Z') { + continue + } + for bit, p := range probes { + if obiseq.SameIUPACNuc(byte(c), p) { + m[c] |= 1 << uint(bit) + } + } + } + return m +} + +// allRegularIUPAC reports whether every base has a non-zero IUPAC mask, in +// which case base compatibility is a plain AND of masks. +func allRegularIUPAC(s []byte) bool { + for _, c := range s { + if iupacMask[c] == 0 { + return false + } + } + return true +} + +// bandOffsets returns the diagonal offsets k = j - i covered by a band of +// half-width g around the diagonals 0..delta, with sequence A (columns j) not +// shorter than sequence B (rows i), clipped to the matrix. +func bandOffsets(lA, lB, g int) (kmin, kmax int) { + return max(-g, -lB), min(lA-lB+g, lA) +} + +// BandedLCS computes the score of the longest common subsequence of two +// nucleotide sequences, and the length of the shortest alignment achieving +// it (match: +1, mismatch and gap: 0; bases are compared with the IUPAC +// rules of obiseq.SameIUPACNuc). The distance is alilen - lcs. +// +// If maxError >= 0, the result is the optimal (lcs, alilen) when its distance +// is at most maxError, and (-1, -1) otherwise. This is exact, not a +// heuristic: the DP is restricted to a band, and a bound on what a path +// leaving the band can score proves when the band was wide enough. maxError < 0 +// means no limit. +// +// An error is returned when a sequence is longer than MaxBandedLCSLength. +func BandedLCS(seqA, seqB []byte, maxError int) (lcs, alilen int, err error) { + if len(seqA) > MaxBandedLCSLength || len(seqB) > MaxBandedLCSLength { + return -1, -1, ErrSequenceTooLong + } + if len(seqA) < len(seqB) { + seqA, seqB = seqB, seqA + } + ws := bandedLCSPool.Get().(*bandedLCSWorkspace) + lcs, alilen = bandedLCS(seqA, seqB, maxError, ws, haveSIMD) + bandedLCSPool.Put(ws) + return lcs, alilen, nil +} + +// bandedLCS expects len(a) >= len(b), both within MaxBandedLCSLength. +// +// With delta = lA - lB and a band of half-width g (offsets -g .. delta+g): +// - an alignment with at most E errors needs at least delta+2h gaps to reach +// the offset delta+h or -h, so it lies in the band of g0 = (E-delta)/2; +// - a path leaving a band of half-width g makes at least g+1 gaps on the B +// side, hence scores at most lB-g-1 matches, and has more than E errors +// as soon as g >= g0. +// +// Let Q be the best path of the band g0. If its distance exceeds E, the +// global optimum cannot be within E either (it would be in the band and be +// Q): the answer is -1. Otherwise, if score(Q) >= lB-g0-1 no outside path can +// beat Q (a tie in score loses on length, the outside path having more errors +// than Q): Q is the global optimum. Else the band is widened once to +// g1 = min(E-delta, lB-score(Q)-1); since score(Q) >= lA-E for any alignment +// within E errors, g1 always certifies the result. +func bandedLCS(a, b []byte, maxError int, ws *bandedLCSWorkspace, useSIMD bool) (int, int) { + lA, lB := len(a), len(b) + delta := lA - lB + if maxError >= 0 && maxError < delta { + return -1, -1 + } + if lB == 0 { + return 0, lA + } + simdOK := useSIMD && allRegularIUPAC(a) && allRegularIUPAC(b) + if simdOK { + encodeBases(ws, a, b) + } + run := func(g int) (int, int) { + kmin, kmax := bandOffsets(lA, lB, g) + var key uint32 + if simdOK { + key = bandedLCSSIMD(lA, lB, kmin, kmax, ws) + } else { + key = bandedLCSScalar(a, b, kmin, kmax, ws) + } + return int(key >> 16), int(lcsLenBase - key&0xFFFF) + } + + if maxError < 0 { + return run(lA) + } + g0 := (maxError - delta) / 2 + lcs, alilen := run(g0) + if alilen-lcs > maxError { + return -1, -1 + } + if lcs >= lB-g0-1 { + return lcs, alilen + } + ws.widened++ + lcs, alilen = run(min(maxError-delta, lB-lcs-1)) + if alilen-lcs > maxError { + return -1, -1 + } + return lcs, alilen +} diff --git a/pkg/obialign/lcsband_nosimd.go b/pkg/obialign/lcsband_nosimd.go new file mode 100644 index 0000000..818a5e2 --- /dev/null +++ b/pkg/obialign/lcsband_nosimd.go @@ -0,0 +1,13 @@ +//go:build !goexperiment.simd + +package obialign + +const haveSIMD = false + +// encodeBases and bandedLCSSIMD are never called without GOEXPERIMENT=simd +// (haveSIMD is false); they exist to keep the dispatcher uniform. +func encodeBases(ws *bandedLCSWorkspace, a, b []byte) {} + +func bandedLCSSIMD(lA, lB, kmin, kmax int, ws *bandedLCSWorkspace) uint32 { + return 0 +} diff --git a/pkg/obialign/lcsband_scalar.go b/pkg/obialign/lcsband_scalar.go new file mode 100644 index 0000000..d6220d0 --- /dev/null +++ b/pkg/obialign/lcsband_scalar.go @@ -0,0 +1,53 @@ +package obialign + +import "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" + +// bandedLCSScalar is the portable reference of the banded DP, row by row. It +// returns the packed cell (lB, lA). len(a) >= len(b) >= 1. +// +// Cells outside the band hold 0 (unavailable). Row i spans the columns +// [i+kmin, i+kmax] clipped to [1, lA]; the cell just left of that span and +// the one just right of it are reset so that the next row never reads a stale +// value. +func bandedLCSScalar(a, b []byte, kmin, kmax int, ws *bandedLCSWorkspace) uint32 { + lA, lB := len(a), len(b) + + prev := growUint32(ws.rows[0], lA+2) + cur := growUint32(ws.rows[1], lA+2) + ws.rows[0], ws.rows[1] = prev, cur + + for j := 0; j <= lA; j++ { + prev[j] = 0 + if j <= kmax { + prev[j] = lcsLenBase - uint32(j) + } + } + prev[lA+1] = 0 + + for i := 1; i <= lB; i++ { + if -i >= kmin { + cur[0] = lcsLenBase - uint32(i) + } else { + cur[0] = 0 + } + jlo := max(1, i+kmin) + jhi := min(lA, i+kmax) + if jlo > 1 { + cur[jlo-1] = 0 + } + bi := b[i-1] + for j := jlo; j <= jhi; j++ { + d := prev[j-1] + if obiseq.SameIUPACNuc(a[j-1], bi) { + d += lcsScoreUnit + } + m := max(d, prev[j], cur[j-1], 1) + cur[j] = m - 1 + } + if jhi < lA { + cur[jhi+1] = 0 + } + prev, cur = cur, prev + } + return prev[lA] +} diff --git a/pkg/obialign/lcsband_simd.go b/pkg/obialign/lcsband_simd.go new file mode 100644 index 0000000..42b4a09 --- /dev/null +++ b/pkg/obialign/lcsband_simd.go @@ -0,0 +1,107 @@ +//go:build goexperiment.simd + +package obialign + +import "simd" + +const haveSIMD = true + +func simdLanes() int { + var probe simd.Uint32s + return probe.Len() +} + +// encodeIUPAC writes in dst the 4-bit IUPAC mask of each base of s, reversed +// if rev is set, followed by zero padding up to len(dst). +func encodeIUPAC(dst []uint32, s []byte, rev bool) { + n := len(s) + for x := range dst { + dst[x] = 0 + } + if rev { + for x, c := range s { + dst[n-1-x] = uint32(iupacMask[c]) + } + return + } + for x, c := range s { + dst[x] = uint32(iupacMask[c]) + } +} + +// encodeBases encodes both sequences once for all the passes of an alignment. +func encodeBases(ws *bandedLCSWorkspace, a, b []byte) { + vl := simdLanes() + ws.aRev = growUint32(ws.aRev, len(a)+vl) + ws.bCode = growUint32(ws.bCode, len(b)+vl) + encodeIUPAC(ws.aRev, a, true) + encodeIUPAC(ws.bCode, b, false) +} + +// bandedLCSSIMD runs the banded DP along antidiagonals t = i + j. Lane i of +// the buffer of antidiagonal t holds the cell (i, t-i), so that +// +// diag = buffer[t-2][i-1] up = buffer[t-1][i-1] left = buffer[t-1][i] +// +// are plain unaligned loads. Along an antidiagonal the offset j-i = t-2i +// must stay in [kmin, kmax], which gives the window [lo, hi] of rows. The +// vector loop may run past hi; the lanes beyond are then reset to 0 +// (unavailable), as is the lane just before the window, because the next +// antidiagonals read them as neighbours. Lane 0 and lane t carry the first +// row and the first column of the matrix. +// +// Requires encodeBases to have been called, lA >= lB >= 1 and kmin <= 0 <= +// kmax. +func bandedLCSSIMD(lA, lB, kmin, kmax int, ws *bandedLCSWorkspace) uint32 { + vl := simdLanes() + aRev, bCode := ws.aRev, ws.bCode + + size := lB + 2*vl + 2 + for r := range ws.rows { + ws.rows[r] = growUint32(ws.rows[r], size) + clear(ws.rows[r]) + } + p2, p1, cur := ws.rows[0], ws.rows[1], ws.rows[2] + + p2[0] = lcsLenBase // cell (0,0) + if kmax >= 1 { + p1[0] = lcsLenBase - 1 // cell (0,1) + } + if kmin <= -1 { + p1[1] = lcsLenBase - 1 // cell (1,0) + } + + zero := simd.BroadcastUint32s(0) + one := simd.BroadcastUint32s(1) + unit := simd.BroadcastUint32s(lcsScoreUnit) + + for t := 2; t <= lA+lB; t++ { + lo := max(1, t-lA, (t-kmax+1)>>1) + hi := min(lB, t-1, (t-kmin)>>1) + + for i := lo; i <= hi; i += vl { + ca := simd.LoadUint32s(aRev[lA-t+i:]) + cb := simd.LoadUint32s(bCode[i-1:]) + inc := unit.Masked(ca.And(cb).NotEqual(zero)) + d := simd.LoadUint32s(p2[i-1:]).Add(inc) + u := simd.LoadUint32s(p1[i-1:]) + l := simd.LoadUint32s(p1[i:]) + d.Max(u).Max(l).Sub(one).Store(cur[i:]) + } + + zero.Store(cur[hi+1:]) + if lo > 1 { + cur[lo-1] = 0 + } + cur[0] = 0 + if t <= lA && t <= kmax { + cur[0] = lcsLenBase - uint32(t) + } + if t <= lB && -t >= kmin { + cur[t] = lcsLenBase - uint32(t) + } + + p2, p1, cur = p1, cur, p2 + } + return p1[lB] +} diff --git a/pkg/obialign/lcsband_test.go b/pkg/obialign/lcsband_test.go new file mode 100644 index 0000000..f957717 --- /dev/null +++ b/pkg/obialign/lcsband_test.go @@ -0,0 +1,288 @@ +package obialign + +import ( + "math/rand" + "testing" + + "git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obiseq" +) + +// bandDP is the oracle: a plain DP over the whole matrix, written +// independently of the production code, in which only the cells with a +// diagonal offset j-i in [kmin, kmax] exist (kmin=-lB, kmax=lA when +// unbounded). It returns the optimal LCS score and the shortest alignment +// length achieving it among the alignments staying in the band. +func bandDP(a, b []byte, kmin, kmax int) (int, int) { + lA, lB := len(a), len(b) + type cell struct { + s, l int + ok bool + } + better := func(x, y cell) bool { + return !y.ok || (x.ok && (x.s > y.s || (x.s == y.s && x.l < y.l))) + } + m := make([][]cell, lB+1) + for i := range m { + m[i] = make([]cell, lA+1) + } + for i := 0; i <= lB; i++ { + for j := 0; j <= lA; j++ { + if k := j - i; k < kmin || k > kmax { + continue + } + if i == 0 || j == 0 { + m[i][j] = cell{0, i + j, true} + continue + } + var best cell + if d := m[i-1][j-1]; d.ok { + best = cell{d.s, d.l + 1, true} + if obiseq.SameIUPACNuc(a[j-1], b[i-1]) { + best.s++ + } + } + if u := m[i-1][j]; u.ok && better(cell{u.s, u.l + 1, true}, best) { + best = cell{u.s, u.l + 1, true} + } + if l := m[i][j-1]; l.ok && better(cell{l.s, l.l + 1, true}, best) { + best = cell{l.s, l.l + 1, true} + } + m[i][j] = best + } + } + return m[lB][lA].s, m[lB][lA].l +} + +// fullDP is the unbanded optimum. +func fullDP(a, b []byte) (int, int) { return bandDP(a, b, -len(b), len(a)) } + +// expected is the contract of BandedLCS: the unbounded optimum when it has at +// most maxError errors, (-1,-1) otherwise. +func expected(a, b []byte, maxError int) (int, int) { + s, l := fullDP(a, b) + if maxError >= 0 && l-s > maxError { + return -1, -1 + } + return s, l +} + +func randomFrom(rng *rand.Rand, alphabet string, n int) []byte { + s := make([]byte, n) + for i := range s { + s[i] = alphabet[rng.Intn(len(alphabet))] + } + return s +} + +func mutateFrom(rng *rand.Rand, alphabet string, s []byte, rate float64) []byte { + out := make([]byte, 0, len(s)+8) + for _, c := range s { + if rng.Float64() >= rate { + out = append(out, c) + continue + } + switch rng.Intn(3) { + case 0: + out = append(out, alphabet[rng.Intn(len(alphabet))]) + case 1: + case 2: + out = append(out, c, alphabet[rng.Intn(len(alphabet))]) + } + } + return out +} + +var testAlphabets = []struct { + name, letters string + regular bool // only IUPAC letters of non-zero mask +}{ + {"acgt", "ACGT", true}, + {"low-complexity", "AC", true}, + {"mixed-case", "ACGTacgt", true}, + {"iupac", "ACGTRYSWKMBDHVNUacgtrywsn", true}, + {"exotic", "ACGT-.EXz*", false}, +} + +// TestBandedLCSDifferential checks the scalar reference against the oracle, +// and the dispatcher (SIMD kernel when built with GOEXPERIMENT=simd) against +// the scalar reference, for random pairs of unequal lengths and many bands. +func TestBandedLCSDifferential(t *testing.T) { + rng := rand.New(rand.NewSource(7)) + ws := new(bandedLCSWorkspace) + for _, alpha := range testAlphabets { + for it := 0; it < 3000; it++ { + lA := 1 + rng.Intn(90) + a := randomFrom(rng, alpha.letters, lA) + var b []byte + if rng.Intn(3) > 0 { + b = mutateFrom(rng, alpha.letters, a, 0.02+rng.Float64()*0.5) + } else { + b = randomFrom(rng, alpha.letters, 1+rng.Intn(90)) + } + if len(b) == 0 { + continue + } + if len(a) < len(b) { + a, b = b, a + } + delta := len(a) - len(b) + for _, maxError := range []int{-1, 0, 1, 2, 3, 5, 8, 13, 21, 40, 200, delta, delta + 1, delta + 2, delta + 3} { + ws, wl := expected(a, b, maxError) + ss, sl := bandedLCS(a, b, maxError, new(bandedLCSWorkspace), false) + if ss != ws || sl != wl { + t.Fatalf("%s scalar maxError=%d a=%q b=%q: got (%d,%d), oracle (%d,%d)", + alpha.name, maxError, a, b, ss, sl, ws, wl) + } + ds, dl := bandedLCS(a, b, maxError, new(bandedLCSWorkspace), haveSIMD) + if ds != ss || dl != sl { + t.Fatalf("%s dispatch(simd=%v) maxError=%d a=%q b=%q: got (%d,%d), scalar (%d,%d)", + alpha.name, haveSIMD, maxError, a, b, ds, dl, ss, sl) + } + } + } + } + _ = ws +} + +// TestBandedLCSLongPairs covers lengths well above the vector width, with +// reused workspaces (stale buffers must never leak between calls). +func TestBandedLCSLongPairs(t *testing.T) { + rng := rand.New(rand.NewSource(11)) + ws := new(bandedLCSWorkspace) + for it := 0; it < 60; it++ { + a := randomFrom(rng, "ACGT", 150+rng.Intn(500)) + b := mutateFrom(rng, "ACGT", a, 0.01+rng.Float64()*0.2) + if len(a) < len(b) { + a, b = b, a + } + for _, maxError := range []int{-1, 4, 20, 60, 150} { + wantS, wantL := expected(a, b, maxError) + gotS, gotL := bandedLCS(a, b, maxError, ws, haveSIMD) + if gotS != wantS || gotL != wantL { + t.Fatalf("maxError=%d lens=%d/%d: got (%d,%d), oracle (%d,%d)", + maxError, len(a), len(b), gotS, gotL, wantS, wantL) + } + } + } +} + +// TestBandedLCSMatchesLegacy compares with the historical implementation +// (FastLCSEGFScoreByte, not end-gap-free) wherever that one is exact: no +// limit, or a result within the limit. +func TestBandedLCSMatchesLegacy(t *testing.T) { + rng := rand.New(rand.NewSource(3)) + for it := 0; it < 4000; it++ { + a := randomFrom(rng, "ACGTRYN", 1+rng.Intn(120)) + b := mutateFrom(rng, "ACGTRYN", a, rng.Float64()*0.4) + if len(b) == 0 { + continue + } + for _, maxError := range []int{-1, 2, 7, 25} { + ls, ll, err := BandedLCS(a, b, maxError) + if err != nil { + t.Fatal(err) + } + if ls < 0 { + continue + } + os, ol, _ := FastLCSEGFScoreByte(a, b, maxError, false, nil) + if os != ls || ol != ll { + t.Fatalf("maxError=%d a=%q b=%q: BandedLCS (%d,%d), legacy (%d,%d)", + maxError, a, b, ls, ll, os, ol) + } + } + } +} + +func TestBandedLCSEdgeCases(t *testing.T) { + cases := []struct { + a, b string + maxError int + s, l int + }{ + {"", "", -1, 0, 0}, + {"ACGT", "", -1, 0, 4}, + {"ACGT", "", 3, -1, -1}, + {"ACGT", "", 4, 0, 4}, + {"ACGT", "ACGT", 0, 4, 4}, + {"ACGT", "ACGA", 0, -1, -1}, + {"ACGT", "ACGA", 1, 3, 4}, + {"ACGTT", "ACGT", 0, -1, -1}, // length difference above the limit + {"ACGTT", "ACGT", 1, 4, 5}, + } + for _, c := range cases { + s, l, err := BandedLCS([]byte(c.a), []byte(c.b), c.maxError) + if err != nil || s != c.s || l != c.l { + t.Errorf("%q/%q maxError=%d: got (%d,%d,%v), want (%d,%d)", c.a, c.b, c.maxError, s, l, err, c.s, c.l) + } + // the arguments are symmetric + s, l, _ = BandedLCS([]byte(c.b), []byte(c.a), c.maxError) + if s != c.s || l != c.l { + t.Errorf("swapped %q/%q maxError=%d: got (%d,%d), want (%d,%d)", c.b, c.a, c.maxError, s, l, c.s, c.l) + } + } +} + +func TestBandedLCSTooLong(t *testing.T) { + long := make([]byte, MaxBandedLCSLength+1) + for i := range long { + long[i] = 'A' + } + if _, _, err := BandedLCS(long, []byte("ACGT"), -1); err != ErrSequenceTooLong { + t.Errorf("A too long: err = %v", err) + } + if _, _, err := BandedLCS([]byte("ACGT"), long, 5); err != ErrSequenceTooLong { + t.Errorf("B too long: err = %v", err) + } + // the maximum length itself is accepted + ok := long[:MaxBandedLCSLength] + if s, l, err := BandedLCS(ok, ok, 0); err != nil || s != len(ok) || l != len(ok) { + t.Errorf("max length: got (%d,%d,%v)", s, l, err) + } +} + +// TestBandedLCSSecondPass checks the exactness over many bands, also on +// low-complexity sequences where the best LCS is not the closest alignment, +// and reports how often the band had to be widened (never more than once). +func TestBandedLCSSecondPass(t *testing.T) { + rng := rand.New(rand.NewSource(5)) + for _, alpha := range testAlphabets[:4] { + for _, rate := range []float64{0.02, 0.1, 0.3} { + ws := new(bandedLCSWorkspace) + var total int + for it := 0; it < 1500; it++ { + a := randomFrom(rng, alpha.letters, 20+rng.Intn(120)) + b := mutateFrom(rng, alpha.letters, a, rate) + if len(b) == 0 { + continue + } + if len(a) < len(b) { + a, b = b, a + } + for _, maxError := range []int{1, 3, 6, 10, 20, 40} { + want, wantL := expected(a, b, maxError) + got, gotL := bandedLCS(a, b, maxError, ws, haveSIMD) + if got != want || gotL != wantL { + t.Fatalf("%s maxError=%d a=%q b=%q: got (%d,%d), want (%d,%d)", + alpha.name, maxError, a, b, got, gotL, want, wantL) + } + total++ + } + } + t.Logf("%-15s mutation %.2f: %d bounded calls, %d widened (%.1f%%)", + alpha.name, rate, total, ws.widened, 100*float64(ws.widened)/float64(total)) + } + } +} + +// TestFastLCSScoreWrapper checks the BioSequence entry point. +func TestFastLCSScoreWrapper(t *testing.T) { + sa := obiseq.NewBioSequence("a", []byte("ACGTACGTAC"), "") + sb := obiseq.NewBioSequence("b", []byte("ACGTTCGTAC"), "") + if s, l := FastLCSScore(sa, sb, -1, nil); s != 9 || l != 10 { + t.Errorf("got (%d,%d), want (9,10)", s, l) + } + if s, l := FastLCSScore(sa, sb, 0, nil); s != -1 || l != -1 { + t.Errorf("got (%d,%d), want (-1,-1)", s, l) + } +}