bug fix in cdirect - handle case where different indices have same (x,y)

darcs-hash:20070828034729-c8de0-af82fa8b6d60721808387791951cddc46e154067.gz

diff --git a/cdirect/cdirect.c b/cdirect/cdirect.c
index 210907ffdd963fc78f84a8c7e68feb7665846197..8906446f3372db0e491517fd7af2a0787c28a823 100644 (file)
--- a/cdirect/cdirect.c
+++ b/cdirect/cdirect.c
@@ -164,7 +164,7 @@ static nlopt_result divide_rect(int *N, int *Na, int idiv, params *p)
          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)))
+              if (!(node = rb_tree_find_exact(&p->rtree, idiv)))
                    return NLOPT_FAILURE;
               w[isort[i]] *= THIRD;
               r[L*idiv + 1] = rect_diameter(n, w, p);
@@ -194,7 +194,7 @@ static nlopt_result divide_rect(int *N, int *Na, int idiv, params *p)
               i = imax; /* trisect longest side */
          ALLOC_RECTS(n, Na, r, (*N)+2);
           p->rects = r; c = r + L*idiv + 2; w = c + n;
-         if (!(node = rb_tree_find(&p->rtree, idiv)))
+         if (!(node = rb_tree_find_exact(&p->rtree, idiv)))
               return NLOPT_FAILURE;
          w[i] *= THIRD;
          r[L*idiv + 1] = rect_diameter(n, w, p);
@@ -218,6 +218,8 @@ static nlopt_result divide_rect(int *N, int *Na, int idiv, params *p)
 /* Find the lower convex hull of a set of points (xy[s*i+1], xy[s*i]), where
    0 <= i < N and s >= 2.
 
+   Unlike standard convex hulls, we allow redundant points on the hull.
+
    The return value is the number of points in the hull, with indices
    stored in ihull.  ihull should point to arrays of length >= N.
    rb_tree should be a red-black tree of indices (keys == i) sorted
@@ -265,7 +267,7 @@ static int convex_hull(int N, double *xy, int s, int *ihull, rb_tree *t)
               int t1 = ihull[nhull - 1], t2 = ihull[nhull - 2];
               /* cross product (t1-t2) x (k-t2) > 0 for a left turn: */
               if ((xy[s*t1+1]-xy[s*t2+1]) * (xy[s*k]-xy[s*t2])
-                  - (xy[s*t1]-xy[s*t2]) * (xy[s*k+1]-xy[s*t2+1]) > 0)
+                  - (xy[s*t1]-xy[s*t2]) * (xy[s*k+1]-xy[s*t2+1]) >= 0)
                    break;
               --nhull;
          }

diff --git a/cdirect/redblack.c b/cdirect/redblack.c
index 7472c732a1a9779b79b896729c73bc27cc9edfb2..3001486e8c0d621155af5d6624b6cab164a698d0 100644 (file)
--- a/cdirect/redblack.c
+++ b/cdirect/redblack.c
@@ -6,7 +6,7 @@
 
 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->nil.c = BLACK; t->nil.l = t->nil.r = t->nil.p = &t->nil; t->nil.k = -1;
      t->root = &t->nil;
      t->N = 0;
      t->Nalloc = 100; /* allocate some space to start with */
@@ -160,17 +160,24 @@ int rb_tree_insert(rb_tree *t, int k)
      return 1;
 }
 
-static int check_node(rb_node *n, int *nblack, rb_node *nil)
+static int check_node(rb_node *n, int *nblack, rb_tree *t)
 {
+     rb_node *nil = &t->nil;
      int nbl, nbr;
+     rb_compare compare = t->compare;
+     void *c_data = t->c_data;
      if (n == nil) { *nblack = 0; return 1; }
      if (n->r != nil && n->r->p != n) return 0;
+     if (n->r != nil && compare(n->r->k, n->k, c_data) < 0)
+         return 0;
      if (n->l != nil && n->l->p != n) return 0;
+     if (n->l != nil && compare(n->l->k, n->k, c_data) > 0)
+         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))) 
+     if (!(check_node(n->r, &nbl, t) && check_node(n->l, &nbr, t))) 
          return 0;
      if (nbl != nbr) return 0;
      *nblack = nbl + n->c == BLACK;
@@ -183,7 +190,7 @@ int rb_tree_check(rb_tree *t)
      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);
+     return check_node(t->root, &nblack, t);
 }
 
 rb_node *rb_tree_find(rb_tree *t, int k)
@@ -200,6 +207,78 @@ rb_node *rb_tree_find(rb_tree *t, int k)
      return NULL;
 }
 
+/* like rb_tree_find, but guarantees that returned node n will have
+   n->k == k (may not be true above if compare(k,k') == 0 for some k != k') */
+rb_node *rb_tree_find_exact(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) break;
+         p = comp <= 0 ? p->l : p->r;
+     }
+     if (p == nil)
+         return NULL;
+     while (p->l != nil && !compare(k, p->l->k, c_data)) p = p->l;
+     if (p->l != nil) p = p->l;
+     do {
+         if (p->k == k) return p;
+         p = rb_tree_succ(t, p);
+     } while (p && compare(p->k, k, c_data) <= 0);
+     return NULL;
+}
+
+/* find greatest point in subtree p that is <= k */
+rb_node *find_le(rb_node *p, int k, rb_tree *t)
+{
+     rb_node *nil = &t->nil;
+     rb_compare compare = t->compare;
+     void *c_data = t->c_data;
+     while (p != nil) {
+         if (compare(p->k, k, c_data) <= 0) { /* p->k <= k */
+              rb_node *r = find_le(p->r, k, t);
+              if (r) return r;
+              else return p;
+         }
+         else /* p->k > k */
+              p = p->l;
+     }
+     return NULL; /* k < everything in subtree */
+}
+
+/* find greatest point in t <= k */
+rb_node *rb_tree_find_le(rb_tree *t, int k)
+{
+     return find_le(t->root, k, t);
+}
+
+/* find least point in subtree p that is > k */
+rb_node *find_gt(rb_node *p, int k, rb_tree *t)
+{
+     rb_node *nil = &t->nil;
+     rb_compare compare = t->compare;
+     void *c_data = t->c_data;
+     while (p != nil) {
+         if (compare(p->k, k, c_data) > 0) { /* p->k > k */
+              rb_node *l = find_gt(p->l, k, t);
+              if (l) return l;
+              else return p;
+         }
+         else /* p->k <= k */
+              p = p->r;
+     }
+     return NULL; /* k >= everything in subtree */
+}
+
+/* find least point in t > k */
+rb_node *rb_tree_find_gt(rb_tree *t, int k)
+{
+     return find_gt(t->root, k, t);
+}
+
 rb_node *rb_tree_min(rb_tree *t)
 {
      rb_node *nil = &t->nil;

diff --git a/cdirect/redblack.h b/cdirect/redblack.h
index 5aa54c76e28e2085649bc310ebe9c98f19b711bb..6aa0831a81d313c63d8841d8c7cde185ce90e9c7 100644 (file)
--- a/cdirect/redblack.h
+++ b/cdirect/redblack.h
@@ -36,6 +36,7 @@ extern void rb_tree_destroy(rb_tree *t);
 extern int rb_tree_insert(rb_tree *t, int k);
 extern int rb_tree_check(rb_tree *t);
 extern rb_node *rb_tree_find(rb_tree *t, int k);
+extern rb_node *rb_tree_find_exact(rb_tree *t, int k);
 extern rb_node *rb_tree_resort(rb_tree *t, rb_node *n);
 extern rb_node *rb_tree_min(rb_tree *t);
 extern rb_node *rb_tree_max(rb_tree *t);

diff --git a/cdirect/redblack_test.c b/cdirect/redblack_test.c
index f69676067f57a0bc2f0f40bf95276b7dcaa42460..7aa1757bc1c36c8c02fe13aa950b125957f00cb6 100644 (file)
--- a/cdirect/redblack_test.c
+++ b/cdirect/redblack_test.c
@@ -42,7 +42,7 @@ int main(int argc, char **argv)
      }
      
      for (i = 0; i < N; ++i)
-         if (!rb_tree_find(&t, k[i])) {
+         if (!rb_tree_find(&t, k[i]) || !rb_tree_find_exact(&t, k[i])) {
               fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
               return 1;
          }
@@ -76,14 +76,14 @@ int main(int argc, char **argv)
      }
      
      for (M = N; M > 0; --M) {
-         int knew;
-         j = rand() % M;
+         int knew = rand() % N; /* random new key */
+         j = rand() % M; /* random original key to replace */
          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]))) {
+         if (!(n = rb_tree_find(&t, k[i])) || !rb_tree_find_exact(&t, k[i])) {
                fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
                return 1;
           }
@@ -109,7 +109,7 @@ int main(int argc, char **argv)
                    if (j-- == 0)
                         break;
          if (i >= N) abort();
-         if (!(n = rb_tree_find(&t, k[i]))) {
+         if (!(n = rb_tree_find(&t, k[i])) || !rb_tree_find_exact(&t, k[i])) {
               fprintf(stderr, "rb_tree_find lost %d!\n", k[i]);
               return 1;
          }