巡回セールスマン問題をPerlで(4)

はじめに

前回に引き続き、巡回セールスマン問題の解法をPythonからPerlに移植します。
今回は、greed1関数内で既存ライブラリ(List::PriorityQueue, Graph::UnionFind)を使っていた部分を、元サイトのロジックどおりに実装し直して、draw_path関数の移植をミスっていたのを修正しました。
結果的に全サンプルデータで元サイトどおりの結果が出るようになったので、今回で終了にしたいと思います。

参考サイト

http://www.geocities.jp/m_hiroi/light/pyalgo64.html
http://www.geocities.jp/m_hiroi/light/pyalgo03.html
http://www.geocities.jp/m_hiroi/light/pyalgo61.html

greed1で使っているモジュール

extlib/lib/perl5/My/Edge.pm

package My::Edge;

sub new {
    my $class = shift;
    my ($p1, $p2, $weight) = @_;
    return bless {
        p1     => $p1,
        p2     => $p2,
        weight => $weight,
    }, $class;
}

1;

extlib/lib/perl5/My/PQueue.pm

package My::PQueue;

use strict;
use warnings;

sub new {
    my $class = shift;
    my ($cmp) = @_;
    my $self = bless {
        buff => [],
        cmp  => $cmp || sub { $_[0] - $_[1] },
    }, $class;
    for my $n (0 .. $#{$self->{buff}}) {
        $self->downheap($n);
    }
    return $self;
}

sub downheap {
    my $self = shift;
    my ($n) = @_;
    my $size = scalar(@{$self->{buff}});
    while () {
        my $c = 2 * $n + 1;
        if ($c >= $size) {
            last;
        }
        if ($c + 1 < $size) {
            if ($self->{cmp}->($self->{buff}->[$c], $self->{buff}->[$c + 1]) > 0) {
                $c += 1;
            }
        }
        if ($self->{cmp}->($self->{buff}->[$n], $self->{buff}->[$c]) <= 0) {
            last;
        }
        my $temp = $self->{buff}->[$n];
        $self->{buff}->[$n] = $self->{buff}->[$c];
        $self->{buff}->[$c] = $temp;
        $n = $c;
    }
}

sub upheap {
    my $self = shift;
    my ($n) = @_;
    while () {
        my $p = ($n - 1) / 2;
        if ($p < 0 || $self->{cmp}->($self->{buff}->[$p], $self->{buff}->[$n]) <= 0) {
            last;
        }
        my $temp = $self->{buff}->[$n];
        $self->{buff}->[$n] = $self->{buff}->[$p];
        $self->{buff}->[$p] = $temp;
        $n = $p;
    }
}

sub push {
    my $self = shift;
    my ($data) = @_;
    push @{$self->{buff}}, $data;
    $self->upheap(scalar(@{$self->{buff}}) - 1);
}

sub pop {
    my $self = shift;
    if (scalar(@{$self->{buff}}) == 0) {
        die;
    }
    my $value = $self->{buff}->[0];
    my $last = pop @{$self->{buff}};
    if (scalar(@{$self->{buff}}) > 0) {
        $self->{buff}->[0] = $last;
        $self->downheap(0);
    }
    return $value;
}

1;

extlib/lib/perl5/My/UnionFind.pm

package My::UnionFind;

use strict;
use warnings;

sub new {
    my $class = shift;
    my ($size) = @_;
    my $table = [ map { -1 } (0 .. ($size - 1)) ];
    return bless {
        table => $table,
    }, $class;
}

sub find {
    my $self = shift;
    my ($x) = @_;
    if ($self->{table}->[$x] < 0) {
        return $x;
    }
    else {
        $self->{table}->[$x] = $self->find($self->{table}->[$x]);
        return $self->{table}->[$x];
    }
}

sub union {
    my $self = shift;
    my ($x, $y) = @_;
    my $s1 = $self->find($x);
    my $s2 = $self->find($y);
    if ($s1 != $s2) {
        if ($self->{table}->[$s1] <= $self->{table}->[$s2]) {
            $self->{table}->[$s1] += $self->{table}->[$s2];
            $self->{table}->[$s2] = $s1;
        }
        else {
            $self->{table}->[$s2] += $self->{table}->[$s1];
            $self->{table}->[$s1] = $s2;
        }
        return 1;
    }
    return 0;
}

1;

実行スクリプト

#!/usr/bin/perl
use strict;
use warnings;
use lib "./extlib/lib/perl5/";

use Time::HiRes qw(gettimeofday tv_interval);
use Data::Dumper;
use My::Edge;
use My::PQueue;
use My::UnionFind;
use Tk;
use List::Util qw(max);

sub read_data {
    my $buff = [];
    while (my $line = <STDIN>) {
        my @vals = split /\s+/, $line;
        push @$buff, \@vals;
    }
    return $buff;
}

sub distance {
    my $ps = shift; # [ [x, y], .. ]
    my $max = $#{$ps};
    my $table = [ map { [0] } @$ps ];
    for my $i (0..$max) {
        for my $j (0..$max) {
            if ($i != $j) {
                my $dx = $ps->[$i]->[0] - $ps->[$j]->[0];
                my $dy = $ps->[$i]->[1] - $ps->[$j]->[1];
                $table->[$i]->[$j] = sqrt($dx ** 2 + $dy ** 2);
            }
        }
    }
    return $table;
}

sub path_length {
    my ($path, $distance_table) = @_;
    my $n = 0;
    for my $i (1..$#{$path}) {
        $n += $distance_table->[$path->[$i - 1]]->[$path->[$i]];
    }
    $n += $distance_table->[$path->[0]]->[$path->[-1]];
    return $n;
}

sub greedy0 {
    my ($size, $distance_table) = @_;
    my @path = (0 .. $size-1);
    for my $i (0 .. $size-2) {
        my $min_len = 1000000;
        my $min_pos = 0;
        for my $j ($i+1 .. $size-1) {
            my $len = $distance_table->[$path[$i]]->[$path[$j]];
            if ($len < $min_len) {
                $min_len = $len;
                $min_pos = $j;
            }
        }
        ($path[$i + 1], $path[$min_pos]) = ($path[$min_pos], $path[$i + 1]);
    }
    return \@path;
}

sub make_edge {
    my ($size, $distance_table) = @_;
    my $edges = My::PQueue->new(sub { $_[0]->{weight} - $_[1]->{weight} });
    for my $i (0 .. $size-2) {
        for my $j ($i+1 .. $size-1) {
            my $priority = $distance_table->[$i]->[$j];
            my $edge = My::Edge->new($i, $j, $priority);
            $edges->push($edge);
        }
    }
    return $edges;
}

sub edge_to_path {
    my ($edges, $size) = @_;
    my $search_edge = sub {
        my $x = shift;
        my @r;
        for my $i (0 .. $size-1) {
            my $edge = $edges->[$i];
            my ($p1, $p2) = ($edge->{p1}, $edge->{p2});
            if ($p1 == $x) {
                push @r, $p2;
            }
            elsif ($p2 == $x) {
                push @r, $p1;
            }
        }
        return @r;
    };
    my @path = map { 0 } (0 .. $size-1);
    for my $i (0 .. $size-2) {
        my ($x, $y) = $search_edge->($path[$i]);
        if ($i == 0) {
            $path[$i + 1] = $x;
            $path[-1] = $y;
        }
        elsif ($path[$i - 1] == $x) {
            $path[$i + 1] = $y;
        }
        else {
            $path[$i + 1] = $x;
        }
    }
    return \@path;
}

sub greedy1 {
    my ($size, $edges) = @_;
    
    my @edge_count = map { 0 } (0 .. $size-1);
    
    my $u = My::UnionFind->new($size);
    
    my $i = 0;
    my @select_edge;
    while ($i < $size) {
        my $edge = $edges->pop();
        my ($p1, $p2) = ($edge->{p1}, $edge->{p2});
        if (
            $edge_count[$p1] < 2 &&
            $edge_count[$p2] < 2 &&
            ( $u->find($p1) != $u->find($p2) || $i == $size -1 )
        ) {
            $u->union($p1, $p2);
            $edge_count[$p1] += 1;
            $edge_count[$p2] += 1;
            push @select_edge, $edge;
            $i += 1;
        }
    }
    return edge_to_path(\@select_edge, $size);
}

sub opt_2 {
    my ($size, $path, $distance_table) = @_;
    my $total = 0;
    while () {
        my $count = 0;
        for my $i (0 .. $size-3) {
            my $i1 = $i + 1;
            for my $j ($i+2 .. $size-1) {
                my $j1;
                if ($j == $size-1) {
                    $j1 = 0;
                }
                else {
                    $j1 = $j + 1;
                }
                if ($i != 0 || $j1 != 0) {
                    my $l1 = $distance_table->[$path->[$i]]->[$path->[$i1]];
                    my $l2 = $distance_table->[$path->[$j]]->[$path->[$j1]];
                    my $l3 = $distance_table->[$path->[$i]]->[$path->[$j]];
                    my $l4 = $distance_table->[$path->[$i1]]->[$path->[$j1]];
                    if ($l1 + $l2 > $l3 + $l4) {
                        my @new_path = @{$path}[$i1 .. $j];
                        @{$path}[$i1 .. $j] = reverse @new_path;
                        $count += 1;
                    }
                }
            }
        }
        my $path_len = path_length($path, $distance_table);
        print "len:$path_len count:$count\n";
        $total += $count;
        if ($count == 0) {
            last;
        }
    }
    return $path, $total;
}

sub opt_or {
    my ($size, $path, $distance_table) = @_;
    my $total = 0;
    while () {
        my $count = 0;
        for my $i (0 .. $size-1) {
            my $i0 = $i - 1;
            my $i1 = $i + 1;
            $i0 = $size - 1 if $i0 < 0;
            $i1 = 0         if $i1 == $size;
            for my $j (0 .. $size-1) {
                my $j1 = $j + 1;
                $j1 = 0 if $j1 == $size;
                if ($j != $i and $j1 != $i) {
                    my $l1 = $distance_table->[$path->[$i0]]->[$path->[$i]];
                    my $l2 = $distance_table->[$path->[$i]]->[$path->[$i1]];
                    my $l3 = $distance_table->[$path->[$j]]->[$path->[$j1]];
                    my $l4 = $distance_table->[$path->[$i0]]->[$path->[$i1]];
                    my $l5 = $distance_table->[$path->[$j]]->[$path->[$i]];
                    my $l6 = $distance_table->[$path->[$i]]->[$path->[$j1]];
                    if ($l1 + $l2 + $l3 > $l4 + $l5 + $l6) {
                        my $p = splice(@$path, $i, 1);
                        if ($i < $j) {
                            splice(@$path, $j, 0, $p);
                        }
                        else {
                            splice(@$path, $j1, 0, $p);
                        }
                        $count += 1;
                    }
                }
            }
        }
        my $path_len = path_length($path, $distance_table);
        print "len:$path_len count:$count\n";
        $total += $count;
        if ($count == 0) {
            last;
        }
    }
    return $path, $total;
}

sub draw {
    my ($path, $point_table) = @_;
    my $top = new MainWindow;

    my $max_x = max map { $point_table->[$_]->[0] } @$path;
    my $max_y = max map { $point_table->[$_]->[1] } @$path;
    my $canvas = $top->Canvas( -width => $max_x + 20, -height => $max_y + 20 )->pack();
    
    draw_path($canvas, $path, $point_table);
 
    MainLoop;
}

sub draw_path {
    my ($canvas, $path, $point_table) = @_;
    my @_path = map { $point_table->[$_] } @$path;
    my ($x0, $y0) = @{$_path[0]};
    for my $i (1 .. $#_path) {
        my ($x1, $y1) = @{$_path[$i]};
        $canvas->createLine($x0, $y0, $x1, $y1);
        $x0 = $x1;
        $y0 = $y1;
    }
    my ($x1, $y1) = @{$_path[0]};
    $canvas->createLine($x0, $y0, $x1, $y1);
    for my $p (@$point_table) {
        my ($x, $y) = ($p->[0], $p->[1]);
        $canvas->createOval($x-4, $y-4, $x+4, $y+4, -fill => 'green' );
    }
}

sub main {

    if ($#ARGV < 1) {
        die "invalid args 1";
    }

    my $start = [gettimeofday];

    # 初期化
    my $point_table = read_data();
    my $point_size = scalar(@$point_table);
    my $distance_table = distance($point_table);

    # 貪欲法で最適化
    my $path;
    if ($ARGV[0] eq "g0") {
        $path = greedy0($point_size, $distance_table);
    }
    elsif ($ARGV[0] eq "g1") {
        my $edges = make_edge($point_size, $distance_table);
        $path = greedy1($point_size, $edges);
    }
    else {
        die "invalid args 2";
    }
    my $path_len = path_length($path, $distance_table);

    my $end1 = [gettimeofday];
    my $spend1 = tv_interval $start, $end1;

    printf "\nlen:%d\ntime1:%s\n", $path_len, $spend1;

    # 2-opt / or-opt で最適化
    if ($ARGV[1] eq "2opt") {
        ($path, undef) = opt_2($point_size, $path, $distance_table);
    }
    elsif ($ARGV[1] eq "oropt") {
        ($path, undef) = opt_or($point_size, $path, $distance_table);
    }
    elsif ($ARGV[1] eq "2or") {
        while () {
            my $flag;
            ($path, undef) = opt_2($point_size, $path, $distance_table);
            ($path, $flag) = opt_or($point_size, $path, $distance_table);
            last if $flag == 0;
        }
    }
    elsif ($ARGV[1] eq "or2") {
        while () {
            my $flag;
            ($path, undef) = opt_or($point_size, $path, $distance_table);
            ($path, $flag) = opt_2($point_size, $path, $distance_table);
            last if $flag == 0;
        }
    }
    else {
        #die "invalid args 3";
    }
    $path_len = path_length($path, $distance_table);

    my $end2 = [gettimeofday];
    my $spend2 = tv_interval $end1, $end2;

    printf "\nlen:%d\ntime2:%s\n\n", $path_len, $spend2;

    draw($path, $point_table);
}


main();

__END__

=head1 USAGE

 perl SCRIPT_NAME [g0 | g1] [2opt | oropt | 2or | or2] < DATA_FILE

=head1 DATA FILE FORMAT

 position1 [space] position2

=head1 DATA FILE EXAMPLE

 20  20
 120 20
 220 20
 70  120
 170 120
 270 120
 20  220
 120 220