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

はじめに

前回に引き続き、巡回セールスマン問題の解法をPythonからPerlに移植します。
今回は、2-opt法とor-opt法の対応です。

参考サイト

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

「2-opt 法のプログラム」と「or-opt 法のプログラム」をPerlで

• Pythonの配列のスライスの記法がよくわからずにはまりました
• 分割統治法は元サイトでそこまで効果が出ていなかったので省略しています
• data0, 1, 2 など小さいデータだと元サイトの結果通りになりましたが、data3 は元サイトより少し距離が長くなりました
• ランダムデータ100件とかでやってみたら全然駄目でした。。

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

use Time::HiRes qw(gettimeofday tv_interval);
use Data::Dumper;
use List::PriorityQueue;
use Graph::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 = new List::PriorityQueue;
    for my $i (0 .. $size-2) {
        for my $j ($i+1 .. $size-1) {
            my $edge = sprintf "%d,%d", $i, $j;
            my $priority = $distance_table->[$i]->[$j];
            $edges->insert($edge, $priority);
        }
    }
    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) = split /,/, $edge;
            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 = Graph::UnionFind->new;
    for my $n (0 .. $size-1) {
        $u->add($n);
    }
    
    my $i = 0;
    my @select_edge;
    while ($i < $size) {
        my $edge = $edges->pop();
        my ($p1, $p2) = split /,/, $edge;
        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 ($x0, $y0) = ($point_table->[$path->[0]]->[0], $point_table->[$path->[0]]->[1]);
    for my $i (1 .. $#{$path}) {
        my $p = $point_table->[$i];
        my ($x1, $y1) = ($p->[0], $p->[1]);
        $canvas->createLine($x0, $y0, $x1, $y1);
        ($x0, $y0) = ($x1, $y1);
    }
    my $p0 = $point_table->[0];
    $canvas->createLine($x0, $y0, $p0->[0], $p0->[1]);
    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 DEPENDENCY LIBRARY INATALL

 cpanm -l ./extlib List::PriorityQueue
 cpanm -l ./extlib Graph::UnionFind

=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