#!/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-d5.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 $NN = $N+1;
my @partner = (0) x $NN;
matrix_distance();

my $Dmin=9e99;
my $iter=0;
my ($D, $rutas, $ropt);
my @p = 1..$N;
my @in;
my ($fault, $faultant, $count);
$faultant=9e99;
do {
    @in = (0) x $NN;
    @ruta=();
    $D=$rutas=$count=0;
#    for (my $n = $N ; $n > 0 ; $n--){
# for $n (1..$N){ # es peor, el orden importa
    $in[0] = 1 + $N * rand() / 2;
    for my $n ( shuffle( @p ) ){
#    towards(random p)
	if (++$count <= $in[0]) { 
	    next if $in[$n]; # scalar(grep($n==$_, @ruta));
#	    if ($n>100 ) { $n=1; }
	    $D += $d[0][$n];
	    goto cut if $D > $Dmin;
	    push @ruta, 0;
	    push @ruta, $n;
	    $in[$n]=1;
	    $rutas++;
	    recursive($n);
	    goto cut if $D > $Dmin;
	}else{
	    $D += 2*$d[0][$n];
#	    $in[$n]=1;
	}
    }
    $fault = scalar(grep ($_==0,@in));
#    if ($fault <= $faultant) {    # && $D <= $Dmin
	if ($D <= $Dmin) {
	    print "$0 $iter: $rutas rutas; distancia $D faltas $fault\n";
	    $ropt=$rutas;
	    @xopt=@x;
	    @yopt=@y;
	    @rutaopt = @ruta;
	    @codeopt=@code;
	    $Dmin=$D;
	    $faultant=$fault;
	}
#    }
cut:    
#      swap();      # this line could be commented out to gain time at a minor cost
#     select ( undef, undef, undef, 0.001 );   # delay for cooling
}until (++$iter > 2_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){
    for my $z (1..$N){
	if ($in[$z]==0){
	    $image -> line(color=>$white, x1=>r($x[0]), x2=>r($xopt[$z]),
		y1=>r($y[0]), y2=>r($yopt[$z]), aa=>0, endp=>1 );
	}
    }
}
print "\n" if ($debug);
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)) {
	# no $in[$z]
#	push @ruta, $z;
	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]){
#		if ( ($x[$i]-$x[0])*($y[$j]-$y[0])-($y[$i]-$y[0])*($x[$j]-$x[0]) >= 0 ){
		if (defined $partner[$j]){
		    if ($partner[$j] != $i){
			$mini[$i] = $d[$i][$j];
			$partner[$i] = $j;
		    }
		}
#		}
                # Graham scan
                # Three points are a counter-clockwise turn if ccw > 0, clockwise if
		# ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
		# gives twice the signed  area of the triangle formed by p1, p2 and p3.
		# function ccw(p1, p2, p3):
		#      return (p2.x - p1.x)*(p3.y - p1.y) - (p2.y - p1.y)*(p3.x - p1.x)
		  
	    }
	}
    }
}
__END__
  
vrp-d solo dist
vrp-d3 only distance
vrp-d5 changed algo; first three routes, adding until cost not reduces; consider local search
 Problema: