Algorithm-Bertsekas
view release on metacpan or search on metacpan
lib/Algorithm/Bertsekas.pm view on Meta::CPAN
} else {
$args{verbose} = 7;
warn "\n *** Could not open <$target> for writing: $! *** \n";
}
}
$maximize_total_benefit = $args{maximize_total_benefit};
my $optimal_benefit = 0;
my %assignment_hash; # assignment: 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 @matrix;
my @matrix_index;
foreach ( @matrix_input ){ # copy the orginal matrix N x M
push @matrix, [ @$_ ];
}
if ( $max_size <= 1 ){ # matrix_input 1 x 1
$assignment_hash{0} = 0;
$output_index[0] = 0;
$matrix_index[0] = 0;
$optimal_benefit = $matrix_input[0]->[0];
}
$need_transpose = 1 if ( $array1_size > $array2_size ); # will always be chosen N <= M
if ( $need_transpose ){
my $transposed = transpose(\@matrix);
@matrix = @$transposed;
}
get_matrix_info( \@matrix, $args{verbose} );
delete_multiple_columns( \@matrix, $args{verbose} ) if ( $min_size >= 2 and $min_size != $max_size );
# epsilon is the stepsize and auction algorithm terminates with a feasible assignment if the problem data are integer and epsilon < 1/min(N,M).
# There is a trade-off between runtime and the chosen stepsize. Using the largest possible increment accelerates the algorithm.
$inicial_price = $args{inicial_price};
$price_object{$_} = $inicial_price for ( 0 .. $max_size - 1 );
#$price_object{$_} = sprintf( "%.0f", rand($max_matrix_value) ) for ( 0 .. $max_size - 1 ); # random values for initial prices
# inicial epsilon value
my $epsilon = int (($max_matrix_value/2) * exp ( -1 * $max_size/$min_size )); # exp (1) = e = 2.71828182845905
$epsilon = $args{inicial_stepsize} if ( defined $args{inicial_stepsize} );
$epsilon = 1/(1+$min_size) if ($epsilon < 1/$min_size);
# (inicial epsilon value) / factor ** k = 1/$min_size --> epsilon * $min_size = factor ** k
# log(epsilon * $min_size) = k * log(factor) --> log(factor) = (log(epsilon * $min_size)) / k --> factor = exp ( (log(epsilon * $min_size)) / k )
#$max_epsilon_scaling = 20;
#my $factor = exp ( (log( $epsilon * (1+$min_size) )) / ($max_epsilon_scaling-1) ); print "\n \$max_epsilon_scaling = $max_epsilon_scaling ; \$factor = $factor \n";
my $factor = 4; # $factor > 1
$max_epsilon_scaling = 2 + int( log( (1+$min_size) * $epsilon )/log($factor) ); # print "\n \$max_epsilon_scaling = $max_epsilon_scaling \n";
my $feasible_assignment_condition = 0;
# The preceding observations suggest the idea of epsilon-scaling, which consists of applying the algorithm several times,
# starting with a large value of epsilon and successively reducing epsilon until it is less than some critical value.
while( $epsilon >= 1/(1+$min_size) and $max_size >= 2 ){
$epsilon_scaling++;
$iter_count_local = 0;
%assignned_object = ();
%assignned_person = ();
%seen_person = ();
$seen_assignned_objects{$_} = 0 for ( 0 .. $max_size - 1 );
while ( (scalar keys %assignned_person) < $max_size ){ # while there is at least one element not assigned.
$iter_count_global++;
$iter_count_local++;
auctionRound( \@matrix, $epsilon, $args{verbose} );
if ( $args{verbose} >= 10 ){
for my $i ( -1 .. $#matrix ) {
if ($i >= 0){ printf $output " %2s [", $i; } else{ printf $output "object"; }
for my $j ( 0 .. $#{$matrix[$i]} ) {
if ($i >= 0){ printf $output (" %${matrix_spaces}.${decimals}f", $matrix[$i]->[$j]); } else{ printf $output (" %${matrix_spaces}.0f", $j); }
if ( defined $assignned_person{$i} and $j == $assignned_person{$i} ){ print $output "**"; } else{ print $output " "; }
}
if ($i >= 0){ print $output "]\n"; } else{ print $output "\n\n"; }
}
}
}
$epsilon = $epsilon / $factor ; # (1/2): smooth convergence
if ( not $feasible_assignment_condition and $epsilon < 1/$min_size ){
$epsilon = 1/(1+$min_size);
$feasible_assignment_condition = 1;
}
}
$epsilon = $factor * $epsilon; # correcting information for printing
my %seeN;
my %seeM;
foreach my $person ( sort { $a <=> $b } keys %assignned_person ){
my $object = $assignned_person{$person};
$matrix_index[$person] = $object;
#print " \$need_transpose = $need_transpose ; \$matrix_index[$person] = $object ; \$index_i = $person ; \$index_j = $object --> $index_correlation{$object} ;";
my $index_i = $need_transpose ? $index_correlation{$object} // $object : $person;
my $index_j = $need_transpose ? $person : $index_correlation{$object} // $object;
$output_index[$index_i] = $index_j;
$seeN{$index_i}++;
$seeM{$index_j}++;
#print " \$output_index[$index_i] = $index_j \n";
next unless ( defined $matrix_input[$index_i] and defined $matrix_input[$index_i]->[$index_j] );
$assignment_hash{ $index_i } = $index_j;
$optimal_benefit += $matrix_input[$index_i]->[$index_j];
lib/Algorithm/Bertsekas.pm view on Meta::CPAN
$price_object{$object} += $bid;
if ( $verbose >= 8 ){
printf $output " --> Assigning to personI = %3s the objectJ = %3s with highestBidForJ = %20.5f and update the price vector ; \$assignned_person{%3s} = %3s ; \$price_object{%3s} = %.5f \n", $person, $object, $bid, $person, $assignned_person{$person...
}
}
if ( $verbose >= 9 )
{
$number_of_assignned_object = scalar keys %assignned_object;
print $output "\n Final: Matrix Size N x M: $min_size x $max_size ; epsilon_scaling = $epsilon_scaling ; Number of Global Iterations = $iter_count_global ; Number of Local Iterations = $iter_count_local ; epsilon = $epsilon ; \$number_of_assignned_...
foreach my $person ( sort { $a <=> $b } keys %assignned_person ){
my $object = $assignned_person{$person};
printf $output " \$assignned_person{%3s} --> object %3s --> \$price_object{%3s} = $price_object{$object} \n", $person, $object, $object;
}
foreach my $object ( sort { $a <=> $b } keys %price_object ){
printf $output " \$price_object{%3s} = $price_object{$object} \n", $object;
}
print $output "\n";
}
}
1; # don't forget to return a true value from the file
__END__
=head1 NAME
Algorithm::Bertsekas - auction algorithm for the assignment problem.
This is a perl implementation for the auction algorithm for the asymmetric (N<=M) assignment problem.
=head1 DESCRIPTION
The assignment problem in the general form can be stated as follows:
"Given N jobs (or persons), M tasks (or objects) and the effectiveness of each job for each task,
the problem is to assign each job to one and only one task in such a way that the measure of
effectiveness is optimised (Maximised or Minimised)."
"Each assignment problem has associated with a table or matrix. Generally, the rows contain the
jobs (or persons) we wish to assign, and the columns comprise the tasks (or objects) we want them
assigned to. The numbers in the table are the costs associated with each particular assignment."
In Auction Algorithm (AA) the N persons iteratively submit the bids to M objects.
The AA take cost Matrix N×M = [aij] as an input and produce assignment as an output.
In the AA persons iteratively submit the bids to the objects which are then reassigned
to the bidders which offer them the best bid.
Another application is to find the (nearest/more distant) neighbors.
The distance between neighbors can be represented by a matrix or a weight function, for example:
1: f(i,j) = abs ($array1[i] - $array2[j])
2: f(i,j) = ($array1[i] - $array2[j]) ** 2
=head1 SYNOPSIS
### --- simple and direct application --- ###
### --- start --- ###
#!/usr/bin/perl
use strict;
use warnings FATAL => 'all';
use diagnostics;
use Algorithm::Bertsekas qw(auction); # To install this modulus: 'cpan Algorithm::Bertsekas' or 'ppm install Algorithm-Bertsekas'
my @array1; my @array2;
my @input_matrix;
my $N = 10;
my $M = 10;
my $range = 1000;
for my $i (1..$N) {
push @array1, sprintf( "%.0f", rand($range) );
}
for my $i (1..$M) {
push @array2, sprintf( "%.0f", rand($range) );
}
print "\n \@array1 = ( ";
for my $value (@array1) { printf("%4.0f ", $value); }
print ")\n";
print " \@array2 = ( ";
for my $value (@array2) { printf("%4.0f ", $value); }
print ")\n";
for my $i ( 0 .. $#array1 ){
my @weight_function;
for my $j ( 0 .. $#array2 ){
#my $weight = sprintf( "%.0f", rand($range) );
my $weight = abs ($array1[$i] - $array2[$j]);
push @weight_function, $weight;
}
push @input_matrix, \@weight_function;
}
print "\n The Nearest Neighbors and the Matrix of the weight function f(i,j) between each element of the two vectors \@array1 and \@array2.";
print "\n The weight function chosen can be the modulus of the difference between two real numbers: f(i,j) = abs (\$array1[i] - \$array2[j]). \n\n \@input_matrix = \n\n ";
print " " x 7;
printf("%4.0f ", $_ ) for @array2;
print "\n\n";
for my $i ( 0 .. $#input_matrix ) {
printf(" %4.0f [ ", $array1[$i] );
for my $j ( 0 .. $#{$input_matrix[$i]} ) {
printf("%4.0f ", $input_matrix[$i]->[$j] );
}
print "]\n";
}
my ( $optimal, $assignment_ref, $output_index_ref ) = auction( matrix_ref => \@input_matrix, maximize_total_benefit => 0, verbose => 10 );
print "\n";
my $sum = 0;
for my $i ( 0 .. $#{$output_index_ref} ){
lib/Algorithm/Bertsekas.pm view on Meta::CPAN
original matrix 10 x 10 with solution:
[ 84 94 75 56 66 95** 39 53 73 4 ]
[ 76** 71 56 49 29 1 40 40 72 72 ]
[ 85 100** 71 23 47 18 82 70 30 71 ]
[ 2 95 71 89 73 73 48 52 90** 51 ]
[ 65 28 77 73 24 28 75 48 8 81**]
[ 25 27 35 89 98 10 99** 3 27 4 ]
[ 58 15 99** 37 92 55 52 82 73 96 ]
[ 11 75 2 1 88** 43 8 28 98 20 ]
[ 52 95 10 38 41 64 20 75** 1 47 ]
[ 50 80 31 90** 10 83 51 55 57 40 ]
Pairs (in ascending order of matrix values):
indexes ( 8, 7 ), matrix value = 75 ; sum of values = 75
indexes ( 1, 0 ), matrix value = 76 ; sum of values = 151
indexes ( 4, 9 ), matrix value = 81 ; sum of values = 232
indexes ( 7, 4 ), matrix value = 88 ; sum of values = 320
indexes ( 3, 8 ), matrix value = 90 ; sum of values = 410
indexes ( 9, 3 ), matrix value = 90 ; sum of values = 500
indexes ( 0, 5 ), matrix value = 95 ; sum of values = 595
indexes ( 5, 6 ), matrix value = 99 ; sum of values = 694
indexes ( 6, 2 ), matrix value = 99 ; sum of values = 793
indexes ( 2, 1 ), matrix value = 100 ; sum of values = 893
=head1 OPTIONS
matrix_ref => \@input_matrix, reference to array: matrix N x M.
maximize_total_benefit => 0, 0: minimize the total benefit ; 1: maximize the total benefit.
inicial_stepsize => 1, auction algorithm terminates with a feasible assignment if the problem data are integer and stepsize < 1/min(N,M).
inicial_price => 0,
verbose => 3, print messages on the screen, level of verbosity, 0: quiet; 1, 2, 3, 4, 5, 8, 9, 10: debug information.
=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.
=head1 OUTPUT
The $output_index_ref is the reference to the output_index array.
The $assignment_ref is the reference to the assignment hash.
The $optimal is the total benefit which can be a minimum or maximum value.
=head1 SEE ALSO
1. Network Optimization: Continuous and Discrete Models (1998).
Dimitri P. Bertsekas
http://web.mit.edu/dimitrib/www/netbook_Full_Book.pdf
2. Towards auction algorithms for large dense assignment problems (2008).
Libor Bus and Pavel Tvrdik
https://pdfs.semanticscholar.org/b759/b8fb205df73c810b483b5be2b1ded62309b4.pdf
3. https://github.com/EvanOman/AuctionAlgorithmCPP/blob/master/auction.cpp
This Perl algorithm started from this C++ implementation.
4. https://en.wikipedia.org/wiki/Assignment_problem
5. https://en.wikipedia.org/wiki/Auction_algorithm
=head1 AUTHOR
Claudio Fernandes de Souza Rodrigues
May 21, 2018
Sao Paulo, Brasil
claudiofsr@yahoo.com
=head1 COPYRIGHT AND LICENSE
Copyright (c) 2018 Claudio Fernandes de Souza Rodrigues. All rights reserved.
This program is free software; you can redistribute it and/or modify
it under the same terms as Perl itself.
=cut
( run in 0.944 second using v1.01-cache-2.11-cpan-cdf2f3d4e48 )