Updated Flute 2.4 to Flute 3.1

* Theoretically slightly better
* Removed the artifical limitation to 150 terminals
This commit is contained in:
Gabriel Gouvine 2015-04-03 16:56:02 +02:00
parent 1aa416e82a
commit e97d4ca06a
21 changed files with 4377 additions and 77 deletions

View File

@ -2,7 +2,7 @@
# include ( ${QT_USE_FILE} )
include_directories ( ${KNIK_SOURCE_DIR}/src
${KNIK_SOURCE_DIR}/src/flute-2.4/src
${KNIK_SOURCE_DIR}/src/flute-3.1/src
${HURRICANE_INCLUDE_DIR}
${CORIOLIS_INCLUDE_DIR}
${UTILITIES_INCLUDE_DIR}
@ -39,13 +39,29 @@
KnikEngine.cpp
GraphicKnikEngine.cpp
)
set ( fluteIncludes flute-2.4/src/knik/flute.h )
set ( fluteCpps flute-2.4/src/flute.cpp )
set ( fluteIncludes flute-3.1/src/knik/flute.h
flute-3.1/src/knik/dl.h
flute-3.1/src/knik/mst2.h
flute-3.1/src/knik/err.h
flute-3.1/src/knik/heap.h
flute-3.1/src/knik/dist.h
flute-3.1/src/knik/global.h
flute-3.1/src/knik/neighbors.h
)
set ( fluteCpps flute-3.1/src/flute.cpp
flute-3.1/src/flute_mst.cpp
flute-3.1/src/dist.cpp
flute-3.1/src/dl.cpp
flute-3.1/src/err.cpp
flute-3.1/src/mst2.cpp
flute-3.1/src/heap.cpp
flute-3.1/src/neighbors.cpp
)
qtX_wrap_cpp ( mocCpps ${mocIncludes} )
add_library ( flute ${fluteCpps} )
set_target_properties ( flute PROPERTIES VERSION 2.4 SOVERSION 2 )
set_target_properties ( flute PROPERTIES VERSION 3.1 SOVERSION 3 )
target_link_libraries ( flute ${HURRICANE_LIBRARIES}
${CORIOLIS_LIBRARIES}
)
@ -71,5 +87,5 @@
install ( FILES ${includes}
${mocIncludes}
${fluteIncludes} DESTINATION include/coriolis2/knik )
install ( FILES flute-2.4/etc/POST9.dat
flute-2.4/etc/POWV9.dat DESTINATION share/coriolis2/flute-2.4 )
install ( FILES flute-3.1/etc/POST9.dat
flute-3.1/etc/POWV9.dat DESTINATION share/coriolis2/flute-3.1 )

View File

@ -39,7 +39,6 @@
#define __USE_MATRIXVERTEX__
#define EPSILON 10e-4
#define FLUTE_LIMIT 150
#define HISTORIC_INC 1.5 // define the increment of historic cost for ripup & reroute
namespace Knik {
@ -1547,13 +1546,6 @@ void Graph::UpdateEstimateCongestion ( bool create )
{
if ( _vertexes_to_route.size() < 2 )
return;
if ( _vertexes_to_route.size() >= FLUTE_LIMIT ) {
if ( create )
cerr << Warning( "Graph::UpdateEstimateCongestion(): %s\n"
" has more than %d vertex/terminals and cannot be handled by FLUTE."
, getString(_working_net).c_str(), FLUTE_LIMIT ) << endl;
return;
}
//cerr << "Running FLUTE for net : " << _working_net << endl;
auto_ptr<FTree> flutetree ( createFluteTree() );

View File

@ -0,0 +1,44 @@
#include "knik/global.h"
/*********************************************************************/
/*
Return the Manhattan distance between two points
*/
long dist(
Point p,
Point q
)
{
long dx, dy;
dx = (p.x) - (q.x);
if( dx < 0 ) dx = -dx;
dy = (p.y) - (q.y);
if( dy < 0 ) dy = -dy;
return dx + dy;
}
/*********************************************************************/
/*
Return the Manhattan distance between two points
*/
long dist2(
Point* p,
Point* q
)
{
long dx, dy;
dx = (p->x) - (q->x);
if( dx < 0 ) dx = -dx;
dy = (p->y) - (q->y);
if( dy < 0 ) dy = -dy;
return dx + dy;
}
/*********************************************************************/
/*********************************************************************/

View File

@ -0,0 +1,161 @@
#include "knik/dl.h"
#include <assert.h>
#include <stdio.h>
dl_t dl_alloc()
{
dl_t dl = (dl_t)malloc(sizeof(dl_s));
if (!dl) {
printf("Out of memory!!\n");
} else {
dl->first = dl->last = 0; dl->count = 0;
}
return dl;
}
void dl_delete(dl_t dl, dl_el *el)
{
if (dl->first == el) {
dl->first = el->next;
}
if (dl->last == el) {
dl->last = el->prev;
}
if (el->next) {
el->next->prev = el->prev;
}
if (el->prev) {
el->prev->next = el->next;
}
free(el); dl->count--;
}
void dl_clear(dl_t dl)
{
dl_el *el, *next;
if (dl->count > 0) {
for (el=dl->first; el; el=next) {
next = el->next;
free(el);
}
}
dl->first = dl->last = 0;
dl->count = 0;
}
void dl_concat(dl_t first_list, dl_t second_list)
{
if (first_list->count <= 0) {
*first_list = *second_list;
} else if (second_list->count > 0) {
first_list->last->next = second_list->first;
second_list->first->prev = first_list->last;
first_list->last = second_list->last;
first_list->count += second_list->count;
}
free(second_list);
}
static void dl_insertion_sort(dl_t dl, size_t el_size,
int(*compar)(void *, void *))
{
char *buf;
void *curr_d, *srch_d;
dl_el *curr, *srch;
if (dl_length(dl) <= 1) {
return;
}
buf = (char*)malloc(el_size);
for (curr=dl->first; curr!=dl->last; curr=curr->next) {
curr_d = (void*)(((dl_el*)curr)+1);
for (srch=dl->last; srch!=curr; srch=srch->prev) {
srch_d = (void*)(((dl_el*)srch)+1);
if (compar(curr_d, srch_d) > 0) {
memcpy((void*)buf, curr_d, el_size);
memcpy(curr_d, srch_d, el_size);
memcpy(srch_d, (void*)buf, el_size);
}
}
}
free(buf);
}
void dl_sort(dl_t dl, size_t el_size, int(*compar)(void *, void *))
{
dl_el *el, *first_head, *second_head;
dl_s first_list, second_list;
void *first_item, *second_item;
int i, len;
if (dl_length(dl) <= 25) {
dl_insertion_sort(dl, el_size, compar);
return;
}
len = dl_length(dl)/2;
for (i=0, el=dl->first; i<len; i++) {
el = el->next;
}
first_list.first = dl->first;
first_list.last = el->prev;
first_list.count = len;
first_list.last->next = 0;
second_list.first = el;
second_list.last = dl->last;
second_list.count = dl_length(dl)-len;
second_list.first->prev = 0;
dl_sort(&first_list, el_size, compar);
dl_sort(&second_list, el_size, compar);
/* in-place merging */
first_head = first_list.first;
second_head = second_list.first;
first_item = (void*)(((dl_el*)first_head)+1);
second_item = (void*)(((dl_el*)second_head)+1);
if (compar(first_item, second_item) <= 0) {
dl->first = el = first_head;
first_head = first_head->next;
} else {
dl->first = el = second_head;
second_head = second_head->next;
}
while (1) {
first_item = (void*)(((dl_el*)first_head)+1);
second_item = (void*)(((dl_el*)second_head)+1);
if (compar(first_item, second_item) <= 0) {
el->next = first_head;
first_head->prev = el;
el = first_head;
first_head = first_head->next;
if (!first_head) {
el->next = second_head;
second_head->prev = el;
dl->last = second_list.last;
break;
}
} else {
el->next = second_head;
second_head->prev = el;
el = second_head;
second_head = second_head->next;
if (!second_head) {
el->next = first_head;
first_head->prev = el;
dl->last = first_list.last;
break;
}
}
}
}

View File

@ -0,0 +1,28 @@
#include <stdio.h>
#include <stdlib.h>
/**************************************************************************/
/*
print error message and continue
*/
void err_msg(
char* msg
)
{
fprintf(stderr, "%s\n", msg);
}
/**************************************************************************/
/*
print error message and exit
*/
void err_exit(
char* msg
)
{
fprintf(stderr, "%s\n", msg);
exit(1);
}

View File

@ -17,11 +17,6 @@ using CRL::AllianceFramework;
#include <math.h>
#include "knik/flute.h"
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define abs(x) ((x)<0?(-x):(x))
#define ADIFF(x,y) ((x)>(y)?(x-y):(y-x)) // Absolute difference
#if D<=7
#define MGROUP 5040/4 // Max. # of groups, 7! = 5040
#define MPOWV 15 // Max. # of POWVs per group
@ -37,13 +32,19 @@ int numgrp[10]={0,0,0,0,6,30,180,1260,10080,90720};
struct csoln
{
unsigned char parent;
unsigned char seg[11]; // add: 0..i, Sub: j..10; seg[i+1]=seg[j-1]=0
unsigned char seg[11]; // Add: 0..i, Sub: j..10; seg[i+1]=seg[j-1]=0
unsigned char rowcol[D-2]; // row = rowcol[]/16, col = rowcol[]%16,
unsigned char neighbor[2*D-2];
};
struct csoln *LUT[D+1][MGROUP]; // storing 4 .. D
int numsoln[D+1][MGROUP];
struct point
{
DTYPE x, y;
int o;
};
void readLUT();
DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc);
DTYPE flutes_wl_LD(int d, DTYPE xs[], DTYPE ys[], int s[]);
@ -56,10 +57,12 @@ FTree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
FTree dmergetree(FTree t1, FTree t2);
FTree hmergetree(FTree t1, FTree t2, int s[]);
FTree vmergetree(FTree t1, FTree t2);
void local_refinement(FTree *tp, int p);
DTYPE wirelength(FTree t);
void printtree(FTree t);
void plottree(FTree t);
void readLUT()
{
unsigned char charnum[256], line[32], *linep, c;
@ -151,15 +154,13 @@ void readLUT()
fclose ( fprt );
}
DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc)
{
DTYPE xs[MAXD], ys[MAXD], minval, l, xu, xl, yu, yl;
int s[MAXD];
int i, j, /*k,*/ minidx;
struct point {
DTYPE x, y;
int o;
} pt[MAXD], *ptp[MAXD], *tmpp;
int i, j, k, minidx;
struct point pt[MAXD], *ptp[MAXD], *tmpp;
if (d==2)
l = ADIFF(x[0], x[1]) + ADIFF(y[0], y[1]);
@ -248,7 +249,7 @@ DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc)
// 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
// i.e., s[i] = s_i as 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
@ -405,13 +406,13 @@ DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
if (lb < 2) lb = 2;
ub=d-1-lb;
// compute scores
// Compute scores
#define AAWL 0.6
#define BBWL 0.3
float CCWL = 7.4/((d+10.)*(d-3.));
float DDWL = 4.8/(d-1);
// compute penalty[]
// Compute penalty[]
dx = CCWL*(xs[d-2]-xs[1]);
dy = CCWL*(ys[d-2]-ys[1]);
for (r = d/2, pnlty = 0; r>=0; r--, pnlty += dx)
@ -422,7 +423,7 @@ DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
// 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[]
// 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];
@ -568,15 +569,36 @@ DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
return minl;
}
static int orderx(const void *a, const void *b)
{
struct point *pa, *pb;
pa = *(struct point**)a;
pb = *(struct point**)b;
if (pa->x < pb->x) return -1;
if (pa->x > pb->x) return 1;
return 0;
}
static int ordery(const void *a, const void *b)
{
struct point *pa, *pb;
pa = *(struct point**)a;
pb = *(struct point**)b;
if (pa->y < pb->y) return -1;
if (pa->y > pb->y) return 1;
return 0;
}
FTree flute(int d, DTYPE x[], DTYPE y[], int acc)
{
DTYPE xs[MAXD], ys[MAXD], minval;
int s[MAXD];
int i, j, /*k,*/ minidx;
struct point {
DTYPE x, y;
int o;
} pt[MAXD], *ptp[MAXD], *tmpp;
DTYPE *xs, *ys, minval;
int *s;
int i, j, k, minidx;
struct point *pt, **ptp, *tmpp;
FTree t;
if (d==2) {
@ -591,27 +613,37 @@ FTree flute(int d, DTYPE x[], DTYPE y[], int acc)
t.branch[1].n = 1;
}
else {
xs = (DTYPE *)malloc(sizeof(DTYPE)*(d));
ys = (DTYPE *)malloc(sizeof(DTYPE)*(d));
s = (int *)malloc(sizeof(int)*(d));
pt = (struct point *)malloc(sizeof(struct point)*(d+1));
ptp = (struct point **)malloc(sizeof(struct point*)*(d+1));
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;
minidx = j;
if (d<200) {
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;
minidx = j;
}
}
tmpp = ptp[i];
ptp[i] = ptp[minidx];
ptp[minidx] = tmpp;
}
tmpp = ptp[i];
ptp[i] = ptp[minidx];
ptp[minidx] = tmpp;
} else {
qsort(ptp, d, sizeof(struct point *), orderx);
}
#if REMOVE_DUPLICATE_PIN==1
ptp[d] = &pt[d];
ptp[d]->x = ptp[d]->y = -999999;
@ -632,31 +664,46 @@ FTree flute(int d, DTYPE x[], DTYPE y[], int acc)
}
// 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;
if (d<200) {
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;
} else {
qsort(ptp, d, sizeof(struct point *), ordery);
for (i=0; i<d; i++) {
ys[i] = ptp[i]->y;
s[i] = ptp[i]->o;
}
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
// i.e., s[i] = s_i as 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
@ -917,13 +964,13 @@ FTree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
if (lb < 2) lb = 2;
ub=d-1-lb;
// compute scores
// Compute scores
#define AA 0.6 // 2.0*BB
#define BB 0.3
float CC = 7.4/((d+10.)*(d-3.));
float DD = 4.8/(d-1);
// compute penalty[]
// Compute penalty[]
dx = CC*(xs[d-2]-xs[1]);
dy = CC*(ys[d-2]-ys[1]);
for (r = d/2, pnlty = 0; r>=2; r--, pnlty += dx)
@ -939,7 +986,7 @@ FTree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
// for (r=0; r<d; r++)
// penalty[r] = v(r)*dx + v(si[r])*dy;
// compute distx[], disty[]
// 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];
@ -1096,9 +1143,22 @@ FTree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc)
}
}
if (BreakInX(bestbp))
#if LOCAL_REFINEMENT==1
if (BreakInX(bestbp)) {
t = hmergetree(bestt1, bestt2, s);
else t = vmergetree(bestt1, bestt2);
local_refinement(&t, si[BreakPt(bestbp)]);
} else {
t = vmergetree(bestt1, bestt2);
local_refinement(&t, BreakPt(bestbp));
}
#else
if (BreakInX(bestbp)) {
t = hmergetree(bestt1, bestt2, s);
} else {
t = vmergetree(bestt1, bestt2);
}
#endif
free(bestt1.branch);
free(bestt2.branch);
@ -1289,6 +1349,108 @@ FTree vmergetree(FTree t1, FTree t2)
return t;
}
void local_refinement(FTree *tp, int p)
{
int d, dd, i, ii, j, prev, curr, next, root;
int SteinerPin[2*MAXD], index[2*MAXD];
DTYPE x[MAXD], xs[D], ys[D];
int ss[D];
FTree tt;
d = tp->deg;
root = tp->branch[p].n;
// Reverse edges to point to root
prev = root;
curr = tp->branch[prev].n;
next = tp->branch[curr].n;
while (curr != next) {
tp->branch[curr].n = prev;
prev = curr;
curr = next;
next = tp->branch[curr].n;
}
tp->branch[curr].n = prev;
tp->branch[root].n = root;
// Find Steiner nodes that are at pins
for (i=d; i<=2*d-3; i++)
SteinerPin[i] = -1;
for (i=0; i<d; i++) {
next = tp->branch[i].n;
if (tp->branch[i].x == tp->branch[next].x &&
tp->branch[i].y == tp->branch[next].y)
SteinerPin[next] = i; // Steiner 'next' at Pin 'i'
}
SteinerPin[root] = p;
// Find pins that are directly connected to root
dd = 0;
for (i=0; i<d; i++) {
curr = tp->branch[i].n;
if (SteinerPin[curr] == i)
curr = tp->branch[curr].n;
while (SteinerPin[curr] < 0)
curr = tp->branch[curr].n;
if (curr == root) {
x[dd] = tp->branch[i].x;
if (SteinerPin[tp->branch[i].n] == i && tp->branch[i].n != root)
index[dd++] = tp->branch[i].n; // Steiner node
else index[dd++] = i; // Pin
}
}
if (4 <= dd && dd <= D) {
// Find Steiner nodes that are directly connected to root
ii=dd;
for (i=0; i<dd; i++) {
curr = tp->branch[index[i]].n;
while (SteinerPin[curr] < 0) {
index[ii++] = curr;
SteinerPin[curr] = INT_MAX;
curr = tp->branch[curr].n;
}
}
index[ii] = root;
for (ii=0; ii<dd; ii++) {
ss[ii] = 0;
for (j=0; j<ii; j++)
if (x[j] < x[ii])
ss[ii]++;
for (j=ii+1; j<dd; j++)
if (x[j] <= x[ii])
ss[ii]++;
xs[ss[ii]] = x[ii];
ys[ii] = tp->branch[index[ii]].y;
}
tt = flutes_LD(dd, xs, ys, ss);
// Find new wirelength
tp->length += tt.length;
for (ii=0; ii<2*dd-3; ii++) {
i = index[ii];
j = tp->branch[i].n;
tp->length -= ADIFF(tp->branch[i].x, tp->branch[j].x)
+ ADIFF(tp->branch[i].y, tp->branch[j].y);
}
// Copy tt into t
for (ii=0; ii<dd; ii++) {
tp->branch[index[ii]].n = index[tt.branch[ii].n];
}
for (; ii<=2*dd-3; ii++) {
tp->branch[index[ii]].x = tt.branch[ii].x;
tp->branch[index[ii]].y = tt.branch[ii].y;
tp->branch[index[ii]].n = index[tt.branch[ii].n];
}
free(tt.branch);
}
return;
}
DTYPE wirelength(FTree t)
{
int i, j;

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,177 @@
/****************************************************************************/
/*
Binary heap routines for use in Prim's algorithm,
with points are numbered from 0 to n-1
*/
#include <stdlib.h>
#include "knik/heap.h"
#include "knik/err.h"
Heap* _heap = (Heap*)NULL;
long _max_heap_size = 0;
long _heap_size = 0;
/****************************************************************************/
/*
*/
void allocate_heap( long n )
{
if( _max_heap_size < n )
{
_heap = (Heap*)realloc( (void*)_heap, (size_t)(n+1)*sizeof(Heap) );
if( ! _heap )
{
err_exit( "Cannot reallocate memory in allocate_heap!" );
}
_max_heap_size = n;
}
}
/****************************************************************************/
/*
*/
void deallocate_heap()
{
_max_heap_size = 0;
if( _heap )
{
free( (void*)_heap );
_heap = (Heap*)NULL;
}
}
/****************************************************************************/
void heap_init( long n )
{
register long p;
allocate_heap( n );
_heap_size = 0;
for( p = 0; p < n; p++ )
{
heap_idx( p ) = 0;
}
} /* END heap_init() */
/****************************************************************************/
void heap_insert(
long p,
long key
)
{
register long k; /* hole in the heap */
register long j; /* parent of the hole */
register long q; /* heap_elt(j) */
heap_key( p ) = key;
if( _heap_size == 0 )
{
_heap_size = 1;
heap_elt( 1 ) = p;
heap_idx( p ) = 1;
return;
}
k = ++ _heap_size;
j = k >> 1; /* k/2 */
while( (j > 0) && (heap_key(q=heap_elt(j)) > key) ) {
heap_elt( k ) = q;
heap_idx( q ) = k;
k = j;
j = k>>1; /* k/2 */
}
/* store p in the position of the hole */
heap_elt( k ) = p;
heap_idx( p ) = k;
} /* END heap_insert() */
/****************************************************************************/
void heap_decrease_key
(
long p,
long new_key
)
{
register long k; /* hole in the heap */
register long j; /* parent of the hole */
register long q; /* heap_elt(j) */
heap_key( p ) = new_key;
k = heap_idx( p );
j = k >> 1; /* k/2 */
if( (j > 0) && (heap_key(q=heap_elt(j)) > new_key) ) { /* change is needed */
do {
heap_elt( k ) = q;
heap_idx( q ) = k;
k = j;
j = k>>1; /* k/2 */
} while( (j > 0) && (heap_key(q=heap_elt(j)) > new_key) );
/* store p in the position of the hole */
heap_elt( k ) = p;
heap_idx( p ) = k;
}
} /* END heap_decrease_key() */
/****************************************************************************/
long heap_delete_min()
{
long min, last;
register long k; /* hole in the heap */
register long j; /* child of the hole */
register long l_key; /* key of last point */
if( _heap_size == 0 ) /* heap is empty */
return( -1 );
min = heap_elt( 1 );
last = heap_elt( _heap_size -- );
l_key = heap_key( last );
k = 1; j = 2;
while( j <= _heap_size ) {
if( heap_key(heap_elt(j)) > heap_key(heap_elt(j+1)) )
j++;
if( heap_key(heap_elt(j)) >= l_key)
break; /* found a position to insert 'last' */
/* else, sift hole down */
heap_elt(k) = heap_elt(j); /* Note that j <= _heap_size */
heap_idx( heap_elt(k) ) = k;
k = j;
j = k << 1;
}
heap_elt( k ) = last;
heap_idx( last ) = k;
heap_idx( min ) = -1; /* mark the point visited */
return( min );
} /* END heap_delete_min() */
/****************************************************************************/

View File

@ -0,0 +1,16 @@
#ifndef _DIST_H_
#define _DIST_H_
#include "global.h"
long dist(
Point p,
Point q
);
long dist2(
Point* p,
Point* q
);
#endif

View File

@ -0,0 +1,180 @@
#ifndef DL_H
#define DL_H
#include <string.h>
#include <stdlib.h>
typedef struct dl_el_s {
struct dl_el_s *prev, *next;
} dl_el;
typedef struct {
dl_el *first, *last;
unsigned int count;
} dl_s, *dl_t;
dl_t dl_alloc(void);
void dl_delete(dl_t dl, dl_el *el);
void dl_clear(dl_t dl);
void dl_concat(dl_t list1, dl_t list2);
void dl_sort(dl_t dl, size_t el_size, int(*compar)(void *, void *));
#define dl_length(dl) (dl)->count
#define dl_empty(dl) ((dl)->count <= 0)
#define dl_data(type, el) \
*(type*)(((dl_el*)(el))+1)
#define dl_data_p(type, el) \
((type*)(((dl_el*)(el))+1))
#define dl_forall(type, dl, data) \
{ \
dl_el *_el, *_next; \
dl_t _curr_dl = (dl); \
for (_el=_curr_dl->first; _el; _el=_next) { \
_next = _el->next; \
(data) = dl_data(type, _el);
#define dl_forall_p(type, dl, data_p) \
{ \
dl_el *_el, *_next; \
dl_t _curr_dl = (dl); \
for (_el=_curr_dl->first; _el; _el=_next) { \
_next = _el->next; \
(data_p) = dl_data_p(type, _el);
#define dl_current() _el
#define dl_delete_current() dl_delete(_curr_dl, _el)
#define dl_endfor \
} \
}
#define dl_forall_reverse(type, dl, data) \
{ \
dl_el *_el, *_next; \
dl_t _curr_dl = (dl); \
for (_el=_curr_dl->last; _el; _el=_next) { \
_next = _el->prev; \
(data) = dl_data(type, _el);
#define dl_forall_reverse_p(type, dl, data_p) \
{ \
dl_el *_el, *_next; \
dl_t _curr_dl = (dl); \
for (_el=_curr_dl->last; _el; _el=_next) { \
_next = _el->prev; \
(data_p) = dl_data_p(type, _el);
#define dl_first(type, dl) \
dl_data(type, (dl)->first)
#define dl_first_element(dl) (dl)->first
#define dl_last(type, dl) \
dl_data(type, (dl)->last)
#define dl_pop_first(type, dl, data) \
{ \
(data) = dl_first(type, dl); \
dl_delete((dl), (dl)->first); \
}
#define dl_pop_last(type, dl, data) \
{ (data) = dl_last(type, dl); dl_delete((dl), (dl)->last); }
#define dl_insert_before(type, dl, element, data) \
{ \
if ((element) == (dl)->first) { \
dl_prepend(type, dl, data); \
} else { \
dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \
if (!_el) { \
printf("Out of memory!!\n"); \
} else { \
memcpy(_el+1, &(data), sizeof(type)); \
_el->prev = (element)->prev; _el->next = (element); \
(element)->prev->next = _el; (element)->prev = _el; \
(dl)->count++; \
} \
} \
}
#define dl_insert_after(type, dl, element, data) \
{ \
if ((element) == (dl)->last) { \
dl_append(type, dl, data); \
} else { \
dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \
if (!_el) { \
printf("Out of memory!!\n"); \
} else { \
memcpy(_el+1, &(data), sizeof(type)); \
_el->next = (element)->next; _el->prev = (element); \
(element)->next->prev = _el; (element)->next = _el; \
(dl)->count++; \
} \
} \
}
#define dl_append(type, dl, data) \
{ \
dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \
if (!_el) { \
printf("Out of memory!!\n"); \
} else { \
memcpy(_el+1, &(data), sizeof(type)); \
_el->next = 0; \
if ((dl)->count <= 0) { \
_el->prev = 0; \
(dl)->first = (dl)->last = _el; \
(dl)->count = 1; \
} else { \
_el->prev = (dl)->last; \
(dl)->last->next = _el; \
(dl)->last = _el; \
(dl)->count++; \
} \
} \
}
#define dl_prepend(type, dl, data) \
{ \
dl_el *_el = (dl_el*) malloc(sizeof(dl_el)+sizeof(type)); \
if (!_el) { \
printf("Out of memory!!\n"); \
} else { \
memcpy(_el+1, &(data), sizeof(type)); \
_el->prev = 0; \
if ((dl)->count <= 0) { \
_el->next = 0; \
(dl)->first = (dl)->last = _el; \
(dl)->count = 1; \
} else { \
_el->next = (dl)->first; \
(dl)->first->prev = _el; \
(dl)->first = _el; \
(dl)->count++; \
} \
} \
}
#define dl_free(dl) \
{ \
dl_clear(dl); free(dl); dl = 0; \
}
#define dl_duplicate(dest, src, type) \
{ \
dest = dl_alloc(); \
type _data_el; \
dl_forall(type, src, _data_el) { \
dl_append(type, dest, _data_el); \
} dl_endfor; \
}
#endif

View File

@ -0,0 +1,12 @@
#ifndef _ERR_H_
#define _ERR_H_
void err_msg(
char* msg
);
void err_exit(
char* msg
);
#endif

View File

@ -1,19 +1,43 @@
#ifndef _KNIK_FLUTE_H
#define _KNIK_FLUTE_H
#define POWVFILE "POWV9.dat" // LUT for POWV (Wirelength Vector)
#define POSTFILE "POST9.dat" // LUT for POST (Steiner FTree)
#define MAXD 150 // max. degree of a net that can be handled
// setting MAXD to more than 150 is not recommended
#define D 9 // LUT is used for d <= D, D <= 9
#define ROUTING 1 // 1 to construct routing, 0 to estimate WL only
#define REMOVE_DUPLICATE_PIN 0 // remove dup. pin for flute_wl() & flute()
/*****************************/
/* User-Defined Parameters */
/*****************************/
#define MAXD 150 // max. degree that can be handled
#define ACCURACY 3 // Default accuracy
#define ROUTING 1 // 1 to construct routing, 0 to estimate WL only
#define LOCAL_REFINEMENT 0 // Suggestion: Set to 1 if ACCURACY >= 5
#define REMOVE_DUPLICATE_PIN 0 // Remove dup. pin for flute_wl() & flute()
#ifndef DTYPE // Data type for distance
#define DTYPE int
#endif
/*****************************/
/* User-Callable Functions */
/*****************************/
// void readLUT();
// DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc);
// DTYPE flutes_wl(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
// FTree flute(int d, DTYPE x[], DTYPE y[], int acc);
// FTree flutes(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
// DTYPE wirelength(FTree t);
// void printtree(FTree t);
// void plottree(FTree t);
/*************************************/
/* Internal Parameters and Functions */
/*************************************/
#define POWVFILE "POWV9.dat" // LUT for POWV (Wirelength Vector)
#define POSTFILE "POST9.dat" // LUT for POST (Steiner FTree)
#define D 9 // LUT is used for d <= D, D <= 9
#define TAU(A) (8+1.3*(A))
#define D1(A) (25+120/((A)*(A))) // flute_mr is used for D1 < d <= D2
#define D2(A) ((A)<=6 ? 500 : 75+5*(A))
typedef struct
{
DTYPE x, y; // starting point of the branch
@ -27,7 +51,7 @@ struct FTree
Branch *branch; // array of tree branches
};
// Major functions
// User-Callable Functions
extern void readLUT();
extern DTYPE flute_wl(int d, DTYPE x[], DTYPE y[], int acc);
//Macro: DTYPE flutes_wl(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
@ -37,13 +61,14 @@ extern DTYPE wirelength(FTree t);
extern void printtree(FTree t);
extern void plottree(FTree t);
// Other useful functions
extern void init_param();
extern DTYPE flutes_wl_LD(int d, DTYPE xs[], DTYPE ys[], int s[]);
extern DTYPE flutes_wl_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
extern DTYPE flutes_wl_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
extern FTree flutes_LD(int d, DTYPE xs[], DTYPE ys[], int s[]);
extern FTree flutes_MD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
extern FTree flutes_HD(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
extern FTree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
#if REMOVE_DUPLICATE_PIN==1
@ -55,12 +80,20 @@ extern FTree flutes_RDP(int d, DTYPE xs[], DTYPE ys[], int s[], int acc);
#endif
#define flutes_wl_ALLD(d, xs, ys, s, acc) flutes_wl_LMD(d, xs, ys, s, acc)
#define flutes_ALLD(d, xs, ys, s, acc) flutes_LMD(d, xs, ys, s, acc)
//#define flutes_ALLD(d, xs, ys, s, acc) (d<=D ? flutes_LD(d, xs, ys, s) : (d<=D2 ? flutes_MD(d, xs, ys, s, acc) : flutes_HD(d, xs, ys, s, acc)))
#define flutes_ALLD(d, xs, ys, s, acc) \
(d<=D ? flutes_LD(d, xs, ys, s) \
: (d<=D1(acc) ? flutes_MD(d, xs, ys, s, acc) \
: flutes_HD(d, xs, ys, s, acc)))
#define flutes_wl_LMD(d, xs, ys, s, acc) \
(d<=D ? flutes_wl_LD(d, xs, ys, s) : flutes_wl_MD(d, xs, ys, s, acc))
#define flutes_LMD(d, xs, ys, s, acc) \
(d<=D ? flutes_LD(d, xs, ys, s) : flutes_MD(d, xs, ys, s, acc))
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
#define abs(x) ((x)<0?(-x):(x))
#define ADIFF(x,y) ((x)>(y)?(x-y):(y-x)) // Absolute difference
#endif

View File

@ -0,0 +1,19 @@
#ifndef _GLOBAL_H_
#define _GLOBAL_H_
#include <stdio.h>
#define TRUE 1
#define FALSE 0
#define MAXLONG 0x7fffffffL
struct point
{
long x, y;
};
typedef struct point Point;
typedef long nn_array[8];
#endif /* _GLOBAL_H_ */

View File

@ -0,0 +1,31 @@
#ifndef _HEAP_H_
#define _HEAP_H_
#include "global.h"
struct heap_info
{
long key;
long idx;
long elt;
};
typedef struct heap_info Heap;
extern Heap* _heap;
#define heap_key( p ) ( _heap[p].key )
#define heap_idx( p ) ( _heap[p].idx )
#define heap_elt( k ) ( _heap[k].elt )
#define in_heap( p ) ( heap_idx(p) > 0 )
#define never_seen( p ) ( heap_idx(p) == 0 )
void allocate_heap( long n );
void deallocate_heap();
void heap_init( long n );
void heap_insert( long p, long key );
void heap_decrease_key( long p, long new_key );
long heap_delete_min();
#endif /* _HEAP_H_ */

View File

@ -0,0 +1,11 @@
#ifndef _MST2_H_
#define _MST2_H_
#include "global.h"
void mst2_package_init( long n );
void mst2_package_done();
void mst2( long n, Point* pt, long* parent );
#endif

View File

@ -0,0 +1,19 @@
#include "global.h"
void allocate_nn_arrays( long n );
void deallocate_nn_arrays();
void brute_force_nearest_neighbors
(
long n,
Point* pt,
nn_array* nn
);
void dq_nearest_neighbors
(
long n,
Point* pt,
nn_array* nn
);

View File

@ -0,0 +1,92 @@
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include "knik/global.h"
#include "knik/neighbors.h"
#include "knik/dist.h"
#include "knik/heap.h"
#include "knik/err.h"
void mst2_package_init( long n )
{
allocate_heap( n );
allocate_nn_arrays( n );
}
/****************************************************************************/
/*
*/
void mst2_package_done()
{
deallocate_heap();
deallocate_nn_arrays();
}
/****************************************************************************/
/*
*/
void mst2
(
long n,
Point* pt,
long* parent
)
{
long i, k, nn1;
long d;
long oct;
long root = 0;
extern nn_array* nn;
// brute_force_nearest_neighbors( n, pt, nn );
dq_nearest_neighbors( n, pt, nn );
/*
Binary heap implementation of Prim's algorithm.
Runs in O(n*log(n)) time since at most 8n edges are considered
*/
heap_init( n );
heap_insert( root, 0 );
parent[root] = root;
for( k = 0; k < n; k++ ) /* n points to be extracted from heap */
{
i = heap_delete_min();
if (i<0) break;
#ifdef DEBUG
assert( i >= 0 );
#endif
/*
pt[i] entered the tree, update heap keys for its neighbors
*/
for( oct = 0; oct < 8; oct++ )
{
nn1 = nn[i][oct];
if( nn1 >= 0 )
{
d = dist( pt[i], pt[nn1] );
if( in_heap(nn1) && (d < heap_key(nn1)) )
{
heap_decrease_key( nn1, d );
parent[nn1] = i;
}
else if( never_seen(nn1) )
{
heap_insert( nn1, d );
parent[nn1] = i;
}
}
}
}
}
/****************************************************************************/
/****************************************************************************/

View File

@ -0,0 +1,527 @@
#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include "knik/global.h"
#include "knik/err.h"
#include "knik/dist.h"
long octant
(
Point from,
Point to
);
static Point* _pt;
/***************************************************************************/
/*
For efficiency purposes auxiliary arrays are allocated as globals
*/
long max_arrays_size = 0;
nn_array* nn = (nn_array*)NULL;
Point* sheared = (Point*)NULL;
long* sorted = (long*)NULL;
long* aux = (long*)NULL;
/***************************************************************************/
/*
resize the auxiliary arrays to fit the specified number of points
*/
void allocate_nn_arrays( long n )
{
if( max_arrays_size < n )
{
nn = (nn_array*)realloc( (void*)nn, (size_t)n*sizeof(nn_array) );
sheared = (Point*)realloc( (void*)sheared, (size_t)n*sizeof(Point) );
sorted = (long*)realloc( (void*)sorted, (size_t)n*sizeof(long) );
aux = (long*)realloc( (void*)aux, (size_t)n*sizeof(long) );
if( !nn || !sheared || !sorted || !aux )
{
err_exit( "Cannot allocate memory in allocate_nn_arrays!" );
}
max_arrays_size = n;
}
}
/***************************************************************************/
/*
free memory used by auxiliary arrays
*/
void deallocate_nn_arrays()
{
max_arrays_size = 0;
if( nn )
{
free( (void*)nn );
nn = (nn_array*)NULL;
}
if( sheared )
{
free( (void*)sheared );
sheared = (Point*)NULL;
}
if( sorted )
{
free( (void*)sorted );
sorted = (long*)NULL;
}
if( aux )
{
free( (void*)aux );
aux = (long*)NULL;
}
}
/***************************************************************************/
/*
comparison function for use in quicksort
*/
static int compare_x
(
const void* i,
const void* j
)
{
/*
points with the same x must appear in increasing order of y
*/
if( sheared[*((long*)i)].x == sheared[*((long*)j)].x)
{
return sheared[*((long*)i)].y - sheared[*((long*)j)].y;
}
else
{
return sheared[*((long*)i)].x - sheared[*((long*)j)].x;
}
}
/***************************************************************************/
/*
Combine step of the Guibas-Stolfi divide-and-conquer NE nearest neighbor
algorithm. For efficiency purposes SW nearest neighbors are computed
at the same time.
*/
void ne_sw_combine
(
long left,
long mid,
long right,
Point* pt,
long* sorted,
long* aux,
long oct,
nn_array* nn
)
{
long i, j, k, y2;
long i1;
long i2;
long best_i2; /* index of current best nearest-neighbor */
long best_dist; /* distance to best nearest-neighbor */
long d;
#ifdef DEBUG
assert( right > mid );
assert( mid > left );
#endif
/*
update north-east nearest neighbors accross the mid-line
*/
i1 = left;
i2 = mid; y2 = pt[ sorted[i2] ].y;
while( (i1 < mid) && (pt[ sorted[i1] ].y >= y2) )
{
i1++;
}
if( i1 < mid )
{
best_i2 = i2;
best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] );
i2++;
while( (i1 < mid) && (i2 < right) )
{
if( pt[ sorted[i1] ].y < pt[ sorted[i2] ].y )
{
d = dist2( pt + sorted[i1], pt + sorted[i2] );
if( d < best_dist )
{
best_i2 = i2;
best_dist = d;
}
i2++;
}
else
{
if( (nn[ sorted[i1] ][oct] == -1) ||
( best_dist < dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) )
)
{
nn[ sorted[i1] ][oct] = sorted[best_i2];
}
i1++;
if( i1 < mid )
{
best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] );
}
}
}
while( i1 < mid )
{
if( (nn[ sorted[i1] ][oct] == -1) ||
( dist2( pt + sorted[i1], pt + sorted[best_i2] ) <
dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) )
)
{
nn[ sorted[i1] ][oct] = sorted[best_i2];
}
i1++;
}
}
/*
repeat for south-west nearest neighbors
*/
oct = (oct + 4) % 8;
i1 = right - 1;
i2 = mid - 1; y2 = pt[ sorted[i2] ].y;
while( (i1 >= mid) && (pt[ sorted[i1] ].y <= y2) )
{
i1--;
}
if( i1 >= mid )
{
best_i2 = i2;
best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] );
i2--;
while( (i1 >= mid) && (i2 >= left) )
{
if( pt[ sorted[i1] ].y > pt[ sorted[i2] ].y )
{
d = dist2( pt + sorted[i1], pt + sorted[i2] );
if( d < best_dist )
{
best_i2 = i2;
best_dist = d;
}
i2--;
}
else
{
if( (nn[ sorted[i1] ][oct] == -1) ||
( best_dist < dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) )
)
{
nn[ sorted[i1] ][oct] = sorted[best_i2];
}
i1--;
if( i1 >= mid )
{
best_dist = dist2( pt + sorted[i1], pt + sorted[best_i2] );
}
}
}
while( i1 >= mid )
{
if( (nn[ sorted[i1] ][oct] == -1) ||
( dist2( pt + sorted[i1], pt + sorted[best_i2] ) <
dist2( pt + sorted[i1], pt + nn[ sorted[i1] ][oct]) )
)
{
nn[ sorted[i1] ][oct] = sorted[best_i2];
}
i1--;
}
}
/*
merge sorted[left..mid-1] with sorted[mid..right-1] by y-coordinate
*/
i = left; /* first unprocessed element in left list */
j = mid; /* first unprocessed element in right list */
k = left; /* first free available slot in output list */
while( (i < mid) && (j < right) )
{
if( pt[ sorted[i] ].y >= pt[ sorted[j] ].y )
{
aux[k++] = sorted[i++];
}
else
{
aux[k++] = sorted[j++];
}
}
/*
copy leftovers
*/
while( i < mid ) { aux[k++] = sorted[i++]; }
while( j < right ) { aux[k++] = sorted[j++]; }
/*
now copy sorted points from 'aux' to 'sorted'
*/
for( i = left; i < right; i++ ) { sorted[i] = aux[i]; }
#if 0
memcpy( (void*)(sorted+left), /* destination */
(void*)(aux+left), /* source */
(size_t)(right-left)*sizeof(long) /* number of bytes */
);
#endif
}
/***************************************************************************/
/*
compute north-east and south-west nearest neighbors for points indexed
by {sorted[left],...,sorted[right-1]}
*/
void ne_sw_nearest_neighbors
(
long left,
long right,
Point* pt,
long* sorted,
long* aux,
long oct,
nn_array* nn
)
{
long mid;
#ifdef DEBUG
assert( right > left );
#endif
if( right == left + 1 )
{
nn[ sorted[left] ][oct] = nn[ sorted[left]][(oct+4) % 8] = -1;
}
else
{
mid = (left + right) / 2;
ne_sw_nearest_neighbors( left, mid, pt, sorted, aux, oct, nn );
ne_sw_nearest_neighbors( mid, right, pt, sorted, aux, oct, nn );
ne_sw_combine( left, mid, right, pt, sorted, aux, oct, nn );
}
}
/***************************************************************************/
/*
Guibas-Stolfi algorithm for computing nearest NE neighbors
*/
void dq_nearest_neighbors
(
long n,
Point* pt,
nn_array* nn
)
{
long i, oct;
void check_nn( long, Point*, nn_array* );
long shear[4][4] = {
{1, -1, 0, 2},
{2, 0, -1, 1},
{1, 1, -2, 0},
{0, 2, -1, -1}
};
_pt = pt;
for( oct = 0; oct < 4; oct++ )
{
for( i = 0; i < n; i++ )
{
sheared[i].x = shear[oct][0]*pt[i].x + shear[oct][1]*pt[i].y;
sheared[i].y = shear[oct][2]*pt[i].x + shear[oct][3]*pt[i].y;
sorted[i] = i;
}
qsort( sorted, n, sizeof(long), compare_x );
ne_sw_nearest_neighbors( 0, n, sheared, sorted, aux, oct, nn );
}
#ifdef DEBUG
check_nn( n, pt, nn );
#endif
}
/***************************************************************************/
/***************************************************************************/
/*
Brute-force nearest-neighbor computation for debugging purposes
*/
/***************************************************************************/
/*
Half-open octants are numbered from 0 to 7 in anti-clockwise order
starting from ( dx >= dy > 0 ).
*/
#define sgn(x) ( x>0 ? 1 : (x < 0 ? -1 : 0) )
long octant
(
Point from,
Point to
)
{
long dx = to.x - from.x;
long dy = to.y - from.y;
long sgn1 = sgn(dx)*sgn(dy);
long sgn2 = sgn(dx+dy)*sgn(dx-dy);
long oct = 0x0;
if( (dy < 0) || ((dy==0) && (dx>0)) ) oct += 4;
if( (sgn1 < 0) || (dy==0) ) oct += 2;
if( (sgn1*sgn2 < 0) || (dy==0) || (dx==0) ) oct += 1;
return oct;
}
/***************************************************************************/
/*
O(n^2) algorithm for computing all nearest neighbors
*/
void brute_force_nearest_neighbors
(
long n,
Point* pt,
nn_array* nn
)
{
long i, j, oct;
long d;
/*
compute nearest neighbors by inspecting all pairs of points
*/
for( i = 0; i < n; i++ )
{
for( oct = 0; oct < 8; oct++ )
{
nn[i][oct] = -1;
}
}
for( i = 0; i < n; i++ )
{
for( j = i+1; j < n; j++ )
{
d = dist(pt[i], pt[j]);
oct = octant( pt[i], pt[j] );
if( ( nn[i][oct] == -1 ) ||
( d < dist(pt[i], pt[ nn[i][oct] ]) )
)
{
nn[i][oct] = j;
}
oct = (oct + 4) % 8;
if( ( nn[j][oct] == -1 ) ||
( d < dist(pt[j], pt[ nn[j][oct] ]) )
)
{
nn[j][oct] = i;
}
}
}
}
/***************************************************************************/
/*
compare nearest neighbors against those computed by brute force
*/
void check_nn
(
long n,
Point* pt,
nn_array* nn
)
{
long i, j, oct;
nn_array* nn1;
nn1 = (nn_array*)calloc( (size_t)n, (size_t)sizeof(nn_array) );
brute_force_nearest_neighbors( n, pt, nn1 );
for( i = 0; i < n; i++ )
{
for( oct = 0; oct < 8; oct++ )
{
if( nn[i][oct] == -1 )
{
assert( nn1[i][oct] == -1 );
}
else
{
assert( nn1[i][oct] != -1 );
if( octant(pt[i], pt[ nn[i][oct] ]) != oct )
{
printf( "WRONG OCTANT!\noct=%ld\n", oct );
printf( "i=%ld, x=%ld, y=%ld\n", i, pt[i].x, pt[i].y );
j = nn[i][oct];
printf( "nn=%ld, x=%ld, y=%ld, dist = %ld\n", j, pt[j].x, pt[j].y,
dist(pt[i], pt[j ]) );
}
// assert( octant(pt[i], pt[ nn[i][oct] ]) == oct );
assert( octant(pt[i], pt[ nn1[i][oct] ]) == oct );
if( dist(pt[i], pt[ nn[i][oct] ]) !=
dist(pt[i], pt[ nn1[i][oct] ]) )
{
printf( "NNs DON'T MATCH!\noct=%ld\n", oct );
printf( "i=%ld, x=%ld, y=%ld\n", i, pt[i].x, pt[i].y );
j = nn[i][oct];
printf( "nn=%ld, x=%ld, y=%ld, dist = %ld\n", j, pt[j].x, pt[j].y,
dist(pt[i], pt[j ]) );
j = nn1[i][oct];
printf( "nn1=%ld, x=%ld, y=%ld, dist = %ld\n", j, pt[j].x, pt[j].y,
dist(pt[i], pt[ j ]) );
}
// assert( dist(pt[i], pt[ nn[i][oct] ]) ==
// dist(pt[i], pt[ nn1[i][oct] ]) );
}
}
}
free( nn1 );
}
/***************************************************************************/
/***************************************************************************/