obitools4/pkg/obisimilarity/scoring.go

package obisimilarity

import (
	"errors"
	"math"

	"git.metabarcoding.org/obitools/obitools4/obitools4/pkg/obialign"
	log "github.com/sirupsen/logrus"
)

// Reference is the length used to normalize a score (or to express a raw
// distance): the alignment length, the longest or the shortest sequence.
type Reference int

const (
	RefAlignment Reference = iota
	RefLongest
	RefShortest
)

// Params describe how a pair of sequences is scored and compared to a threshold.
type Params struct {
	// Threshold: with a normalized similarity an identity (0.95 for 95%), with
	// a normalized distance 1 - identity, with a raw similarity the length of
	// the LCS, with a raw distance the number of differences with the
	// reference length. 0 means no threshold: all the pairs are reported.
	Threshold  float64
	Normalized bool // the score is divided by the reference length
	Distance   bool // the score is a distance, not a similarity
	Reference  Reference
	// Ratio is the maximum ratio between the counts of a sequence and of the
	// one it is compared to (clustering only): a larger ratio gives no match.
	Ratio float64
}

// Scoring scores the pairs of sequences and decides whether they pass the
// threshold. It must not be modified once built.
type Scoring struct {
	threshold float64 // threshold as used for the comparison (1 - t for a normalized distance)
	normalize bool
	lcsmode   bool // true: score is a similarity, false: it is a distance
	ref       Reference
	ratio     float64
	worst     float64 // score of a pair that does not pass
	all       bool    // no threshold: every pair is reported

	// byLength gives, for the pairs whose longest sequence has a length l,
	// the values that depend on l only. It is built when the reference is not
	// the shortest sequence, otherwise those values depend on both lengths.
	byLength []lengthInfo
}

type lengthInfo struct {
	lcsMin    int32 // see lcsMin
	minCommon int32 // see minCommon
}

// NewScoring builds the scoring of the sequences `entries`, which are the
// ones to be compared (the longest of them is needed for a raw distance).
func NewScoring(p Params, entries []Entry) (*Scoring, error) {
	if p.Threshold < 0 || (p.Normalized && p.Threshold > 1.0) {
		return nil, errors.New("the threshold must not be negative, and must be at most 1 when scores are normalized")
	}

	sc := &Scoring{
		normalize: p.Normalized,
		lcsmode:   !p.Distance,
		ref:       p.Reference,
		ratio:     p.Ratio,
	}
	sc.threshold = p.Threshold
	sc.all = p.Threshold == 0
	if !sc.lcsmode && sc.normalize && p.Threshold > 0 {
		sc.threshold = 1.0 - p.Threshold
	}

	lmax := 0
	for i := range entries {
		lmax = max(lmax, entries[i].Length)
	}
	if !sc.similarityLike() {
		sc.worst = float64(lmax)
	}

	if sc.ref != RefShortest {
		sc.byLength = make([]lengthInfo, lmax+1)
		for l := range sc.byLength {
			sc.byLength[l] = lengthInfo{int32(sc.lcsMin(l, l)), int32(sc.minCommon(l, l))}
		}
	}
	return sc, nil
}

// Worst is the score of a pair that does not pass the threshold.
func (sc *Scoring) Worst() float64 { return sc.worst }

// Threshold is the threshold the scores are compared to.
func (sc *Scoring) Threshold() float64 { return sc.threshold }

// SimilarityLike tells if the comparison with the threshold is "score >= threshold".
func (sc *Scoring) SimilarityLike() bool { return sc.similarityLike() }

func (sc *Scoring) similarityLike() bool { return sc.lcsmode || sc.normalize }

// lcsMin is the smallest LCS score that a pair of lengths l1 >= l2 can have
// and still pass the threshold.
func (sc *Scoring) lcsMin(l1, l2 int) int {
	switch {
	case sc.threshold <= 0:
		return 0
	case sc.normalize:
		if sc.ref == RefShortest {
			return int(sc.threshold * float64(l2))
		}
		return int(sc.threshold * float64(l1))
	case sc.lcsmode:
		return int(sc.threshold)
	case sc.ref == RefShortest:
		return int(float64(l2) - sc.threshold)
	default:
		return int(float64(l1) - sc.threshold)
	}
}

// minCommon is the number of 4-mers that two sequences of lengths l1 >= l2
// must share to be possibly similar enough.
func (sc *Scoring) minCommon(l1, l2 int) int {
	if sc.all {
		return 0 // no threshold, no filter
	}
	reflen := l2
	if sc.ref == RefAlignment || sc.ref == RefLongest {
		reflen = l1
	}
	var lcs int
	if sc.normalize {
		lcs = int(math.Ceil(float64(reflen) * sc.threshold))
	} else {
		th := sc.threshold
		if !sc.lcsmode {
			th = float64(reflen) - th
		}
		lcs = int(th)
	}
	return minCommonWords(l1, lcs)
}

// lcsMinFor is lcsMin for sequences of lengths l1 >= l2, from the table when it exists.
func (sc *Scoring) lcsMinFor(l1, l2 int) int {
	if sc.byLength != nil {
		return int(sc.byLength[l1].lcsMin)
	}
	return sc.lcsMin(l1, l2)
}

// LengthCompatible tells if two sequences of these lengths may be similar
// enough: the shortest cannot be shorter than the minimum LCS score of the pair.
func (sc *Scoring) LengthCompatible(la, lb int) bool {
	l1, l2 := max(la, lb), min(la, lb)
	return l2 >= sc.lcsMinFor(l1, l2)
}

// PairScore scores the pair (q, c). ok is false when the pair is known not to
// reach the threshold, in which case no score is given.
//
// common, if not negative, is the sum of the minima of the 4-mer counters of
// the two sequences, already computed by a block of tables (see Index).
func (sc *Scoring) PairScore(q, c *Entry, common int, buf *obialign.BandedLCSBuffer) (score float64, ok bool) {
	// The quotient of two counts is not above 1 when the first is not above
	// the second, which spares the division in the usual case.
	if !(sc.ratio >= 1 && q.Count <= c.Count) && float64(q.Count)/float64(c.Count) > sc.ratio {
		return 0, false
	}

	l1, l2 := q.Length, c.Length
	long, short := q, c
	if l1 < l2 {
		l1, l2 = l2, l1
		long, short = c, q
	}
	var lcsMin, minCommon int
	if sc.byLength != nil {
		info := sc.byLength[l1]
		lcsMin, minCommon = int(info.lcsMin), int(info.minCommon)
	} else {
		lcsMin, minCommon = sc.lcsMin(l1, l2), sc.minCommon(l1, l2)
	}
	if l2 < lcsMin {
		return 0, false
	}
	// The 4-mer filter counts identical words: it does not hold when IUPAC
	// codes can match different letters. It lets short sequences pass.
	if l1 >= 12 && l2 >= 12 && !long.Ambiguous && !short.Ambiguous {
		if common >= 0 {
			common += min(long.Table.Overflow, short.Table.Overflow)
		} else {
			common = commonWords(long.Table, short.Table)
		}
		if common < minCommon {
			return 0, false
		}
	}

	// The length of the alignment is needed only to normalize by it, or to
	// express a distance with it; the score alone is much cheaper.
	var lcs, alilen int
	var err error
	if sc.ref == RefAlignment && (sc.normalize || !sc.lcsmode) {
		lcs, alilen, err = buf.AtLeast(q.Seq, c.Seq, lcsMin)
	} else {
		lcs, err = buf.ScoreAtLeast(q.Seq, c.Seq, lcsMin)
	}
	if err != nil {
		log.Fatalf("Cannot align two sequences: %v", err)
	}
	if lcs < 0 {
		// The best alignment has a score under lcsMin: the pair cannot pass.
		return 0, false
	}

	id := float64(lcs)
	if !sc.lcsmode && !sc.normalize {
		switch sc.ref {
		case RefAlignment:
			id = float64(alilen) - id
		case RefLongest:
			id = float64(l1) - id
		case RefShortest:
			id = float64(l2) - id
		}
	}
	if sc.normalize {
		switch sc.ref {
		case RefAlignment:
			id /= float64(alilen)
		case RefLongest:
			id /= float64(l1)
		case RefShortest:
			id /= float64(l2)
		}
	}
	return id, true
}

// Passes tells if the score s of a pair reaches the threshold.
func (sc *Scoring) Passes(s float64) bool {
	if sc.all {
		return true
	}
	if sc.similarityLike() {
		return s >= sc.threshold
	}
	return s <= sc.threshold
}

// Beats tells if the score s is better than the score best (without any
// threshold).
func (sc *Scoring) Beats(s, best float64) bool {
	if sc.similarityLike() {
		return s > best
	}
	return s < best
}

// SelfScore is the score of a sequence with itself.
func (sc *Scoring) SelfScore(length int) float64 {
	switch {
	case sc.normalize && sc.lcsmode:
		return 1.0
	case !sc.lcsmode:
		return 0.0
	default:
		return float64(length)
	}
}

// Reported is the score to report for a pair whose score is s: a normalized
// distance is computed as a similarity against the converted threshold.
func (sc *Scoring) Reported(s float64) float64 {
	if !sc.lcsmode && sc.normalize {
		return 1.0 - s
	}
	return s
}