#!/usr/bin/perl -w

use strict;
use List::Util qw(shuffle);

srand(999);

my $plot = 1;
my $debug = 1;
my $N = 100;

my (@x,@y,@code, @ruta, @rutaopt, @xopt, @yopt, @codeopt, %track);
my ($image, $white);
my $file = 'image-vrp-d3.png';


$x[0]=0.5;    # warehouse
$y[0]=0.5; 
$code[0]="W";
$track{"W"} = 0;

my $tmp;
for (1..$N){   # points of demand
    push @x, rand;
    push @y, rand;
    $tmp = $code[$_-1];
    $code[$_] = ++$tmp;
    $track{$code[$_]} = $_;
    print "$_ $code[$_] $x[$_] $y[$_]\n" if ($debug);
}

my @d;
my @partner;
matrix_distance();

my $Dmin=9e99;
my $iter=0;
my ($D, $rutas, $ropt);
my @p = 1..$N;
my @in;
do {
    @in = (0) x ($N+1);
    @ruta=();
    $D=$rutas=0;
#    for (my $n = $N ; $n > 0 ; $n--){
# for $n (1..$N){ # es peor, el orden importa
    for my $n ( shuffle( @p ) ){
#    towards(random p)
	next if $in[$n]; # scalar(grep($n==$_, @ruta));
	$D += $d[0][$n];
	goto cut if $D > $Dmin;
	push @ruta, 0;
	push @ruta, $n;
	$in[$n]=1;
	$rutas++;
	recursive($n);    
    }
    if ($D <= $Dmin) {
	print "$0 $iter: $rutas rutas; distancia $D\n";
	$ropt=$rutas;
	@xopt=@x;
	@yopt=@y;
	@rutaopt = @ruta;
	@codeopt=@code;
	$Dmin=$D;
    }
cut:    
      swap();      # this line could be commented out to gain time at a minor cost
    select ( undef, undef, undef, 0.01 );   # delay for cooling
}until (++$iter > 1_000_000);  # enough is enough ($N**3)

if ($plot == 1) { 
    use POSIX;
    use Imager;
    use Imager::Fill;
    
#    $file = 'image-vrp-d3.png';
    $image=Imager->new(xsize=>1000,ysize=>1000,channels=>3);
    
#    $img_x_max = $image->getwidth;
#    $img_y_max = $image->getheight;
    
#    $x0 = $img_x_max / 2;  # X = 0
#    $y0 = $img_y_max / 2;  # Y = 0
    
    $white  = Imager::Color->new( red=>255, green=>255, blue=>255 );
    
    # Draw the coordinate axises and tick marks
#    $image->line(color=>$white, x1=>$x0, x2=>$x0, y1=>0, y2=>$img_y_max, aa=>0, endp=>1 );
    
#    $image->line(color=>$white, x1=>0, x2=>$img_x_max, y1=>$y0, y2=>$y0, aa=>0, endp=>1 );

} 

print "$0 rutas $ropt, distancia total $Dmin\n" if ($debug);
my ($i,$j);
$j = shift @rutaopt;
while (@rutaopt){
    $i = $j;
    $j = shift @rutaopt;
    print "$track{$code[$i]}, " if ($debug);
    print "\n" if ($debug && $track{$code[$i]}==0);
    if ($plot == 1) { 
	$image->line(color=>$white, x1=>r($xopt[$i]), x2=>r($xopt[$j]),
	    y1=>r($yopt[$i]), y2=>r($yopt[$j]), aa=>0, endp=>1 );
    }
    
}
if ($plot == 1){
    $image->write(file => $file) or die "Cannot write: ",$image->errstr;
}

# --------------------------------

sub r {
    my $p = shift;
    return sprintf("%d", 1000*$p);
}

sub recursive {
    my $z = shift;
    return if $z==0;
#    return if ($in[$z]==1); # scalar(grep($z==$_,@ruta)) && $z!=$ruta[-1]);
    if ($in[$partner[$z]] == 1){                 # $partner[$z]==$_,@ruta)) {
	push @ruta, 0;
	$D+= $d[$z][0];
	return;
    }else{
	$D += $d[$z][$partner[$z]];
	push @ruta , $partner[$z];
	$in[$partner[$z]]=1;
	recursive($partner[$z]);
    }
}

sub swap {
    my $t1 = 1 + int ($N*rand);
    my $t2 = 1 + int ($N*rand);
#    @d[$t1,$t2] = @d[$t2,$t1];
#    my $tmp;
#    for my $i (0..$N){
#	$tmp=$d[$t1][$i];
#	$d[$t1][$i] = $d[$t2][$i];
#	$d[$t2][$i] = $tmp;
#    }
#    ($partner[$t1],$partner[$t2]) = ($partner[$t2],$partner[$t1]);
    ($x[$t1],$x[$t2]) = ($x[$t2],$x[$t1]);
    ($y[$t1],$y[$t2]) = ($y[$t2],$y[$t1]);
    ($code[$t1],$code[$t2]) = ($code[$t2],$code[$t1]);
    matrix_distance();
}

sub matrix_distance {
    for my $i (0..$N){
	for my $j (0..$i){
	    $d[$i][$j] = sqrt(($x[$i]-$x[$j])**2+($y[$i]-$y[$j])**2);
	    $d[$j][$i] = $d[$i][$j];
	}
    }
    my @mini;
    for my $i (1..$N){
	$mini[$i] = 9e99;
	for my $j (1..$N){
	    next if $i==$j;
	    if ($d[$i][$j] < $mini[$i]){
		$mini[$i] = $d[$i][$j];
		$partner[$i] = $j;
	    }
	}
    }
}
__END__
  
vrp-d solo dist
vrp-d3 only distance