Xplace_for_ICCAD/thirdparty/flute_mp/flute.cpp

1379 lines
37 KiB
C++

#include "flute.hpp"
namespace Flute {
#if DEGREE <= 7
#define MGROUP 5040 / 4 // Max. # of groups, 7! = 5040
#define MPOWV 15 // Max. # of POWVs per group
#elif DEGREE == 8
#define MGROUP 40320 / 4 // Max. # of groups, 8! = 40320
#define MPOWV 33 // Max. # of POWVs per group
#elif DEGREE == 9
#define MGROUP 362880 / 4 // Max. # of groups, 9! = 362880
#define MPOWV 79 // Max. # of POWVs per group
#endif
int numgrp[10] = {0, 0, 0, 0, 6, 30, 180, 1260, 10080, 90720};
struct Point {
DType x, y;
int o;
};
struct csoln {
unsigned char parent;
unsigned char seg[12]; // Add: 0..i, Sub: j..11; seg[i+1]=seg[j-1]=0
unsigned char row[DEGREE - 2], col[DEGREE - 2];
unsigned char neighbor[2 * DEGREE - 2];
};
struct csoln* LUT[DEGREE + 1][MGROUP]; // storing 4 .. D
int numsoln[DEGREE + 1][MGROUP];
Tree dmergetree(Tree t1, Tree t2);
Tree hmergetree(Tree t1, Tree t2, int* s);
Tree vmergetree(Tree t1, Tree t2);
void readLUT(const char* powv, const char* post) {
static bool LUTread = false;
if (LUTread) return;
LUTread = true;
FILE *fpwv, *fprt;
struct csoln* p;
int d, i, j, k, kk, ns, nn, ne;
unsigned char line[99], *linep, c;
unsigned char charnum[256];
unsigned char divd[256], modd[256], div16[256], mod16[256], dsq;
for (i = 0; i <= 255; i++) {
if ('0' <= i && i <= '9')
charnum[i] = i - '0';
else if (i >= 'A')
charnum[i] = i - 'A' + 10;
else // if (i=='$' || i=='\n' || ... )
charnum[i] = 0;
div16[i] = i / 16;
mod16[i] = i % 16;
}
fpwv = fopen(powv, "r");
if (fpwv == nullptr) {
printf("Error in opening %s\n", powv);
exit(1);
}
#if FLUTEROUTING == 1
fprt = fopen(post, "r");
if (fprt == nullptr) {
printf("Error in opening %s\n", post);
exit(1);
}
#endif
for (d = 4; d <= DEGREE; d++) {
for (i = 0; i <= 255; i++) {
divd[i] = i / d;
modd[i] = i % d;
}
dsq = d * d;
// Return value of fscanf/fread functions.
int ret = 0;
ret = fscanf(fpwv, "d=%d\n", &d);
#if FLUTEROUTING == 1
ret = fscanf(fprt, "d=%d\n", &d);
#endif
for (k = 0; k < numgrp[d]; k++) {
ns = (int)charnum[fgetc(fpwv)];
if (ns == 0) { // same as some previous group
ret = fscanf(fpwv, "%d\n", &kk);
numsoln[d][k] = numsoln[d][kk];
LUT[d][k] = LUT[d][kk];
} else {
fgetc(fpwv); // '\n'
numsoln[d][k] = ns;
p = (struct csoln*)malloc(ns * sizeof(struct csoln));
LUT[d][k] = p;
for (i = 1; i <= ns; i++) {
linep = (unsigned char*)fgets((char*)line, 99, fpwv);
p->parent = charnum[*(linep++)];
j = 0;
while ((p->seg[j++] = charnum[*(linep++)]) != 0)
;
j = 10;
while ((p->seg[j--] = charnum[*(linep++)]) != 0)
;
#if FLUTEROUTING == 1
nn = 2 * d - 2;
ne = 2 * d - 3;
ret = fread(line, 1, d - 2, fprt);
linep = line;
for (j = d; j < nn; j++) {
c = *(linep++);
if (c >= dsq) {
c -= dsq;
p->neighbor[divd[c]] = j;
ne--;
}
p->row[j - d] = divd[c];
p->col[j - d] = modd[c];
p->neighbor[j] = j; // initialized
}
ret = fread(line, 1, ne + 1, fprt);
linep = line; // last char = \n
for (j = 0; j < ne; j++) {
c = *(linep++);
p->neighbor[div16[c]] = mod16[c];
}
#endif
p++;
}
}
}
}
}
DType flute_wl(int d, DType* x, DType* y, int acc) {
auto allocateSize = std::max(MAXD, d + 1);
DType* xs = (DType*)malloc(sizeof(DType) * allocateSize);
DType* ys = (DType*)malloc(sizeof(DType) * allocateSize);
int* s = (int*)malloc(sizeof(int) * allocateSize);
Point* pt = (Point*)malloc(sizeof(Point) * allocateSize);
/* replaced the selection sort with stl stable_sort. mckim
Point* *ptp = (Point**) malloc(sizeof(Point*)*allocateSize); */
Point* tmpp;
DType minval, l, xu, xl, yu, yl;
int i, j, minidx;
if (d == 2)
l = std::abs(x[0] - x[1]) + std::abs(y[0] - y[1]);
else if (d == 3) {
if (x[0] > x[1]) {
xu = std::max(x[0], x[2]);
xl = std::min(x[1], x[2]);
} else {
xu = std::max(x[1], x[2]);
xl = std::min(x[0], x[2]);
}
if (y[0] > y[1]) {
yu = std::max(y[0], y[2]);
yl = std::min(y[1], y[2]);
} else {
yu = std::max(y[1], y[2]);
yl = std::min(y[0], y[2]);
}
l = (xu - xl) + (yu - yl);
} else {
for (i = 0; i < d; i++) {
pt[i].x = x[i];
pt[i].y = y[i];
}
// sort x
std::stable_sort(&pt[0], &pt[d],
[](const auto& pa, const auto& pb) {
return pa.x < pb.x;
});
#if REMOVE_DUPLICATE_PIN == 1
ptp[d] = &pt[d];
ptp[d]->x = ptp[d]->y = -999999;
j = 0;
for (i = 0; i < d; i++) {
for (k = i + 1; ptp[k]->x == ptp[i]->x; k++)
if (ptp[k]->y == ptp[i]->y) // pins k and i are the same
break;
if (ptp[k]->x != ptp[i]->x)
ptp[j++] = ptp[i];
}
d = j;
#endif
for (i = 0; i < d; i++) {
xs[i] = pt[i].x; // ptp[i]->x;
pt[i].o = i;
}
// sort y to find s[]
std::stable_sort(&pt[0], &pt[d],
[](const auto& pa, const auto& pb) {
return pa.y < pb.y;
});
for (i = 0; i < d; i++) {
ys[i] = pt[i].y;
s[i] = pt[i].o;
}
l = flutes_wl(d, xs, ys, s, acc);
}
free(xs);
free(ys);
free(s);
free(pt);
return l;
}
// xs[] and ys[] are coords in x and y in sorted order
// s[] is a list of nodes in increasing y direction
// if nodes are indexed in the order of increasing x coord
// i.e., s[i] = s_i in defined in paper
// The Points are (xs[s[i]], ys[i]) for i=0..d-1
// or (xs[i], ys[si[i]]) for i=0..d-1
DType flutes_wl_RDP(int d, DType* xs, DType* ys, int* s, int acc) {
int i, j, ss;
for (i = 0; i < d - 1; i++) {
if (xs[s[i]] == xs[s[i + 1]] && ys[i] == ys[i + 1]) {
if (s[i] < s[i + 1])
ss = s[i + 1];
else {
ss = s[i];
s[i] = s[i + 1];
}
for (j = i + 2; j < d; j++) {
ys[j - 1] = ys[j];
s[j - 1] = s[j];
}
for (j = ss + 1; j < d; j++)
xs[j - 1] = xs[j];
for (j = 0; j <= d - 2; j++)
if (s[j] > ss) s[j]--;
i--;
d--;
}
}
return flutes_wl_ALLD(d, xs, ys, s, acc);
}
// For low-degree, i.e., 2 <= d <= D
DType flutes_wl_LD(int d, DType* xs, DType* ys, int* s) {
int k, pi, i, j;
struct csoln* rlist;
DType dd[2 * DEGREE - 2]; // 0..D-2 for v, D-1..2*D-3 for h
DType minl, sum, l[MPOWV + 1];
if (d <= 3)
minl = xs[d - 1] - xs[0] + ys[d - 1] - ys[0];
else {
k = 0;
if (s[0] < s[2]) k++;
if (s[1] < s[2]) k++;
for (i = 3; i <= d - 1; i++) { // p0=0 always, skip i=1 for symmetry
pi = s[i];
for (j = d - 1; j > i; j--)
if (s[j] < s[i])
pi--;
k = pi + (i + 1) * k;
}
if (k < numgrp[d]) // no horizontal flip
for (i = 1; i <= d - 3; i++) {
dd[i] = ys[i + 1] - ys[i];
dd[d - 1 + i] = xs[i + 1] - xs[i];
}
else {
k = 2 * numgrp[d] - 1 - k;
for (i = 1; i <= d - 3; i++) {
dd[i] = ys[i + 1] - ys[i];
dd[d - 1 + i] = xs[d - 1 - i] - xs[d - 2 - i];
}
}
minl = l[0] = xs[d - 1] - xs[0] + ys[d - 1] - ys[0];
rlist = LUT[d][k];
for (i = 0; rlist->seg[i] > 0; i++)
minl += dd[rlist->seg[i]];
l[1] = minl;
j = 2;
while (j <= numsoln[d][k]) {
rlist++;
sum = l[rlist->parent];
for (i = 0; rlist->seg[i] > 0; i++)
sum += dd[rlist->seg[i]];
for (i = 10; rlist->seg[i] > 0; i--)
sum -= dd[rlist->seg[i]];
minl = std::min(minl, sum);
l[j++] = sum;
}
}
return minl;
}
// For medium-degree, i.e., D+1 <= d
DType flutes_wl_MD(int d, DType* xs, DType* ys, int* s, int acc) {
auto allocateSize = std::max(MAXD, d + 1);
DType* x1 = (DType*)malloc(sizeof(DType) * allocateSize);
DType* x2 = (DType*)malloc(sizeof(DType) * allocateSize);
DType* y1 = (DType*)malloc(sizeof(DType) * allocateSize);
DType* y2 = (DType*)malloc(sizeof(DType) * allocateSize);
int* si = (int*)malloc(sizeof(int) * allocateSize);
int* s1 = (int*)malloc(sizeof(int) * allocateSize);
int* s2 = (int*)malloc(sizeof(int) * allocateSize);
float* score = (float*)malloc(sizeof(float) * 2 * allocateSize);
float* penalty = (float*)malloc(sizeof(float) * allocateSize);
DType* distx = (DType*)malloc(sizeof(DType) * allocateSize);
DType* disty = (DType*)malloc(sizeof(DType) * allocateSize);
// DType x1[MAXD], x2[MAXD], y1[MAXD], y2[MAXD];
// int si[MAXD], s1[MAXD], s2[MAXD];
// float score[2*MAXD], penalty[MAXD];
// DType distx[MAXD], disty[MAXD];
float pnlty, dx, dy;
DType ll, minl, extral;
int i, r, p, maxbp, nbp, bp, ub, lb, n1, n2, ms, newacc;
int mins, maxs, minsi, maxsi;
DType xydiff;
DType retVal;
if (s[0] < s[d - 1]) {
ms = std::max(s[0], s[1]);
for (i = 2; i <= ms; i++)
ms = std::max(ms, s[i]);
if (ms <= d - 3) {
for (i = 0; i <= ms; i++) {
x1[i] = xs[i];
y1[i] = ys[i];
s1[i] = s[i];
}
x1[ms + 1] = xs[ms];
y1[ms + 1] = ys[ms];
s1[ms + 1] = ms + 1;
s2[0] = 0;
for (i = 1; i <= d - 1 - ms; i++)
s2[i] = s[i + ms] - ms;
retVal = 0;
retVal += flutes_wl_LMD(ms + 2, x1, y1, s1, acc);
retVal += flutes_wl_LMD(d - ms, xs + ms, ys + ms, s2, acc);
free(x1);
free(x2);
free(y1);
free(y2);
free(si);
free(s1);
free(s2);
free(score);
free(penalty);
free(distx);
free(disty);
return retVal;
}
} else { // (s[0] > s[d-1])
ms = std::min(s[0], s[1]);
for (i = 2; i <= d - 1 - ms; i++)
ms = std::min(ms, s[i]);
if (ms >= 2) {
x1[0] = xs[ms];
y1[0] = ys[0];
s1[0] = s[0] - ms + 1;
for (i = 1; i <= d - 1 - ms; i++) {
x1[i] = xs[i + ms - 1];
y1[i] = ys[i];
s1[i] = s[i] - ms + 1;
}
x1[d - ms] = xs[d - 1];
y1[d - ms] = ys[d - 1 - ms];
s1[d - ms] = 0;
s2[0] = ms;
for (i = 1; i <= ms; i++)
s2[i] = s[i + d - 1 - ms];
retVal = 0;
retVal += flutes_wl_LMD(d + 1 - ms, x1, y1, s1, acc);
retVal += flutes_wl_LMD(ms + 1, xs, ys + d - 1 - ms, s2, acc);
free(x1);
free(x2);
free(y1);
free(y2);
free(si);
free(s1);
free(s2);
free(score);
free(penalty);
free(distx);
free(disty);
return retVal;
}
}
// Find inverse si[] of s[]
for (r = 0; r < d; r++)
si[s[r]] = r;
// Determine breaking directions and positions dp[]
lb = 2;
ub = d - 1 - lb;
// Compute scores
#define AA 0.6
#define BB 0.3
#define CC 0.32
#define DD 4.8
float DDD = DD / (d - 1);
// Compute penalty[]
dx = CC * (xs[d - 2] - xs[1]) / (d - 3);
dy = CC * (ys[d - 2] - ys[1]) / (d - 3);
for (r = d / 2, pnlty = 0; r >= 0; r--, pnlty += dx)
penalty[r] = pnlty, penalty[d - 1 - r] = pnlty;
for (r = d / 2 - 1, pnlty = dy; r >= 0; r--, pnlty += dy)
penalty[s[r]] += pnlty, penalty[s[d - 1 - r]] += pnlty;
// #define CC 0.16
// for (r=0; r<d; r++)
// penalty[r] = abs(d-1-r-r)*dx + abs(d-1-si[r]-si[r])*dy;
// Compute distx[], disty[]
xydiff = (xs[d - 1] - xs[0]) - (ys[d - 1] - ys[0]);
if (s[0] < s[1])
mins = s[0], maxs = s[1];
else
mins = s[1], maxs = s[0];
if (si[0] < si[1])
minsi = si[0], maxsi = si[1];
else
minsi = si[1], maxsi = si[0];
for (r = 2; r <= ub; r++) {
if (s[r] < mins)
mins = s[r];
else if (s[r] > maxs)
maxs = s[r];
distx[r] = xs[maxs] - xs[mins];
if (si[r] < minsi)
minsi = si[r];
else if (si[r] > maxsi)
maxsi = si[r];
disty[r] = ys[maxsi] - ys[minsi] + xydiff;
}
if (s[d - 2] < s[d - 1])
mins = s[d - 2], maxs = s[d - 1];
else
mins = s[d - 1], maxs = s[d - 2];
if (si[d - 2] < si[d - 1])
minsi = si[d - 2], maxsi = si[d - 1];
else
minsi = si[d - 1], maxsi = si[d - 2];
for (r = d - 3; r >= lb; r--) {
if (s[r] < mins)
mins = s[r];
else if (s[r] > maxs)
maxs = s[r];
distx[r] += xs[maxs] - xs[mins];
if (si[r] < minsi)
minsi = si[r];
else if (si[r] > maxsi)
maxsi = si[r];
disty[r] += ys[maxsi] - ys[minsi];
}
nbp = 0;
for (r = lb; r <= ub; r++) {
if (si[r] == 0 || si[r] == d - 1)
score[nbp] = (xs[r + 1] - xs[r - 1]) - penalty[r] - AA * (ys[d - 2] - ys[1]) - DDD * disty[r];
else
score[nbp] = (xs[r + 1] - xs[r - 1]) - penalty[r] - BB * (ys[si[r] + 1] - ys[si[r] - 1]) - DDD * disty[r];
nbp++;
if (s[r] == 0 || s[r] == d - 1)
score[nbp] = (ys[r + 1] - ys[r - 1]) - penalty[s[r]] - AA * (xs[d - 2] - xs[1]) - DDD * distx[r];
else
score[nbp] = (ys[r + 1] - ys[r - 1]) - penalty[s[r]] - BB * (xs[s[r] + 1] - xs[s[r] - 1]) - DDD * distx[r];
nbp++;
}
if (acc <= 3)
newacc = 1;
else {
newacc = acc / 2;
if (acc >= nbp) acc = nbp - 1;
}
minl = (DType)INT_MAX;
for (i = 0; i < acc; i++) {
maxbp = 0;
for (bp = 1; bp < nbp; bp++)
if (score[maxbp] < score[bp]) maxbp = bp;
score[maxbp] = -9e9f;
#define BreakPt(bp) ((bp) / 2 + lb)
#define BreakInX(bp) ((bp) % 2 == 0)
p = BreakPt(maxbp);
// Breaking in p
if (BreakInX(maxbp)) { // break in x
n1 = n2 = 0;
for (r = 0; r < d; r++) {
if (s[r] < p) {
s1[n1] = s[r];
y1[n1] = ys[r];
n1++;
} else if (s[r] > p) {
s2[n2] = s[r] - p;
y2[n2] = ys[r];
n2++;
} else { // if (s[r] == p) i.e., r = si[p]
s1[n1] = p;
s2[n2] = 0;
if (r == d - 1 || r == d - 2) {
y1[n1] = y2[n2] = ys[r - 1];
extral = ys[r] - ys[r - 1];
}
if (r == 0 || r == 1) {
y1[n1] = y2[n2] = ys[r + 1];
extral = ys[r + 1] - ys[r];
} else {
y1[n1] = y2[n2] = ys[r];
extral = 0;
}
n1++;
n2++;
}
}
ll = extral + flutes_wl_LMD(p + 1, xs, y1, s1, newacc) + flutes_wl_LMD(d - p, xs + p, y2, s2, newacc);
} else { // if (!BreakInX(maxbp))
n1 = n2 = 0;
for (r = 0; r < d; r++) {
if (si[r] < p) {
s1[si[r]] = n1;
x1[n1] = xs[r];
n1++;
} else if (si[r] > p) {
s2[si[r] - p] = n2;
x2[n2] = xs[r];
n2++;
} else { // if (si[r] == p) i.e., r = s[p]
s1[p] = n1;
s2[0] = n2;
if (r == d - 1 || r == d - 2) {
x1[n1] = x2[n2] = xs[r - 1];
extral = xs[r] - xs[r - 1];
}
if (r == 0 || r == 1) {
x1[n1] = x2[n2] = xs[r + 1];
extral = xs[r + 1] - xs[r];
} else {
x1[n1] = x2[n2] = xs[r];
extral = 0;
}
n1++;
n2++;
}
}
ll = extral + flutes_wl_LMD(p + 1, x1, ys, s1, newacc) + flutes_wl_LMD(d - p, x2, ys + p, s2, newacc);
}
if (minl > ll) minl = ll;
}
free(x1);
free(x2);
free(y1);
free(y2);
free(si);
free(s1);
free(s2);
free(score);
free(penalty);
free(distx);
free(disty);
return minl;
}
Tree flute(int d, DType* x, DType* y, int acc) {
auto allocateSize = std::max(MAXD, d + 1);
DType* xs = (DType*)malloc(sizeof(DType) * allocateSize);
DType* ys = (DType*)malloc(sizeof(DType) * allocateSize);
int* s = (int*)malloc(sizeof(int) * allocateSize);
Point* pt = (Point*)malloc(sizeof(Point) * allocateSize);
Point** ptp = (Point**)malloc(sizeof(Point*) * allocateSize);
Point* tmpp;
DType minval;
int i, j, minidx;
Tree t;
if (d == 2) {
t.deg = 2;
t.length = std::abs(x[0] - x[1]) + std::abs(y[0] - y[1]);
t.branch = (Branch*)malloc(2 * sizeof(Branch));
t.branch[0].x = x[0];
t.branch[0].y = y[0];
t.branch[0].n = 1;
t.branch[1].x = x[1];
t.branch[1].y = y[1];
t.branch[1].n = 1;
} else {
for (i = 0; i < d; i++) {
pt[i].x = x[i];
pt[i].y = y[i];
ptp[i] = &pt[i];
}
// sort x
for (i = 0; i < d - 1; i++) {
minval = ptp[i]->x;
minidx = i;
for (j = i + 1; j < d; j++) {
if (minval > ptp[j]->x || (minval == ptp[j]->x && ptp[minidx]->y < ptp[j]->y)) {
minval = ptp[j]->x;
minidx = j;
}
}
tmpp = ptp[i];
ptp[i] = ptp[minidx];
ptp[minidx] = tmpp;
}
#if REMOVE_DUPLICATE_PIN == 1
ptp[d] = &pt[d];
ptp[d]->x = ptp[d]->y = -999999;
j = 0;
for (i = 0; i < d; i++) {
for (k = i + 1; ptp[k]->x == ptp[i]->x; k++)
if (ptp[k]->y == ptp[i]->y) // pins k and i are the same
break;
if (ptp[k]->x != ptp[i]->x)
ptp[j++] = ptp[i];
}
d = j;
#endif
for (i = 0; i < d; i++) {
xs[i] = ptp[i]->x;
ptp[i]->o = i;
}
// sort y to find s[]
for (i = 0; i < d - 1; i++) {
minval = ptp[i]->y;
minidx = i;
for (j = i + 1; j < d; j++) {
if (minval > ptp[j]->y) {
minval = ptp[j]->y;
minidx = j;
}
}
ys[i] = ptp[minidx]->y;
s[i] = ptp[minidx]->o;
ptp[minidx] = ptp[i];
}
ys[d - 1] = ptp[d - 1]->y;
s[d - 1] = ptp[d - 1]->o;
t = flutes(d, xs, ys, s, acc);
}
free(xs);
free(ys);
free(s);
free(pt);
free(ptp);
return t;
}
// xs[] and ys[] are coords in x and y in sorted order
// s[] is a list of nodes in increasing y direction
// if nodes are indexed in the order of increasing x coord
// i.e., s[i] = s_i in defined in paper
// The Points are (xs[s[i]], ys[i]) for i=0..d-1
// or (xs[i], ys[si[i]]) for i=0..d-1
Tree flutes_RDP(int d, DType* xs, DType* ys, int* s, int acc) {
int i, j, ss;
for (i = 0; i < d - 1; i++) {
if (xs[s[i]] == xs[s[i + 1]] && ys[i] == ys[i + 1]) {
if (s[i] < s[i + 1])
ss = s[i + 1];
else {
ss = s[i];
s[i] = s[i + 1];
}
for (j = i + 2; j < d; j++) {
ys[j - 1] = ys[j];
s[j - 1] = s[j];
}
for (j = ss + 1; j < d; j++)
xs[j - 1] = xs[j];
for (j = 0; j <= d - 2; j++)
if (s[j] > ss) s[j]--;
i--;
d--;
}
}
return flutes_ALLD(d, xs, ys, s, acc);
}
// For low-degree, i.e., 2 <= d <= D
Tree flutes_LD(int d, DType* xs, DType* ys, int* s) {
int k, pi, i, j;
struct csoln *rlist, *bestrlist;
DType dd[2 * DEGREE - 2]; // 0..D-2 for v, D-1..2*D-3 for h
DType minl, sum;
DType l[MPOWV + 1];
int hflip;
Tree t;
t.deg = d;
t.branch = (Branch*)malloc((2 * d - 2) * sizeof(Branch));
if (d == 2) {
minl = xs[1] - xs[0] + ys[1] - ys[0];
t.branch[0].x = xs[s[0]];
t.branch[0].y = ys[0];
t.branch[0].n = 1;
t.branch[1].x = xs[s[1]];
t.branch[1].y = ys[1];
t.branch[1].n = 1;
} else if (d == 3) {
minl = xs[2] - xs[0] + ys[2] - ys[0];
t.branch[0].x = xs[s[0]];
t.branch[0].y = ys[0];
t.branch[0].n = 3;
t.branch[1].x = xs[s[1]];
t.branch[1].y = ys[1];
t.branch[1].n = 3;
t.branch[2].x = xs[s[2]];
t.branch[2].y = ys[2];
t.branch[2].n = 3;
t.branch[3].x = xs[1];
t.branch[3].y = ys[1];
t.branch[3].n = 3;
} else {
k = 0;
if (s[0] < s[2]) k++;
if (s[1] < s[2]) k++;
for (i = 3; i <= d - 1; i++) { // p0=0 always, skip i=1 for symmetry
pi = s[i];
for (j = d - 1; j > i; j--)
if (s[j] < s[i])
pi--;
k = pi + (i + 1) * k;
}
if (k < numgrp[d]) { // no horizontal flip
hflip = 0;
for (i = 1; i <= d - 3; i++) {
dd[i] = ys[i + 1] - ys[i];
dd[d - 1 + i] = xs[i + 1] - xs[i];
}
} else {
hflip = 1;
k = 2 * numgrp[d] - 1 - k;
for (i = 1; i <= d - 3; i++) {
dd[i] = ys[i + 1] - ys[i];
dd[d - 1 + i] = xs[d - 1 - i] - xs[d - 2 - i];
}
}
minl = l[0] = xs[d - 1] - xs[0] + ys[d - 1] - ys[0];
rlist = LUT[d][k];
for (i = 0; rlist->seg[i] > 0; i++)
minl += dd[rlist->seg[i]];
bestrlist = rlist;
l[1] = minl;
j = 2;
while (j <= numsoln[d][k]) {
rlist++;
sum = l[rlist->parent];
for (i = 0; rlist->seg[i] > 0; i++)
sum += dd[rlist->seg[i]];
for (i = 10; rlist->seg[i] > 0; i--)
sum -= dd[rlist->seg[i]];
if (sum < minl) {
minl = sum;
bestrlist = rlist;
}
l[j++] = sum;
}
t.branch[0].x = xs[s[0]];
t.branch[0].y = ys[0];
t.branch[1].x = xs[s[1]];
t.branch[1].y = ys[1];
for (i = 2; i < d - 2; i++) {
t.branch[i].x = xs[s[i]];
t.branch[i].y = ys[i];
t.branch[i].n = bestrlist->neighbor[i];
}
t.branch[d - 2].x = xs[s[d - 2]];
t.branch[d - 2].y = ys[d - 2];
t.branch[d - 1].x = xs[s[d - 1]];
t.branch[d - 1].y = ys[d - 1];
if (hflip) {
if (s[1] < s[0]) {
t.branch[0].n = bestrlist->neighbor[1];
t.branch[1].n = bestrlist->neighbor[0];
} else {
t.branch[0].n = bestrlist->neighbor[0];
t.branch[1].n = bestrlist->neighbor[1];
}
if (s[d - 1] < s[d - 2]) {
t.branch[d - 2].n = bestrlist->neighbor[d - 1];
t.branch[d - 1].n = bestrlist->neighbor[d - 2];
} else {
t.branch[d - 2].n = bestrlist->neighbor[d - 2];
t.branch[d - 1].n = bestrlist->neighbor[d - 1];
}
for (i = d; i < 2 * d - 2; i++) {
t.branch[i].x = xs[d - 1 - bestrlist->col[i - d]];
t.branch[i].y = ys[bestrlist->row[i - d]];
t.branch[i].n = bestrlist->neighbor[i];
}
} else { // !hflip
if (s[0] < s[1]) {
t.branch[0].n = bestrlist->neighbor[1];
t.branch[1].n = bestrlist->neighbor[0];
} else {
t.branch[0].n = bestrlist->neighbor[0];
t.branch[1].n = bestrlist->neighbor[1];
}
if (s[d - 2] < s[d - 1]) {
t.branch[d - 2].n = bestrlist->neighbor[d - 1];
t.branch[d - 1].n = bestrlist->neighbor[d - 2];
} else {
t.branch[d - 2].n = bestrlist->neighbor[d - 2];
t.branch[d - 1].n = bestrlist->neighbor[d - 1];
}
for (i = d; i < 2 * d - 2; i++) {
t.branch[i].x = xs[bestrlist->col[i - d]];
t.branch[i].y = ys[bestrlist->row[i - d]];
t.branch[i].n = bestrlist->neighbor[i];
}
}
}
t.length = minl;
return t;
}
// For medium-degree, i.e., D+1 <= d <= D2
Tree flutes_MD(int d, DType* xs, DType* ys, int* s, int acc) {
auto allocateSize = std::max(MAXD, d + 1);
DType* x1 = (DType*)malloc(sizeof(DType) * allocateSize);
DType* x2 = (DType*)malloc(sizeof(DType) * allocateSize);
DType* y1 = (DType*)malloc(sizeof(DType) * allocateSize);
DType* y2 = (DType*)malloc(sizeof(DType) * allocateSize);
int* si = (int*)malloc(sizeof(int) * allocateSize);
int* s1 = (int*)malloc(sizeof(int) * allocateSize);
int* s2 = (int*)malloc(sizeof(int) * allocateSize);
float* score = (float*)malloc(sizeof(float) * 2 * allocateSize);
float* penalty = (float*)malloc(sizeof(float) * allocateSize);
DType* distx = (DType*)malloc(sizeof(DType) * allocateSize);
DType* disty = (DType*)malloc(sizeof(DType) * allocateSize);
// DType x1[MAXD], x2[MAXD], y1[MAXD], y2[MAXD];
// int si[MAXD], s1[MAXD], s2[MAXD];
// float score[2*MAXD], penalty[MAXD];
// DType distx[MAXD], disty[MAXD];
float pnlty, dx, dy;
DType ll, minl, coord1, coord2;
int i, r, p, maxbp, bestbp, bp, nbp, ub, lb, n1, n2, nn1, nn2, ms, newacc;
Tree t, t1, t2, bestt1, bestt2;
int mins, maxs, minsi, maxsi;
DType xydiff;
for (i = 0; i < allocateSize; ++i) {
distx[i] = disty[i] = 0;
}
if (s[0] < s[d - 1]) {
ms = std::max(s[0], s[1]);
for (i = 2; i <= ms; i++)
ms = std::max(ms, s[i]);
if (ms <= d - 3) {
for (i = 0; i <= ms; i++) {
x1[i] = xs[i];
y1[i] = ys[i];
s1[i] = s[i];
}
x1[ms + 1] = xs[ms];
y1[ms + 1] = ys[ms];
s1[ms + 1] = ms + 1;
s2[0] = 0;
for (i = 1; i <= d - 1 - ms; i++)
s2[i] = s[i + ms] - ms;
t1 = flutes_LMD(ms + 2, x1, y1, s1, acc);
t2 = flutes_LMD(d - ms, xs + ms, ys + ms, s2, acc);
t = dmergetree(t1, t2);
free(t1.branch);
free(t2.branch);
free(x1);
free(x2);
free(y1);
free(y2);
free(si);
free(s1);
free(s2);
free(score);
free(penalty);
free(distx);
free(disty);
return t;
}
} else { // (s[0] > s[d-1])
ms = std::min(s[0], s[1]);
for (i = 2; i <= d - 1 - ms; i++)
ms = std::min(ms, s[i]);
if (ms >= 2) {
x1[0] = xs[ms];
y1[0] = ys[0];
s1[0] = s[0] - ms + 1;
for (i = 1; i <= d - 1 - ms; i++) {
x1[i] = xs[i + ms - 1];
y1[i] = ys[i];
s1[i] = s[i] - ms + 1;
}
x1[d - ms] = xs[d - 1];
y1[d - ms] = ys[d - 1 - ms];
s1[d - ms] = 0;
s2[0] = ms;
for (i = 1; i <= ms; i++)
s2[i] = s[i + d - 1 - ms];
t1 = flutes_LMD(d + 1 - ms, x1, y1, s1, acc);
t2 = flutes_LMD(ms + 1, xs, ys + d - 1 - ms, s2, acc);
t = dmergetree(t1, t2);
free(t1.branch);
free(t2.branch);
free(x1);
free(x2);
free(y1);
free(y2);
free(si);
free(s1);
free(s2);
free(score);
free(penalty);
free(distx);
free(disty);
return t;
}
}
// Find inverse si[] of s[]
for (r = 0; r < d; r++)
si[s[r]] = r;
// Determine breaking directions and positions dp[]
lb = std::max((d - acc) / 5, 2);
ub = d - 1 - lb;
// Compute scores
#define AA 0.6 // 2.0*BB
#define BB 0.3
// #define CCC std::max(0.425-0.005*d-0.015*acc, 0.1)
// #define CCC std::max(0.43-0.005*d-0.01*acc, 0.1)
#define CCC std::max(0.41 - 0.005 * d, 0.1)
#define DD 4.8
float DDD = DD / (d - 1);
// Compute penalty[]
dx = CCC * (xs[d - 2] - xs[1]) / (d - 3);
dy = CCC * (ys[d - 2] - ys[1]) / (d - 3);
for (r = d / 2, pnlty = 0; r >= 2; r--, pnlty += dx)
penalty[r] = pnlty, penalty[d - 1 - r] = pnlty;
penalty[1] = pnlty, penalty[d - 2] = pnlty;
penalty[0] = pnlty, penalty[d - 1] = pnlty;
for (r = d / 2 - 1, pnlty = dy; r >= 2; r--, pnlty += dy)
penalty[s[r]] += pnlty, penalty[s[d - 1 - r]] += pnlty;
penalty[s[1]] += pnlty, penalty[s[d - 2]] += pnlty;
penalty[s[0]] += pnlty, penalty[s[d - 1]] += pnlty;
// #define CC 0.16
// #define v(r) ((r==0||r==1||r==d-2||r==d-1) ? d-3 : abs(d-1-r-r))
// for (r=0; r<d; r++)
// penalty[r] = v(r)*dx + v(si[r])*dy;
// Compute distx[], disty[]
xydiff = (xs[d - 1] - xs[0]) - (ys[d - 1] - ys[0]);
if (s[0] < s[1])
mins = s[0], maxs = s[1];
else
mins = s[1], maxs = s[0];
if (si[0] < si[1])
minsi = si[0], maxsi = si[1];
else
minsi = si[1], maxsi = si[0];
for (r = 2; r <= ub; r++) {
if (s[r] < mins)
mins = s[r];
else if (s[r] > maxs)
maxs = s[r];
distx[r] = xs[maxs] - xs[mins];
if (si[r] < minsi)
minsi = si[r];
else if (si[r] > maxsi)
maxsi = si[r];
disty[r] = ys[maxsi] - ys[minsi] + xydiff;
}
if (s[d - 2] < s[d - 1])
mins = s[d - 2], maxs = s[d - 1];
else
mins = s[d - 1], maxs = s[d - 2];
if (si[d - 2] < si[d - 1])
minsi = si[d - 2], maxsi = si[d - 1];
else
minsi = si[d - 1], maxsi = si[d - 2];
for (r = d - 3; r >= lb; r--) {
if (s[r] < mins)
mins = s[r];
else if (s[r] > maxs)
maxs = s[r];
distx[r] += xs[maxs] - xs[mins];
if (si[r] < minsi)
minsi = si[r];
else if (si[r] > maxsi)
maxsi = si[r];
disty[r] += ys[maxsi] - ys[minsi];
}
nbp = 0;
for (r = lb; r <= ub; r++) {
if (si[r] <= 1)
score[nbp] = (xs[r + 1] - xs[r - 1]) - penalty[r] - AA * (ys[2] - ys[1]) - DDD * disty[r];
else if (si[r] >= d - 2)
score[nbp] = (xs[r + 1] - xs[r - 1]) - penalty[r] - AA * (ys[d - 2] - ys[d - 3]) - DDD * disty[r];
else
score[nbp] = (xs[r + 1] - xs[r - 1]) - penalty[r] - BB * (ys[si[r] + 1] - ys[si[r] - 1]) - DDD * disty[r];
nbp++;
if (s[r] <= 1)
score[nbp] = (ys[r + 1] - ys[r - 1]) - penalty[s[r]] - AA * (xs[2] - xs[1]) - DDD * distx[r];
else if (s[r] >= d - 2)
score[nbp] = (ys[r + 1] - ys[r - 1]) - penalty[s[r]] - AA * (xs[d - 2] - xs[d - 3]) - DDD * distx[r];
else
score[nbp] = (ys[r + 1] - ys[r - 1]) - penalty[s[r]] - BB * (xs[s[r] + 1] - xs[s[r] - 1]) - DDD * distx[r];
nbp++;
}
if (acc <= 3)
newacc = 1;
else {
newacc = acc / 2;
if (acc >= nbp) acc = nbp - 1;
}
minl = (DType)INT_MAX;
bestt1.branch = bestt2.branch = nullptr;
for (i = 0; i < acc; i++) {
maxbp = 0;
for (bp = 1; bp < nbp; bp++)
if (score[maxbp] < score[bp]) maxbp = bp;
score[maxbp] = -9e9f;
#define BreakPt(bp) ((bp) / 2 + lb)
#define BreakInX(bp) ((bp) % 2 == 0)
p = BreakPt(maxbp);
// Breaking in p
if (BreakInX(maxbp)) { // break in x
n1 = n2 = 0;
for (r = 0; r < d; r++) {
if (s[r] < p) {
s1[n1] = s[r];
y1[n1] = ys[r];
n1++;
} else if (s[r] > p) {
s2[n2] = s[r] - p;
y2[n2] = ys[r];
n2++;
} else { // if (s[r] == p) i.e., r = si[p]
s1[n1] = p;
s2[n2] = 0;
y1[n1] = y2[n2] = ys[r];
nn1 = n1;
nn2 = n2;
n1++;
n2++;
}
}
t1 = flutes_LMD(p + 1, xs, y1, s1, newacc);
t2 = flutes_LMD(d - p, xs + p, y2, s2, newacc);
ll = t1.length + t2.length;
coord1 = t1.branch[t1.branch[nn1].n].y;
coord2 = t2.branch[t2.branch[nn2].n].y;
if (t2.branch[nn2].y > std::max(coord1, coord2))
ll -= t2.branch[nn2].y - std::max(coord1, coord2);
else if (t2.branch[nn2].y < std::min(coord1, coord2))
ll -= std::min(coord1, coord2) - t2.branch[nn2].y;
} else { // if (!BreakInX(maxbp))
n1 = n2 = 0;
for (r = 0; r < d; r++) {
if (si[r] < p) {
s1[si[r]] = n1;
x1[n1] = xs[r];
n1++;
} else if (si[r] > p) {
s2[si[r] - p] = n2;
x2[n2] = xs[r];
n2++;
} else { // if (si[r] == p) i.e., r = s[p]
s1[p] = n1;
s2[0] = n2;
x1[n1] = x2[n2] = xs[r];
n1++;
n2++;
}
}
t1 = flutes_LMD(p + 1, x1, ys, s1, newacc);
t2 = flutes_LMD(d - p, x2, ys + p, s2, newacc);
ll = t1.length + t2.length;
coord1 = t1.branch[t1.branch[p].n].x;
coord2 = t2.branch[t2.branch[0].n].x;
if (t2.branch[0].x > std::max(coord1, coord2))
ll -= t2.branch[0].x - std::max(coord1, coord2);
else if (t2.branch[0].x < std::min(coord1, coord2))
ll -= std::min(coord1, coord2) - t2.branch[0].x;
}
if (minl > ll) {
minl = ll;
free(bestt1.branch);
free(bestt2.branch);
bestt1 = t1;
bestt2 = t2;
bestbp = maxbp;
} else {
free(t1.branch);
free(t2.branch);
}
}
if (BreakInX(bestbp))
t = hmergetree(bestt1, bestt2, s);
else
t = vmergetree(bestt1, bestt2);
free(bestt1.branch);
free(bestt2.branch);
free(x1);
free(x2);
free(y1);
free(y2);
free(si);
free(s1);
free(s2);
free(score);
free(penalty);
free(distx);
free(disty);
return t;
}
Tree dmergetree(Tree t1, Tree t2) {
int i, d, prev, curr, next, offset1, offset2;
Tree t;
t.deg = d = t1.deg + t2.deg - 2;
t.length = t1.length + t2.length;
t.branch = (Branch*)malloc((2 * d - 2) * sizeof(Branch));
offset1 = t2.deg - 2;
offset2 = 2 * t1.deg - 4;
for (i = 0; i <= t1.deg - 2; i++) {
t.branch[i].x = t1.branch[i].x;
t.branch[i].y = t1.branch[i].y;
t.branch[i].n = t1.branch[i].n + offset1;
}
for (i = t1.deg - 1; i <= d - 1; i++) {
t.branch[i].x = t2.branch[i - t1.deg + 2].x;
t.branch[i].y = t2.branch[i - t1.deg + 2].y;
t.branch[i].n = t2.branch[i - t1.deg + 2].n + offset2;
}
for (i = d; i <= d + t1.deg - 3; i++) {
t.branch[i].x = t1.branch[i - offset1].x;
t.branch[i].y = t1.branch[i - offset1].y;
t.branch[i].n = t1.branch[i - offset1].n + offset1;
}
for (i = d + t1.deg - 2; i <= 2 * d - 3; i++) {
t.branch[i].x = t2.branch[i - offset2].x;
t.branch[i].y = t2.branch[i - offset2].y;
t.branch[i].n = t2.branch[i - offset2].n + offset2;
}
prev = t2.branch[0].n + offset2;
curr = t1.branch[t1.deg - 1].n + offset1;
next = t.branch[curr].n;
while (curr != next) {
t.branch[curr].n = prev;
prev = curr;
curr = next;
next = t.branch[curr].n;
}
t.branch[curr].n = prev;
return t;
}
Tree hmergetree(Tree t1, Tree t2, int* s) {
int i, prev, curr, next, extra, offset1, offset2;
int p, ii, n1, n2, nn1, nn2;
DType coord1, coord2;
Tree t;
t.deg = t1.deg + t2.deg - 1;
t.length = t1.length + t2.length;
t.branch = (Branch*)malloc((2 * t.deg - 2) * sizeof(Branch));
offset1 = t2.deg - 1;
offset2 = 2 * t1.deg - 3;
p = t1.deg - 1;
n1 = n2 = 0;
for (i = 0; i < t.deg; i++) {
if (s[i] < p) {
t.branch[i].x = t1.branch[n1].x;
t.branch[i].y = t1.branch[n1].y;
t.branch[i].n = t1.branch[n1].n + offset1;
n1++;
} else if (s[i] > p) {
t.branch[i].x = t2.branch[n2].x;
t.branch[i].y = t2.branch[n2].y;
t.branch[i].n = t2.branch[n2].n + offset2;
n2++;
} else {
t.branch[i].x = t2.branch[n2].x;
t.branch[i].y = t2.branch[n2].y;
t.branch[i].n = t2.branch[n2].n + offset2;
nn1 = n1;
nn2 = n2;
ii = i;
n1++;
n2++;
}
}
for (i = t.deg; i <= t.deg + t1.deg - 3; i++) {
t.branch[i].x = t1.branch[i - offset1].x;
t.branch[i].y = t1.branch[i - offset1].y;
t.branch[i].n = t1.branch[i - offset1].n + offset1;
}
for (i = t.deg + t1.deg - 2; i <= 2 * t.deg - 4; i++) {
t.branch[i].x = t2.branch[i - offset2].x;
t.branch[i].y = t2.branch[i - offset2].y;
t.branch[i].n = t2.branch[i - offset2].n + offset2;
}
extra = 2 * t.deg - 3;
coord1 = t1.branch[t1.branch[nn1].n].y;
coord2 = t2.branch[t2.branch[nn2].n].y;
if (t2.branch[nn2].y > std::max(coord1, coord2)) {
t.branch[extra].y = std::max(coord1, coord2);
t.length -= t2.branch[nn2].y - t.branch[extra].y;
} else if (t2.branch[nn2].y < std::min(coord1, coord2)) {
t.branch[extra].y = std::min(coord1, coord2);
t.length -= t.branch[extra].y - t2.branch[nn2].y;
} else
t.branch[extra].y = t2.branch[nn2].y;
t.branch[extra].x = t2.branch[nn2].x;
t.branch[extra].n = t.branch[ii].n;
t.branch[ii].n = extra;
prev = extra;
curr = t1.branch[nn1].n + offset1;
next = t.branch[curr].n;
while (curr != next) {
t.branch[curr].n = prev;
prev = curr;
curr = next;
next = t.branch[curr].n;
}
t.branch[curr].n = prev;
return t;
}
Tree vmergetree(Tree t1, Tree t2) {
int i, prev, curr, next, extra, offset1, offset2;
DType coord1, coord2;
Tree t;
t.deg = t1.deg + t2.deg - 1;
t.length = t1.length + t2.length;
t.branch = (Branch*)malloc((2 * t.deg - 2) * sizeof(Branch));
offset1 = t2.deg - 1;
offset2 = 2 * t1.deg - 3;
for (i = 0; i <= t1.deg - 2; i++) {
t.branch[i].x = t1.branch[i].x;
t.branch[i].y = t1.branch[i].y;
t.branch[i].n = t1.branch[i].n + offset1;
}
for (i = t1.deg - 1; i <= t.deg - 1; i++) {
t.branch[i].x = t2.branch[i - t1.deg + 1].x;
t.branch[i].y = t2.branch[i - t1.deg + 1].y;
t.branch[i].n = t2.branch[i - t1.deg + 1].n + offset2;
}
for (i = t.deg; i <= t.deg + t1.deg - 3; i++) {
t.branch[i].x = t1.branch[i - offset1].x;
t.branch[i].y = t1.branch[i - offset1].y;
t.branch[i].n = t1.branch[i - offset1].n + offset1;
}
for (i = t.deg + t1.deg - 2; i <= 2 * t.deg - 4; i++) {
t.branch[i].x = t2.branch[i - offset2].x;
t.branch[i].y = t2.branch[i - offset2].y;
t.branch[i].n = t2.branch[i - offset2].n + offset2;
}
extra = 2 * t.deg - 3;
coord1 = t1.branch[t1.branch[t1.deg - 1].n].x;
coord2 = t2.branch[t2.branch[0].n].x;
if (t2.branch[0].x > std::max(coord1, coord2)) {
t.branch[extra].x = std::max(coord1, coord2);
t.length -= t2.branch[0].x - t.branch[extra].x;
} else if (t2.branch[0].x < std::min(coord1, coord2)) {
t.branch[extra].x = std::min(coord1, coord2);
t.length -= t.branch[extra].x - t2.branch[0].x;
} else
t.branch[extra].x = t2.branch[0].x;
t.branch[extra].y = t2.branch[0].y;
t.branch[extra].n = t.branch[t1.deg - 1].n;
t.branch[t1.deg - 1].n = extra;
prev = extra;
curr = t1.branch[t1.deg - 1].n + offset1;
next = t.branch[curr].n;
while (curr != next) {
t.branch[curr].n = prev;
prev = curr;
curr = next;
next = t.branch[curr].n;
}
t.branch[curr].n = prev;
return t;
}
DType wirelength(Tree t) {
int i, j;
DType l = 0;
for (i = 0; i < 2 * t.deg - 2; i++) {
j = t.branch[i].n;
l += std::abs(t.branch[i].x - t.branch[j].x) +
std::abs(t.branch[i].y - t.branch[j].y);
}
return l;
}
void printtree(Tree t) {
int i;
for (i = 0; i < t.deg; i++)
printf(" %-2d: x=%4g y=%4g e=%d\n",
i, (float)t.branch[i].x, (float)t.branch[i].y, t.branch[i].n);
for (i = t.deg; i < 2 * t.deg - 2; i++)
printf("s%-2d: x=%4g y=%4g e=%d\n",
i, (float)t.branch[i].x, (float)t.branch[i].y, t.branch[i].n);
printf("\n");
}
} // namespace flute