Mercurial > octave-nkf
changeset 8814:de16ebeef93d
improve lookup, provide Array<T>::lookup
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Thu, 19 Feb 2009 15:19:59 +0100 |
| parents | 70d06ed27c08 |
| children | af907aeedbf4 |
| files | liboctave/Array.cc liboctave/Array.h liboctave/ChangeLog liboctave/Makefile.in liboctave/oct-lookup.h liboctave/oct-sort.cc liboctave/oct-sort.h src/ChangeLog src/DLD-FUNCTIONS/lookup.cc |
| diffstat | 9 files changed, 511 insertions(+), 394 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/Array.cc Thu Feb 19 07:34:15 2009 -0500 +++ b/liboctave/Array.cc Thu Feb 19 15:19:59 2009 +0100 @@ -2383,6 +2383,67 @@ } +// Do a binary lookup in a sorted array. +template <class T> +octave_idx_type +Array<T>::lookup (const T& value, sortmode mode) const +{ + octave_idx_type n = numel (); + octave_sort<T> lsort; + + if (mode == UNSORTED) + { + // auto-detect mode + if (n > 1 && lsort.descending_compare (elem (0), elem (n-1))) + mode = DESCENDING; + else + mode = ASCENDING; + } + + lsort.set_compare (mode); + + return lsort.lookup (data (), n, value); +} + +// Ditto, but for an array of values, specializing on long runs. +// Adds optional offset to all indices. +template <class T> +Array<octave_idx_type> +Array<T>::lookup (const Array<T>& values, sortmode mode, + bool linf, bool rinf) const +{ + octave_idx_type n = numel (); + octave_sort<T> lsort; + Array<octave_idx_type> idx (values.dims ()); + + if (mode == UNSORTED) + { + // auto-detect mode + if (n > 1 && lsort.descending_compare (elem (0), elem (n-1))) + mode = DESCENDING; + else + mode = ASCENDING; + } + + lsort.set_compare (mode); + + // set offset and shift size. + octave_idx_type offset = 0; + + if (linf && n > 0) + { + offset++; + n--; + } + if (rinf && n > 0) + n--; + + lsort.lookup (data () + offset, n, values.data (), values.numel (), + idx.fortran_vec (), offset); + + return idx; +} + #define INSTANTIATE_ARRAY_SORT(T) template class octave_sort<T>; #define NO_INSTANTIATE_ARRAY_SORT(T) \ @@ -2409,6 +2470,13 @@ template <> sortmode \ Array<T>::is_sorted_rows (sortmode) const \ { return UNSORTED; } \ + \ +template <> octave_idx_type \ +Array<T>::lookup (const T&, sortmode) const \ +{ return 0; } \ +template <> Array<octave_idx_type> \ +Array<T>::lookup (const Array<T>&, sortmode, bool, bool) const \ +{ return Array<octave_idx_type> (); } \
--- a/liboctave/Array.h Thu Feb 19 07:34:15 2009 -0500 +++ b/liboctave/Array.h Thu Feb 19 15:19:59 2009 +0100 @@ -567,6 +567,18 @@ // Ordering is auto-detected or can be specified. sortmode is_sorted_rows (sortmode mode = UNSORTED) const; + // Do a binary lookup in a sorted array. + // Mode can be specified or is auto-detected by comparing 1st and last element. + octave_idx_type lookup (const T& value, sortmode mode = UNSORTED) const; + + // Ditto, but for an array of values, specializing on long runs. + // If linf is true, the leftmost interval is extended to infinity + // (indices will be >= 1). + // If rinf is true, the rightmost interval is extended to infinity + // (indices will be <= length ()-1). + Array<octave_idx_type> lookup (const Array<T>& values, sortmode mode = UNSORTED, + bool linf = false, bool rinf = false) const; + Array<T> diag (octave_idx_type k = 0) const; template <class U, class F>
--- a/liboctave/ChangeLog Thu Feb 19 07:34:15 2009 -0500 +++ b/liboctave/ChangeLog Thu Feb 19 15:19:59 2009 +0100 @@ -1,3 +1,16 @@ +2009-02-19 Jaroslav Hajek <highegg@gmail.com> + + * oct-types.h (sortmode): Move enum here. + * oct-sort.h (octave_sort<T>::ms): Declare as pointer. + (octave_sort<T>::lookup): New overloaded method. + * oct-sort.cc: Reflect change to ms. + (octave_sort<T>::lookup): New overloaded method. + (out_of_range_pred): New helper class. + (out_of_range): New helper function. + * oct-lookup.h: Remove file. + * Array.cc (Array<T>::lookup): New overloaded method. + * Array.h: Declare it. + 2009-02-18 John W. Eaton <jwe@octave.org> * dbleQR.cc (QR::init, QR::form): Cast int to octave_idx_type in
--- a/liboctave/Makefile.in Thu Feb 19 07:34:15 2009 -0500 +++ b/liboctave/Makefile.in Thu Feb 19 15:19:59 2009 +0100 @@ -88,7 +88,7 @@ lo-ieee.h lo-mappers.h lo-math.h lo-specfun.h lo-sysdep.h \ lo-traits.h lo-utils.h mach-info.h md5.h oct-alloc.h oct-cmplx.h \ oct-env.h oct-fftw.h oct-getopt.h oct-group.h oct-inttypes.h \ - oct-locbuf.h oct-lookup.h oct-md5.h oct-mutex.h oct-norm.h \ + oct-locbuf.h oct-md5.h oct-mutex.h oct-norm.h \ oct-passwd.h oct-rand.h oct-rl-edit.h oct-rl-hist.h oct-shlib.h \ oct-sort.h oct-spparms.h oct-syscalls.h oct-sparse.h oct-time.h \ oct-uname.h pathlen.h pathsearch.h prog-args.h \
--- a/liboctave/oct-lookup.h Thu Feb 19 07:34:15 2009 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,175 +0,0 @@ -/* - -Copyright (C) 2008 VZLU Prague, a.s., Czech Republic - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -// Author: Jaroslav Hajek <highegg@gmail.com> - -#if !defined (octave_oct_lookup) -#define octave_oct_lookup 1 - -#include <algorithm> -#include <functional> - -#include "oct-types.h" - -// a simple binary lookup -template<typename T, typename bpred> -inline octave_idx_type -bin_lookup (const T *table, octave_idx_type size, - const T& val, - bpred comp) -{ - return std::upper_bound (table, table + size, val, comp) - table; -} - -// version using < operator -template<typename T> -inline octave_idx_type -bin_lookup (const T *table, octave_idx_type size, - const T& val) -{ - return std::upper_bound (table, table + size, val) - table; -} - -// a unary functor that checks whether a value is outside [a,b) range -template<class T, class bpred> -class out_range : public std::unary_function<T, bool> -{ -public: - out_range (const T& aa, const T& bb, const bpred& ccomp) - : a(aa), b(bb), comp(ccomp) { } - - bool operator() (const T& x) { return comp (x, a) || ! comp (x, b); } - -private: - T a; - T b; - - bpred comp; -}; - -// conveniently constructs the above functor -// NOTE: with SGI extensions, this can be written as -// compose2 (logical_and<bool>(), -// bind2nd (less<T>(), a), -// not1 (bind2nd (less<T>(), b))) -template<class T, class bpred> -inline out_range<T, bpred> -chk_out_range (const T& a, const T& b, bpred comp) -{ - return out_range<T, bpred> (a, b, comp); -} - -template<typename T, typename bpred> -inline void -seq_lookup (const T *table, octave_idx_type offset, octave_idx_type size, - const T *vals, octave_idx_type nvals, - octave_idx_type *idx, bpred comp) -{ - const T *begin = table + offset; - - if (size == 0) - // the trivial case of empty table - std::fill_n (idx, nvals, offset); - else - { - const T *vcur = vals; - const T *vend = vals + nvals; - - const T *cur = begin; - const T *end = begin + size; - - while (vcur < vend) - { - // determine the enclosing interval for next value, trying - // ++cur as a special case; - if (cur == end || comp (*vcur, *cur)) - cur = std::upper_bound (begin, cur, *vcur, comp); - else - { - ++cur; - if (cur < end && ! comp (*vcur, *cur)) - cur = std::upper_bound (cur + 1, end, *vcur, comp); - } - - // store index of the current interval. - *(idx++) = (cur - table); - ++vcur; - - // find first value not in current subrange - const T *vnew; - if (cur < end) - if (cur > begin) - // inner interval - vnew = std::find_if (vcur, vend, - chk_out_range (*(cur-1), *cur, comp)); - - else - // special case: lowermost range (-Inf, min) - vnew = std::find_if (vcur, vend, - std::not1 (std::bind2nd (comp, *cur))); - else - // special case: uppermost range [max, Inf) - vnew = std::find_if (vcur, vend, - std::bind2nd (comp, *(cur-1))); - - // store index of the current interval. - std::fill_n (idx, vnew - vcur, cur - table); - idx += (vnew - vcur); - vcur = vnew; - - } - } -} - -// overload using < operator -template<typename T, typename bpred> -inline void -seq_lookup (const T *table, octave_idx_type offset, octave_idx_type size, - const T *vals, octave_idx_type nvals, - octave_idx_type *idx) -{ - seq_lookup (table, offset, size, vals, nvals, idx, std::less<T>()); -} - -// helper functions - determine whether an array is descending -template<typename T> -inline bool -is_descending (const T *table, octave_idx_type size) -{ - return size > 1 && table[size-1] < table[0]; -} - -template<typename T, typename bpred> -inline bool -is_descending (const T *table, octave_idx_type size, - bpred comp) -{ - return size > 1 && comp (table[size-1], table[0]); -} - -#endif - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/liboctave/oct-sort.cc Thu Feb 19 07:34:15 2009 -0500 +++ b/liboctave/oct-sort.cc Thu Feb 19 15:19:59 2009 +0100 @@ -96,6 +96,7 @@ #include <cassert> #include <algorithm> +#include <functional> #include <cstring> #include <stack> @@ -105,15 +106,22 @@ #include "oct-locbuf.h" template <class T> -octave_sort<T>::octave_sort (void) : compare (ascending_compare) +octave_sort<T>::octave_sort (void) : + compare (ascending_compare), ms (0) { - merge_init (); } template <class T> -octave_sort<T>::octave_sort (compare_fcn_type comp) : compare (comp) +octave_sort<T>::octave_sort (compare_fcn_type comp) + : compare (comp), ms (0) { - merge_init (); +} + +template <class T> +octave_sort<T>::~octave_sort () +{ + merge_freemem (); + delete ms; } template <class T> @@ -476,11 +484,12 @@ void octave_sort<T>::merge_init (void) { - ms.a = 0; - ms.ia = 0; - ms.alloced = 0; - ms.n = 0; - ms.min_gallop = MIN_GALLOP; + if (! ms) ms = new MergeState; + ms->a = 0; + ms->ia = 0; + ms->alloced = 0; + ms->n = 0; + ms->min_gallop = MIN_GALLOP; } /* Free all the temp memory owned by the MergeState. This must be called @@ -491,10 +500,13 @@ void octave_sort<T>::merge_freemem (void) { - delete [] ms.a; - delete [] ms.ia; - ms.alloced = 0; - ms.a = 0; + if (ms) + { + delete [] ms->a; + delete [] ms->ia; + ms->alloced = 0; + ms->a = 0; + } } static inline octave_idx_type @@ -541,7 +553,7 @@ int octave_sort<T>::merge_getmem (octave_idx_type need) { - if (need <= ms.alloced) + if (need <= ms->alloced) return 0; need = roundupsize (need); @@ -549,10 +561,10 @@ * we don't care what's in the block. */ merge_freemem (); - ms.a = new T[need]; - if (ms.a) + ms->a = new T[need]; + if (ms->a) { - ms.alloced = need; + ms->alloced = need; return 0; } merge_freemem (); /* reset to sane state */ @@ -564,7 +576,7 @@ int octave_sort<T>::merge_getmemi (octave_idx_type need) { - if (need <= ms.alloced && ms.a && ms.ia) + if (need <= ms->alloced && ms->a && ms->ia) return 0; need = roundupsize (need); @@ -572,11 +584,11 @@ * we don't care what's in the block. */ merge_freemem (); - ms.a = new T[need]; - ms.ia = new octave_idx_type[need]; - if (ms.a && ms.ia) + ms->a = new T[need]; + ms->ia = new octave_idx_type[need]; + if (ms->a && ms->ia) { - ms.alloced = need; + ms->alloced = need; return 0; } merge_freemem (); /* reset to sane state */ @@ -600,13 +612,13 @@ octave_idx_type k; T *dest; int result = -1; /* guilty until proved innocent */ - octave_idx_type min_gallop = ms.min_gallop; + octave_idx_type min_gallop = ms->min_gallop; if (merge_getmem (na) < 0) return -1; - std::copy (pa, pa + na, ms.a); + std::copy (pa, pa + na, ms->a); dest = pa; - pa = ms.a; + pa = ms->a; *dest++ = *pb++; --nb; @@ -662,7 +674,7 @@ do { min_gallop -= min_gallop > 1; - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; k = gallop_right (*pb, pa, na, 0, comp); acount = k; if (k) @@ -706,7 +718,7 @@ while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP); ++min_gallop; /* penalize it for leaving galloping mode */ - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; } Succeed: @@ -736,14 +748,14 @@ T *dest; octave_idx_type *idest; int result = -1; /* guilty until proved innocent */ - octave_idx_type min_gallop = ms.min_gallop; + octave_idx_type min_gallop = ms->min_gallop; if (merge_getmemi (na) < 0) return -1; - std::copy (pa, pa + na, ms.a); - std::copy (ipa, ipa + na, ms.ia); + std::copy (pa, pa + na, ms->a); + std::copy (ipa, ipa + na, ms->ia); dest = pa; idest = ipa; - pa = ms.a; ipa = ms.ia; + pa = ms->a; ipa = ms->ia; *dest++ = *pb++; *idest++ = *ipb++; --nb; @@ -796,7 +808,7 @@ do { min_gallop -= min_gallop > 1; - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; k = gallop_right (*pb, pa, na, 0, comp); acount = k; if (k) @@ -842,7 +854,7 @@ while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP); ++min_gallop; /* penalize it for leaving galloping mode */ - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; } Succeed: @@ -883,15 +895,15 @@ T *dest; int result = -1; /* guilty until proved innocent */ T *basea, *baseb; - octave_idx_type min_gallop = ms.min_gallop; + octave_idx_type min_gallop = ms->min_gallop; if (merge_getmem (nb) < 0) return -1; dest = pb + nb - 1; - std::copy (pb, pb + nb, ms.a); + std::copy (pb, pb + nb, ms->a); basea = pa; - baseb = ms.a; - pb = ms.a + nb - 1; + baseb = ms->a; + pb = ms->a + nb - 1; pa += na - 1; *dest-- = *pa--; @@ -944,7 +956,7 @@ do { min_gallop -= min_gallop > 1; - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; k = gallop_right (*pb, basea, na, na-1, comp); if (k < 0) goto Fail; @@ -989,7 +1001,7 @@ goto Succeed; } while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP); ++min_gallop; /* penalize it for leaving galloping mode */ - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; } Succeed: @@ -1022,17 +1034,17 @@ int result = -1; /* guilty until proved innocent */ T *basea, *baseb; octave_idx_type *ibasea, *ibaseb; - octave_idx_type min_gallop = ms.min_gallop; + octave_idx_type min_gallop = ms->min_gallop; if (merge_getmemi (nb) < 0) return -1; dest = pb + nb - 1; idest = ipb + nb - 1; - std::copy (pb, pb + nb, ms.a); - std::copy (ipb, ipb + nb, ms.ia); + std::copy (pb, pb + nb, ms->a); + std::copy (ipb, ipb + nb, ms->ia); basea = pa; ibasea = ipa; - baseb = ms.a; ibaseb = ms.ia; - pb = ms.a + nb - 1; ipb = ms.ia + nb - 1; + baseb = ms->a; ibaseb = ms->ia; + pb = ms->a + nb - 1; ipb = ms->ia + nb - 1; pa += na - 1; ipa += na - 1; *dest-- = *pa--; *idest-- = *ipa--; @@ -1085,7 +1097,7 @@ do { min_gallop -= min_gallop > 1; - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; k = gallop_right (*pb, basea, na, na-1, comp); if (k < 0) goto Fail; @@ -1132,7 +1144,7 @@ goto Succeed; } while (acount >= MIN_GALLOP || bcount >= MIN_GALLOP); ++min_gallop; /* penalize it for leaving galloping mode */ - ms.min_gallop = min_gallop; + ms->min_gallop = min_gallop; } Succeed: @@ -1169,19 +1181,19 @@ octave_idx_type na, nb; octave_idx_type k; - pa = data + ms.pending[i].base; - na = ms.pending[i].len; - pb = data + ms.pending[i+1].base; - nb = ms.pending[i+1].len; + pa = data + ms->pending[i].base; + na = ms->pending[i].len; + pb = data + ms->pending[i+1].base; + nb = ms->pending[i+1].len; /* Record the length of the combined runs; if i is the 3rd-last * run now, also slide over the last run (which isn't involved * in this merge). The current run i+1 goes away in any case. */ - ms.pending[i].len = na + nb; - if (i == ms.n - 3) - ms.pending[i+1] = ms.pending[i+2]; - --ms.n; + ms->pending[i].len = na + nb; + if (i == ms->n - 3) + ms->pending[i+1] = ms->pending[i+2]; + ms->n--; /* Where does b start in a? Elements in a before that can be * ignored (already in place). @@ -1221,21 +1233,21 @@ octave_idx_type na, nb; octave_idx_type k; - pa = data + ms.pending[i].base; - ipa = idx + ms.pending[i].base; - na = ms.pending[i].len; - pb = data + ms.pending[i+1].base; - ipb = idx + ms.pending[i+1].base; - nb = ms.pending[i+1].len; + pa = data + ms->pending[i].base; + ipa = idx + ms->pending[i].base; + na = ms->pending[i].len; + pb = data + ms->pending[i+1].base; + ipb = idx + ms->pending[i+1].base; + nb = ms->pending[i+1].len; /* Record the length of the combined runs; if i is the 3rd-last * run now, also slide over the last run (which isn't involved * in this merge). The current run i+1 goes away in any case. */ - ms.pending[i].len = na + nb; - if (i == ms.n - 3) - ms.pending[i+1] = ms.pending[i+2]; - --ms.n; + ms->pending[i].len = na + nb; + if (i == ms->n - 3) + ms->pending[i+1] = ms->pending[i+2]; + ms->n--; /* Where does b start in a? Elements in a before that can be * ignored (already in place). @@ -1279,11 +1291,11 @@ int octave_sort<T>::merge_collapse (T *data, Comp comp) { - struct s_slice *p = ms.pending; + struct s_slice *p = ms->pending; - while (ms.n > 1) + while (ms->n > 1) { - octave_idx_type n = ms.n - 2; + octave_idx_type n = ms->n - 2; if (n > 0 && p[n-1].len <= p[n].len + p[n+1].len) { if (p[n-1].len < p[n+1].len) @@ -1308,11 +1320,11 @@ int octave_sort<T>::merge_collapse (T *data, octave_idx_type *idx, Comp comp) { - struct s_slice *p = ms.pending; + struct s_slice *p = ms->pending; - while (ms.n > 1) + while (ms->n > 1) { - octave_idx_type n = ms.n - 2; + octave_idx_type n = ms->n - 2; if (n > 0 && p[n-1].len <= p[n].len + p[n+1].len) { if (p[n-1].len < p[n+1].len) @@ -1342,11 +1354,11 @@ int octave_sort<T>::merge_force_collapse (T *data, Comp comp) { - struct s_slice *p = ms.pending; + struct s_slice *p = ms->pending; - while (ms.n > 1) + while (ms->n > 1) { - octave_idx_type n = ms.n - 2; + octave_idx_type n = ms->n - 2; if (n > 0 && p[n-1].len < p[n+1].len) --n; if (merge_at (n, data, comp) < 0) @@ -1361,11 +1373,11 @@ int octave_sort<T>::merge_force_collapse (T *data, octave_idx_type *idx, Comp comp) { - struct s_slice *p = ms.pending; + struct s_slice *p = ms->pending; - while (ms.n > 1) + while (ms->n > 1) { - octave_idx_type n = ms.n - 2; + octave_idx_type n = ms->n - 2; if (n > 0 && p[n-1].len < p[n+1].len) --n; if (merge_at (n, data, idx, comp) < 0) @@ -1406,8 +1418,13 @@ octave_sort<T>::sort (T *data, octave_idx_type nel, Comp comp) { /* Re-initialize the Mergestate as this might be the second time called */ - ms.n = 0; - ms.min_gallop = MIN_GALLOP; + if (ms) + { + ms->n = 0; + ms->min_gallop = MIN_GALLOP; + } + else + merge_init (); if (nel > 1) { @@ -1437,10 +1454,10 @@ n = force; } /* Push run onto pending-runs stack, and maybe merge. */ - assert (ms.n < MAX_MERGE_PENDING); - ms.pending[ms.n].base = lo; - ms.pending[ms.n].len = n; - ++ms.n; + assert (ms->n < MAX_MERGE_PENDING); + ms->pending[ms->n].base = lo; + ms->pending[ms->n].len = n; + ms->n++; if (merge_collapse (data, comp) < 0) goto fail; /* Advance to find next run. */ @@ -1462,10 +1479,6 @@ octave_sort<T>::sort (T *data, octave_idx_type *idx, octave_idx_type nel, Comp comp) { - /* Re-initialize the Mergestate as this might be the second time called */ - ms.n = 0; - ms.min_gallop = MIN_GALLOP; - if (nel > 1) { octave_idx_type nremaining = nel; @@ -1497,10 +1510,10 @@ n = force; } /* Push run onto pending-runs stack, and maybe merge. */ - assert (ms.n < MAX_MERGE_PENDING); - ms.pending[ms.n].base = lo; - ms.pending[ms.n].len = n; - ++ms.n; + assert (ms->n < MAX_MERGE_PENDING); + ms->pending[ms->n].base = lo; + ms->pending[ms->n].len = n; + ms->n++; if (merge_collapse (data, idx, comp) < 0) goto fail; /* Advance to find next run. */ @@ -1520,7 +1533,17 @@ void octave_sort<T>::sort (T *data, octave_idx_type nel) { + /* Re-initialize the Mergestate as this might be the second time called */ + if (ms) + { + ms->n = 0; + ms->min_gallop = MIN_GALLOP; + } + else + merge_init (); + merge_getmem (1024); + #ifdef INLINE_ASCENDING_SORT if (compare == ascending_compare) sort (data, nel, std::less<T> ()); @@ -1539,7 +1562,17 @@ void octave_sort<T>::sort (T *data, octave_idx_type *idx, octave_idx_type nel) { + /* Re-initialize the Mergestate as this might be the second time called */ + if (ms) + { + ms->n = 0; + ms->min_gallop = MIN_GALLOP; + } + else + merge_init (); + merge_getmemi (1024); + #ifdef INLINE_ASCENDING_SORT if (compare == ascending_compare) sort (data, idx, nel, std::less<T> ()); @@ -1761,6 +1794,188 @@ } +template <class T> template <class Comp> +octave_idx_type +octave_sort<T>::lookup (const T *data, octave_idx_type nel, + const T& value, Comp comp) +{ + return std::upper_bound (data, data + nel, value, comp) - data; +} + +template <class T> +octave_idx_type +octave_sort<T>::lookup (const T *data, octave_idx_type nel, + const T& value) +{ + octave_idx_type retval = 0; + +#ifdef INLINE_ASCENDING_SORT + if (compare == ascending_compare) + retval = lookup (data, nel, value, std::less<T> ()); + else +#endif +#ifdef INLINE_DESCENDING_SORT + if (compare == descending_compare) + retval = lookup (data, nel, value, std::greater<T> ()); + else +#endif + if (compare) + retval = lookup (data, nel, value, std::ptr_fun (compare)); + + return retval; +} + +// a unary functor that checks whether a value is outside [a,b) range +template<class T, class Comp> +class out_of_range_pred : public std::unary_function<T, bool> +{ +public: + out_of_range_pred (const T& aa, const T& bb, Comp c) + : a (aa), b (bb), comp (c) { } + bool operator () (const T& x) { return comp (x, a) || ! comp (x, b); } + +private: + T a, b; + Comp comp; +}; + +// a unary functor that checks whether a value is < a +template<class T, class Comp> +class less_than_pred : public std::unary_function<T, bool> +{ + typedef typename ref_param<T>::type param_type; +public: + less_than_pred (param_type aa, Comp c) + : a (aa), comp (c) { } + bool operator () (const T& x) { return comp (x, a); } + +private: + T a; + Comp comp; +}; + +// a unary functor that checks whether a value is >= a +template<class T, class Comp> +class greater_or_equal_pred : public std::unary_function<T, bool> +{ +public: + greater_or_equal_pred (const T& aa, Comp c) + : a (aa), comp (c) { } + bool operator () (const T& x) { return ! comp (x, a); } + +private: + T a; + Comp comp; +}; + +// conveniently constructs the above functors. +// NOTE: with SGI extensions, this one can be written as +// compose2 (logical_and<bool>(), bind2nd (less<T>(), a), +// not1 (bind2nd (less<T>(), b))) +template<class T, class Comp> +inline out_of_range_pred<T, Comp> +out_of_range (const T& a, + const T& b, Comp comp) +{ + return out_of_range_pred<T, Comp> (a, b, comp); +} + +// Note: these could be written as +// std::not1 (std::bind2nd (comp, *cur)) +// and +// std::bind2nd (comp, *(cur-1))); +// but that doesn't work for functions with reference parameters in g++ 4.3. +template<class T, class Comp> +inline less_than_pred<T, Comp> +less_than (const T& a, Comp comp) +{ + return less_than_pred<T, Comp> (a, comp); +} +template<class T, class Comp> +inline greater_or_equal_pred<T, Comp> +greater_or_equal (const T& a, Comp comp) +{ + return greater_or_equal_pred<T, Comp> (a, comp); +} + + +template <class T> template <class Comp> +void +octave_sort<T>::lookup (const T *data, octave_idx_type nel, + const T *values, octave_idx_type nvalues, + octave_idx_type *idx, octave_idx_type offset, Comp comp) +{ + if (nel == 0) + // the trivial case of empty table + std::fill_n (idx, nvalues, offset); + else + { + const T *vcur = values; + const T *vend = values + nvalues; + + const T *cur = data; + const T *end = data + nel; + + while (vcur != vend) + { + // determine the enclosing interval for next value, trying + // ++cur as a special case; + if (cur == end || comp (*vcur, *cur)) + cur = std::upper_bound (data, cur, *vcur, comp); + else + { + ++cur; + if (cur != end && ! comp (*vcur, *cur)) + cur = std::upper_bound (cur + 1, end, *vcur, comp); + } + + octave_idx_type vidx = cur - data + offset; + // store index of the current interval. + *(idx++) = vidx; + ++vcur; + + // find first value not in current subrange + const T *vnew; + if (cur != end) + if (cur != data) + // inner interval + vnew = std::find_if (vcur, vend, + out_of_range (*(cur-1), *cur, comp)); + else + // special case: lowermost range (-Inf, min) + vnew = std::find_if (vcur, vend, greater_or_equal (*cur, comp)); + else + // special case: uppermost range [max, Inf) + vnew = std::find_if (vcur, vend, less_than (*(cur-1), comp)); + + // store index of the current interval. + std::fill_n (idx, vnew - vcur, vidx); + idx += (vnew - vcur); + vcur = vnew; + } + } +} + +template <class T> +void +octave_sort<T>::lookup (const T *data, octave_idx_type nel, + const T* values, octave_idx_type nvalues, + octave_idx_type *idx, octave_idx_type offset) +{ +#ifdef INLINE_ASCENDING_SORT + if (compare == ascending_compare) + lookup (data, nel, values, nvalues, idx, offset, std::less<T> ()); + else +#endif +#ifdef INLINE_DESCENDING_SORT + if (compare == descending_compare) + lookup (data, nel, values, nvalues, idx, offset, std::greater<T> ()); + else +#endif + if (compare) + lookup (data, nel, values, nvalues, idx, offset, std::ptr_fun (compare)); +} + template <class T> bool octave_sort<T>::ascending_compare (typename ref_param<T>::type x,
--- a/liboctave/oct-sort.h Thu Feb 19 07:34:15 2009 -0500 +++ b/liboctave/oct-sort.h Thu Feb 19 15:19:59 2009 +0100 @@ -114,7 +114,7 @@ octave_sort (compare_fcn_type); - ~octave_sort (void) { merge_freemem (); } + ~octave_sort (void); void set_compare (compare_fcn_type comp) { compare = comp; } @@ -138,6 +138,16 @@ bool is_sorted_rows (const T *data, octave_idx_type rows, octave_idx_type cols); + // Do a binary lookup in a sorted array. + octave_idx_type lookup (const T *data, octave_idx_type nel, + const T& value); + + // Ditto, but for an array of values, specializing on long runs. + // Adds offset to all indices. + void lookup (const T *data, octave_idx_type nel, + const T* values, octave_idx_type nvalues, + octave_idx_type *idx, octave_idx_type offset = 0); + static bool ascending_compare (typename ref_param<T>::type, typename ref_param<T>::type); @@ -187,7 +197,7 @@ compare_fcn_type compare; - MergeState ms; + MergeState *ms; template <class Comp> @@ -277,6 +287,15 @@ bool is_sorted_rows (const T *data, octave_idx_type rows, octave_idx_type cols, Comp comp); + template <class Comp> + octave_idx_type lookup (const T *data, octave_idx_type nel, + const T& value, Comp comp); + + template <class Comp> + void lookup (const T *data, octave_idx_type nel, + const T* values, octave_idx_type nvalues, + octave_idx_type *idx, octave_idx_type offset, Comp comp); + }; template <class T>
--- a/src/ChangeLog Thu Feb 19 07:34:15 2009 -0500 +++ b/src/ChangeLog Thu Feb 19 15:19:59 2009 +0100 @@ -1,3 +1,8 @@ +2009-02-19 Jaroslav Hajek <highegg@gmail.com> + + * DLD-FUNCTIONS/lookup.cc (Flookup): Use Array<T>::lookup if possible. + Do not compare octave_values directly. Properly check for iscellstr. + 2009-02-19 John W. Eaton <jwe@octave.org> * data.cc, graphics.cc, help.cc, lex.l, load-path.cc, parse.y:
--- a/src/DLD-FUNCTIONS/lookup.cc Thu Feb 19 07:34:15 2009 -0500 +++ b/src/DLD-FUNCTIONS/lookup.cc Thu Feb 19 15:19:59 2009 +0100 @@ -32,7 +32,6 @@ #include "dNDArray.h" #include "CNDArray.h" -#include "oct-lookup.h" #include "Cell.h" #include "defun-dld.h" @@ -49,20 +48,6 @@ || str.find (std::toupper (c)) != std::string::npos); } -// normal ascending comparator -static bool -ov_str_less (const octave_value& a, const octave_value& b) -{ - return a.string_value () < b.string_value (); -} - -// normal descending comparator -static bool -ov_str_greater (const octave_value& a, const octave_value& b) -{ - return a.string_value () > b.string_value (); -} - // case-insensitive character comparison functors struct icmp_char_lt : public std::binary_function<char, char, bool> { @@ -76,30 +61,47 @@ { return std::toupper (x) > std::toupper (y); } }; +// FIXME: maybe these should go elsewhere? // case-insensitive ascending comparator static bool -ov_stri_less (const octave_value& a, const octave_value& b) +stri_comp_lt (const std::string& a, const std::string& b) { - std::string as = a.string_value (); - std::string bs = b.string_value (); - - return std::lexicographical_compare (as.begin (), as.end (), - bs.begin (), bs.end (), + return std::lexicographical_compare (a.begin (), a.end (), + b.begin (), b.end (), icmp_char_lt()); } // case-insensitive descending comparator static bool -ov_stri_greater (const octave_value& a, const octave_value& b) +stri_comp_gt (const std::string& a, const std::string& b) { - std::string as = a.string_value (); - std::string bs = b.string_value (); - - return std::lexicographical_compare (as.begin (), as.end (), - bs.begin (), bs.end (), + return std::lexicographical_compare (a.begin (), a.end (), + b.begin (), b.end (), icmp_char_gt()); } +template <class T> +inline sortmode +get_sort_mode (const Array<T>& array, + typename octave_sort<T>::compare_fcn_type desc_comp + = octave_sort<T>::descending_compare) +{ + octave_idx_type n = array.numel (); + if (n > 1 && desc_comp (array (0), array (n-1))) + return DESCENDING; + else + return ASCENDING; +} + +// FIXME: perhaps there should be octave_value::lookup? +// The question is, how should it behave w.r.t. the second argument's type. +// We'd need a dispatch on two arguments. Hmmm... + +#define INT_ARRAY_LOOKUP(TYPE) \ + (table.is_ ## TYPE ## _type () && y.is_ ## TYPE ## _type ()) \ + idx = table.TYPE ## _array_value ().lookup (y.TYPE ## _array_value (), \ + UNSORTED, left_inf, right_inf); + DEFUN_DLD (lookup, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {@var{idx} =} lookup (@var{table}, @var{y}, @var{opt})\n\ @@ -118,14 +120,15 @@ \n\ The algorithm used by lookup is standard binary search, with optimizations\n\ to speed up the case of partially ordered arrays (dense downsampling).\n\ -In particular, looking up a single entry is of binary complexity.\n\ +In particular, looking up a single entry is of logarithmic complexity\n\ +(unless a conversion occurs due to non-numeric or unequal types).\n\ \n\ -@var{table} and @var{y} can also be a cell array of strings\n\ +@var{table} and @var{y} can also be cell arrays of strings\n\ (or @var{y} can be a single string). In this case, string lookup\n\ is performed using lexicographical comparison.\n\ +\n\ If @var{opts} is specified, it shall be a string with letters indicating\n\ additional options.\n\ -\n\ For numeric lookup, 'l' in @var{opts} indicates that\n\ the leftmost subinterval shall be extended to infinity (i.e. all indices\n\ at least 1), and 'r' indicates that the rightmost subinterval shall be\n\ @@ -144,12 +147,12 @@ return retval; } - octave_value argtable = args(0), argy = args(1); - if (argtable.ndims () > 2 || (argtable.columns () > 1 && argtable.rows () > 1)) + octave_value table = args(0), y = args(1); + if (table.ndims () > 2 || (table.columns () > 1 && table.rows () > 1)) warning ("lookup: table is not a vector"); - bool num_case = argtable.is_numeric_type () && argy.is_numeric_type (); - bool str_case = argtable.is_cell () && (argy.is_cell () || argy.is_string ()); + bool num_case = table.is_numeric_type () && y.is_numeric_type (); + bool str_case = table.is_cellstr () && (y.is_string () || y.is_cellstr ()); if (num_case) { @@ -163,74 +166,42 @@ right_inf = contains_char (opt, 'r'); } - // in the case of a complex array, absolute values will be used for compatibility + // In the case of a complex array, absolute values will be used for compatibility // (though it's not too meaningful). - ArrayN<octave_idx_type> idx; - - if (argtable.is_single_type () || argy.is_single_type ()) - { - FloatNDArray table = (argtable.is_complex_type ()) - ? argtable.float_complex_array_value ().abs () - : argtable.float_array_value (); + + if (table.is_complex_type ()) + table = table.abs (); - FloatNDArray y = (argy.is_complex_type ()) - ? argy.float_complex_array_value ().abs () - : argy.float_array_value (); - - idx = ArrayN<octave_idx_type> (y.dims ()); + if (y.is_complex_type ()) + y = y.abs (); - // determine whether the array is descending. - bool desc = is_descending (table.data (), table.length ()); - octave_idx_type offset = left_inf ? 1 : 0; - octave_idx_type size = table.length () - offset - (right_inf ? 1 : 0); - if (size < 0) - size = 0; + Array<octave_idx_type> idx; - if (desc) - seq_lookup (table.data (), offset, size, - y.data (), y.length (), idx.fortran_vec (), - std::greater<float> ()); - else - seq_lookup (table.data (), offset, size, - y.data (), y.length (), idx.fortran_vec (), - std::less<float> ()); - } + // PS: I learned this from data.cc + if INT_ARRAY_LOOKUP (int8) + else if INT_ARRAY_LOOKUP (int16) + else if INT_ARRAY_LOOKUP (int32) + else if INT_ARRAY_LOOKUP (int64) + else if INT_ARRAY_LOOKUP (uint8) + else if INT_ARRAY_LOOKUP (uint16) + else if INT_ARRAY_LOOKUP (uint32) + else if INT_ARRAY_LOOKUP (uint64) + else if (table.is_single_type () || y.is_single_type ()) + idx = table.float_array_value ().lookup (y.float_array_value (), + UNSORTED, left_inf, right_inf); else - { - NDArray table = (argtable.is_complex_type ()) - ? argtable.complex_array_value ().abs () - : argtable.array_value (); - - NDArray y = (argy.is_complex_type ()) - ? argy.complex_array_value ().abs () - : argy.array_value (); - - idx = ArrayN<octave_idx_type> (y.dims ()); - - // determine whether the array is descending. - bool desc = is_descending (table.data (), table.length ()); - octave_idx_type offset = left_inf ? 1 : 0; - octave_idx_type size = table.length () - offset - (right_inf ? 1 : 0); - if (size < 0) - size = 0; - - if (desc) - seq_lookup (table.data (), offset, size, - y.data (), y.length (), idx.fortran_vec (), - std::greater<double> ()); - else - seq_lookup (table.data (), offset, size, - y.data (), y.length (), idx.fortran_vec (), - std::less<double> ()); - } + idx = table.array_value ().lookup (y.array_value (), + UNSORTED, left_inf, right_inf); retval(0) = NDArray (idx); + } else if (str_case) { - Cell table = argtable.cell_value (); + Array<std::string> str_table = table.cellstr_value (); - bool (*ov_str_comp) (const octave_value&, const octave_value&); + // Here we'll use octave_sort directly to avoid converting the array + // for case-insensitive comparison. bool icase = false; @@ -241,51 +212,40 @@ icase = contains_char (opt, 'i'); } + sortmode mode = (icase ? get_sort_mode (str_table, stri_comp_gt) + : get_sort_mode (str_table)); + + bool (*str_comp) (const std::string&, const std::string&); + // pick the correct comparator - if (icase) - { - if (is_descending (table.data (), table.length (), ov_stri_less)) - ov_str_comp = ov_stri_greater; - else - ov_str_comp = ov_stri_less; - } + if (mode == DESCENDING) + str_comp = icase ? stri_comp_gt : octave_sort<std::string>::descending_compare; else - { - if (is_descending (table.data (), table.length (), ov_str_less)) - ov_str_comp = ov_str_greater; - else - ov_str_comp = ov_str_less; - } + str_comp = icase ? stri_comp_lt : octave_sort<std::string>::ascending_compare; - - // query just the first cell to verify it's a string - if (table.is_empty () || table(0).is_string ()) + octave_sort<std::string> lsort (str_comp); + if (y.is_cellstr ()) { - if (argy.is_cell ()) - { - Cell y = argy.cell_value (); - ArrayN<octave_idx_type> idx (y.dims ()); + Array<std::string> str_y = y.cellstr_value (); - + Array<octave_idx_type> idx (str_y.dims ()); - for (int i = 0; i < y.numel (); i++) - idx(i) = bin_lookup (table.data (), table.length (), y(i), - std::ptr_fun (ov_str_comp)); + lsort.lookup (str_table.data (), str_table.nelem (), str_y.data (), + str_y.nelem (), idx.fortran_vec ()); - retval(0) = NDArray (idx); - } - else - { - octave_idx_type idx; + retval(0) = NDArray (idx); + } + else if (y.is_string ()) + { + std::string str_y = y.string_value (); - idx = bin_lookup (table.data (), table.length (), argy, - std::ptr_fun (ov_str_comp)); + octave_idx_type idx; - retval(0) = static_cast<double> (idx); - } + lsort.lookup (str_table.data (), str_table.nelem (), &str_y, + 1, &idx); + + retval(0) = idx; } - else - error("lookup: table is not a cell array of strings."); } else print_usage ();