#!/usr/bin/perl

use 5.040;

# RDEA stands for 'randomized DEA', similar to 'random forests' applied 
# to DMUs.                 (c) Jesús Lozano Mosterín, 2024.

my $version = 0.003;

# bug hunting vie 27 dic 2024 22:52:18 CET

# STEP 1: read the data file

my @datafile;
if (@ARGV){
    open my $in, '<', $ARGV[0] or die $!;
    @datafile = <$in>;
    close $in;
}else{
    say "Randomized DEA, RDEA";
    say "";
    say "Usage: [perl] $0 <datafile>";
    say "";
    say "- where in the first row are number of inputs <space> number of outputs";
    say "- rest of the rows contains data for each DMU, values of inputs and outputs";
    say "- the separator between data is the 'space'";
    say "- comments in a row are avilable if the row begins with #";
    say "";
    say "version $version";
    exit 4;
}

my $header = 0;
my (@input, @output);
my ($nin, $nout, $rest, $ndmu);

for (@datafile){
    next if /^#/;
    chomp;
    if ($header == 0){
	($nin, $nout, $rest) = split /\s+/, $_;
	sured( $nin );
	sured( $nout );
	$header = 1;
	next;
    }else{
	$ndmu++;
    }	
    my @list = split /\s+/, $_;
    for my $i (0 .. $nin-1){
	$input[$ndmu][$i] = shift @list;
	sured( $input[$ndmu][$i] );
    }
    for my $i (0 .. $nout-1){
	$output[$ndmu][$i] = shift @list;
	sured( $output[$ndmu][$i] );
    }
}

# STEP 2: generate a set of weigths and shuffle that a number of times

my $nloop = 30*60;
my $vueltas = (factorial_bound($nin+$nout) < 100_000) ? factorial_bound($nin+$nout) : 100_000 ;

# STEP 3: Las restricciones normales (DEA - Charnes, Cooper, Rhodes) son:
# - suma denominador == 1  --> Con RDEA, se podría normalizar ambos, pero creo que con solo cumplir lo siguiente vale
# Es decir, el denominador=1 lo puso la transformación de Charnes para convertir program fraccionario en LP. No es el caso.
# - numerador <= denominador  --> no problem
# - testar para todos los dmus  --> no problem

# STEP 4: probar un número bajo de permutaciones de un set 100_000
# - cambiar de set 30*60 veces

# STEP 5:
# - sub evaluate { se debe almacenar las veces que una dmu obtiene el premio MAX, pero el rdo. sería un ranking }
# - sub evaluateplus { si se pone un threshold normal, como 0.001, o un redondeo, se puede poner a la dmu en una 
# categoría (lista) de eficiencia. Lo contrario no existiría. }

# STEP 6: 
# - Mirar cuáles son eficientes según evaluateplus
# - Listar en ranking y porcentual ordenadas todas las DMUs

my @RDEA;
my @effs;

for (1 .. $nloop) {
    my @set;
    for (1 .. $nin+$nout){
	# problema básico: rand() devuelve en intervalo [0,1) y lo que se
	# necesita es (0,1] ¡justo lo contrario! mersenne-twister ofrece opciones
	# sobre el mismo generador... Otra opción podría ser tomar (1-rand)
	# Vale, aprobada la moción.
    
	push @set, 1-rand;      # sprintf ("%.03f", 1-rand);  No, to not have 0.000
    }
    
    
    for (1 .. $vueltas){
	
	my @weigths = jlm_shuffle(@set);
    
	# AQUI VA LO IMPORTANTE
	
	my $RDEA_MAX = 0;
	my $candidate;
	my $tmp;
	
	for my $d (1 .. $ndmu){
	    
	    my ($numerador, $denominador);
	    for my $o ( 0 .. $nout-1 ){
		$numerador += $weigths[$o] * $output[$d][$o];
	    }
	    for my $i ($nout .. $nin + $nout -1){
		$denominador += $weigths[$i] * $input[$d][$i-$nout];
	    }
	    if ($numerador > $denominador){
		undef $candidate;
		last;
	    }elsif ($denominador > 0) {
		my $tmp = ($numerador / $denominador) / ($nloop * $vueltas);
		if ($tmp > $RDEA_MAX){
		    $RDEA_MAX = $tmp;
		    $candidate = $d;
		}
	    }  	
	}
	
	if ($candidate){
	    $RDEA[$candidate] += $RDEA_MAX;
	    push @effs, $candidate if defined $candidate;
	}else{
	    # huh ?
	}
    }
}
# FINAL REPORT
my @counter;
for my $i (@effs){
    if ($i){
	$counter[$i]++;
	$counter[0]++;
    }
}
for my $j (1 .. $ndmu) {
    unless (defined $counter[$j]) {
	$counter[$j] = 0;
    }
    unless ( defined $RDEA[$j] ) {
	$RDEA[$j] = 0;
    }
}

for my $i (sort { $counter[$b] <=> $counter[$a] || $RDEA[$b] <=> $RDEA[$a] } 1..$ndmu){
    if ($i){
	say "DMU $i)    counter $counter[$i]    corrected percent ", 100*($counter[$i]/$counter[0] - 1/$ndmu), "%    index = $RDEA[$i]";
    }
}


exit 3;


sub sured {
    my $n = shift;
    return if (0+$n eq $n);
    die "Not numeric -> $n";
}

sub factorial_bound {
    my $n = shift;
    return 100_111 if $n >= 100_000;
    return 1 if $n <= 1;
    return $n * factorial_bound($n-1);
}
   
sub jlm_shuffle {
    my @deck = @_;
    my @deck2;
    my $n = scalar(@deck);
    while ($n){
	my $t = int $n*rand;
	push @deck2, splice(@deck,$t,1);
	--$n;
    }
    return @deck2;
}