#include "nlopt-util.h"
#include "nlopt.h"
#include "cdirect.h"
+#include "redblack.h"
#include "config.h"
#define MIN(a,b) ((a) < (b) ? (a) : (b))
/***************************************************************************/
/* basic data structure:
*
- * a hyper-rectangle is stored as an array of length 2n+2, where [0]
- * is the value of the function at the center, [1] is the "size"
- * measure of the rectangle, [2..n+1] are the coordinates of the
- * center, and [n+2..2n+1] are the widths of the sides.
+ * a hyper-rectangle is stored as an array of length L = 2n+2, where [0]
+ * is the value (f) of the function at the center, [1] is the "size"
+ * measure (d) of the rectangle, [2..n+1] are the coordinates of the
+ * center (c), and [n+2..2n+1] are the widths of the sides (w).
*
* a list of rectangles is just an array of N hyper-rectangles
- * stored as an N x (2n+1) in row-major order. Generally,
+ * stored as an N x L in row-major order. Generally,
* we allocate more than we need, allocating Na hyper-rectangles.
*
* n > 0 always
/* parameters of the search algorithm and various information that
needs to be passed around */
typedef struct {
+ int n; /* dimension */
+ int L; /* RECT_LEN(n) */
+ double *rects; /* the hyper-rectangles */
double magic_eps; /* Jones' epsilon parameter (1e-4 is recommended) */
int which_diam; /* which measure of hyper-rectangle diam to use:
0 = Jones, 1 = Gablonsky */
2: Jones Encyc. Opt.: pick random longest side */
const double *lb, *ub;
nlopt_stopping *stop; /* stopping criteria */
- int n;
nlopt_func f; void *f_data;
double *work; /* workspace, of length >= 2*n */
int *iwork, iwork_len; /* workspace, of length iwork_len >= n */
double fmin, *xmin; /* minimum so far */
+
+ /* red-black tree of hyperrect indices, sorted by (d,f) in
+ lexographical order */
+ rb_tree rtree;
} params;
/***************************************************************************/
#define EQUAL_SIDE_TOL 5e-2 /* tolerance to equate side sizes */
-/* divide rectangle idiv in the list rects */
-static nlopt_result divide_rect(int *N, int *Na, double **rects, int idiv,
- params *p)
+/* divide rectangle idiv in the list p->rects */
+static nlopt_result divide_rect(int *N, int *Na, int idiv, params *p)
{
int i;
const const int n = p->n;
- const int L = RECT_LEN(n);
- double *r = *rects;
+ const int L = p->L;
+ double *r = p->rects;
double *c = r + L*idiv + 2; /* center of rect to divide */
double *w = c + n; /* widths of rect to divide */
double wmax = w[0];
int imax = 0, nlongest = 0;
+ rb_node *node;
for (i = 1; i < n; ++i)
if (w[i] > wmax)
}
sort_fv(n, fv, isort);
ALLOC_RECTS(n, Na, r, (*N)+2*nlongest);
- *rects = r; c = r + L*idiv + 2; w = c + n;
+ p->rects = r; c = r + L*idiv + 2; w = c + n;
for (i = 0; i < nlongest; ++i) {
int k;
+ if (!(node = rb_tree_find(&p->rtree, idiv)))
+ return NLOPT_FAILURE;
w[isort[i]] *= THIRD;
r[L*idiv + 1] = rect_diameter(n, w, p);
+ rb_tree_resort(&p->rtree, node);
for (k = 0; k <= 1; ++k) {
memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
r[L*(*N) + 2 + isort[i]] += w[isort[i]] * (2*k-1);
r[L*(*N)] = fv[2*isort[i]+k];
+ if (!rb_tree_insert(&p->rtree, *N))
+ return NLOPT_FAILURE;
++(*N);
}
}
else
i = imax; /* trisect longest side */
ALLOC_RECTS(n, Na, r, (*N)+2);
- *rects = r; c = r + L*idiv + 2; w = c + n;
+ p->rects = r; c = r + L*idiv + 2; w = c + n;
+ if (!(node = rb_tree_find(&p->rtree, idiv)))
+ return NLOPT_FAILURE;
w[i] *= THIRD;
r[L*idiv + 1] = rect_diameter(n, w, p);
+ rb_tree_resort(&p->rtree, node);
for (k = 0; k <= 1; ++k) {
memcpy(r + L*(*N) + 1, c-1, sizeof(double) * (2*n+1));
r[L*(*N) + 2 + i] += w[i] * (2*k-1);
FUNCTION_EVAL(r[L*(*N)], r + L*(*N) + 2, p);
+ if (!rb_tree_insert(&p->rtree, *N))
+ return NLOPT_FAILURE;
++(*N);
}
}
/* O(N log N) convex hull algorithm, used later to find the potentially
optimal points */
-/* sort ihull by xy in lexographic order by x,y */
-static int s_qsort = 1; static double *xy_qsort = 0;
-static int sort_xy_compare(const void *a_, const void *b_)
-{
- int a = *((const int *) a_), b = *((const int *) b_);
- double xa = xy_qsort[a*s_qsort+1], xb = xy_qsort[b*s_qsort+1];
- if (xa < xb) return -1;
- else if (xb < xa) return +1;
- else {
- double ya = xy_qsort[a*s_qsort], yb = xy_qsort[b*s_qsort];
- if (ya < yb) return -1;
- else if (ya > yb) return +1;
- else return 0;
- }
-}
-static void sort_xy(int N, double *xy, int s, int *isort)
-{
- int i;
-
- for (i = 0; i < N; ++i) isort[i] = i;
- s_qsort = s; xy_qsort = xy;
- qsort(isort, (unsigned) N, sizeof(int), sort_xy_compare);
- xy_qsort = 0;
-}
-
/* Find the lower convex hull of a set of points (xy[s*i+1], xy[s*i]), where
0 <= i < N and s >= 2.
The return value is the number of points in the hull, with indices
- stored in ihull. ihull and is should point to arrays of length >= N.
-
- Note that we don't allow redundant points along the same line in the
- hull, similar to Gablonsky's version of DIRECT and differing from
- Jones'. */
-static int convex_hull(int N, double *xy, int s, int *ihull, int *is)
+ stored in ihull. ihull should point to arrays of length >= N.
+ rb_tree should be a red-black tree of indices (keys == i) sorted
+ in lexographic order by (xy[s*i+1], xy[s*i]).
+*/
+static int convex_hull(int N, double *xy, int s, int *ihull, rb_tree *t)
{
- int minmin; /* first index (smallest y) with min x */
- int minmax; /* last index (largest y) with min x */
- int maxmin; /* first index (smallest y) with max x */
- int maxmax; /* last index (largest y) with max x */
- int i, nhull = 0;
+ int nhull;
double minslope;
- double xmin, xmax;
+ double xmin, xmax, yminmin, ymaxmin;
+ rb_node *n, *nmax;
+ if (N == 0) return 0;
+
/* Monotone chain algorithm [Andrew, 1979]. */
- sort_xy(N, xy, s, is);
+ n = rb_tree_min(t);
+ nmax = rb_tree_max(t);
- xmin = xy[s*is[minmin=0]+1]; xmax = xy[s*is[maxmax=N-1]+1];
+ xmin = xy[s*(n->k)+1];
+ yminmin = xy[s*(n->k)];
+ xmax = xy[s*(nmax->k)+1];
- if (xmin == xmax) { /* degenerate case */
- ihull[nhull++] = is[minmin];
- return nhull;
- }
+ ihull[nhull = 1] = n->k;
+ if (xmin == xmax) return nhull;
- for (minmax = minmin; minmax+1 < N && xy[s*is[minmax+1]+1]==xmin;
- ++minmax);
- for (maxmin = maxmax; maxmin-1>=0 && xy[s*is[maxmin-1]+1]==xmax;
- --maxmin);
+ /* set nmax = min mode with x == xmax */
+ while (xy[s*(nmax->k)+1] == xmax)
+ nmax = rb_tree_pred(t, nmax); /* non-NULL since xmin != xmax */
+ nmax = rb_tree_succ(t, nmax);
- minslope = (xy[s*is[maxmin]] - xy[s*is[minmin]]) / (xmax - xmin);
+ ymaxmin = xy[s*(nmax->k)];
+ minslope = (ymaxmin - yminmin) / (xmax - xmin);
- ihull[nhull++] = is[minmin];
- for (i = minmax + 1; i < maxmin; ++i) {
- int k = is[i];
- if (xy[s*k] > xy[s*is[minmin]] + (xy[s*k+1] - xmin) * minslope)
+ /* set n = first node with x != xmin */
+ while (xy[s*(n->k)+1] == xmin)
+ n = rb_tree_succ(t, n); /* non-NULL since xmin != xmax */
+
+ for (; n != nmax; n = rb_tree_succ(t, n)) {
+ int k = n->k;
+ if (xy[s*k] > yminmin + (xy[s*k+1] - xmin) * minslope)
continue;
/* remove points until we are making a "left turn" to k */
while (nhull > 1) {
}
ihull[nhull++] = k;
}
- ihull[nhull++] = is[maxmin];
+ ihull[nhull++] = nmax->k;
return nhull;
}
return 1;
}
-static nlopt_result divide_good_rects(int *N, int *Na, double **rects,
- params *p)
+static nlopt_result divide_good_rects(int *N, int *Na, params *p)
{
const int n = p->n;
- const int L = RECT_LEN(n);
+ const int L = p->L;
int *ihull, nhull, i, xtol_reached = 1, divided_some = 0;
- double *r = *rects;
+ double *r = p->rects;
double magic_eps = p->magic_eps;
- if (p->iwork_len < n + 2*(*N)) {
- p->iwork_len = p->iwork_len + n + 2*(*N);
+ if (p->iwork_len < *N) {
+ p->iwork_len = p->iwork_len + *N;
p->iwork = (int *) realloc(p->iwork, sizeof(int) * p->iwork_len);
if (!p->iwork)
return NLOPT_OUT_OF_MEMORY;
}
ihull = p->iwork;
- nhull = convex_hull(*N, r, L, ihull, ihull + *N);
+ nhull = convex_hull(*N, r, L, ihull, &p->rtree);
divisions:
for (i = 0; i < nhull; ++i) {
double K1 = -HUGE_VAL, K2 = -HUGE_VAL, K;
/* "potentially optimal" rectangle, so subdivide */
divided_some = 1;
nlopt_result ret;
- ret = divide_rect(N, Na, rects, ihull[i], p);
- r = *rects; /* may have grown */
+ ret = divide_rect(N, Na, ihull[i], p);
+ r = p->rects; /* may have grown */
if (ret != NLOPT_SUCCESS) return ret;
xtol_reached = xtol_reached && small(r + L*ihull[i] + 2+n, p);
}
for (i = 1; i < *N; ++i)
if (r[L*i+1] > wmax)
wmax = r[L*(imax=i)+1];
- return divide_rect(N, Na, rects, imax, p);
+ return divide_rect(N, Na, imax, p);
}
}
return xtol_reached ? NLOPT_XTOL_REACHED : NLOPT_SUCCESS;
/***************************************************************************/
+/* lexographic sort order (d,f) of hyper-rects, for red-black tree */
+static int hyperrect_compare(int i, int j, void *p_)
+{
+ params *p = (params *) p_;
+ int L = p->L;
+ double *r = p->rects;
+ double di = r[i*L+1], dj = r[j*L+1], fi, fj;
+ if (di < dj) return -1;
+ if (dj < di) return +1;
+ fi = r[i*L]; fj = r[j*L];
+ if (fi < fj) return -1;
+ if (fj < fi) return +1;
+ return 0;
+}
+
+/***************************************************************************/
+
nlopt_result cdirect_unscaled(int n, nlopt_func f, void *f_data,
const double *lb, const double *ub,
double *x,
double magic_eps, int which_alg)
{
params p;
- double *rects;
int Na = 100, N = 1, i, x_center = 1;
nlopt_result ret = NLOPT_OUT_OF_MEMORY;
p.lb = lb; p.ub = ub;
p.stop = stop;
p.n = n;
+ p.L = RECT_LEN(n);
p.f = f;
p.f_data = f_data;
p.xmin = x;
p.fmin = f(n, x, NULL, f_data); stop->nevals++;
p.work = 0;
p.iwork = 0;
- rects = 0;
+ p.rects = 0;
+
+ if (!rb_tree_init(&p.rtree, hyperrect_compare, &p)) goto done;
p.work = (double *) malloc(sizeof(double) * 2*n);
if (!p.work) goto done;
- rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
- if (!rects) goto done;
- p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = 2*Na + n));
+ p.rects = (double *) malloc(sizeof(double) * Na * RECT_LEN(n));
+ if (!p.rects) goto done;
+ p.iwork = (int *) malloc(sizeof(int) * (p.iwork_len = Na + n));
if (!p.iwork) goto done;
for (i = 0; i < n; ++i) {
- rects[2+i] = 0.5 * (lb[i] + ub[i]);
+ p.rects[2+i] = 0.5 * (lb[i] + ub[i]);
x_center = x_center
- && (fabs(rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
- rects[2+n+i] = ub[i] - lb[i];
+ && (fabs(p.rects[2+i]-x[i]) < 1e-13*(1+fabs(x[i])));
+ p.rects[2+n+i] = ub[i] - lb[i];
}
- rects[1] = rect_diameter(n, rects+2+n, &p);
+ p.rects[1] = rect_diameter(n, p.rects+2+n, &p);
if (x_center)
- rects[0] = p.fmin; /* avoid computing f(center) twice */
+ p.rects[0] = p.fmin; /* avoid computing f(center) twice */
else
- rects[0] = function_eval(rects+2, &p);
+ p.rects[0] = function_eval(p.rects+2, &p);
+ if (!rb_tree_insert(&p.rtree, 0)) {
+ ret = NLOPT_FAILURE;
+ goto done;
+ }
- ret = divide_rect(&N, &Na, &rects, 0, &p);
+ ret = divide_rect(&N, &Na, 0, &p);
if (ret != NLOPT_SUCCESS) goto done;
while (1) {
double fmin0 = p.fmin;
- ret = divide_good_rects(&N, &Na, &rects, &p);
+ ret = divide_good_rects(&N, &Na, &p);
if (ret != NLOPT_SUCCESS) goto done;
if (nlopt_stop_f(p.stop, p.fmin, fmin0)) {
ret = NLOPT_FTOL_REACHED;
}
done:
+ rb_tree_destroy(&p.rtree);
free(p.iwork);
- free(rects);
+ free(p.rects);
free(p.work);
*fmin = p.fmin;
--- /dev/null
+/* simple implementation of red-black trees optimized for use with DIRECT */
+
+#include <stddef.h>
+#include <stdlib.h>
+#include "redblack.h"
+
+int rb_tree_init(rb_tree *t, rb_compare compare, void *c_data) {
+ t->compare = compare; t->c_data = c_data;
+ t->nil.c = BLACK; t->nil.l = t->nil.r = t->nil.p = &t->nil;
+ t->root = &t->nil;
+ t->N = 0;
+ t->Nalloc = 100; /* allocate some space to start with */
+ t->nodes = (rb_node*) malloc(sizeof(rb_node) * t->Nalloc);
+ return t->nodes != NULL;
+}
+
+void rb_tree_destroy(rb_tree *t)
+{
+ t->root = 0; t->N = 0; t->Nalloc = 0;
+ free(t->nodes); t->nodes = 0;
+}
+
+/* in our application, we can optimize memory allocation because
+ we never delete two nodes in a row (we always add a node after
+ deleting)... or rather, we never delete but the value of
+ the key sometimes changes. ... this means we can just
+ allocate a linear, exponentially growing stack (nodes) of
+ nodes, and don't have to worry about holes in the stack ...
+ otherwise, alloc1 should be replaced by an implementation that
+ malloc's each node separately */
+static rb_node *alloc1(rb_tree *t, int k)
+{
+ rb_node *nil = &t->nil;
+ rb_node *n;
+ if (t->Nalloc == t->N) { /* grow allocation */
+ rb_node *old_nodes = t->nodes;
+ ptrdiff_t change;
+ int i;
+ t->Nalloc = 2*t->Nalloc + 1;
+ t->nodes = (rb_node*) realloc(t->nodes, sizeof(rb_node) * t->Nalloc);
+ if (!t->nodes) return NULL;
+ change = t->nodes - old_nodes;
+ if (t->root != nil) t->root += change;
+ for (i = 0; i < t->N; ++i) { /* shift all pointers, ugh */
+ if (t->nodes[i].p != nil) t->nodes[i].p += change;
+ if (t->nodes[i].r != nil) t->nodes[i].r += change;
+ if (t->nodes[i].l != nil) t->nodes[i].l += change;
+ }
+ }
+ n = t->nodes + t->N++;
+ n->k = k;
+ n->p = n->l = n->r = nil;
+ return n;
+}
+
+static void rotate_left(rb_node *p, rb_tree *t)
+{
+ rb_node *nil = &t->nil;
+ rb_node *n = p->r; /* must be non-NULL */
+ p->r = n->l;
+ n->l = p;
+ if (p->p != nil) {
+ if (p == p->p->l) p->p->l = n;
+ else p->p->r = n;
+ }
+ else
+ t->root = n;
+ n->p = p->p;
+ p->p = n;
+ if (p->r != nil) p->r->p = p;
+}
+
+static void rotate_right(rb_node *p, rb_tree *t)
+{
+ rb_node *nil = &t->nil;
+ rb_node *n = p->l; /* must be non-NULL */
+ p->l = n->r;
+ n->r = p;
+ if (p->p != nil) {
+ if (p == p->p->l) p->p->l = n;
+ else p->p->r = n;
+ }
+ else
+ t->root = n;
+ n->p = p->p;
+ p->p = n;
+ if (p->l != nil) p->l->p = p;
+}
+
+static void insert_node(rb_tree *t, rb_node *n)
+{
+ rb_node *nil = &t->nil;
+ rb_compare compare = t->compare;
+ void *c_data = t->c_data;
+ int k = n->k;
+ rb_node *p = t->root;
+ n->c = RED;
+ if (p == nil) {
+ t->root = n;
+ n->c = BLACK;
+ return;
+ }
+ /* insert (RED) node into tree */
+ while (1) {
+ if (compare(k, p->k, c_data) <= 0) { /* k <= p->k */
+ if (p->l != nil)
+ p = p->l;
+ else {
+ p->l = n;
+ n->p = p;
+ break;
+ }
+ }
+ else {
+ if (p->r != nil)
+ p = p->r;
+ else {
+ p->r = n;
+ n->p = p;
+ break;
+ }
+ }
+ }
+ fixtree:
+ if (n->p->c == RED) { /* red cannot have red child */
+ rb_node *u = p == p->p->l ? p->p->r : p->p->l;
+ if (u != nil && u->c == RED) {
+ p->c = u->c = BLACK;
+ n = p->p;
+ if ((p = n->p) != nil) {
+ n->c = RED;
+ goto fixtree;
+ }
+ }
+ else {
+ if (n == p->r && p == p->p->l) {
+ rotate_left(p, t);
+ p = n; n = n->l;
+ }
+ else if (n == p->l && p == p->p->r) {
+ rotate_right(p, t);
+ p = n; n = n->r;
+ }
+ p->c = BLACK;
+ p->p->c = RED;
+ if (n == p->l && p == p->p->l)
+ rotate_right(p->p, t);
+ else if (n == p->r && p == p->p->r)
+ rotate_left(p->p, t);
+ }
+
+ }
+}
+
+int rb_tree_insert(rb_tree *t, int k)
+{
+ rb_node *n = alloc1(t, k);
+ if (!n) return 0;
+ insert_node(t, n);
+ return 1;
+}
+
+static int check_node(rb_node *n, int *nblack, rb_node *nil)
+{
+ int nbl, nbr;
+ if (n == nil) { *nblack = 0; return 1; }
+ if (n->r != nil && n->r->p != n) return 0;
+ if (n->l != nil && n->l->p != n) return 0;
+ if (n->c == RED) {
+ if (n->r != nil && n->r->c == RED) return 0;
+ if (n->l != nil && n->l->c == RED) return 0;
+ }
+ if (!(check_node(n->r, &nbl, nil) && check_node(n->l, &nbr, nil)))
+ return 0;
+ if (nbl != nbr) return 0;
+ *nblack = nbl + n->c == BLACK;
+ return 1;
+}
+int rb_tree_check(rb_tree *t)
+{
+ rb_node *nil = &t->nil;
+ int nblack;
+ if (nil->c != BLACK) return 0;
+ if (t->root == nil) return 1;
+ if (t->root->c != BLACK) return 0;
+ return check_node(t->root, &nblack, nil);
+}
+
+rb_node *rb_tree_find(rb_tree *t, int k)
+{
+ rb_node *nil = &t->nil;
+ rb_compare compare = t->compare;
+ void *c_data = t->c_data;
+ rb_node *p = t->root;
+ while (p != nil) {
+ int comp = compare(k, p->k, c_data);
+ if (!comp) return p;
+ p = comp <= 0 ? p->l : p->r;
+ }
+ return NULL;
+}
+
+rb_node *rb_tree_min(rb_tree *t)
+{
+ rb_node *nil = &t->nil;
+ rb_node *n = t->root;
+ while (n != nil && n->l != nil)
+ n = n->l;
+ return(n == nil ? NULL : n);
+}
+
+rb_node *rb_tree_max(rb_tree *t)
+{
+ rb_node *nil = &t->nil;
+ rb_node *n = t->root;
+ while (n != nil && n->r != nil)
+ n = n->r;
+ return(n == nil ? NULL : n);
+}
+
+rb_node *rb_tree_succ(rb_tree *t, rb_node *n)
+{
+ rb_node *nil = &t->nil;
+ if (n->r == nil) {
+ rb_node *prev;
+ do {
+ prev = n;
+ n = n->p;
+ } while (prev == n->r && n != nil);
+ return n == nil ? NULL : n;
+ }
+ else {
+ n = n->r;
+ while (n->l != nil)
+ n = n->l;
+ return n;
+ }
+}
+
+rb_node *rb_tree_pred(rb_tree *t, rb_node *n)
+{
+ rb_node *nil = &t->nil;
+ if (n->l == nil) {
+ rb_node *prev;
+ do {
+ prev = n;
+ n = n->p;
+ } while (prev == n->l && n != nil);
+ return n == nil ? NULL : n;
+ }
+ else {
+ n = n->l;
+ while (n->r != nil)
+ n = n->r;
+ return n;
+ }
+}
+
+rb_node *rb_tree_remove(rb_tree *t, rb_node *n)
+{
+ rb_node *nil = &t->nil;
+ rb_node *m;
+ if (n->l != nil && n->r != nil) {
+ rb_node *lmax = n->l;
+ while (lmax->r != nil) lmax = lmax->r;
+ n->k = lmax->k;
+ n = lmax;
+ }
+ m = n->l != nil? n->l : n->r;
+ if (n->p != nil) {
+ if (n->p->r == n) n->p->r = m;
+ else n->p->l = m;
+ }
+ else
+ t->root = m;
+ m->p = n->p;
+ if (n->c == BLACK) {
+ if (m->c == RED)
+ m->c = BLACK;
+ else {
+ deleteblack:
+ if (m->p != nil) {
+ rb_node *s = m == m->p->l ? m->p->r : m->p->l;
+ if (s->c == RED) {
+ m->p->c = RED;
+ s->c = BLACK;
+ if (m == m->p->l) rotate_left(m->p, t);
+ else rotate_right(m->p, t);
+ s = m == m->p->l ? m->p->r : m->p->l;
+ }
+ if (m->p->c == BLACK && s->c == BLACK
+ && s->l->c == BLACK && s->r->c == BLACK) {
+ if (s != nil) s->c = RED;
+ m = m->p;
+ goto deleteblack;
+ }
+ else if (m->p->c == RED && s->c == BLACK &&
+ s->l->c == BLACK && s->r->c == BLACK) {
+ if (s != nil) s->c = RED;
+ m->p->c = BLACK;
+ }
+ else {
+ if (m == m->p->l && s->c == BLACK &&
+ s->l->c == RED && s->r->c == BLACK) {
+ s->c = RED;
+ s->l->c = BLACK;
+ rotate_right(s, t);
+ s = m == m->p->l ? m->p->r : m->p->l;
+ }
+ else if (m == m->p->r && s->c == BLACK &&
+ s->r->c == RED && s->l->c == BLACK) {
+ s->c = RED;
+ s->r->c = BLACK;
+ rotate_left(s, t);
+ s = m == m->p->l ? m->p->r : m->p->l;
+ }
+ s->c = m->p->c;
+ m->p->c = BLACK;
+ if (m == m->p->l) {
+ s->r->c = BLACK;
+ rotate_left(m->p, t);
+ }
+ else {
+ s->l->c = BLACK;
+ rotate_right(m->p, t);
+ }
+ }
+ }
+ }
+ }
+ return n; /* the node that was deleted may be different from initial n */
+}
+
+rb_node *rb_tree_resort(rb_tree *t, rb_node *n)
+{
+ int k = n->k;
+ n = rb_tree_remove(t, n);
+ n->p = n->l = n->r = &t->nil;
+ n->k = k; /* n may have changed during remove */
+ insert_node(t, n);
+ return n;
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <time.h>
+#include "redblack.h"
+
+static int comp(int k1, int k2, void *dummy)
+{
+ (void) dummy;
+ return k1 - k2;
+}
+
+int main(int argc, char **argv)
+{
+ int N, M;
+ int *k;
+ rb_tree t;
+ rb_node *n;
+ int i, j;
+
+ if (argc != 2) {
+ fprintf(stderr, "Usage: redblack_test Ntest\n");
+ return 1;
+ }
+
+ N = atoi(argv[1]);
+ k = (int *) malloc(N * sizeof(int));
+ if (!rb_tree_init(&t, comp, NULL)) {
+ fprintf(stderr, "error in rb_tree_init\n");
+ return 1;
+ }
+
+ srand((unsigned) time(NULL));
+ for (i = 0; i < N; ++i) {
+ if (!rb_tree_insert(&t, k[i] = rand() % N)) {
+ fprintf(stderr, "error in rb_tree_insert\n");
+ return 1;
+ }
+ if (!rb_tree_check(&t)) {
+ fprintf(stderr, "rb_tree_check_failed after insert!\n");
+ return 1;
+ }
+ }
+
+ for (i = 0; i < N; ++i)
+ if (!rb_tree_find(&t, k[i])) {
+ fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
+ return 1;
+ }
+
+ n = rb_tree_min(&t);
+ for (i = 0; i < N; ++i) {
+ if (!n) {
+ fprintf(stderr, "not enough successors %d\n!", i);
+ return 1;
+ }
+ printf("%d: %d\n", i, n->k);
+ n = rb_tree_succ(&t, n);
+ }
+ if (n) {
+ fprintf(stderr, "too many successors!\n");
+ return 1;
+ }
+
+ n = rb_tree_max(&t);
+ for (i = 0; i < N; ++i) {
+ if (!n) {
+ fprintf(stderr, "not enough predecessors %d\n!", i);
+ return 1;
+ }
+ printf("%d: %d\n", i, n->k);
+ n = rb_tree_pred(&t, n);
+ }
+ if (n) {
+ fprintf(stderr, "too many predecessors!\n");
+ return 1;
+ }
+
+ for (M = N; M > 0; --M) {
+ int knew;
+ j = rand() % M;
+ for (i = 0; i < N; ++i)
+ if (k[i] >= 0)
+ if (j-- == 0)
+ break;
+ if (i >= N) abort();
+ if (!(n = rb_tree_find(&t, k[i]))) {
+ fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
+ return 1;
+ }
+ n->k = knew;
+ if (!rb_tree_resort(&t, n)) {
+ fprintf(stderr, "error in rb_tree_resort\n");
+ return 1;
+ }
+ if (!rb_tree_check(&t)) {
+ fprintf(stderr, "rb_tree_check_failed after change!\n");
+ return 1;
+ }
+ k[i] = -1 - knew;
+ }
+
+ for (i = 0; i < N; ++i)
+ k[i] = -1 - k[i];
+
+ for (M = N; M > 0; --M) {
+ j = rand() % M;
+ for (i = 0; i < N; ++i)
+ if (k[i] >= 0)
+ if (j-- == 0)
+ break;
+ if (i >= N) abort();
+ if (!(n = rb_tree_find(&t, k[i]))) {
+ fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
+ return 1;
+ }
+ rb_tree_remove(&t, n);
+ if (!rb_tree_check(&t)) {
+ fprintf(stderr, "rb_tree_check_failed after remove!\n");
+ return 1;
+ }
+ k[i] = -1 - k[i];
+ }
+
+ rb_tree_destroy(&t);
+ free(k);
+ return 0;
+}
~ianmdlvl / nlopt.git / commitdiff