#!/usr/bin/perl -w

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

# srand(999);

my $plot = 1;
my $debug = 1;
my $N = 400;
my $N3 = 100_000 + $N**3; # limit of iterations
my $subN = $N/3;      # partición de N de enrutar
my $index = 9e99;

my (@x,@y,@code, @ruta, @rutaopt, @xopt, @yopt, @codeopt, %track);
my ($image, $white);
my $file = 'image-vrp-9.0.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 + int $subN * rand();
    for my $n ( jlm( @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]= ++$rutas;
#	    $rutas++;
	    recursive($n);
	    goto cut if $D > $Dmin;
	}else{
	    $D += 2*$d[0][$n];
#	    $in[$n]=0;     No sé porqué está comentado esto. Puede ser que no puedan estar en el mismo $in 2 rutas. 0 ya es por defecto.
	}
    }
#    if ($fault <= $faultant) {    # && $D <= $Dmin
    $fault = scalar(grep ($_==0,@in));
    if ($D * $fault <= $index) {
	$index = $D * $fault;
	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 > $N3);  # 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 );
    my $red    = Imager::Color->new( red=>255, green=>0, blue=>0 );
    my $green  = Imager::Color->new( red=>0, green=>255, blue=>0 );
    my $blue   = Imager::Color->new( red=>0, green=>0, 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, $c);
$c = 0;
$j = shift @rutaopt;
while (@rutaopt){
    $i = $j;
    $j = shift @rutaopt;
    if ($debug){
	print "$track{$code[$i]}, ";
	if ($track{$code[$i]}==0){ $c++; }
    }
    print "\n" if ($debug && $track{$code[$i]}==0 && $c==2);
    $c = 0;
    if ($plot == 1) { 
	$image->line(color=>$green, 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;
}

exit 1;

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

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]] > 0){                 # $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]]=$rutas;
	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)
		  
	    }
	}
    }
}

sub jlm {
    my @deck = @_;
    my @deck2;
    my $n = scalar(@deck);
    while ($n){
	my $t = int $n*rand;
	push @deck2, splice(@deck,$t,1);
	--$n;
    }
    return @deck2;
}


__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: 
VRP-8 revigorizing