4
\$\begingroup\$

The Hungarian algorithm solves the assignment problem, and I'm looking for any suggestions regarding improvement on my implementation (also coding style). It is based on the wikipedia entry and references therein. I have to use it repeatedly in calculations having cost matrices with anywhere between 15 and 40 rows or columns, note the cost matrix is not generally square.

I've got a working C version that will run in about 45 sec on my laptop with -O3 compiler options for a 15 x 15 cost matrix. The algorithm moves between steps 4 and 6 multiple times before convergence.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <assert.h>
#include <stdbool.h>
#include <time.h>

#define TOL 0.000001

typedef struct idx_t { int row;
              int col;
              int nrows;  // rows - n
              int ncols;  // cols - m
              int path_ct;
              double min;
              double* dm;
              int* M;
              int** path;
              int* path_row;
              int* path_col;
              int* row_cover;
              int* col_cover;
              int path_count;
              int path_row_0;
              int path_col_0;
           } idx_t;

double get_random() { return ((double)rand() / (double)RAND_MAX); }



struct idx_t* initIdxs( double** C, int n, int m )
{
    size_t i, j;
    struct idx_t* idxs = malloc( sizeof( struct idx_t ));
    assert( idxs != NULL );

    idxs->nrows = n;  // rows
    idxs->ncols = m;  // columns
    idxs->min = DBL_MAX;
    idxs->dm = ( double* )calloc(( m * n ), sizeof( double ));
    idxs->M = ( int* ) calloc( (n * m), sizeof( int ));
    idxs->path = ( int** ) calloc(( m + n ), sizeof( int* ));
    idxs->row_cover = ( int* ) calloc( n, sizeof( int ));
    idxs->col_cover = ( int* ) calloc( m, sizeof( int ));
    idxs->path_row = ( int* ) calloc( n, sizeof( int ));
    idxs->path_col = ( int* ) calloc( m, sizeof( int ));

    for ( i = 0; i < m + n; i++ )
         idxs->path[i] = ( int* ) calloc( 2, sizeof( int ));
    for ( i = 0; i < n; i++ ){
        for ( j = 0; j < m; j++ ){
            idxs->dm[i*n+j] = C[i][j];
        }
    }
    return idxs;
}

// For each row of the cost matrix, find the smallest element and subtract
// it from every element in its row. Go to Step 2.
int step_one( struct idx_t* idxs ) {
    double min_in_row;
    size_t r, c;
    size_t n = idxs->nrows;

    for ( r = 0; r < idxs->nrows; r++ ){
        min_in_row = idxs->dm[r*n];
        for ( c = 0; c < idxs->ncols; c++ ){
             if ( idxs->dm[r*n+c] < min_in_row )
                 min_in_row = idxs->dm[r*n+c];
        }
        for ( c = 0; c < idxs->ncols; c++ )
            idxs->dm[r*n+c] -= min_in_row;
    }
    return 2;  // step two
}

//Find a zero (Z) in the resulting matrix.  If there is no starred
//zero in its row or column, star Z. Repeat for each element, then go to Step 3.
int step_two( struct idx_t* idxs ){
    size_t r, c;
    size_t n = idxs->nrows;
    size_t m = idxs->ncols;
    size_t L = n < m ? n : m;
    int i;
    int* z;
if ( L == m )
    z = &idxs->col_cover[L];
else
    z = &idxs->row_cover[L];


for ( int i = 0; i < L; i++ ){
    if ( idxs->row_cover[r] == 0 && idxs->col_cover[c] == 0 ){
        if ( idxs->dm[r*n+c] < TOL ){
            idxs->M[r*n+c] = 1;
            idxs->dm[r*n+c] = 0.0;
            idxs->row_cover[r] = 1;
            idxs->col_cover[c] = 1;
        }
    }
}

for ( i = 0; i < L; i++ ){
    idxs->row_cover[i] = 0;
    idxs->col_cover[i] = 0;
}
for ( i = 0; i < m-L; i++ )
    z[i] = 0;

return 3; 
}

//Cover each column containing a starred zero.  If all columns are covered,
//these  describe a complete set of unique assignments. In such
//case, we're DONE. Otherwise, go to Step 4.
  int step_three( struct idx_t* idxs ){
    size_t r, c;
    size_t n = idxs->nrows;
    int col_count = 0;
    int step = 0;

    for ( r = 0; r < idxs->nrows; r++ ) {
        for ( c = 0; c < idxs->ncols; c++ ) {
            if ( idxs->M[r*n+c] == 1 ) {
                idxs->col_cover[c] = 1;
                col_count += 1;
            }
        }
    }

    if ( col_count == idxs->ncols || col_count == idxs->nrows ) {
         step = 7;
    } else {
         step = 4;
    }
    return step;
}

void find_zero( struct idx_t* idxs ){
     size_t r = 0;
     size_t c = 0;
     size_t n = idxs->nrows;
     bool done = false;
     idxs->row = -1;
     idxs->col = -1;

     while ( !done ){
        c = 0;
        while ( true ){
            if ( idxs->row_cover[r] == 0 && idxs->col_cover[c] == 0 && idxs->dm[r*n+c] < TOL ){
                 idxs->row = r;
                 idxs->col = c;
                 idxs->dm[r*n+c] = 0;
                 done = true;
            }
            c += 1;
            if ( c >= idxs->ncols || done )
                break;
        }
        r += 1;
        if ( r >= idxs->nrows )
            done = true;
    }
}
bool star_in_row( struct idx_t* idxs ){
    bool tmp = false;
    size_t n = idxs->nrows;
    size_t c = 0;

    for ( c = 0; c < idxs->ncols; c++ ){
        if ( idxs->M[idxs->row*n+c] == 1 ){
            tmp = true;
            idxs->col = c;
        }
    }
    return tmp;
}
void find_star_in_row( struct idx_t* idxs ){
    size_t c = 0;
    idxs->col = -1;
    size_t n = idxs->nrows;

    for ( c = 0; c < idxs->ncols; c++ ){
         if ( idxs->M[idxs->row*n+c] == 1 )
             idxs->col = c;
    }
}
//Find a noncovered zero and prime it.  If there is no starred zero
//in the row containing it, go to 5.  Otherwise,
//cover this row and uncover the column containing the starred zero.
//Continue until there are no uncovered zeros left.
//Save the smallest uncovered value and go to Step 6.
int step_four( struct idx_t* idxs ){
    int step = 0;
    size_t n = idxs->nrows;
    bool done = false;
    idxs->row = -1;
    idxs->col = -1;

    while ( !done ){
        find_zero( idxs );
        if ( idxs->row == -1 ){
            done = true;
            step = 6;
        } else {
            idxs->M[idxs->row*n+idxs->col] = 2;
            if ( star_in_row( idxs )){
                idxs->row_cover[idxs->row] = 1;
                idxs->col_cover[idxs->col] = 0;
            } else {
                done = true;
                step = 5;
                idxs->path_row_0 = idxs->row;
                idxs->path_col_0 = idxs->col;
            }
         }
    }
    return step;
}
void find_star_in_col( struct idx_t* idxs ){
    size_t i;
    size_t col = idxs->path[idxs->path_ct][1];
    idxs->row = -1;
    size_t n = idxs->nrows;

    for ( i = 0; i < idxs->nrows; i++ ){
        if ( idxs->M[i*n+col] == 1 )
            idxs->row = i;
    }
}

void find_prime_in_row( struct idx_t* idxs ){
    size_t j;
    size_t row = idxs->path[idxs->path_ct][0];
    size_t n = idxs->nrows;

    for ( j = 0; j < idxs->ncols; j++ ){
        if ( idxs->M[row*n+j] == 2 )
            idxs->col = j;
    }
}

void augment_path( struct idx_t* idxs ){
    size_t n = idxs->nrows;
    for ( int p = 0; p < idxs->path_ct + 1; p++ )
        if ( idxs->M[idxs->path[p][0]*n+idxs->path[p][1]] == 1 ){
            idxs->M[idxs->path[p][0]*n+idxs->path[p][1]] = 0;
        } else {
            idxs->M[idxs->path[p][0]*n+idxs->path[p][1]] = 1;
        }
}
void clear_covers( struct idx_t* idxs ){
    size_t r, c;
    size_t n = idxs->nrows;
    size_t m = idxs->ncols;
    size_t L = n < m ? n : m;
    int i;
    int* z;
    if ( L == m )
        z = &idxs->col_cover[L];
    else
        z = &idxs->row_cover[L];
    for ( i = 0; i < L; i++ ){
        idxs->row_cover[i] = 0;
        idxs->col_cover[i] = 0;
    }
    for ( i = 0; i < m-L; i++ )
        z[i] = 0;
}
void erase_primes( struct idx_t* idxs ){
    size_t r, c;
    size_t n = idxs->nrows;

    for ( r = 0; r < idxs->nrows; r++ ){
        for ( c = 0; c < idxs->ncols; c++ )
            if ( idxs->M[r*n+c] == 2 )
                idxs->M[r*n+c] = 0;
    }
}
//Construct a series of alternating primed and starred zeros.
//Z0 is the uncovered primed zero found in Step 4.  Z1 is
//the starred zero in the column of Z0 (if any). Z2 is the primed zero
//in the row of Z1. Continue until the series
//ends at a primed zero that has no starred zero in its column.
//Clear each starred zero of the series, star each primed zero of the series.
//Then erase all primes and uncover every line in the matrix. Return to Step 3.
int step_five( struct idx_t* idxs ) {
    bool done = false;
    int tmp;

    idxs->path_ct = 0;
    idxs->row = -1;
    idxs->col = -1;
    idxs->path[0][0] = idxs->path_row_0;
    idxs->path[0][1] = idxs->path_col_0;

    while ( !done ){
        find_star_in_col( idxs );
        if ( idxs->row > -1 ){
            idxs->path_ct += 1;
            tmp = idxs->path_ct - 1;
            idxs->path[idxs->path_ct][0] = idxs->row;
            idxs->path[idxs->path_ct][1] = idxs->path[tmp][1];
        } else {
            done = true;
        }
        if ( !done ){
            find_prime_in_row( idxs );
            idxs->path_ct += 1;
            tmp = idxs->path_ct - 1;
            idxs->path[idxs->path_ct][0] = idxs->path[tmp][0];
            idxs->path[idxs->path_ct][1] = idxs->col;
        }
    }
    augment_path( idxs );
    clear_covers( idxs );
    erase_primes( idxs );
    return 3;
}

void min_val( struct idx_t* idxs ) {
    size_t r, c;
    size_t n = idxs->nrows;

    for ( r = 0; r < idxs->nrows; r++ ){
        if ( idxs->row_cover[r] == 0 ) {
            for ( c = 0; c < idxs->ncols; c++ ){
                if ( idxs->col_cover[c] == 0 ){
                    if ( idxs->min > idxs->dm[r*n+c] )
                        idxs->min = idxs->dm[r*n+c];
                }
            }
        }
    }
}
//Add the value found in Step 4 to every element of each covered row, and subtract
//it from every element of each uncovered column.  Return to Step 4 
int step_six( struct idx_t* idxs ) {
     size_t r, c;
     size_t n = idxs->nrows;
     min_val( idxs );

    for ( r = 0; r < idxs->nrows; r++ ){
        for ( c = 0; c < idxs->ncols; c++ ){
            if ( idxs->row_cover[r] == 1 )
                idxs->dm[r*n+c] += idxs->min;
            if( idxs->col_cover[c] == 0 )
                idxs->dm[r*n+c] -= idxs->min;
        }
    }
    return 4;
}
double calc_cost( double** C, struct idx_t* idxs ){
    size_t i, j;
    size_t n = idxs->nrows;
    double cost = 0;

    for ( i = 0; i < idxs->nrows; i++ ){
        for ( j = 0; j < idxs->ncols; j++ )
             if ( idxs->M[i*n+j] == 1 )
                cost += C[i][j];
    }
    return cost;
}
void deleteMem(struct idx_t* idxs ){
    size_t i;
    for ( i = 0; i < idxs->ncols + idxs->nrows; i++ )
        free( idxs->path[i] );

    free( idxs->path_row );
    free( idxs->path_col );
    free( idxs->row_cover );
    free( idxs->col_cover );
    free( idxs->path);
    free( idxs->dm );
    free( idxs->M );
    free( idxs );
     return;
 }
double assignment( double** C, int m, int n ){
    bool done = false;
    unsigned int start;
    double diff;
    double msec;
    double sec;
    double cps = 1000000.0;
    struct idx_t* idxs;
    idxs = initIdxs( C, m, n );
    int step = 1;
   double cost = 0;

    while ( !done ){
        switch ( step ){
            case 1:
                // start = clock();
                step = step_one( idxs );
                // diff = clock() - start;
                // msec = diff * (1000.0 / cps);
                // sec = diff / cps;
                // printf("Step 1: Time taken %f seconds %f milliseconds\n", sec, msec);
                break;
            case 2:
                // start = clock();
                step = step_two( idxs );
                // diff = clock() - start;
                // msec = diff * (1000.0 / cps);
                // sec = diff / cps;
                // printf("Step 2: Time taken %f seconds %f milliseconds\n", sec, msec);
                break;
            case 3:
                // start = clock();
                step = step_three( idxs );
                // diff = clock() - start;
                // msec = diff * (1000.0 / cps);
                // sec = diff / cps;
                // printf("Step 3: Time taken %f seconds %f milliseconds\n", sec, msec);
                break;
            case 4:
                // start = clock();
                step = step_four( idxs );
                // diff = clock() - start;
                // msec = diff * (1000.0 / cps);
                // sec = diff / cps;
                // printf("Step 4: Time taken %f seconds %f milliseconds\n", sec, msec);
                break;
            case 5:
                // start = clock();
                step = step_five( idxs );
                // diff = clock() - start;
                // msec = diff * (1000.0 / cps);
                // sec = diff / cps;
                // printf("Step 5: Time taken %f seconds %f milliseconds\n", sec, msec);
                break;
            case 6:
                // start = clock();
                step = step_six( idxs );
                // diff = clock() - start;
                // msec = diff * (1000.0 / cps);
                // sec = diff / cps;
                // printf("Step 6: Time taken %f seconds %f milliseconds\n", sec, msec);
                break;
            case 7:
                // start = clock();
                done = true;
                cost = calc_cost( C, idxs );
                // diff = clock() - start;
                // msec = diff * (1000.0 / cps);
                // sec = diff / cps;
                // printf("Step 7: Time taken %f seconds %f milliseconds\n", sec, msec);
                break;
        }
    }
    deleteMem( idxs );
    return cost;
}

int main(){
    size_t i, j;
    unsigned int start;
    double diff;
    double msec;
    double sec;
    double cps = 1000000.0;
    int n = 15;
    double c = 0;
    double** C = malloc( n  * sizeof( double* ));
    double* tmp = malloc(sizeof( double ) * n * n );

    for( i = 0; i < n; i++ )
         C[i] = &(tmp[i * n]);

    for( i = 0; i < n; i++ ){
        for( j = 0; j < n; j++ )
            C[i][j] = (i + 1) * (j + 1) + get_random();
    }

    start = clock();
    c = assignment(C, n, n);
    printf("\nIn c, cost = %10.4lf\n", c);
    diff = clock() - start;
    msec = diff * (1000.0 / cps);
    sec = diff / cps;
    printf("Total time: %f seconds %f milliseconds\n", sec, msec);

    free( C );
    free( tmp );

    return 0;
}
\$\endgroup\$
  • \$\begingroup\$ What range of values do you expect with get_random()? should it include 1.0? \$\endgroup\$ – chux Sep 28 '18 at 3:33
1
\$\begingroup\$

I'm looking for any suggestions regarding improvement on my implementation (also coding style).

Just some basic things.

No need for casting.

See Do I cast the result of malloc?

// idxs->dm = ( double* )calloc(( m * n ), sizeof( double ));
idxs->dm = calloc(( m * n ), sizeof( double ));

Use the size of the referenced object.

sizeof( double ) obliges the right type used with idxs->dm. By using sizeof idxs->dm[0] coordinating the type is not necessary. It is easier to code right the first time, easier to review and maintain.

// idxs->dm = calloc(( m * n ), sizeof( double ));
idxs->dm = calloc( m * n, sizeof idxs->dm[0]);

*alloc() can fail

Run time allocations can fail. Robust code check for that.

idxs->dm = calloc( m * n, sizeof idxs->dm[0]);
if (idxs->dm == NULL) Handle_Error_Somehow();

Recommend clearing

As struct idx_t has so many members, I'd recommend clear the allocation first via calloc() ormemset().

// struct idx_t* idxs = malloc( sizeof( struct idx_t ));
struct idx_t* idxs = calloc(1, sizeof *idxs);

Avoid casual int size_t mixing

size_t is unsigned and int is signed. With all the array indexing and sizing, best o use one type. Usually it is best to use size_t for all sizing and indexing as that type is neither too wide nor too narrow. Many standard functions employing a size parameter use size_t.

Be careful to avoid indexing before 0. Can't do that with size_t.

Unused variables

double diff, msec, sec, cps are not used.

\$\endgroup\$
  • \$\begingroup\$ Thanks for the style suggestions. Do you have nay ideas about improving the running time? Currently the version in scipy beats mine. \$\endgroup\$ – Bob Dobalina Oct 9 '18 at 19:23
  • \$\begingroup\$ I did not find the code (with its 60+ warnings) well partitioned nor description clear enough to begin a performance improvement suggestion. IMO first try smaller tests that need 2-5 seconds and improve that. Posting your exact output would help, mine was In c, cost = 688.2540 Total time: 0.073594 seconds 73.594000 milliseconds yet took many more seconds than to 73.5ms to run. \$\endgroup\$ – chux Oct 9 '18 at 19:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.