File Coverage

blib/lib/Algorithm/NeedlemanWunsch.pm

Criterion	Covered	Total	%
statement	346	360	96.1
branch	113	138	81.8
condition	11	15	73.3
subroutine	26	26	100.0
pod	5	5	100.0
total	501	544	92.1

package Algorithm::NeedlemanWunsch;

use warnings;
use strict;

use List::Util qw(max);
use Carp;

our $VERSION = '0.04';

my $from_diag = 1;
my $from_up = 2;
my $from_left = 4;
my $from_diag_idx = 0;
my $from_up_idx = 1;
my $from_left_idx = 2;

sub _curry_callback {
    my ($univ_cb, $spec_name) = @_;

    my $cb;
    if ($spec_name eq 'align') {
        $cb = sub {
	    my $arg = { align => [ @_ ] };
	    my $rv = &$univ_cb($arg);
	    croak "select_align callback returned invalid selection $rv."
	        unless $rv eq 'align';
	};
    } else {
        $cb = sub {
	    my $arg = { $spec_name => $_[0] };
	    my $rv = &$univ_cb($arg);
	    croak "select_align callback returned invalid selection $rv."
	        unless $rv eq $spec_name;
	};
    }

    return $cb;
}

sub _canonicalize_callbacks {
    my $cb;
    if (@_) {
        $cb = $_[0];
    } else {
        $cb = { };
    }

    if (exists($cb->{select_align})) {
        my @cn = qw(align shift_a shift_b);
	foreach (@cn) {
	    if (!exists($cb->{$_})) {
	        $cb->{$_} = _curry_callback($cb->{select_align}, $_);
	    }
	}
    }

    return $cb;
}

sub new {
    my $class = shift;
    my $score_sub = shift;

    my $self = { score_sub => $score_sub, local => 0 };
    if (@_) {
        $self->{gap_penalty} = $_[0];
    }

    return bless $self, $class;
}

sub local {
    my $self = shift;

    if (@_) {
        $self->{local} = $_[0];
    }

    return $self->{local};
}

sub gap_open_penalty {
    my $self = shift;

    if (@_) {
        $self->{gap_open_penalty} = $_[0];
    }

    return $self->{gap_open_penalty};
}

sub gap_extend_penalty {
    my $self = shift;

    if (@_) {
        $self->{gap_extend_penalty} = $_[0];
    }

    return $self->{gap_extend_penalty};
}

sub align {
    my $self = shift;

    my $a = shift;
    my $b = shift;

    $self->{callbacks} = _canonicalize_callbacks(@_);

    if (!exists($self->{gap_open_penalty})) {
	if (exists($self->{gap_extend_penalty})) {
	    croak "gap_open_penalty must be defined together with gap_extend_penalty";
	}

        if (!exists($self->{gap_penalty})) {
	    $self->{gap_penalty} = &{$self->{score_sub}}();
	}

	return $self->_align_basic($a, $b);
    } else {
	if (!exists($self->{gap_extend_penalty})) {
	    croak "gap_extend_penalty must be defined together with gap_open_penalty";
	}

	if ($self->{gap_open_penalty} >= $self->{gap_extend_penalty}) {
	    croak "gap_open_penalty must be smaller than gap_extend_penalty";
	}

	return $self->_align_affine($a, $b);
    }
}

sub _align_basic {
    my $self = shift;
    my $a = shift;
    my $b = shift;

    my $A = [ [ 0 ] ];
    my $D = [ [ 0 ] ];
    my $m = scalar(@$b);
    my $n = scalar(@$a);

    my $score_diag = sub {
        my ($i, $j) = @_;

	$D->[$i - 1]->[$j - 1] +
	    &{$self->{score_sub}}($a->[$j - 1], $b->[$i - 1]);
    };

    my $score_up = sub {
	my ($i, $j) = @_;

	$D->[$i - 1]->[$j] + $self->{gap_penalty};
    };

    my $score_left;
    if (!$self->{local}) {
	$score_left = sub {
	    my ($i, $j) = @_;

	    $D->[$i]->[$j - 1] + $self->{gap_penalty};
	};
    } else {
	$score_left = sub {
	    my ($i, $j) = @_;

	    ($i < $m) ?
		$D->[$i]->[$j - 1] + $self->{gap_penalty} :
		$D->[$i]->[$j - 1];
	};
    }

    # order must correspond to $from_* constants
    my @subproblems = ( $score_diag, $score_up, $score_left );

    my $j = 1;
    while ($j <= $n) {
        $A->[0]->[$j] = $from_left;
	++$j;
    }

    if (!$self->{local}) {
	$j = 1;
	while ($j <= $n) {
	    $D->[0]->[$j] = $j * $self->{gap_penalty};
	    ++$j;
	}
    } else {
        $j = 1;
	while ($j <= $n) {
	    $D->[0]->[$j] = 0;
	    ++$j;
        }
    }

    my $i = 1;
    while ($i <= $m) {
        $A->[$i]->[0] = $from_up;
	++$i;
    }

    $i = 1;
    while ($i <= $m) {
	$D->[$i]->[0] = $i * $self->{gap_penalty};
	++$i;
    }

    $i = 1;
    while ($i <= $m) {
	$j = 1;
	while ($j <= $n) {
	    my @scores = map { &$_($i, $j); } @subproblems;
	    my $d = max(@scores);

	    my $a = 0;
	    my $from = 1;
	    my $k = 0;
	    while ($k < scalar(@scores)) {
		if ($scores[$k] == $d) {
		  $a |= $from;
		}

		$from *= 2;
		++$k;
	    }

	    $A->[$i]->[$j] = $a;
	    $D->[$i]->[$j] = $d;

	    # my $x = join ', ', @scores;
	    # warn "$i, $j: $x -> ", $D->[$i]->[$j], "\n";
	    ++$j;
	}

	++$i;
    }

    $i = $m;
    $j = $n;
    while (($i > 0) || ($j > 0)) {
        my $a = $A->[$i]->[$j];
	my @alt;
	if ($a & $from_diag) {
	    die "internal error" unless ($i > 0) && ($j > 0);
	    push @alt, [ $i - 1, $j - 1 ];
	}

	if ($a & $from_up) {
	    die "internal error" unless ($i > 0);
	    push @alt, [ $i - 1, $j ];
	}

	if ($a & $from_left) {
	    die "internal error" unless ($j > 0);
	    push @alt, [ $i, $j - 1];
	}

	if (!@alt) {
	    die "internal error";
	}

	my $cur = [ $i, $j ];
	my $move;
	if (@alt == 1) {
	    $move = $self->_simple_trace_back($cur, $alt[0],
					      $self->{callbacks});
	} else {
	    $move = $self->_trace_back($cur, \@alt);
	}

	if ($move eq 'align') {
	    --$i;
	    --$j;
	} elsif ($move eq 'shift_a') {
	    --$j;
	} elsif ($move eq 'shift_b') {
	    --$i;
	} else {
	    die "internal error";
	}
    }

    return $D->[$m]->[$n];
}

sub _align_affine {
    my $self = shift;
    my $a = shift;
    my $b = shift;

    my @D = ([ [ 0 ] ], [ [ 0 ] ], [ [ 0 ] ]); # indexed by $from_*_idx
    my $m = scalar(@$b);
    my $n = scalar(@$a);

    my $score_diag = sub {
        my ($i, $j) = @_;

	my @base = map { $_->[$i - 1]->[$j - 1]; } @D;
	my $base = max(@base);
	$base + &{$self->{score_sub}}($a->[$j - 1], $b->[$i - 1]);
    };

    my $score_up = sub {
	my ($i, $j) = @_;

	my @base = map { $_->[$i - 1]->[$j]; } @D;
	$base[$from_diag_idx] += $self->{gap_open_penalty};
	$base[$from_up_idx] += $self->{gap_extend_penalty};
	$base[$from_left_idx] += $self->{gap_open_penalty};
	max(@base);
    };

    my $score_left;
    if (!$self->{local}) {
	$score_left = sub {
	    my ($i, $j) = @_;

	    my @base = map { $_->[$i]->[$j - 1]; } @D;
	    $base[$from_diag_idx] += $self->{gap_open_penalty};
	    $base[$from_up_idx] += $self->{gap_open_penalty};
	    $base[$from_left_idx] += $self->{gap_extend_penalty};

	    max(@base);
	};
    } else {
	$score_left = sub {
	    my ($i, $j) = @_;

	    my @base = map { $_->[$i]->[$j - 1]; } @D;
	    if ($i < $m) {
		$base[$from_diag_idx] += $self->{gap_open_penalty};
		$base[$from_up_idx] += $self->{gap_open_penalty};
		$base[$from_left_idx] += $self->{gap_extend_penalty};
	    }

	    max(@base);
	};
    }

    my $j;
    if (!$self->{local}) {
	$j = 1;
	while ($j <= $n) {
	    foreach (@D) {
		$_->[0]->[$j] = $self->{gap_open_penalty} +
		    ($j - 1) * $self->{gap_extend_penalty};
	    }

	    ++$j;
	}
    } else {
	$j = 1;
	while ($j <= $n) {
	    foreach (@D) {
		$_->[0]->[$j] = 0;
	    }

	    ++$j;
	}
    }

    my $i = 1;
    while ($i <= $m) {
	foreach (@D) {
	    $_->[$i]->[0] = $self->{gap_open_penalty} +
		($i - 1) * $self->{gap_extend_penalty};
	}

	++$i;
    }

    # order must correspond to $from_* constants
    my @subproblems = ( $score_diag, $score_up, $score_left );

    $i = 1;
    while ($i <= $m) {
	$j = 1;
	while ($j <= $n) {
	    my $k = 0;
	    while ($k < 3) { # scalar(@D), scalar(@subproblems)
		$D[$k]->[$i]->[$j] = &{$subproblems[$k]}($i, $j);
		++$k;
	    }

	    # my $x = join ', ', map { $_->[$i]->[$j]; } @D;
	    # warn "$i, $j: $x\n";

	    ++$j;
	}

	++$i;
    }

    # like $score_up
    my @delta_up = ( $self->{gap_open_penalty}, $self->{gap_extend_penalty},
		     $self->{gap_open_penalty} );

    # like $score_left
    my @delta_left = ( $self->{gap_open_penalty}, $self->{gap_open_penalty},
		       $self->{gap_extend_penalty} );

    my @no_delta = (0, 0, 0);

    my @score = map { $_->[$m]->[$n]; } @D;
    my $res = max(@score);

    my $arrow = 0;
    my $flag = 1;
    my $idx = 0;
    while ($idx < 3) { # scalar(@score)
	if ($score[$idx] == $res) {
	    $arrow |= $flag;
	}

	$flag *= 2;
	++$idx;
    }

    $i = $m;
    $j = $n;
    while (($i > 0) || ($j > 0)) {
	my @alt;
	if ($arrow & $from_diag) {
	    die "internal error" unless ($i > 0) && ($j > 0);
	    push @alt, [ $i - 1, $j - 1 ];
	}

	if ($arrow & $from_up) {
	    die "internal error" unless ($i > 0);
	    push @alt, [ $i - 1, $j ];
	}

	if ($arrow & $from_left) {
	    die "internal error" unless ($j > 0);
	    push @alt, [ $i, $j - 1];
	}

	if (!@alt) {
	    die "internal error";
	}

	# my $x = join ', ', map { "[ " . $_->[0] . ", " . $_->[1] . " ]"; } @alt;
	# warn "$i, $j: $x\n";

	my $cur = [ $i, $j ];
	my $move;
	if (@alt == 1) {
	    $move = $self->_simple_trace_back($cur, $alt[0],
					      $self->{callbacks});
	} else {
	    $move = $self->_trace_back($cur, \@alt);
	}

	if ($move eq 'align') {
	    --$i;
	    --$j;

	    @score = map { $_->[$i]->[$j]; } @D;
	    if ($i == 0) {
		$arrow = $from_left;
	    } elsif ($j == 0) {
	        $arrow = $from_up;
	    } else {
		my $d = max(@score);
		$arrow = 0;
		$flag = 1;
		$idx = 0;
		while ($idx < 3) { # scalar(@score)
		    if ($score[$idx] == $d) {
			$arrow |= $flag;
		    }

		    $flag *= 2;
		    ++$idx;
		}
	    }
	} elsif ($move eq 'shift_a') {
	    --$j;

	    my @base = map { $_->[$i]->[$j] } @D;
	    my $delta;
	    if ($self->{local} && ($i == $m)) {
		$delta = \@no_delta;
	    } else {
		$delta = \@delta_left;
	    }

	    $arrow = $self->_retread($score[$from_left_idx],  $i, $j,
                \@base, $delta);
	    @score = @base;
	} elsif ($move eq 'shift_b') {
	    --$i;

	    my @base = map { $_->[$i]->[$j] } @D;
	    $arrow = $self->_retread($score[$from_up_idx], $i, $j,
		\@base, \@delta_up);
	    @score = @base;
	} else {
	    die "internal error";
	}
    }

    return $res;
}

sub _retread {
    my ($self, $to_score, $i, $j, $base, $delta) = @_;

    if ($i == 0) {
	return $from_left;
    } elsif ($j == 0) {
	return $from_up;
    }

    my $a = 0;
    my $flag = 1;
    my $idx = 0;
    while ($idx < 3) {
	if ($base->[$idx] + $delta->[$idx] == $to_score) {
	    $a |= $flag;
	}

	$flag *= 2;
	++$idx;
    }

    return $a;
}

sub _trace_back {
    my ($self, $cur, $sources) = @_;

    my $arg = { };
    foreach my $next (@$sources) {
        my $m = $self->_simple_trace_back($cur, $next, { });
	if ($m eq 'align') {
	    $arg->{align} = [ $cur->[1] - 1, $cur->[0] - 1 ];
	} elsif ($m eq 'shift_a') {
	    $arg->{shift_a} = $cur->[1] - 1;
	} elsif ($m eq 'shift_b') {
	    $arg->{shift_b} = $cur->[0] - 1;
	} else {
	    die "internal error";
	}
    }

    my $move;
    my $cb = $self->{callbacks};
    if (exists($cb->{select_align})) {
        $move = &{$cb->{select_align}}($arg);
	if (!exists($arg->{$move})) {
	    die "select_align callback returned invalid selection $move.";
	}
    } else {
        my @cn = qw(align shift_a shift_b);
	foreach my $m (@cn) {
	    if (exists($arg->{$m})) {
	        $move = $m;
		last;
	    }
	}

	if (!$move) {
	    die "internal error";
	}

	if (exists($cb->{$move})) {
	    if ($move eq 'align') {
	        &{$cb->{align}}(@{$arg->{align}});
	    } else {
	        &{$cb->{$move}}($arg->{$move});
	    }
	}
    }

    return $move;
}

sub _simple_trace_back {
    my ($self, $cur, $next, $cb) = @_;

    if ($next->[0] == $cur->[0] - 1) {
        if ($next->[1] == $cur->[1] - 1) {
	    if (exists($cb->{align})) {
	        &{$cb->{align}}($next->[1], $next->[0]);
	    }

	    return 'align';
	} else {
	    if ($next->[1] != $cur->[1]) {
	        die "internal error";
	    }

	    if (exists($cb->{shift_b})) {
	        &{$cb->{shift_b}}($cur->[0] - 1);
	    }

	    return 'shift_b';
	}
    } else {
        if ($next->[0] != $cur->[0]) {
	    die "internal error";
	}

	if ($next->[1] != $cur->[1] - 1) {
	    die "internal error";
	}

	if (exists($cb->{shift_a})) {
	    &{$cb->{shift_a}}($cur->[1] - 1);
	}

	return 'shift_a';
    }
}

1;

__END__