github.com/biogo/biogo@v1.0.4/align/nw_affine_type.got

// Copyright ©2011-2012 The bíogo Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.

package align

import (
	"github.com/biogo/biogo/alphabet"
	"github.com/biogo/biogo/feat"

	"fmt"
	"os"
	"text/tabwriter"
)

//line nw_affine_type.got:17
func drawNWAffineTableType(rSeq, qSeq Type, index alphabet.Index, table [][3]int, a NWAffine) {
	tw := tabwriter.NewWriter(os.Stdout, 0, 0, 0, ' ', tabwriter.AlignRight|tabwriter.Debug)
	fmt.Printf("rSeq: %s\n", rSeq)
	fmt.Printf("qSeq: %s\n", qSeq)
	for l := 0; l < 3; l++ {
		fmt.Fprintf(tw, "%c\tqSeq\t", "MUL"[l])
		for _, l := range qSeq {
			fmt.Fprintf(tw, "%c\t", l)
		}
		fmt.Fprintln(tw)

		r, c := rSeq.Len()+1, qSeq.Len()+1
		fmt.Fprint(tw, "rSeq\t")
		for i := 0; i < r; i++ {
			if i != 0 {
				fmt.Fprintf(tw, "%c\t", rSeq[i-1])
			}

			for j := 0; j < c; j++ {
				p := pointerNWAffineType(rSeq, qSeq, i, j, l, table, index, a, c)
				var vi interface{}
				if vi = table[i*c+j][l]; vi == minInt {
					vi = "-Inf"
				}
				if p != "" {
					fmt.Fprintf(tw, "%s % 4v\t", p, vi)
				} else {
					fmt.Fprintf(tw, "%v\t", vi)
				}
			}
			fmt.Fprintln(tw)
		}
		fmt.Fprintln(tw)
	}
	tw.Flush()
}

func pointerNWAffineType(rSeq, qSeq Type, i, j, l int, table [][3]int, index alphabet.Index, a NWAffine, c int) string {
	switch {
	case i == 0 && j == 0:
		return ""
	case i == 0:
		if j == 1 {
			return "⬅ m"
		}
		return "⬅ l"
	case j == 0:
		if i == 1 {
			return "⬆ m"
		}
		return "⬆ u"
	}
	rVal := index[rSeq[i-1]]
	qVal := index[qSeq[j-1]]
	if rVal < 0 || qVal < 0 {
		return ""
	}
	switch p := i*c + j; table[p][l] {
	case table[p-c][up] + a.Matrix[rVal][gap]:
		return "⬆ u"
	case table[p-1][left] + a.Matrix[gap][qVal]:
		return "⬅ l"

	case table[p-c][diag] + a.GapOpen + a.Matrix[rVal][gap]:
		return "⬆ m"
	case table[p-1][diag] + a.GapOpen + a.Matrix[gap][qVal]:
		return "⬅ m"

	case table[p-c-1][up] + a.Matrix[rVal][qVal]:
		return "⬉ u"
	case table[p-c-1][left] + a.Matrix[rVal][qVal]:
		return "⬉ l"
	case table[p-c-1][diag] + a.Matrix[rVal][qVal]:
		return "⬉ m"
	default:
		return [3]string{"", "⬆ u", "⬅ l"}[l]
	}
}

func (a NWAffine) alignType(rSeq, qSeq Type, alpha alphabet.Alphabet) ([]feat.Pair, error) {
	let := len(a.Matrix)
	if let < alpha.Len() {
		return nil, ErrMatrixWrongSize{Size: let, Len: alpha.Len()}
	}
	la := make([]int, 0, let*let)
	for _, row := range a.Matrix {
		if len(row) != let {
			return nil, ErrMatrixNotSquare
		}
		la = append(la, row...)
	}

	index := alpha.LetterIndex()
	r, c := rSeq.Len()+1, qSeq.Len()+1
	table := make([][3]int, r*c)
	table[0] = [3]int{
		diag: 0,
		up:   minInt,
		left: minInt,
	}
	table[1] = [3]int{
		diag: minInt,
		up:   minInt,
		left: a.GapOpen + la[index[qSeq[0]]],
	}
	for j := range table[2:c] {
		table[j+2] = [3]int{
			diag: minInt,
			up:   minInt,
			left: table[j+1][left] + la[index[qSeq[j+1]]],
		}
	}
	table[c] = [3]int{
		diag: minInt,
		up:   a.GapOpen + la[index[rSeq[0]]*let],
		left: minInt,
	}
	for i := 2; i < r; i++ {
		table[i*c] = [3]int{
			diag: minInt,
			up:   table[(i-1)*c][up] + la[index[rSeq[i-1]]*let],
			left: minInt,
		}
	}

	for i := 1; i < r; i++ {
		for j := 1; j < c; j++ {
			var (
				rVal = index[rSeq[i-1]]
				qVal = index[qSeq[j-1]]
			)
			if rVal < 0 {
				return nil, fmt.Errorf("align: illegal letter %q at position %d in rSeq", rSeq[i-1], i-1)
			}
			if qVal < 0 {
				return nil, fmt.Errorf("align: illegal letter %q at position %d in qSeq", qSeq[j-1], j-1)
			}
			p := i*c + j

			diagScore := table[p-c-1][diag]
			upScore := table[p-c-1][up]
			leftScore := table[p-c-1][left]

			table[p][diag] = max3(diagScore, upScore, leftScore) + la[rVal*let+qVal]

			table[p][up] = max2(
				add(table[p-c][diag], a.GapOpen+la[rVal*let]),
				add(table[p-c][up], la[rVal*let]),
			)

			table[p][left] = max2(
				add(table[p-1][diag], a.GapOpen+la[qVal]),
				add(table[p-1][left], la[qVal]),
			)
		}
	}
	if debugNeedleAffine {
		drawNWAffineTableType(rSeq, qSeq, index, table, a)
	}

	var (
		aln   []feat.Pair
		layer int
	)
	score, last := 0, diag
	i, j := r-1, c-1
	t := table[i*c+j]
	best := t[0]
	for i, s := range t[1:] {
		if s > best {
			best, layer = s, i+1
		}
	}
	maxI, maxJ := i, j
	for i > 0 && j > 0 {
		var (
			rVal = index[rSeq[i-1]]
			qVal = index[qSeq[j-1]]
		)
		switch p := i*c + j; table[p][layer] {
		case table[p-c][up] + la[rVal*let]:
			if last != up && p != len(table)-1 {
				aln = append(aln, &featPair{
					a:     feature{start: i, end: maxI},
					b:     feature{start: j, end: maxJ},
					score: score,
				})
				maxI, maxJ = i, j
				score = 0
			}
			score += table[p][layer] - table[p-c][up]
			i--
			layer = up
			last = up
		case table[p-1][left] + la[qVal]:
			if last != left && p != len(table)-1 {
				aln = append(aln, &featPair{
					a:     feature{start: i, end: maxI},
					b:     feature{start: j, end: maxJ},
					score: score,
				})
				maxI, maxJ = i, j
				score = 0
			}
			score += table[p][layer] - table[p-1][left]
			j--
			layer = left
			last = left
		case table[p-c][diag] + a.GapOpen + la[rVal*let]:
			if last != up && p != len(table)-1 {
				aln = append(aln, &featPair{
					a:     feature{start: i, end: maxI},
					b:     feature{start: j, end: maxJ},
					score: score,
				})
				maxI, maxJ = i, j
				score = 0
			}
			score += table[p][layer] - table[p-c][diag]
			i--
			layer = diag
			last = up
		case table[p-1][diag] + a.GapOpen + la[qVal]:
			if last != left && p != len(table)-1 {
				aln = append(aln, &featPair{
					a:     feature{start: i, end: maxI},
					b:     feature{start: j, end: maxJ},
					score: score,
				})
				maxI, maxJ = i, j
				score = 0
			}
			score += table[p][layer] - table[p-1][diag]
			j--
			layer = diag
			last = left
		case table[p-c-1][up] + la[rVal*let+qVal]:
			if last != diag {
				aln = append(aln, &featPair{
					a:     feature{start: i, end: maxI},
					b:     feature{start: j, end: maxJ},
					score: score,
				})
				maxI, maxJ = i, j
				score = 0
			}
			score += table[p][layer] - table[p-c-1][up]
			i--
			j--
			layer = up
			last = diag
		case table[p-c-1][left] + la[rVal*let+qVal]:
			if last != diag {
				aln = append(aln, &featPair{
					a:     feature{start: i, end: maxI},
					b:     feature{start: j, end: maxJ},
					score: score,
				})
				maxI, maxJ = i, j
				score = 0
			}
			score += table[p][layer] - table[p-c-1][left]
			i--
			j--
			layer = left
			last = diag
		case table[p-c-1][diag] + la[rVal*let+qVal]:
			if last != diag {
				aln = append(aln, &featPair{
					a:     feature{start: i, end: maxI},
					b:     feature{start: j, end: maxJ},
					score: score,
				})
				maxI, maxJ = i, j
				score = 0
			}
			score += table[p][layer] - table[p-c-1][diag]
			i--
			j--
			layer = diag
			last = diag

		default:
			panic(fmt.Sprintf("align: nw affine internal error: no path at row: %d col:%d layer:%s\n", i, j, "mul"[layer:layer+1]))
		}
	}

	aln = append(aln, &featPair{
		a:     feature{start: i, end: maxI},
		b:     feature{start: j, end: maxJ},
		score: score,
	})
	if i != j {
		switch 0 {
		case i:
			last = left
		case j:
			last = up
		default:
			panic(fmt.Sprintf("align: nw affine internal error: unexpected final segment: row: %d col: %d", i, j))
		}
		aln = append(aln, &featPair{
			a:     feature{start: 0, end: i},
			b:     feature{start: 0, end: j},
			score: table[i*c+j][last],
		})
	}

	for i, j := 0, len(aln)-1; i < j; i, j = i+1, j-1 {
		aln[i], aln[j] = aln[j], aln[i]
	}

	return aln, nil
}