22 functions analysed · 19 typed · 0 untypeable · 3 timed out · 10 entry points (10 typed) · checker unknown
int ap_assignment(AP *p, int *res)
{
int i;
if(p->s == NULL)
ap_hungarian(p);
for(i = 0; i < p->n; i++)
res[i] = p->s[i];
return p->n;
}AP *ap_create_problem(double *t, int n)
{
int i,j;
AP *p;
p = (AP*) malloc(sizeof(AP));
if(p == NULL)
return NULL;
p->n = n;
p->C = (double **) malloc((n + 1) * sizeof(double *));
p->c = (double **) malloc((n + 1) * sizeof(double *));
if(p->C == NULL || p->c == NULL)
return NULL;
for(i = 1; i <= n; i++){
p->C[i] = (double *) calloc(n + 1, sizeof(double));
p->c[i] = (double *) calloc(n + 1, sizeof(double));
if(p->C[i] == NULL || p->c[i] == NULL)
return NULL;
}
for(i = 1; i <= n; i++)
for( j = 1; j <= n; j++){
p->C[i][j] = t[n*(j - 1) + i - 1];
p->c[i][j] = t[n*(j - 1) + i - 1];
}
p->cost = 0;
p->s = NULL;
p->f = NULL;
return p;
}void ap_free(AP *p)
{
int i;
free(p->s);
free(p->f);
for(i = 1; i <= p->n; i++){
free(p->C[i]);
free(p->c[i]);
}
free(p->C);
free(p->c);
free(p);
}clue_dissimilarity_count_inversions(double *x, double *y, Sint *n,
double *count)
{
Sint i, j;
for(i = 0; i < *n; i++)
for(j = 0; j < *n; j++)
if((clue_sign(x[i] - x[j]) * clue_sign(y[i] - y[j])) < 0)
(*count)++;
}clue_sign(double x)
{
if(x == 0) return(0);
return((x > 0) ? 1 : -1);
}double **clue_vector_to_square_matrix(double *x, Sint n)
{
double **data, *val;
Sint i, j;
data = (double **) R_alloc(n, sizeof(double));
for(i = 0; i < n; i++) {
data[i] = (double *) R_alloc(n, sizeof(double));
val = x + i;
for(j = 0; j < n; j++, val += n)
data[i][j] = *val;
}
return(data);
}deviation_from_additivity(double *x, int *n, double *v, int *max)
{
double **D, p, delta, A, B, C;
int i, j, k, l;
D = clue_vector_to_square_matrix(x, *n);
p = 0;
for(i = 0; i < *n - 3; i++)
for(j = i + 1; j < *n - 2; j++)
for(k = j + 1; k < *n - 1; k++)
for(l = k + 1; l < *n; l++) {
A = D[i][j] + D[k][l];
B = D[i][k] + D[j][l];
C = D[i][l] + D[j][k];
if((A <= B) && (A <= C))
delta = (C - B);
else if(B <= C)
delta = (A - C);
else
delta = (B - A);
if(*max)
p = fmax2(p, fabs(delta));
else
p += delta * delta;
}
*v = p;
}deviation_from_additivity_gradient(double *x, int *n, double *out)
{
double **D, **G, A, B, C, delta;
int i, j, k, l;
D = clue_vector_to_square_matrix(x, *n);
G = clue_vector_to_square_matrix(out, *n);
for(i = 0; i < *n - 3; i++)
for(j = i + 1; j < *n - 2; j++)
for(k = j + 1; k < *n - 1; k++)
for(l = k + 1; l < *n; l++) {
A = D[i][j] + D[k][l];
B = D[i][k] + D[j][l];
C = D[i][l] + D[j][k];
if((A <= B) && (A <= C)) {
delta = 2 * (B - C);
G[i][l] -= delta;
G[j][k] -= delta;
G[i][k] += delta;
G[j][l] += delta;
}
else if(B <= C) {
delta = 2 * (C - A);
G[i][l] += delta;
G[j][k] += delta;
G[i][j] -= delta;
G[k][l] -= delta;
}
else {
delta = 2 * (A - B);
G[i][k] -= delta;
G[j][l] -= delta;
G[i][j] += delta;
G[k][l] += delta;
}
}
for(i = 0; i < *n; i++)
for(j = 0; j < *n; j++)
*out++ = G[i][j];
}deviation_from_ultrametricity(double *x, int *n, double *v, int *max)
{
double **D, p, delta, A, B, C;
int i, j, k;
D = clue_vector_to_square_matrix(x, *n);
p = 0;
for(i = 0; i < *n - 2; i++)
for(j = i + 1; j < *n - 1; j++) {
A = D[i][j];
for(k = j + 1; k < *n; k++) {
B = D[i][k];
C = D[j][k];
if((A <= B) && (A <= C))
delta = C - B;
else if(B <= C)
delta = A - C;
else
delta = B - A;
if(*max)
p = fmax2(p, fabs(delta));
else
p += delta * delta;
}
}
*v = p;
}deviation_from_ultrametricity_gradient(double *x, int *n, double *out)
{
double **D, **G, A, B, C, delta;
int i, j, k;
D = clue_vector_to_square_matrix(x, *n);
G = clue_vector_to_square_matrix(out, *n);
for(i = 0; i < *n - 2; i++)
for(j = i + 1; j < *n - 1; j++) {
A = D[i][j];
for(k = j + 1; k < *n; k++) {
B = D[i][k];
C = D[j][k];
if((A <= B) && (A <= C)) {
delta = 2 * (B - C);
G[i][k] += delta;
G[j][k] -= delta;
}
else if(B <= C) {
delta = 2 * (C - A);
G[j][k] += delta;
G[i][j] -= delta;
}
else {
delta = 2 * (A - B);
G[i][j] += delta;
G[i][k] -= delta;
}
}
}
for(i = 0; i < *n; i++)
for(j = 0; j < *n; j++)
*out++ = G[i][j];
}isort3(int *i, int *j, int *k)
{
iwork3[0] = *i;
iwork3[1] = *j;
iwork3[2] = *k;
R_isort(iwork3, 3);
*i = iwork3[0];
*j = iwork3[1];
*k = iwork3[2];
}isort4(int *i, int *j, int *k, int *l)
{
iwork4[0] = *i;
iwork4[1] = *j;
iwork4[2] = *k;
iwork4[3] = *l;
R_isort(iwork4, 4);
*i = iwork4[0];
*j = iwork4[1];
*k = iwork4[2];
*l = iwork4[3];
}ls_fit_addtree_by_iterative_projection(double *d, int *n, int *order,
int *maxiter, int *iter,
double *tol, int *verbose)
{
double A, B, C, **D, DQ, delta;
int i, i1, j, j1, k, k1, l, l1;
D = clue_vector_to_square_matrix(d, *n);
for(*iter = 0; *iter < *maxiter; (*iter)++) {
delta = 0;
if(*verbose) Rprintf("Iteration: %d, ", *iter);
for(i1 = 0; i1 < *n - 3; i1++)
for(j1 = i1 + 1; j1 < *n - 2; j1++)
for(k1 = j1 + 1; k1 < *n - 1; k1++)
for(l1 = k1 + 1; l1 < *n; l1++) {
i = order[i1];
j = order[j1];
k = order[k1];
l = order[l1];
isort4(&i, &j, &k, &l);
A = D[i][j] + D[k][l];
B = D[i][k] + D[j][l];
C = D[i][l] + D[j][k];
if((A <= B) && (A <= C)) {
DQ = (C - B) / 4;
D[i][l] -= DQ;
D[j][k] -= DQ;
D[i][k] += DQ;
D[j][l] += DQ;
delta += fabs(C - B);
}
else if(B <= C) {
DQ = (A - C) / 4;
D[i][l] += DQ;
D[j][k] += DQ;
D[i][j] -= DQ;
D[k][l] -= DQ;
delta += fabs(A - C);
}
else {
DQ = (B - A) / 4;
D[i][k] -= DQ;
D[j][l] -= DQ;
D[i][j] += DQ;
D[k][l] += DQ;
delta += fabs(B - A);
}
}
if(*verbose) Rprintf("change: %f\n", delta);
if(delta < *tol)
break;
}
for(i = 0; i < *n - 1; i++)
for(j = i + 1; j < *n; j++)
D[j][i] = D[i][j];
/* And now write results back.
Could make this more efficient, of course ...
*/
for(j = 0; j < *n; j++)
for(i = 0; i < *n; i++)
*d++ = D[i][j];
}ls_fit_addtree_by_iterative_reduction(double *d, int *n, int *order,
int *maxiter, int *iter,
double *tol, int *verbose)
{
/* Once we have ls_fit_ultrametric_by_iterative_reduction() we can
always do this as well ...
See page 67f in Barthelemy and Guenoche.
*/
double A, B, C, **D, DQ, delta, tmp, N3;
int i, i1, j, j1, k, k1, l, l1;
D = clue_vector_to_square_matrix(d, *n);
/* And initialize the upper half of D ("work array") to 0.
(Yes, this could be done more efficiently by just propagating the
veclh dist representation.)
*/
for(i = 0; i < *n - 1; i++)
for(j = i + 1; j < *n; j++)
D[i][j] = 0;
N3 = (*n - 2) * (*n - 3) / 2;
for(*iter = 0; *iter < *maxiter; (*iter)++) {
if(*verbose) Rprintf("Iteration: %d, ", *iter);
for(i1 = 0; i1 < *n - 3; i1++)
for(j1 = i1 + 1; j1 < *n - 2; j1++)
for(k1 = j1 + 1; k1 < *n - 1; k1++)
for(l1 = k1 + 1; l1 < *n; l1++) {
i = order[i1];
j = order[j1];
k = order[k1];
l = order[l1];
isort4(&i, &j, &k, &l);
A = D[j][i] + D[l][k];
B = D[k][i] + D[l][j];
C = D[l][i] + D[k][j];
if((A <= B) && (A <= C)) {
/* Case 1: 5090 */
DQ = (C - B) / 4;
D[i][l] -= DQ;
D[j][k] -= DQ;
D[i][k] += DQ;
D[j][l] += DQ;
}
else if(B <= C) {
/* Case 2: 5120 */
DQ = (A - C) / 4;
D[i][l] += DQ;
D[j][k] += DQ;
D[i][j] -= DQ;
D[k][l] -= DQ;
}
else {
/* Case 3: 5150 */
DQ = (B - A) / 4;
D[i][k] -= DQ;
D[j][l] -= DQ;
D[i][j] += DQ;
D[k][l] += DQ;
}
}
delta = 0;
for(i = 0; i < *n - 1; i++)
for(j = i + 1; j < *n; j++) {
tmp = D[i][j] / N3;
D[j][i] += tmp;
D[i][j] = 0;
delta += fabs(tmp);
}
if(*verbose) Rprintf("change: %f\n", delta);
if(delta < *tol)
break;
}
/* And now write results back.
Could make this more efficient, of course ...
*/
for(j = 0; j < *n; j++)
for(i = 0; i < *n; i++)
*d++ = D[i][j];
}ls_fit_ultrametric_by_iterative_projection(double *d, int *n, int *order,
int *maxiter, int *iter,
double *tol, int *verbose)
{
double A, B, C, **D, delta;
int i, i1, j, j1, k, k1;
D = clue_vector_to_square_matrix(d, *n);
for(*iter = 0; *iter < *maxiter; (*iter)++) {
if(*verbose) Rprintf("Iteration: %d, ", *iter);
delta = 0;
for(i1 = 0; i1 < *n - 2; i1++)
for(j1 = i1 + 1; j1 < *n - 1; j1++)
for(k1 = j1 + 1; k1 < *n; k1++) {
i = order[i1];
j = order[j1];
k = order[k1];
isort3(&i, &j, &k);
A = D[i][j];
B = D[i][k];
C = D[j][k];
if((A <= B) && (A <= C)) {
D[i][k] = D[j][k] = (B + C) / 2;
delta += fabs(B - C);
}
else if(B <= C) {
D[i][j] = D[j][k] = (C + A) / 2;
delta += fabs(C - A);
}
else {
D[i][j] = D[i][k] = (A + B) / 2;
delta += fabs(A - B);
}
}
if(*verbose) Rprintf("change: %f\n", delta);
if(delta < *tol)
break;
}
for(i = 0; i < *n - 1; i++)
for(j = i + 1; j < *n; j++)
D[j][i] = D[i][j];
/* And now write results back.
Could make this more efficient, of course ...
*/
for(j = 0; j < *n; j++)
for(i = 0; i < *n; i++)
*d++ = D[i][j];
}ls_fit_ultrametric_by_iterative_reduction(double *d, int *n, int *order,
int *maxiter, int *iter,
double *tol, int *verbose)
{
double A, B, C, **D, DQ, delta, tmp;
int i, i1, j, j1, k, k1, N3;
D = clue_vector_to_square_matrix(d, *n);
/* And initialize the upper half of D ("work array") to 0.
(Yes, this could be done more efficiently by just propagating the
veclh dist representation.)
*/
for(i = 0; i < *n - 1; i++)
for(j = i + 1; j < *n; j++)
D[i][j] = 0;
N3 = (*n - 2);
for(*iter = 0; *iter < *maxiter; (*iter)++) {
if(*verbose) Rprintf("Iteration: %d, ", *iter);
for(i1 = 0; i1 < *n - 2; i1++)
for(j1 = i1 + 1; j1 < *n - 1; j1++)
for(k1 = j1 + 1; k1 < *n; k1++) {
i = order[i1];
j = order[j1];
k = order[k1];
isort3(&i, &j, &k);
A = D[j][i];
B = D[k][i];
C = D[k][j];
/*
<FIXME>
B & G have a divisor of 2 for case 1 and 4 for
cases 2 and 3 ... clearly, we should use the same
in all cases, but should it be 2 or 4?
</FIXME>
*/
if((A <= B) && (A <= C)) {
/* Case 1: 5080 */
DQ = (C - B) / 2;
D[i][k] += DQ;
D[j][k] -= DQ;
}
else if(B <= C) {
/* Case 2: 5100 */
DQ = (C - A) / 2;
D[i][j] += DQ;
D[j][k] -= DQ;
}
else {
/* Case 3: 5120 */
DQ = (B - A) / 2;
D[i][j] += DQ;
D[i][k] -= DQ;
}
}
delta = 0;
for(i = 0; i < *n - 1; i++)
for(j = i + 1; j < *n; j++) {
tmp = D[i][j] / N3;
D[j][i] += tmp;
D[i][j] = 0;
delta += fabs(tmp);
}
if(*verbose) Rprintf("change: %f\n", delta);
if(delta < *tol)
break;
}
/* And now write results back.
Could make this more efficient, of course ...
*/
for(j = 0; j < *n; j++)
for(i = 0; i < *n; i++)
*d++ = D[i][j];
}void preassign(AP *p)
{
int i, j, min, r, c, n, count;
int *ri, *ci, *rz, *cz;
n = p->n;
p->na = 0;
/* row and column markers */
ri = calloc(1 + n, sizeof(int));
ci = calloc(1 + n, sizeof(int));
/* row and column counts of zeroes */
rz = calloc(1 + n, sizeof(int));
cz = calloc(1 + n, sizeof(int));
for(i = 1; i <= n; i++){
count = 0;
for(j = 1; j <= n; j++)
if(p->c[i][j] == 0)
++count;
rz[i] = count;
}
for(i = 1; i <= n; i++){
count = 0;
for(j = 1; j <= n; j++)
if(p->c[j][i] == 0)
++count;
cz[i] = count;
}
while(TRUE){
/* find unassigned row with least number of zeroes > 0 */
min = INT_MAX;
r = 0;
for(i = 1; i <= n; i++)
if(rz[i] > 0 && rz[i] < min && ri[i] == UNASSIGNED){
min = rz[i];
r = i;
}
/* check if we are done */
if(r == 0)
break;
/* find unassigned column in row r with least number of zeroes */
c = 0;
min = INT_MAX;
for(i = 1; i <= n; i++)
if(p->c[r][i] == 0 && cz[i] < min && ci[i] == UNASSIGNED){
min = cz[i];
c = i;
}
if(c){
++p->na;
p->s[r] = c;
p->f[c] = r;
ri[r] = ASSIGNED;
ci[c] = ASSIGNED;
/* adjust zero counts */
cz[c] = 0;
for(i = 1; i <= n; i++)
if(p->c[i][c] == 0)
--rz[i];
}
}
/* free memory */
free(ri);
free(ci);
free(rz);
free(cz);
}void reduce(AP *p, int *ri, int *ci)
{
int i, j, n;
double min;
n = p->n;
/* find minimum in uncovered c-matrix */
min = DBL_MAX;
for(i = 1; i <= n; i++)
for(j = 1; j <= n; j++)
if(ri[i] == UNCOVERED && ci[j] == UNCOVERED){
if(p->c[i][j] < min)
min = p->c[i][j];
}
/* subtract min from each uncovered element and add it to each element */
/* which is covered twice */
for(i = 1; i <= n; i++)
for(j = 1; j <= n; j++){
if(ri[i] == UNCOVERED && ci[j] == UNCOVERED)
p->c[i][j]-= min;
if(ri[i] == COVERED && ci[j] == COVERED)
p->c[i][j]+= min;
}
}solve_LSAP(double *c, Sint *n, Sint *p)
{
AP *ap;
ap = ap_create_problem(c, *n);
ap_hungarian(ap);
ap_assignment(ap, p);
ap_free(ap);
}void ap_hungarian(AP *p)
{
int n; /* size of problem */
int *ri; /* covered rows */
int *ci; /* covered columns */
time_t start, end; /* timer */
int i, j, ok;
start = time(0);
n = p->n;
p->runs = 0;
/* allocate memory */
p->s = calloc(1 + n, sizeof(int));
p->f = calloc(1 + n, sizeof(int));
ri = calloc(1 + n, sizeof(int));
ci = calloc(1 + n, sizeof(int));
if(ri == NULL || ci == NULL || p->s == NULL || p->f == NULL)
error("ap_hungarian: could not allocate memory!");
preprocess(p);
preassign(p);
while(p->na < n){
if(REDUCE == cover(p, ri, ci))
reduce(p, ri, ci);
++p->runs;
}
end = time(0);
p->rtime = end - start;
/* check if assignment is a permutation of (1..n) */
for(i = 1; i <= n; i++){
ok = 0;
for(j = 1; j <= n; j++)
if(p->s[j] == i)
++ok;
if(ok != 1)
error("ap_hungarian: error in assigment, is not a permutation!");
}
/* calculate cost of assignment */
p->cost = 0;
for(i = 1; i <= n; i++)
p->cost+= p->C[i][p->s[i]];
/* reset result back to base-0 indexing */
for(i = 1; i <= n; i++)
p->s[i - 1] = p->s[i] - 1;
/* free memory */
free(ri);
free(ci);
}int cover(AP *p, int *ri, int *ci)
{
int *mr, i, r;
int n;
n = p->n;
mr = calloc(1 + p->n, sizeof(int));
/* reset cover indices */
for(i = 1; i <= n; i++){
if(p->s[i] == UNASSIGNED){
ri[i] = UNCOVERED;
mr[i] = MARKED;
}
else
ri[i] = COVERED;
ci[i] = UNCOVERED;
}
while(TRUE){
/* find marked row */
r = 0;
for(i = 1; i <= n; i++)
if(mr[i] == MARKED){
r = i;
break;
}
if(r == 0)
break;
for(i = 1; i <= n; i++)
if(p->c[r][i] == 0 && ci[i] == UNCOVERED){
if(p->f[i]){
ri[p->f[i]] = UNCOVERED;
mr[p->f[i]] = MARKED;
ci[i] = COVERED;
}else{
if(p->s[r] == UNASSIGNED)
++p->na;
p->f[p->s[r]] = 0;
p->f[i] = r;
p->s[r] = i;
free(mr);
return NOREDUCE;
}
}
mr[r] = UNMARKED;
}
free(mr);
return REDUCE;
}void preprocess(AP *p)
{
int i, j, n;
double min;
n = p->n;
/* subtract column minima in each row */
for(i = 1; i <= n; i++){
min = p->c[i][1];
for(j = 2; j <= n; j++)
if(p->c[i][j] < min)
min = p->c[i][j];
for(j = 1; j <= n; j++)
p->c[i][j]-= min;
}
/* subtract row minima in each column */
for(i = 1; i <= n; i++){
min = p->c[1][i];
for(j = 2; j <= n; j++)
if(p->c[j][i] < min)
min = p->c[j][i];
for(j = 1; j <= n; j++)
p->c[j][i]-= min;
}
}