#!/usr/bin/perl 

use 5.038.2;

my $version = 0.06;

my $header = "Constraint programming schema ad hoc v.$version ". scalar localtime() . "\n";

my $DEBUG = 0;

my $fileout = "$0_$version.out";

# Se trata de ensayar modelos lineales o aproximables a lineales donde
# se mire la factibilidad. Lo que es la lista exahustiva de soluciones 
# el lp_solve 5.5.2.14 no te la va a dar, aunque tengo que mirarlo

# No way. Hay que simular. Qué raro, en casa no toma una segunda taza.

my (@B, @I, @x, @c, $ITER);

my $caso = "";        # cadena donde se deposita el resultado parcial 

my $n =  40;           # 0j0 que igual superamos a una constante positiva de Plank
                      # Ya :-) pues va a ser que no. Hay que asegurarse de que 
                      # estamos en espacio suficiente y que varias restricciones
                      # no conduzcan a lo mismo, que, como se explica al final
                      # no es el tipo de algoritmo
                      # v.6 pruebo a reducir las variables para que haya variedad

my $T = 1_000_000;    # v.6

my $ratio = 0.15;     # parte de T (cantidad) a explotar con coordinación desigual
                      # cantidad y fondos separados 

# Lo que hay que randomizar para saber las opciones son los pares $x[i] $c[i],
# no obstante $x[i] es una porción de $T mientras que @c no está acotado, por
# lo que habría que introducir un $F del total de @c, que a su vez se ponga en
# relación con una parte de SUM $B[i] * $x[i] / $I[i]

# Al principio pensé que esto debería ser típicamente monetario y por tanto típicaamente,
# número alto, pero si hablamos de caudales o cantidades volúmetricas :-)
my $F = 999;    # es otro orden de magnitud, otra cosa.
                  # @x pueden ser unidades de cantidad, @c de dinero

# my $FIJO = 10;    # una constante
# my $fijo = 0;     # pero no tan fija para que no haya microvariaciones

my %h;          # version 5. La cosa está avanzada pero topamos con 'mala aleatoriedad' 
                # o determinismo de la formulación. Vamos a hacer como que importa.  

my $strlist;  # la de la info final

# Empezamos:

for my $i (1 .. $n){
    $I[$i] = $n / $i;  # Nótese que la inversión no es aleatoria
}
$I[0] = 1;     # to trap nefarious known error 

my $contador = 0;      # puede tener interés propio. De momento no.
my $uniq = 0;

# my $R1 = 0.5;            # Radios
# my $R2 = $n * 0.6;
# my $R3 = $n * 0.5;

my $sum;

while (1){
    $B[0] = 0;
    $x[0] = 0;
    $c[0] = 0;
#    $sum = 0;
#    $fijo = $FIJO + (rand() - 0.5);   # microvariación positiva o negativa 
    for my $i (1 .. $n){
	$B[$i] = 10_000 * rand();                   # simple de sobra
	$x[$i] = rand() * ($ratio * $T)/ $n;       # caudal
#	$x[$i] = rand();                   # hiperesfera, conviene que sea positiva
	$c[$i] = rand() * ($F / $n);                      # fondo
#	$sum += $x[$i] ** 2;
    }

    normalizarx();
    
    if ( ! bounds() ){

	$contador++;
	
#	prep_to_digest();
	
	add_to_repport();
	
    }
    
    last if ++$ITER > 1_000_000_000; 

    print "\r$ITER $contador $uniq ";
}

repport();

say "Contador = $contador soluciones";
say "Contador de soluciones únicas = $uniq";
if ($uniq == 1){
    say "Solución única: puede ser que el problema esté formulado incorrectamente";
    say "para este tipo de algoritmo de búsqueda bajo varias restricciones";
}

exit 3;





sub bounds {                   # SUBRUTINA MAS IMPORTANTE: restricciones
    my $ok = 0;
    
#    my $sum1 = 0;
#    my $sum2 = 0;
#    my $sum3 = 0;
    my $sum = 0;
    for (my $i = 1; $i <= $n; $i++){
	return 1 unless ( $x[$i] * ( $B[$i] / $I[$i] ) - $c[$i] >= 0); # &&
#	    $x[$i-1] * $B[$i-1] / $I[$i-1] - $c[$i-1] + $cplus[$i-1] - $fijo <0);   #  &&
#	    $x[$i-2] * $B[$i-2] / $I[$i-2] - $c[$i-2] + $cplus[$i-2] - $fijo <0 );
	    
#	return 1 unless $c[$i] >= 0;                                # con la gorra para este $n
#	$sum1 += $x[$i] * $x[$i] / $i ;
#	$sum2 += $x[$i] * 0.5;
#	$sum3 += $x[$i] * 0.6;
	$sum += $c[$i];   #  - $fijo;
#	$sum2 += $x[$i] * $B[$i]; 
    }
#    return 1 unless ($sum1 <= $R1);
    
#    return 1 unless $sum2 <= $R2;
#    return 2 unless $sum3 <= $R3;
      
#    my $sum = 0;                     # se quita porque ya está normalizado suma 1
#    $sum += $c[$_] for (1 .. $n);
    return 1 if $sum >= $F;      # no ratio for $F

#    return 1 if $sum2 < $sum;
    
#    $x[0] = 0;
#    for my $i (1 .. $n){     # chequeos "cualquier i"
#	return 1 unless $x[$i] > 0;  # puede que haya una posibilidad con x nulo
    
#       IMPORTANTE: PUEDE INTERESAR LA SIGUIENTE LINEA    
#	return 1 unless $x[$i] >= $x[$i-1];         # es un filtro de resultados, realmente
                                                    # 0j0: resultó ser muy difícil, o sea que venga, que tal y
                                                    # ya no es necesario por prep_to_digest()	                                               
#    }
    
    # no se me ocurre nada más todavía
    
    return $ok;
}


sub normalizarx {                                    # SUM x_i = 1 es lo que tiene más sentido 
    my $sum = 0;
    for my $i (1 .. $n){
	$sum += $x[$i];
    }
    for my $i (1 .. $n){
	$x[$i] /= $sum;
    }
}


sub prep_to_digest {                                 #  NO USADO: ordena según x[] todo y nada más
    my @idx = 0 .. $n;
    @idx = sort { $x[$a] <=> $x[$b] || $B[$a] <=> $B[$b] } @idx;
    for my $i (0 .. $n){
	my $j = $idx[$i];
	($x[$j], $x[$i]) = ($x[$i], $x[$j]);
	($B[$j], $B[$i]) = ($B[$i], $B[$j]);
	($c[$j], $c[$i]) = ($c[$i], $c[$j]);	
    }    
}

#    prep_to_digest();         puede haber problema con la posición de esto


sub add_to_repport {
    $caso = "NSOL $contador   lineas:   B  X  C\n"; 
    $caso .= "B @B";
    $caso .= "\n";
    $caso .= "X @x";
    $caso .= "\n";
    $caso .= "C @c";
    $caso .= "\n\n";

    $h{ $caso }++;            

    if ( $h{ $caso } == 1 ){
	
	$uniq++;
	$strlist .= $caso; 
	
    }
}


sub repport {        
    open my $fh, '>', $fileout or die $!;
    say $fh $header;
    say $fh $strlist;
    say $fh "Contador = $contador soluciones\n";
    say $fh "Contador de soluciones únicas = $uniq\n";
    close $fh;
}

__END__
  
Esqueleto-herramienta
  poyaque - Domingo
  porsiaca - /home/src/common/

v.5 tiene problema con la definición del problema y 
  las variables, con lo que el ejemplo a escoger para
  probar sería mejor uno aparentemente simple y 
  arbiratriamente largo:
  
Una Hiperesfera (conjunto convexo no lineal):El interior de una hiperesfera centrada en el origen con radio \(R\)
  se define por una única inecuación no lineal:\(x_{1}^{2}+x_{2}^{2}+\dots +x_{n}^{2}\le R^{2}\) 

  TACHAN! a por la versión 6