package Algorithm::Bertsekas;

use strict;
use warnings;

require Exporter;
our @ISA = qw(Exporter);
our @EXPORT = qw( auction );
our $VERSION = '0.01';

#Variables global to the package

my %assignement;  # assignement: a hash representing edges in the mapping, as in the Algorithm::Kuhn::Munkres.
my @output_index; # output_index: an array giving the number of the value assigned, as in the Algorithm::Munkres.
my $optimal;
my $maximize_total_benefit;
my $matrix_spaces;
my $precision;

sub auction {
   my %args = ( matrix_ref             => undef,     # reference to array: matrix N x M
                stepsize               => undef,     # epsilon is the stepsize, by default epsilon = 1 / ( min(N,M) + 1 )
                maximize_total_benefit => 1,         # 0: minimize_total_benefit ; 1: maximize_total_benefit
				verbose                => 1,         # print messages on the screen
				print_precision        => 2,         # decimal places (the number of digits following the decimal point)
                @_,                                  # argument pair list goes here
	          );
			  
   my @matrix = @{$args{matrix_ref}};                # Input: Reference to the input matrix (MxN)
   $precision = $args{print_precision};
   $maximize_total_benefit = $args{maximize_total_benefit}; 
   
   my %g;
   %g = get_edges( \@matrix );   
   %g = mkbg( \%g );
   #print Dumper \%g;   

   my $size_array1 = $#matrix + 1;
   my $size_array2 = $#{$matrix[0]} + 1;
   
   my $min_size = $size_array1 < $size_array2 ? $size_array1 : $size_array2;
   my $max_size = $size_array1 > $size_array2 ? $size_array1 : $size_array2;

   # eps is the stepsize, if eps < 1/min(n1,n2) then the algorithm is optimal.
   my $epsilon = 1 / ( $min_size + 1 );
   $epsilon = $args{stepsize} if ( defined $args{stepsize} );
   $optimal = bidding( \%g, $epsilon );
   
   my $max_matrix_value = $g{'max_matrix_value'};
   my $total_benefit = $maximize_total_benefit ? $optimal : $min_size * $max_matrix_value - $optimal ;   	
   
   my %seen1;
   my %seen2;
   for my $i ( 0 .. $#{$g{'edges'}} ) {  
   	  my $array1_index = $g{'b'} ? $g{'edges'}[$i]{'v'} : $g{'edges'}[$i]{'u'};	  
	  my $array2_index = $g{'b'} ? $g{'edges'}[$i]{'u'} : $g{'edges'}[$i]{'v'};
      my $weight = $g{'edges'}[$i]{'w'};
	    
      $output_index[$array1_index] = $array2_index;
	  $seen1{$array1_index}++;
	  $seen2{$array2_index}++;
   }
   for my $i ( 0 .. $max_size - 1 ) {
   for my $j ( 0 .. $max_size - 1 ) {
      next if ($seen1{$i} or $seen2{$j});  
      $output_index[$i] = $j;
	  $seen1{$i}++;
	  $seen2{$j}++;
	  last;
   }}
   
   for my $i ( 0 .. $#{$g{'edges'}} ) {
	  my $array1_index = $g{'b'} ? $g{'edges'}[$i]{'v'} : $g{'edges'}[$i]{'u'};	  
	  my $array2_index = $g{'b'} ? $g{'edges'}[$i]{'u'} : $g{'edges'}[$i]{'v'};
      my $weight = $g{'edges'}[$i]{'w'};
	  $assignement{ $array1_index } = $array2_index;
   }
   
   if ($args{verbose}){
      
      printf("\nNumber of left nodes: %u\n",  $size_array1 );
      printf("Number of right nodes: %u\n", $size_array2 );
      printf("Number of edges: %u\n", $size_array1 * $size_array2 );
      printf( $maximize_total_benefit ? "Objective: maximize the total benefit\n" : "Objective: minimize the total benefit\n" );
	  print "\nSolution:\n";
      printf( "Optimal assignment: sum of values = %.${precision}f for stepsize %s= %.4f \n\n", $total_benefit, defined $args{stepsize} ? '' : "= 1 / ( $min_size + 1 ) " , $epsilon );
   
      print "Maximum index size    = [";
      for my $i ( 0 .. $#output_index ) {
         printf("%${matrix_spaces}d ", $i);
      }
      print "]\n";

      print "\@output_index indexes = [";
      for my $index (@output_index) {
         printf("%${matrix_spaces}d ", $index);
      }
      print "]\n";
	  
      print "\@output_index values  = [";
	  
      for my $i ( 0 .. $max_size - 1 ){
         my $j = $output_index[$i];
         my $output = ( defined $matrix[$i] and defined $matrix[$i]->[$j] ) ? sprintf( "%${matrix_spaces}.${precision}f ", $matrix[$i]->[$j] ) : ' ' x ($matrix_spaces+1) ; # %8s 
         print $output;
      }
      print "]\n\n";
	  			   			   
      for my $i ( 0 .. $#matrix ) {
         print " [";
         for my $j ( 0 .. $#{$matrix[$i]} ) {
            printf(" %${matrix_spaces}.${precision}f", $matrix[$i]->[$j] );
            if ( $j == $output_index[$i] ){ print "**"; } else{ print "  "; } # * indica a solução válida
         }
         print "]\n";
      }	  
   }  
   
   return ( $optimal, \%assignement, \@output_index ) ;
}

sub get_edges {
   my $matrix_ref = shift;
   my @matrix = @{$matrix_ref};
   my %g;
   my @edges;
   my $max_matrix_value;   
   
   for my $i ( 0 .. $#matrix ) {
   for my $j ( 0 .. $#{$matrix[$i]} ) {
      $max_matrix_value = $matrix[$i]->[$j] if ( (not defined $max_matrix_value) || ($matrix[$i]->[$j] > $max_matrix_value) );
   }}
   
   $matrix_spaces = length( $max_matrix_value );   
   
   for my $i ( 0 .. $#matrix ) {
   for my $j ( 0 .. $#{$matrix[$i]} ) {

	  my $weight = $maximize_total_benefit ? $matrix[$i]->[$j] : $max_matrix_value - $matrix[$i]->[$j] ;

      my %hash;
      $hash{'u'} = $i;
      $hash{'v'} = $j;
      $hash{'w'} = $weight;

      push @edges, \%hash;   
      #printf ("array1[%3d] = %${matrix_spaces}.${precision}f ; array2[%3d] = %${matrix_spaces}.${precision}f ; matrix = %${matrix_spaces}.${precision}f ; get max weight = %${matrix_spaces}.${precision}f \n", $i, $array1[$i], $j, $array2[$j], $matrix[$i]->[$j], $weight );
   }}

   $g{'edges'} = \@edges; # $#edges + 1 : number of edges
   
   my $size_array1 = $#matrix + 1;
   my $size_array2 = $#{$matrix[0]} + 1;

   $g{'n1'} = $size_array1 + 1; # n1 number of left  nodes + 1
   $g{'n2'} = $size_array2 + 1; # n2 number of right nodes + 1
   $g{'max_matrix_value'} = $max_matrix_value;
   
   return %g;
}

sub mkbg { # Building a special sparse graph structure
    my $hash_ref = shift;
    my %g = %{$hash_ref};
	my @edges = @{$g{'edges'}};
	my @d;

	if ( $g{'n1'} <= $g{'n2'} ){	
	    $g{'b'} = 0;  # bool b; b=1 iff n1>n2		
		for my $i ( 0 .. $#edges ) {
			my $u = $g{'edges'}[$i]{'u'};
			$d[$u]++;
		}		
		$g{'cd'}[0]=0; # cd:  cumulative degree: (start with 0) length=dim+1
		for my $i ( 1 .. $g{'n1'} ) {		
		    $d[$i-1] = 0 if ( not defined $d[$i-1] );
			$g{'cd'}[$i] = $g{'cd'}[$i-1] + $d[$i-1];
			$d[$i-1] = 0;
		}
		for my $i ( 0 .. $#edges ) {
			my $u = $g{'edges'}[$i]{'u'};
			$g{'adj'}[ $g{'cd'}[$u] + $d[$u]   ] = $g{'edges'}[$i]{'v'}; # adj: list of neighbors
            $g{'w'  }[ $g{'cd'}[$u] + $d[$u]++ ] = $g{'edges'}[$i]{'w'}; # w:   weight of the edges
		}				
	}
    else{
		$g{'b'} = 1;  # bool b; b=1 iff n1>n2
		for my $i ( 0 .. $#edges ) {
			my $v = $g{'edges'}[$i]{'v'};
			$d[$v]++;
		}		
		$g{'cd'}[0]=0;
		for my $i ( 1 .. $g{'n2'} ) {
		    $d[$i-1] = 0 if ( not defined $d[$i-1] );
			$g{'cd'}[$i] = $g{'cd'}[$i-1] + $d[$i-1];
			$d[$i-1] = 0;
		}
		for my $i ( 0 .. $#edges ) {
			my $u = $g{'edges'}[$i]{'v'};
			$g{'adj'}[ $g{'cd'}[$u] + $d[$u]   ] = $g{'edges'}[$i]{'u'};
            $g{'w'  }[ $g{'cd'}[$u] + $d[$u]++ ] = $g{'edges'}[$i]{'w'};		
		}		
		my $temp = $g{'n1'};
		$g{'n1'} = $g{'n2'};
		$g{'n2'} = $temp;
	}
	
	return %g;
}

sub bidding { # bidding towards convergence
    my ( $hash_ref, $eps ) = @_;

    my %g = %{$hash_ref};
    my ( $u, $v, $w, $j, $k, $res );
	my ( $max, $max2, $wmax );
	
    my @p;  # p[i]=price of object i
    my @a;  # a[target]=-1 if not assigned = source if assigned
    my @aw; # aw[target]=? if not assigned = weigth of edge source-target if assigned

    for my $v ( 0 .. $g{'n2'} - 1 ) {
       $a[$v] = -1;
    }
	
	my $n_set = $g{'n1'};
	my @l_set;

	for my $u ( 0 .. $n_set - 1 ) {
       $l_set[$u] = $u;
    }
	
    $res = 0;	
	
	while ($n_set>0){
		$u   = $l_set[$n_set-1];		
		$max = 0;
		
		for my $i ( $g{'cd'}[$u] .. $g{'cd'}[$u+1] - 1 ){
			$v = $g{'adj'}[$i];
			$w = $g{'w'}[$i];
            $p[$v] = 0 if ( not defined $p[$v] );
			if ( $w - $p[$v] > $max  ){
				$max  = $w - $p[$v];
				$wmax = $w;
				$k    = $v;
			}
		}
		
		if ( $max == 0 ){
			$n_set--;
			next;
		}
				
		$max2=0;
		
		for my $i ( $g{'cd'}[$u] .. $g{'cd'}[$u+1] - 1 ){
			$v = $g{'adj'}[$i];			
			if ( $v != $k ){
				$w = $g{'w'}[$i];				
				if ( $w - $p[$v] > $max2 ){
					$max2 = $w - $p[$v];
				}			
			}
		}

		$p[$k] += $max - $max2 + $eps;

		if ( $a[$k] != -1 ){
			$l_set[$n_set-1] = $a[$k];
			$res -= $aw[$k];
		}
		else{
			$n_set--;
		}
		
		$a[$k]  = $u;
		$aw[$k] = $wmax;
		$res   += $wmax;
	}

	$j = 0;	
	for my $i ( 0 .. $g{'n2'} - 1 ) {
		if ( $a[$i] != -1 ){
			$g{'edges'}[$j]{'u'} = $a[$i];			
			$g{'edges'}[$j]{'v'} = $i;			
			$g{'edges'}[$j++]{'w'} = $aw[$i];
		}
	}

	splice @{$g{'edges'}}, $j if ( @{$g{'edges'}} > 1 );
	return $res;
}

1;  # don't forget to return a true value from the file


__END__

=head1 NAME

    Algorithm::Bertsekas - This is an efficient perl implementation of the Auction algorithm 
	for the assignement problem such as detailed in: D.P. Bertsekas, A distributed algorithm 
	for the assignment problem,Laboratory for Information and Decision Systems Working 
	Paper (M.I.T., March 1979).

    It should easily scale to hundreds of millions of nodes and/or edges if the data is not too adverserial.		

=head1 SYNOPSIS

 use Algorithm::Bertsekas qw(auction);

 my @array1 = ( 6.70, 3.08, 2.98, 0.53, 9.47, 4.51, 1.83 );
 my @array2 = ( 8.46, 0.31, 1.73, 9.33, );
 
 my $min_size = $#array1 < $#array2 ? $#array1 : $#array2;
 my $max_size = $#array1 > $#array2 ? $#array1 : $#array2;
 
 Find the nearest neighbors given a matrix represented by the weight function f (i, j) between 
 each element of the two vectors @array1 and @array2 or two lists. The weight function chosen can 
 be the modulus of the difference between two real numbers: f(i,j) = abs ($array1[i] - $array2[j]).
 
 In other words, choice of pairs such that the sum of the weights corresponding to the pairs is minimal/maximal.

 for my $i ( 0 .. $#array1 ){
   my @weight_function;		   
   for my $j ( 0 .. $#array2 ){
      my $weight = sprintf( "%.4f", abs ($array1[$i] - $array2[$j]) );			  
      push @weight_function, $weight;
   }
   push @input_matrix, \@weight_function;
 } 
 
        8.46 0.31 1.73 9.33

 6.70 [ 1.76 6.39 4.97 2.63 ]
 3.08 [ 5.38 2.77 1.35 6.25 ]
 2.98 [ 5.48 2.67 1.25 6.35 ]
 0.53 [ 7.93 0.22 1.20 8.80 ]
 9.47 [ 1.01 9.16 7.74 0.14 ]
 4.51 [ 3.95 4.20 2.78 4.82 ]
 1.83 [ 6.63 1.52 0.10 7.50 ]

 @input_matrix = (
     [ 1.76, 6.39, 4.97, 2.63 ],
     [ 5.38, 2.77, 1.35, 6.25 ],
     [ 5.48, 2.67, 1.25, 6.35 ],
     [ 7.93, 0.22, 1.20, 8.80 ],
     [ 1.01, 9.16, 7.74, 0.14 ],
     [ 3.95, 4.20, 2.78, 4.82 ],
     [ 6.63, 1.52, 0.10, 7.50 ],	 
	 );

 my ( $optimal, $assignement_ref, $output_index_ref ) = auction( matrix_ref => \@input_matrix, maximize_total_benefit => 0, verbose => 1, precision => 2 );
	 
 Number of left nodes: 7
 Number of right nodes: 4
 Number of edges: 28
 Objective: minimize the total benefit

 Solution:
 Optimal assignment: sum of values = 2.22 for stepsize = 1 / ( 4 + 1 ) = 0.2000

 Maximum index size    = [   0    1    2    3    4    5    6 ]
 @output_index indexes = [   0    4    5    1    3    6    2 ]
 @output_index values  = [1.76           0.22 0.14      0.10 ]

 [ 1.76** 6.39   4.97   2.63  ]
 [ 5.38   2.77   1.35   6.25  ]
 [ 5.48   2.67   1.25   6.35  ]
 [ 7.93   0.22** 1.20   8.80  ]
 [ 1.01   9.16   7.74   0.14**]
 [ 3.95   4.20   2.78   4.82  ]
 [ 6.63   1.52   0.10** 7.50  ]

 Pairs (in ascending order of weight function values):
   indexes (   6,   2 ), values ( 1.83, 1.73 ), matrix_value = 0.10 ; sum_matrix_value = 0.10
   indexes (   4,   3 ), values ( 9.47, 9.33 ), matrix_value = 0.14 ; sum_matrix_value = 0.24
   indexes (   3,   1 ), values ( 0.53, 0.31 ), matrix_value = 0.22 ; sum_matrix_value = 0.46
   indexes (   0,   0 ), values ( 6.70, 8.46 ), matrix_value = 1.76 ; sum_matrix_value = 2.22
   indexes (   1,   4 ), values ( 3.08,      ), matrix_value =      ; sum_matrix_value = 2.22
   indexes (   2,   5 ), values ( 2.98,      ), matrix_value =      ; sum_matrix_value = 2.22
   indexes (   5,   6 ), values ( 4.51,      ), matrix_value =      ; sum_matrix_value = 2.22

   
 foreach my $array1_index ( sort { $a <=> $b } keys %{$assignement_ref} ){     
   my $i = $array1_index;
   my $j = $assignement_ref->{$array1_index};   
   ...
 }

 Auction Algorithm, assignement hash --> $i =   0 e $value1 = 6.70; $j =   0 e $value2 = 8.46 ; difference = 1.76 ; sum = 1.76
 Auction Algorithm, assignement hash --> $i =   3 e $value1 = 0.53; $j =   1 e $value2 = 0.31 ; difference = 0.22 ; sum = 1.98
 Auction Algorithm, assignement hash --> $i =   4 e $value1 = 9.47; $j =   3 e $value2 = 9.33 ; difference = 0.14 ; sum = 2.12
 Auction Algorithm, assignement hash --> $i =   6 e $value1 = 1.83; $j =   2 e $value2 = 1.73 ; difference = 0.10 ; sum = 2.22
 
 for my $i (0 .. $max_size){
   my $j = $output_index_ref->[$i];
   ...
 }

 Auction Algorithm, output_index array --> $i =   0 e $value1 = 6.70; $j =   0 e $value2 = 8.46 ; difference = 1.76 ; sum = 1.76
 Auction Algorithm, output_index array --> $i =   1 e $value1 = 3.08; $j =   4 e $value2 =      ; difference =      ; sum = 1.76
 Auction Algorithm, output_index array --> $i =   2 e $value1 = 2.98; $j =   5 e $value2 =      ; difference =      ; sum = 1.76
 Auction Algorithm, output_index array --> $i =   3 e $value1 = 0.53; $j =   1 e $value2 = 0.31 ; difference = 0.22 ; sum = 1.98
 Auction Algorithm, output_index array --> $i =   4 e $value1 = 9.47; $j =   3 e $value2 = 9.33 ; difference = 0.14 ; sum = 2.12
 Auction Algorithm, output_index array --> $i =   5 e $value1 = 4.51; $j =   6 e $value2 =      ; difference =      ; sum = 2.12
 Auction Algorithm, output_index array --> $i =   6 e $value1 = 1.83; $j =   2 e $value2 = 1.73 ; difference = 0.10 ; sum = 2.22
 
=head1 OPTIONS
 
 matrix_ref => \@input_matrix     reference to array: matrix N x M
 stepsize => 0.01                 by default epsilon = 1 / ( min(N,M) + 1 )
 maximize_total_benefit => 0      0: minimize the total benefit ; 1: maximize the total benefit
 verbose  => 1                    print messages on the screen
 precision => 2,                  print messages on the screen: decimal places (the number of digits following the decimal point)

=head1 EXPORT

    "auction" function by default.

=head1 INPUT

    The input matrix should be in a two dimensional array (array of array) 
	and the 'auction' subroutine expects a reference to this array. 

    The $output_index_ref is the reference to the output array.

=head1 OUTPUT

    The $output_index_ref is the reference to the output_index array.
	The $assignement_ref  is the reference to the assignement hash.
	The $optimal is the total benefit which can be a minimum or maximum value.
	

=head1 SEE ALSO

    1. Dimitri P. Bertsekas, A distributed algorithm for the assignment problem, 
	   Laboratory for Information and Decision Systems Working Paper (M.I.T., March 1979).
	   
	2. Dimitri P. Bertsekas, Auction algorithms for network flow problems: 
	   A tutorial introduction - Comput Optim Applic (1992). 
	   https://doi.org/10.1007/BF00247653
	
	3. Maximilien Danisch, Efficient C implementation of the auction algorithm.
	   https://github.com/maxdan94/auction
	   This perl algorithm was adapted from the C algorithm.


=head1 AUTHOR

    Claudio Fernandes de Souza Rodrigues
	February 2018
	Sao Paulo, Brasil
	claudiofsr@yahoo.com

=head1 COPYRIGHT AND LICENSE

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.

=cut