**
** Parts of the used code was originally provided by the
** "Stanford GraphGase", but I made changes to this code.
- ** As asked by the copyright node of the "Stanford GraphGase",
+ ** As asked by the copyright node of the "Stanford GraphGase",
** I hereby proclaim that this file are *NOT* part of the
** "Stanford GraphGase" distrubition!
**
#define verbose (0)
typedef struct {
- int num_rows;
- int num_cols;
- double** cost;
- int** assignment;
+ int num_rows;
+ int num_cols;
+ double** cost;
+ int** assignment;
} hungarian_problem_t;
int hungarian_imax(int a, int b) {
- return (a<b)?b:a;
+ return (a<b)?b:a;
}
/** This method initialize the hungarian_problem structure and init
* the cost matrices (missing lines or columns are filled with 0).
* It returns the size of the quadratic(!) assignment matrix. **/
int hungarian_init(hungarian_problem_t* p, double** cost_matrix, int rows, int cols,
- int mode)
+ int mode)
{
- int i,j;
- double max_cost = 0.;
- int org_cols = cols,
- org_rows = rows;
-
- // is the number of cols not equal to number of rows ?
- // if yes, expand with 0-cols / 0-cols
- rows = hungarian_imax(cols, rows);
- cols = rows;
-
- p->num_rows = rows;
- p->num_cols = cols;
-
- p->cost = (double**)calloc(rows,sizeof(double*));
- p->assignment = (int**)calloc(rows,sizeof(int*));
-
-
for(i=0; i<p->num_rows; i++) {
- p->cost[i] = (double*)calloc(cols,sizeof(double));
- p->assignment[i] = (int*)calloc(cols,sizeof(int));
-
for(j=0; j<p->num_cols; j++) {
- p->cost[i][j] = (i < org_rows && j < org_cols) ? cost_matrix[i][j] : 0.;
- p->assignment[i][j] = 0;
-
- if (max_cost < p->cost[i][j])
- max_cost = p->cost[i][j];
- }
- }
-
- if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
-
for(i=0; i<p->num_rows; i++) {
-
for(j=0; j<p->num_cols; j++)
- p->cost[i][j] = max_cost - p->cost[i][j];
- }
- }
- else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
- // nothing to do
- }
-// else
-// fprintf(stderr,"%s: unknown mode. Mode was set to HUNGARIAN_MODE_MINIMIZE_COST !\n", __FUNCTION__);
-
- return rows;
+ int i,j;
+ double max_cost = 0.;
+ int org_cols = cols,
+ org_rows = rows;
+
+ // is the number of cols not equal to number of rows ?
+ // if yes, expand with 0-cols / 0-cols
+ rows = hungarian_imax(cols, rows);
+ cols = rows;
+
+ p->num_rows = rows;
+ p->num_cols = cols;
+
+ p->cost = (double**)calloc(rows,sizeof(double*));
+ p->assignment = (int**)calloc(rows,sizeof(int*));
+
+ for(i=0; i<p->num_rows; i++) {
+ p->cost[i] = (double*)calloc(cols,sizeof(double));
+ p->assignment[i] = (int*)calloc(cols,sizeof(int));
+ for(j=0; j<p->num_cols; j++) {
+ p->cost[i][j] = (i < org_rows && j < org_cols) ? cost_matrix[i][j] : 0.;
+ p->assignment[i][j] = 0;
+
+ if (max_cost < p->cost[i][j])
+ max_cost = p->cost[i][j];
+ }
+ }
+
+ if (mode == HUNGARIAN_MODE_MAXIMIZE_UTIL) {
+ for(i=0; i<p->num_rows; i++) {
+ for(j=0; j<p->num_cols; j++)
+ p->cost[i][j] = max_cost - p->cost[i][j];
+ }
+ }
+ else if (mode == HUNGARIAN_MODE_MINIMIZE_COST) {
+ // nothing to do
+ }
+// else
+// fprintf(stderr,"%s: unknown mode. Mode was set to HUNGARIAN_MODE_MINIMIZE_COST !\n", __FUNCTION__);
+
+ return rows;
}
/** Free the memory allocated by init. **/
void hungarian_free(hungarian_problem_t* p) {
- int i;
-
for(i=0; i<p->num_rows; i++) {
- free(p->cost[i]);
- free(p->assignment[i]);
- }
- free(p->cost);
- free(p->assignment);
- p->cost = NULL;
- p->assignment = NULL;
+ int i;
+ for(i=0; i<p->num_rows; i++) {
+ free(p->cost[i]);
+ free(p->assignment[i]);
+ }
+ free(p->cost);
+ free(p->assignment);
+ p->cost = NULL;
+ p->assignment = NULL;
}
/** This method computes the optimal assignment. **/
void hungarian_solve(hungarian_problem_t* p)
{
- int i, j, m, n, k, l, t, q, unmatched;
- double cost, s;
- int* col_mate;
- int* row_mate;
- int* parent_row;
- int* unchosen_row;
- double* row_dec;
- double* col_inc;
- double* slack;
- int* slack_row;
-
- cost = 0.;
- m =p->num_rows;
- n =p->num_cols;
-
- col_mate = (int*)calloc(p->num_rows,sizeof(int));
- unchosen_row = (int*)calloc(p->num_rows,sizeof(int));
- row_dec = (double*)calloc(p->num_rows,sizeof(double));
- slack_row = (int*)calloc(p->num_rows,sizeof(int));
-
- row_mate = (int*)calloc(p->num_cols,sizeof(int));
- parent_row = (int*)calloc(p->num_cols,sizeof(int));
- col_inc = (double*)calloc(p->num_cols,sizeof(double));
- slack = (double*)calloc(p->num_cols,sizeof(double));
-
-
for (i=0; i<p->num_rows; i++) {
- col_mate[i]=0;
- unchosen_row[i]=0;
- row_dec[i]=0.;
- slack_row[i]=0;
- }
-
for (j=0; j<p->num_cols; j++) {
- row_mate[j]=0;
- parent_row[j] = 0;
- col_inc[j]=0.;
- slack[j]=0.;
- }
-
-
for (i=0; i<p->num_rows; ++i)
-
for (j=0; j<p->num_cols; ++j)
- p->assignment[i][j]=HUNGARIAN_NOT_ASSIGNED;
-
- // Begin subtract column minima in order to start with lots of zeroes 12
- for (l=0; l<n; l++)
- {
- s=p->cost[0][l];
- for (k=1; k<m; k++)
- if (p->cost[k][l]<s)
- s=p->cost[k][l];
- cost+=s;
- if (s!=0.)
- for (k=0; k<m; k++)
- p->cost[k][l]-=s;
- }
- // End subtract column minima in order to start with lots of zeroes 12
-
- // Begin initial state 16
- t=0;
- for (l=0; l<n; l++)
- {
- row_mate[l]= -1;
- parent_row[l]= -1;
- col_inc[l]=0.;
- slack[l]=INF;
- }
- for (k=0; k<m; k++)
- {
- s=p->cost[k][0];
- for (l=1; l<n; l++)
- if (p->cost[k][l]<s)
- s=p->cost[k][l];
- row_dec[k]=s;
- for (l=0; l<n; l++)
- if (s==p->cost[k][l] && row_mate[l]<0)
- {
- col_mate[k]=l;
- row_mate[l]=k;
- goto row_done;
- }
- col_mate[k]= -1;
- unchosen_row[t++]=k;
+ int i, j, m, n, k, l, t, q, unmatched;
+ double cost, s;
+ int* col_mate;
+ int* row_mate;
+ int* parent_row;
+ int* unchosen_row;
+ double* row_dec;
+ double* col_inc;
+ double* slack;
+ int* slack_row;
+
+ cost = 0.;
+ m =p->num_rows;
+ n =p->num_cols;
+
+ col_mate = (int*)calloc(p->num_rows,sizeof(int));
+ unchosen_row = (int*)calloc(p->num_rows,sizeof(int));
+ row_dec = (double*)calloc(p->num_rows,sizeof(double));
+ slack_row = (int*)calloc(p->num_rows,sizeof(int));
+
+ row_mate = (int*)calloc(p->num_cols,sizeof(int));
+ parent_row = (int*)calloc(p->num_cols,sizeof(int));
+ col_inc = (double*)calloc(p->num_cols,sizeof(double));
+ slack = (double*)calloc(p->num_cols,sizeof(double));
+
+ for (i=0; i<p->num_rows; i++) {
+ col_mate[i]=0;
+ unchosen_row[i]=0;
+ row_dec[i]=0.;
+ slack_row[i]=0;
+ }
+ for (j=0; j<p->num_cols; j++) {
+ row_mate[j]=0;
+ parent_row[j] = 0;
+ col_inc[j]=0.;
+ slack[j]=0.;
+ }
+
+ for (i=0; i<p->num_rows; ++i)
+ for (j=0; j<p->num_cols; ++j)
+ p->assignment[i][j]=HUNGARIAN_NOT_ASSIGNED;
+
+ // Begin subtract column minima in order to start with lots of zeroes 12
+ for (l=0; l<n; l++)
+ {
+ s=p->cost[0][l];
+ for (k=1; k<m; k++)
+ if (p->cost[k][l]<s)
+ s=p->cost[k][l];
+ cost+=s;
+ if (s!=0.)
+ for (k=0; k<m; k++)
+ p->cost[k][l]-=s;
+ }
+ // End subtract column minima in order to start with lots of zeroes 12
+
+ // Begin initial state 16
+ t=0;
+ for (l=0; l<n; l++)
+ {
+ row_mate[l]= -1;
+ parent_row[l]= -1;
+ col_inc[l]=0.;
+ slack[l]=INF;
+ }
+ for (k=0; k<m; k++)
+ {
+ s=p->cost[k][0];
+ for (l=1; l<n; l++)
+ if (p->cost[k][l]<s)
+ s=p->cost[k][l];
+ row_dec[k]=s;
+ for (l=0; l<n; l++)
+ if (s==p->cost[k][l] && row_mate[l]<0)
+ {
+ col_mate[k]=l;
+ row_mate[l]=k;
+ goto row_done;
+ }
+ col_mate[k]= -1;
+ unchosen_row[t++]=k;
row_done:
- ;
- }
- // End initial state 16
-
- // Begin Hungarian algorithm 18
- if (t==0)
- goto done;
- unmatched=t;
- while (1)
- {
- q=0;
- while (1)
- {
- while (q<t)
- {
- // Begin explore node q of the forest 19
- {
- k=unchosen_row[q];
- s=row_dec[k];
- for (l=0; l<n; l++)
- if (slack[l] > 0.)
- {
- double del=p->cost[k][l]-s+col_inc[l];
- if (del<slack[l])
- {
- if (del==0.)
- {
- if (row_mate[l]<0)
- goto breakthru;
- slack[l]=0.;
- parent_row[l]=k;
- unchosen_row[t++]=row_mate[l];
- }
- else
- {
- slack[l]=del;
- slack_row[l]=k;
- }
- }
- }
- }
- // End explore node q of the forest 19
- q++;
- }
-
- // Begin introduce a new zero into the matrix 21
- s=INF;
- for (l=0; l<n; l++)
- if (slack[l]>0. && slack[l]<s)
- s=slack[l];
- for (q=0; q<t; q++)
- row_dec[unchosen_row[q]]+=s;
- for (l=0; l<n; l++)
- if (slack[l]>0.)
- {
- slack[l]-=s;
- if (slack[l]==0.)
- {
- // Begin look at a new zero 22
- k=slack_row[l];
- if (row_mate[l]<0)
- {
- for (j=l+1; j<n; j++)
- if (slack[j]==0.)
- col_inc[j]+=s;
- goto breakthru;
- }
- else
- {
- parent_row[l]=k;
- unchosen_row[t++]=row_mate[l];
- }
- // End look at a new zero 22
- }
- }
- else
- col_inc[l]+=s;
- // End introduce a new zero into the matrix 21
- }
+ ;
+ }
+ // End initial state 16
+
+ // Begin Hungarian algorithm 18
+ if (t==0)
+ goto done;
+ unmatched=t;
+ while (1)
+ {
+ q=0;
+ while (1)
+ {
+ while (q<t)
+ {
+ // Begin explore node q of the forest 19
+ {
+ k=unchosen_row[q];
+ s=row_dec[k];
+ for (l=0; l<n; l++)
+ if (slack[l] > 0.)
+ {
+ double del=p->cost[k][l]-s+col_inc[l];
+ if (del<slack[l])
+ {
+ if (del==0.)
+ {
+ if (row_mate[l]<0)
+ goto breakthru;
+ slack[l]=0.;
+ parent_row[l]=k;
+ unchosen_row[t++]=row_mate[l];
+ }
+ else
+ {
+ slack[l]=del;
+ slack_row[l]=k;
+ }
+ }
+ }
+ }
+ // End explore node q of the forest 19
+ q++;
+ }
+
+ // Begin introduce a new zero into the matrix 21
+ s=INF;
+ for (l=0; l<n; l++)
+ if (slack[l]>0. && slack[l]<s)
+ s=slack[l];
+ for (q=0; q<t; q++)
+ row_dec[unchosen_row[q]]+=s;
+ for (l=0; l<n; l++)
+ if (slack[l]>0.)
+ {
+ slack[l]-=s;
+ if (slack[l]==0.)
+ {
+ // Begin look at a new zero 22
+ k=slack_row[l];
+ if (row_mate[l]<0)
+ {
+ for (j=l+1; j<n; j++)
+ if (slack[j]==0.)
+ col_inc[j]+=s;
+ goto breakthru;
+ }
+ else
+ {
+ parent_row[l]=k;
+ unchosen_row[t++]=row_mate[l];
+ }
+ // End look at a new zero 22
+ }
+ }
+ else
+ col_inc[l]+=s;
+ // End introduce a new zero into the matrix 21
+ }
breakthru:
- // Begin update the matching 20
- while (1)
- {
- j=col_mate[k];
- col_mate[k]=l;
- row_mate[l]=k;
- if (j<0)
- break;
- k=parent_row[j];
- l=j;
- }
- // End update the matching 20
- if (--unmatched==0)
- goto done;
- // Begin get ready for another stage 17
- t=0;
- for (l=0; l<n; l++)
- {
- parent_row[l]= -1;
- slack[l]=INF;
- }
- for (k=0; k<m; k++)
- if (col_mate[k]<0)
- unchosen_row[t++]=k;
- // End get ready for another stage 17
- }
+ // Begin update the matching 20
+ while (1)
+ {
+ j=col_mate[k];
+ col_mate[k]=l;
+ row_mate[l]=k;
+ if (j<0)
+ break;
+ k=parent_row[j];
+ l=j;
+ }
+ // End update the matching 20
+ if (--unmatched==0)
+ goto done;
+ // Begin get ready for another stage 17
+ t=0;
+ for (l=0; l<n; l++)
+ {
+ parent_row[l]= -1;
+ slack[l]=INF;
+ }
+ for (k=0; k<m; k++)
+ if (col_mate[k]<0)
+ unchosen_row[t++]=k;
+ // End get ready for another stage 17
+ }
done:
-/* // Begin doublecheck the solution 23
- for (k=0; k<m; k++)
- for (l=0; l<n; l++)
- if (p->cost[k][l]<row_dec[k]-col_inc[l])
- exit(0);
- for (k=0; k<m; k++)
- {
- l=col_mate[k];
- if (l<0 || p->cost[k][l]!=row_dec[k]-col_inc[l])
- exit(0);
- }
- k=0;
- for (l=0; l<n; l++)
- if (col_inc[l])
- k++;
- if (k>m)
- exit(0);
- // End doublecheck the solution 23
-*/ // End Hungarian algorithm 18
-
- for (i=0; i<m; ++i)
- {
- p->assignment[i][col_mate[i]]=HUNGARIAN_ASSIGNED;
- /*TRACE("%d - %d\n", i, col_mate[i]);*/
- }
- for (k=0; k<m; ++k)
- {
- for (l=0; l<n; ++l)
- {
- /*TRACE("%d ",p->cost[k][l]-row_dec[k]+col_inc[l]);*/
- p->cost[k][l]=p->cost[k][l]-row_dec[k]+col_inc[l];
- }
- /*TRACE("\n");*/
- }
- for (i=0; i<m; i++)
- cost+=row_dec[i];
- for (i=0; i<n; i++)
- cost-=col_inc[i];
-
- free(slack);
- free(col_inc);
- free(parent_row);
- free(row_mate);
- free(slack_row);
- free(row_dec);
- free(unchosen_row);
- free(col_mate);
+/* // Begin doublecheck the solution 23
+ for (k=0; k<m; k++)
+ for (l=0; l<n; l++)
+ if (p->cost[k][l]<row_dec[k]-col_inc[l])
+ exit(0);
+ for (k=0; k<m; k++)
+ {
+ l=col_mate[k];
+ if (l<0 || p->cost[k][l]!=row_dec[k]-col_inc[l])
+ exit(0);
+ }
+ k=0;
+ for (l=0; l<n; l++)
+ if (col_inc[l])
+ k++;
+ if (k>m)
+ exit(0);
+ // End doublecheck the solution 23
+*/ // End Hungarian algorithm 18
+
+ for (i=0; i<m; ++i)
+ {
+ p->assignment[i][col_mate[i]]=HUNGARIAN_ASSIGNED;
+ /*TRACE("%d - %d\n", i, col_mate[i]);*/
+ }
+ for (k=0; k<m; ++k)
+ {
+ for (l=0; l<n; ++l)
+ {
+ /*TRACE("%d ",p->cost[k][l]-row_dec[k]+col_inc[l]);*/
+ p->cost[k][l]=p->cost[k][l]-row_dec[k]+col_inc[l];
+ }
+ /*TRACE("\n");*/
+ }
+ for (i=0; i<m; i++)
+ cost+=row_dec[i];
+ for (i=0; i<n; i++)
+ cost-=col_inc[i];
+
+ free(slack);
+ free(col_inc);
+ free(parent_row);
+ free(row_mate);
+ free(slack_row);
+ free(row_dec);
+ free(unchosen_row);
+ free(col_mate);
}
// Get the optimal assignment, by calling hungarian_solve above; "distances" in columns
void hungarianAlgorithm(double* distances, int* pn, int* assignment)
{
- int n = *pn;
- double** mat = array_to_matrix(distances, n, n);
+ int n = *pn;
+ double** mat = array_to_matrix(distances, n, n);
hungarian_problem_t p;
- hungarian_init(&p, mat, n, n, HUNGARIAN_MODE_MINIMIZE_COST);
+ hungarian_init(&p, mat, n, n, HUNGARIAN_MODE_MINIMIZE_COST);
- hungarian_solve(&p);
+ hungarian_solve(&p);
- // Copy the optimal assignment in a R-friendly structure
- for (int i=0; i < n; i++) {
- for (int j=0; j < n; j++) {
- if (p.assignment[i][j])
- assignment[i] = j+1; //add 1 since we work with R
- }
- }
+ // Copy the optimal assignment in a R-friendly structure
+ for (int i=0; i < n; i++) {
+ for (int j=0; j < n; j++) {
+ if (p.assignment[i][j])
+ assignment[i] = j+1; //add 1 since we work with R
+ }
+ }
hungarian_free(&p);
- for (int i=0; i<n; i++)
- free(mat[i]);
- free(mat);
+ for (int i=0; i<n; i++)
+ free(mat[i]);
+ free(mat);
}
projects / morpheus.git / blobdiff