#!/usr/bin/perl -w

use 5.034;

# (c) Jesús Lozano Mosterín, 2023

my $N = 1000;
my (@x, @y, @z, @c);
for (1..$N){
    $x[$_]=rand();
    $y[$_]=rand();
    $z[$_]= 1000*rand();   # peso de cada punto: demanda, habitantes, etc.
    $c[$_]=0;  # al principio, no pertenecen
}   

my $cluster =0;    # counter of clusters
my ($creat, $absort, $fusion) = (0,0,0);
my ($uno,$otro);
my $iteracion = 0;
my $count=0;   # contador de puntos asimilados por el sistema clustérido???
my $lastant = 9e99;

while (1){
    my $min = 9e99;
    my $dist = 0;
    my $maxdist = 0;
    $uno = $otro = 0;
    for my $i (1..$N+$cluster){
	next unless ($c[$i] == 0 || ($i > $N && $c[$i] > 0));
	for my $j ($i+1..$N+$cluster){
	    next unless ( $c[$j] == 0 || ($j>$N && $c[$j] > 0));
	    
	    $maxdist = 0;
	    if ($i > $N && $j > $N){                 # to depreceate fusion:
		for my $k (1 .. $N){                 # distancia máxima entre todos los puntos
		    if ($c[$k] == $i){
			for my $r (1 .. $N){
			    if ($c[$r] ==  $j){
				$dist = sqrt ( ($x[$r]-$x[$k] )**2 + ($y[$r]-$y[$k])**2 );
				$maxdist = $dist if ($dist > $maxdist);
			    }
			}
		    }
		}
		$dist = $maxdist;
	    }else{
		$dist = sqrt ( ($x[$i]-$x[$j] )**2 + ($y[$i]-$y[$j])**2 );
	    }		
	    
	    if ($dist < $min){ 
		$min = $dist;
		$uno = $i;
		$otro = $j; 
	    }
	}
    }
    my $sum = 0;
    if ($uno<=$N && $otro<=$N){    ###                        CREACION
	$cluster++;
	$c[$uno] = $c[$otro] = $N+$cluster;
	$x[$N+$cluster]= ($x[$uno]*$z[$uno]+$x[$otro]*$z[$otro])/($z[$uno]+$z[$otro]);
	$y[$N+$cluster]= ($y[$uno]*$z[$uno]+$y[$otro]*$z[$otro])/($z[$uno]+$z[$otro]);
	$z[$N+$cluster] = $z[$uno] + $z[$otro];
	$c[$N+$cluster] = $N+$cluster;   # esta línea es dudosa pero el centroide NO pretenece a ningún cluster porque no es un punto real
	$count += 2;
	$creat++;
    }elsif (($uno<=$N && $otro>$N) || ($uno>$N && $otro<=$N)){         ###                         ABSORCION
	my ($main, $aux);
	if ($uno<=$N){
	    $c[$uno] = $otro;
	    $main = $otro;
	    $aux = $uno;
	}else{
	    $c[$otro]= $uno;
	    $main = $uno;
	    $aux = $otro;
	}
	$x[$main] = $y[$main]= $z[$main] = $sum = 0;
	for my $i (1..$N){
	    if ($c[$i]==$main){
		$x[$main] += $x[$i]*$z[$i];
		$y[$main] += $y[$i]*$z[$i];
		$z[$main] += $z[$i];
		$sum += $z[$i];
	    }
	}
	$x[$main]/=$sum;
	$y[$main]/=$sum;
	$count++;
	$absort++;
    }elsif ($uno>$N && $otro>$N){     ###                                FUSION
	if ($uno>$otro){ 
	    ($uno,$otro) = ($otro,$uno); # de los dos, se elige el menor y se descarta el mayor
	}
	$c[$otro] = -1;   # cluster borrado y no hace falta más; los números irán salteados
#	delete $x[$otro];
#	delete $y[$otro];
	$x[$uno]= $y[$uno]= $z[$uno] = $sum = 0;
	for my $i (1..$N){
	    if ($c[$i]==$otro){ 
		$c[$i]=$uno; 
	    }
	    if ($c[$i]==$uno){
		$x[$uno] += $x[$i]*$z[$i];
		$y[$uno] += $y[$i]*$z[$i];
		$z[$uno] += $z[$i];
		$sum += $z[$i];		
	    }
	}
	$x[$uno]/=$sum;
	$y[$uno]/=$sum;
	$fusion++;
    }else{
	die "Error del programa";
    }

    die "Error en asimilados" if $count > $N;
    
    ++$iteracion;
    say "=" x 80, "";
    say "Iteration $iteracion";
    for my $i ($N+1..$N+$cluster){
	next if ($c[$i] == -1);
	say "\t\tcluster $i, centroid x=$x[$i] y=$y[$i]   z=$z[$i]"; 
	print "Points: ";
	for my $j (1..$N){
	    if ($c[$j]==$i && $c[$j] >= $N){ print "$j "; }
	}
	say "";       # ,"-" x 70,"\n";
    }
    say "Asimilados $count puntos";
    say "Creation: $creat   Absortion: $absort   Fusion: $fusion";
    print "Outliers: ";
    for my $i (1..$N){
	print "$i " if ($c[$i]==0);
    }
    say "\n";
    my $counter = 0;
    for ($N+1..$N+$cluster){
	next if $c[$_] < 0;
	if ($c[$_] > 0) { $counter++; }
    }
    last if ($count == $N && $counter == 1);
}

exit 2;

__END__
SHA2-224(cluster3)= f08f52c25d65271d8f12dfb523d2390497a58503a50e1d94c49bd91d