#!/usr/bin/perl

use 5.038.2;

my $version = 0.09;

# v.8 -> los resultados son sólo  los resultados @x y lo demás sería datos en un mundo de datos
# v.9 -> mirar $max y $min de %h de soluciones

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

my $DEBUG = 0;

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

my $DEC_CUT = 1;   # cosas de v.8 (mirar abajo)

my $ITERATIONS = 1_000_000_000;

# 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, @x, @c, $ITER);

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

my $n =  10;           # 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
                      # como existe 0, de 0 a 11 son 12 variables, con poco interés la primera
                      # 0j0: malos resultados, bajando a 10

my $T = 10000000.00;    # v.8 -> ESTO NO SON ITERACIONES, es el caudal

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 $I;

# 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.  

# 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 = 0;
my $max = -9e99;
my $nmax = 0;

while (1){
    $B[0] = 0;
    $I = $n * $F * rand();    # arbitrario, pero no tanto $I > $F
    $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] = 100_000 * rand();                   # simple de sobra
	# the above line with decimal cut needed v.8
	$x[$i] = rand() * ($ratio * $T)/ $n;        # caudal
#	$x[$i] = rand();                   # hiperesfera, conviene que sea positiva
	$c[$i] = rand() * ($F / $n);                      # fondo
    }

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

	$contador++;
	
#	ordenarx();                # esto es muy optativo, casi etéreo   
                                   # se mantiene sólo porque se quieren soluciones repetidas	
	add_to_repport();
	
    }
    
    last if ++$ITER > $ITERATIONS; 

    print "\r$ITER  nsol $contador  uniq $uniq  reps ", $contador - $uniq;    # pantalla
}

repport();

say "\nContador = $contador soluciones alcanzadas";
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 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;
	$x[$i] = sprintf "%.02f", $x[$i] if $DEC_CUT;       # 0j0: el problema de la suma que no suma
    }
}

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++ ){
#	if ($i <= $n/5){                             # es normal perder dinero al principio
#	    return 1 unless ( ($x[$i] * $B[$i] - $I * ($n-$i) / $n - $c[$i] <= 0) ); # ||
#	    $sum -= $c[$i];   #  - $fijo;
#	    ($x[$i-1] * $B[$i] / $I[$i] - $c[$i]  >= 0 &&            # proporción de lo anterior
#	$c[$i-1] >= $c[$i]) );   #  &&     # esta restricción supone que se prefiere estabilidad en @x
	    
#	    ($x[$i] * $B[$i-2] / $I[$i-2] - $c[$i] >= 0) );      #        supuesto: que también vale [i-1] y vale [i-2] ; y su contraparte @c
	    
#	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];
#	}else{
	    return 1 unless ( ( $x[$i] * $B[$i] - $I * ($n-$i) / $n - $c[$i] > 0 ) );    # PARTE POSITIVA
	    $sum += $c[$i];
	    # version 8
	    if ($i % 3 == 0 && $i > 2){
		return 1 unless ( $x[$i] + $x[$i-1] + $x[$i-2] <= (1+$i)/($n-2) );   # aquí hay un potencial problema de excesiva restr 
#	    }elsif ($i+1 % 3 == 0){
		
#	    }elsif ($i+2 % 3 == 0){
		
	    }
#	}
#	$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 || $sum < 0);      # no ratio for $F
    
    return 1 if $sum3 != 1;                   # the mother of all restrictions
    
#    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 ordenarx {                                 #  NO USADO: ordena según x[] todo y nada más
                                               # el problema es que al haber > 1 restricción y empates... 
    my @idx = 0 .. $n;                   
                                           
    @idx = sort { $x[$a] <=> $x[$b] || $B[$a] <=> $B[$b] } @idx;
    for my $i (0 .. $n){
	my $j = shift ( @idx );          # no me convence: 1 por 2 y 2 por 3 y 3 por 1 (puede haber una especie de bucle impropio) 
	$x[$i] = $x[$j];     # no sawps allowed = ($x[$j], $x[$i]);
	$B[$i] = $B[$j];     # no swats allowed = ($B[$j], $B[$i]);  # no se imprime
	$c[$i] = $c[$j];     # swingers welcome = ($c[$j], $c[$i]);  # tampoco. Se supone que habría datos	    
    }
}

sub add_to_repport {
#    $caso = "NSOL $contador   linea de @x base 0";
    $caso = "";  # PORQUE UN CONTADOR HACE LA SOLUCION UNICA AUNQUE NO LO SEA
#    $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++;	
    }elsif ( $h{ $caso } > $max ) {
	$nmax = 1;
	$max = $h{ $caso } ;
	# esto es lo que se va a reportar, la pequeña cantidad repetida, no lo grande 
    }elsif ( $h{ $caso } == $max ) {
	$nmax++;
    }
}


sub repport {        
    open my $fh, '>', $fileout or die $!;
    say $fh $header;
    say $fh "En función de un número de iteraciones $ITERATIONS\n";
    say $fh "SOLUCIONES de máxima ocurrencia";
    say $fh "-------------------------------";
    my $k = 0;
    for my $parraf (keys %h){
	if ( $h{ $parraf } == $max ){
	    $k++;
	    say $fh "\nN_printed_SOL = $k";
	    say $fh "Contador de ocurrencias = $h{ $parraf }";
	    say $fh $parraf;
	}
    }
    if ($max > 2){
	say $fh "\nOtras soluciones frecuentes";
	say $fh "---------------------------";
	say $fh "Hay $nmax soluciones con $max ocurrencias, pero";
	say $fh "hay ", $contador - $uniq - $nmax, " que son de";
	say $fh "aparición múltiple y se listan por interés.-";
	for my $parr (keys %h){
	    if ( $h{ $parr } < $max && $h{$parr} > 1){
		say $fh "Contador de ocurrencias = $h{ $parr }";
		say $fh $parr;
	    }
	}
    }
    say $fh "Contador = $contador soluciones\n";
    say $fh "Contador de soluciones únicas = $uniq\n";
    say $fh "Número de soluciones repetidas = ", $contador - $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

La v.7 es mala, no consigue resultados repetidos
  
La v.8 lo consigue simplificando. 

v.9 No quiero laureles, pero 0.01**10 = 1e-20
  Supera a un epsilon normal -> wikipedia