#!/usr/bin/perl

use 5.034;

# $|++;
 
my $file = 'image-vrp-16.40.png';

my $magic1 = 200; # 99:1 times the opt is miniswap, other completely new
# my $magic2 = 0.001; # halt 1 millisecond
my $magic2 = 0; # No delay
my $magic3 = 10_000; # every this iterations, a read of the last optimus is done 
my $magic4 = 20_000; # changes the cuts in the route
my $N3 = 900_000_000; # limit of iterations
my $NRUT = 8;  # hardcoded

my $N = 400;
my @ruta = 0..$N;
my @rutaopt = @ruta;
my @cortes = 0 .. $NRUT-1;
my @corteopt = @cortes;
if ($ARGV[0] =~ /r/i && -e "rutaopt.out"){
    system "cp rutaopt.out rutaopt.cop";
    system "cp corteopt.out corteopt.cop" if (-e "corteopt.out");
    open my $IN, "<", "rutaopt.out" or die $!;
    @rutaopt = <$IN>;
    close $IN;
    open my $IN2, '<', "corteopt.out" or die $!;
    @corteopt = <$IN2>;
    close $IN2;
    chomp @corteopt;
    chomp @rutaopt;
    unshift (@rutaopt, 0) if $rutaopt[0] != 0;
    @ruta = @rutaopt;
    @cortes = @corteopt;
#    goto UNO;
}

my $plot = 1;
my $debug = 1;
my $t0 = time();

my (@x,@y,@z,@code, @codeopt);

$x[0]=5;    # warehouse
$y[0]=5; 
$z[0]=0;
$code[0]="W";
srand(999_999);      # Importante para comparar con GTP-2023.pl
my $tmp;
for (1..$N){   # points of demand
    push @x, 10*rand;
    push @y, 10*rand;
    push @z, int rand(10) +1;  # La restricción de cantidad no está en vigor porque altera lo que pueden ser buenas rutas en distancia mínima
    $tmp = $code[$_-1];
    $code[$_] = ++$tmp;
    say "$_ $code[$_] $x[$_] $y[$_] demand $z[$_]" if ($debug);
}
srand();

my @d;
matrix_distance();
my $iter = 0;
my ($D, $rutas, $Z, $ZZ);

my ($Dmin, $Zmax) = (9e99, 100);

while (++$iter < $N3) {
    print "\r$iter";
    $ruta[0]=0;
    $Z = $D = 0;
    my $nant = 0;

    @ruta = miniswap();
    @ruta = jlm(@ruta) if ($iter % $magic1 == 0);
    
    if ($iter == 1 || $iter % $magic3 == 0){
	@ruta = @rutaopt;
	@cortes = @corteopt;
    }
    
    if ($iter % $magic4 == 0){
	@cortes=();
	push @cortes, 0;
	my %h;
	my $p;
	for (1 .. $NRUT-2) {
	    do {
		$p = 1+int(rand()*$N);
		$h{$p}++;
	    }until ($h{$p} == 1);
	    push @cortes , $p;
	}
	@cortes = sort  { $a <=> $b } @cortes;
    }

    my $j = 0;
    my $new;
    my $nant2 =0;
    my $c = 0;
    for my $n (@ruta){
	if ($c == $cortes[$j]) {
	    $D += $d[$n][0];
	    $j++;
	    $new = 1;
	}elsif ($new == 1){
	    $D += $d[0][$n];
	    $new = 0;
	}else{
	    $D += $d[$nant2][$n];
	}
	goto cut if ($D > $Dmin);
	$c++;
	$Z += $z[$n];
	$nant2 = $n;
    }
    $D += $d[$ruta[-1]][0];
    
    if ( $D < $Dmin ){
	if ( $Z >= $ZZ){      # Repetition is possible, but in order to not repeat swaps, fixed srand() is disconnected ASAP
	    print " $0   ", scalar(localtime()),"   rutas $NRUT   distancia $D   demanda $Z   distancia/points = ", $D/$N;
	    my $t1 = time();
	    say "   Elapsed Time Expected ", (($t1-$t0) * $N3 / $iter)/3600 - ($t1-$t0)/3600 , " h.";
	    @rutaopt = @ruta;
	    @corteopt = @cortes;                                 # CON USO
	    open my $OUT, ">", "rutaopt.out" or die $!;
	    say $OUT $_ for (@rutaopt);
	    close $OUT;
	    open my $OUT2, '>', "corteopt.out" or die $!;
	    say $OUT2 $_ for (@corteopt);
	    close $OUT2;
	    @codeopt=@code;
	    $Dmin = $D;
	    $ZZ = $Z;
	}
    }
cut: 
    select ( undef, undef, undef, $magic2 ) unless ($magic2 <= 0);   # delay for cooling
}

use POSIX;
use Imager;
use Imager::Fill;
    
my $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
    
my $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 );
my $orange = Imager::Color->new( red=>255, green=>128, blue=>0 );

    # 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 );

# say "$0 rutas $ropt,  distancia total $Dmin,  demanda total $ZZ";
my $k = 0;
my $jj = 0;
my $kount = 0;
if ($plot == 1) {
    $image -> line(color=>$white, x1=>r($x[$rutaopt[1]]), x2=>r($x[0]), y1=>r($y[$rutaopt[1]]), y2=>r($y[0]), aa=>0, endp=>1 );
    for my $j (0 .. scalar(@rutaopt)-1) {
	next if $j==0;
#	my $i = $j-1;
	$k++;
	$jj = shift @cortes if $k >= $jj;
	if ($j == $jj){ 
#	    if (defined $rutaopt[$j+1] && $rutaopt[$j+1] == 0){
#		$image -> line(color=>$red, x1 => r($xopt[$i]), x2 => r($xopt[$j]), y1=>r($yopt[$i]), y2=>r($yopt[$j]), aa=>0, endp=>1 );             # CENSURADO: SOLO RUTAS
#		$k++;
#	    }else{
		$image -> line(color=>$white, x1=>r($x[0]), x2=>r($x[$rutaopt[$j]]), y1=>r($y[0]), y2=>r($y[$rutaopt[$j]]), aa=>0, endp=>1 );                # CENSURADO: SOLO RUTAS
#	    }
	}else{
	    $image -> line(color=>$orange, x1=>r($x[$rutaopt[$j-1]]), x2=>r($x[$rutaopt[$j]]), y1=>r($y[$rutaopt[$j-1]]), y2=>r($y[$rutaopt[$j]]), aa=>0, endp=>1 );
	    $kount++;
	}	
    }
    $image -> line(color=>$white, x1=>r($x[$rutaopt[-1]]), x2=>r($x[0]), y1=>r($y[$rutaopt[-1]]), y2=>r($y[0]), aa=>0, endp=>1 );
        
    say "\n$NRUT rutas\n";
    say "$kount tramos naranja\n";
    $image->write(file => $file) or die "Cannot write: ",$image->errstr;
}
say "";
say "@rutaopt";
say "";
say "@corteopt";

exit 2;

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

sub r {                      # EN USO
    my $p = shift;
    return sprintf("%d", 100*$p);                     # como x e y van de 0 a casi 10, en vez de multiplicar por 1000, es por 100 
}

sub miniswap {
    my ($s1, $s2);
    my $n = scalar(@ruta);
    $s1 = 1+int (rand()* $n);
    $s2 = 1+int (rand()* $n);
    @ruta[$s1, $s2] = @ruta[$s2, $s1];
#    @p = $rr -> next;                                # ahora debería llamarse maxiswap
    return @ruta;                               # if $iter > $iter2 + 100_000;
}

sub cross {      # it compares s1 to s2 with s3 to/from 0                     SIN USO
    my ($s1, $s2) = @_;
    return -1 if ($s1 == $s2 || $s1 == 0 || $s2 == 0); # no pueden ser 0
    my $s3 = 1+int(rand()*$N);
                # 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)
    return 1 if ( ($x[$s2]-$x[$s1]) * ($y[$s3] - $y[$s1]) - ($y[$s2] - $y[$s1]) * ($x[$s3] - $x[$s1])  < 0);
#    return 1 if ( ($x[$s1] - $x[$s2] > $x[$s3] - $x[0] && $y[$s1] - $y[$s2] < $y[$s3] - $y[0]) || ($x[$s1] - $x[$s2] < $x[$s3] - $x[0] && $y[$s1] - $y[$s2] > $y[$s3] - $y[0]) );
    return 0;
}	

sub matrix_distance {
    for my $i (0..$N){
	for my $j (0..$N){	    
	    $d[$i][$j] = sqrt(($x[$i]-$x[$j])**2+($y[$i]-$y[$j])**2);
#	    $d[$j][$i] = $d[$i][$j];
	}
    }
}

sub jlm {                   # USADO, NO se prefiere cambios pequeños
    my ($tmp, @deck) = @_;
    my @deck2;
    my $n = scalar(@deck);
    while ($n){
	my $t = int $n*rand;
	push @deck2, splice(@deck,$t,1);
	--$n;
    }
    unshift @deck2, $tmp;    # $tmp is the first 0
    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
mié 09 ago 2023 01:05:45 CEST 
VRPvsGPT tiny data against GPT 3.5 proposal  
VRP-10 messing out
VRP-12 simplify
VRP-16 magic parameters fine tunned 

Example: 
0 0 400 83 255 261 76 86 331 15 203 256 94 89 57 127 170 213 274 178 392 8 72 249 269 323 357 145 299 58 312 55 251 43 132 104 51 366 23 308 324 320 260 222 396 82 35 237 62 56 243 379 70 
113 332 327 106 152 103 228 351 264 153 373 286 117 322 350 78 151 176 154 192 208 349 139 116 365 295 181 341 95 18 212 27 325 187 290 369 298 338 189 114 305 273 0 314 335 161 107 166 202 309 231 
234 302 136 53 311 33 336 319 25 182 168 355 310 146 37 120 272 376 198 9 16 173 267 77 64 193 150 360 81 45 88 381 29 63 67 124 297 278 32 156 383 368 121 97 119 353 91 147 245 270 102 14 291 393 
247 339 248 7 229 361 211 364 109 367 225 190 277 194 100 380 271 144 48 232 30 135 214 356 252 285 201 372 92 141 244 0 4 279 246 74 219 118 241 184 217 79 206 268 209 317 395 224 6 22 226 196 164 
2 253 155 87 54 342 276 292 172 382 262 385 216 266 263 257 167 160 134 3 105 250 31 344 34 174 122 282 128 315 197 300 52 235 159 101 20 287 99 384 5 13 137 133 47 328 333 165 220 337 177 389 199 
80 175 73 284 191 239 296 386 371 65 66 275 254 186 148 363 347 218 28 0 85 10 345 236 343 68 24 60 179 281 294 143 138 130 195 348 1 215 200 96 180 75 183 330 205 111 293 115 259 42 288 59 221 304 
158 242 17 352 207 233 108 188 140 238 318 280 50 321 399 98 316 26 157 390 230 41 289 142 49 397 19 110 84 40 307 123 61 346 149 377 185 358 313 326 112 391 204 265 39 334 38 44 11 258 354 227 126 
329 378 131 394 12 301 0 125 370 240 69 169 93 163 171 387 46 303 359 71 90 162 223 374 340 398 362 210 129 306 283 21 388 36 375 0

# ----

V.16
0   236 111 163 384 215 88 146 117 301 217 77 285 102 34 330 395 2 81 22 267 201 287 228 335 342 396 147 140 194 106 155  103 226 229 24 124 195 347 399 12  323 209 99 259 167 177 291 292 363 144
  316 90 1 181 183 374 203 121 159 302 348 184 387 224 32 158 37 339  293 382 284 63 378 233 38 274 182 247 107 235 241 196 165 208 365 87 54 222 398 249 377 131 80 359 108 114 238 391 199 21 351
  66 104 369 366 153 41 289 232 168 49 31 133 308 320 52 281 295  3 251 9 202 151 92 154 286 173  164 279 118 7 242 29 189 62 116 400 95 96 83 221 357 180 160 346 394 370 150 51 362 18 191 76 122
  174 349 187 20 245 123 19 157 341 223 379 257 327 16 67 204 105 142 10 42 300 13 218 30 263 48 256 385 172 214 373 93 101 389 115 390 305 120 192 186 178 15 231 79 333 129 225 227 125 60 64 138
  136 85 149 312 243 255 298 171 152 294 213 40 127 297 230 368 65 361 262 354 326 244 309 239 35 190 84 314 185 288 139 70 372 328 28 310 375 156 43 188 91 324 134 397 306 205 200 110 46 248 322
  343 78 392 109 246 112 216 98 290 14 75 280 376 47 388 307 162 36 296 276 265 50 371 61 268 332 258 130 33 212 272 220 176 337 311 71 166 113 45 269 303 57 100 250 219 69 355 381 240 25 271 119
  23 211 278 126 197 283 360 74 8 393 143 5 97 82 329 148 353 59 344 206 234 356 73 207 179 264 282 261 319 132 135 26 350 338 260 4 94 252 89  11 170 352 386 39 198 336 318 193 367 321 334 364
  145 266 380 299 86 44 277 72 340 237 273 345 358 175 254  68 253 270 210 53 331 17 141 161 56 317 383 275 325 55 27 169 315 128 137 58 6 304 313
  
0 34 44 73 124 134 385

# ----





0  323 85 163 312 201 267 357 261 320 200 59 319 260 4 11 95 81 255 33 147 107 34 149 253 68 182 228 335 194 40 106 155  44 211 23 195 24 124 347 179 91  236 209 18 150 362 167 138 1 363 144 137 66
  302 291 177 136 313 304 292 104 376 257 348 224 32 37 339 393  293 41 153 63 284 378 219 314 43 70 372 298 241 162 9 143 74 97 174 203 121 249 108 80 359 131 398 114 220 101 275 90 369 222 317 199 38
  185 288 146 392 168 301 31 133 308 344 52 281 295   183 251 165 151 208 286 173 349  329 142 67 10 53 189 191 289 79 272 342 83 96 334 364 175 210 242 29 181 51 88 99 370 160 64 164 16 202 196 20 184
  327 157 373 341 223 379 300 387 158 204 360 283 118 54 377 280 13 218 30 48 57 256 385 27 214 93 19 176 389 231 390 305 354 178 186 61 15 355 274 394 227 225 232 352 102 76 62 270 3 111 235 287 116 330
  2 139 328 240 152 84 45 297 368 65 192 120 262 166 309 326 115 239 250 71 233 382 22 243 171 294 213 28 258 310 375 156 126 324 345 134 306 109 180 117 110 49 248 46 322 343 78 197 188 271 25 226 303
  14 75 325 337 391 388 307 187 36 296 276 371 265 50 361 268 332 140 130 400 247 89 47 245 123 238 366 244 113 230 127 269 100 263 290 212 395 103 112 278 273 216 277 12 264 282 367 217 39 105 8 5 365
  92 154 87 285 318 148 353 135 205 206 234 73 356 399 229 98 396 69 252 94 132 26 350 82 7 122 60 6 190 35 259 346 77 386 336 198 279 338 170 384 221 215 299 380 207 145 266 86 119 381 340 72 237 246
  321 397 358 254  333 193 129 125 21 351 17 331 141 56 161 42 383 311 55 172 169 315 128 58 316 159 374
  
0 34 44 73 124 134 385