/*

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 <address@hidden>


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "NDArray.h"

#include "defun-dld.h"
#include "error.h"
#include "gripes.h"
#include "oct-obj.h"


// change these to select appropriate output index type.
// TODO: if we want to return integer indices someday,
// it should be octave_idx_type. The octave_value class
// can not, however, hold an ArrayN<octave_idx_type> yet.

// typedef octave_idx_type ret_idx_type;
// typedef ???? ret_idx_array;

typedef double ret_idx_type;
typedef NDArray ret_idx_array;


// simple function templates & specialization is used
// to avoid writing two versions of the algorithm

// generic comparison
template<bool rev>
inline bool dcomp (const double& a, const double& b);

template<bool rev>
inline bool dcompe (const double& a, const double& b);

// normal comparison (increasing table)
template<>
inline bool dcomp<false>(const double& a, const double& b)
{ 
  return a < b; 
}
// normal comparison with equality
template<>
inline bool dcompe<false>(const double& a, const double& b)
{ 
  return a <= b; 
}
// reverse comparison (decreasing table)
template<>
inline bool dcomp<true>(const double& a, const double& b)
{ 
  return a > b; 
}
// reverse comparison with equality
template<>
inline bool dcompe<true>(const double& a, const double& b)
{ 
  return a >= b; 
}


// the heart of the lookup. The [begin,end) range is assumed
// sorted so that dcompe<rev>(*iter,*(iter+1)) is true.
//
// outline of the algorithm:
// 1. start with full range
// 2. determine interval enclosing the next value
// 3. fetch values while they stay inside the interval
//    if escaped to the left, go to 4. if escaped to the right,
//    go to 5. If processed all values, return.
// 4. enclose the next value by binary bracketing to the left.
//    if just one step was taken (simple shift), go to 3.
//    otherwise, go to 2.
// 5. enclose the next value by binary bracketing to the right.
//    if just one step was taken (simple shift), go to 3.
//    otherwise, go to 2.

template<bool rev>
static void 
do_seq_lookup (const double *begin, const double *end,
               const double *vals, size_t nvals, ret_idx_type *idx)
{
  if (!nvals) 
    return;

  // these are "global" pointers (i.e. used by all four blocks)
  const double *first, *last;

  // the trivial case of empty array. Store all zeros and return.
  if (begin == end) 
    {
      for (;nvals;--nvals) 
        *(idx++) = 0;
      return;
    }

  // initialize 
  last = end;
  first = begin;

  // perform a binary search in [first,last)
bin_search:    
    {
      size_t dist = last - first, half;

      while (dist > 0) 
        {
          half = dist >> 1;
          last = first + half;
          if (dcomp<rev> (*last, *vals)) 
            {
              first = last + 1;
              dist -= (half + 1);
            }
          else
            dist = half;
        }
      last = first;
    }

  // fetch values as long as they stay in the current interval.
  // once they get out, jump to the appropriate bracketing block.
store_cur:
    {
      size_t cur;
      cur = last - begin;

      if (last == begin) 
        {
          // leftmost interval (can't escape backwards)
          while (nvals) 
            {
              if (dcompe<rev> (*last, *vals)) 
                goto search_forwards;

              *(idx++) = cur;
              nvals--; vals++;
            }
          return;
        } 
      else if (last == end) 
        {
          // rightmost interval (can't escape forwards)
          first = last - 1;
          while (nvals) 
            {
              if (dcomp<rev> (*vals, *first)) 
                goto search_backwards;

              *(idx++) = cur;
              nvals--; vals++;
            }
          return;
        } 
      else 
        {
          // inner interval (can escape both ways)
          first = last - 1;
          while (nvals) 
            {
              if (dcompe<rev> (*last, *vals)) 
                goto search_forwards;

              if (dcomp<rev> (*vals, *first)) 
                goto search_backwards;

              *(idx++) = cur;
              nvals--; vals++;
            }
          return;
        }

    }

  // the current value got to the right of current interval.
  // perform binary bracketing to the right
search_forwards:
    {
      size_t dist = 1, max = end - last;
      while (true) 
        {
          first = last;
          if (dist < max) 
            last += dist;
          else 
            {
              last = end;
              break;
            }
          if (dcomp<rev> (*vals, *last)) 
            break;
          max -= dist;
          dist <<= 1;
        }

      if (!(--dist))
        goto store_cur; // a shortcut. If dist was 1, bin_search is not needed.
      else
        goto bin_search;
    }

  // the current value got to the left of current interval.
  // perform binary bracketing to the left.
search_backwards:
    {
      size_t dist = 1, max = first - begin;
      while (true) 
        {
          last = first;
          if (dist < max) 
            first -= dist;
          else 
            {
              first = begin;
              break;
            }
          if (dcompe<rev> (*first, *vals)) 
            break;
          max -= dist;
          dist <<= 1;
        }

      if (!(--dist))
        goto store_cur; // a shortcut. If dist was 1, bin_search is not needed.
      else
        goto bin_search;
    }

}

// this determines the orientation of the table and calls the appropriate
// instance of do_seq_lookup

static void seq_lookup (const double *table, size_t ntable,
                        const double *vals, size_t nvals, ret_idx_type *idx)
{
  if (ntable > 1 && table[0] > table[ntable-1]) 
    do_seq_lookup<true>(table, table+ntable, vals, nvals, idx);
  else
    do_seq_lookup<false>(table, table+ntable, vals, nvals, idx);
}

DEFUN_DLD (lookup, args, nargout, "\
-*- texinfo -*-\n\
@deftypefn {Function File} address@hidden =} lookup (@var{table}, @var{y})\n\
Lookup values in a sorted table.  Usually used as a prelude to\n\
interpolation.\n\
\n\
If table is strictly increasing and @code{idx = lookup (table, y)}, then\n\
@code{table(idx(i)) <= y(i) < table(idx(i+1))} for all @code{y(i)}\n\
within the table.  If @code{y(i)} is before the table, then\n\
@code{idx(i)} is 0. If @code{y(i)} is after the table then\n\
@code{idx(i)} is @code{table(n)}.\n\
\n\
If the table is strictly decreasing, then the tests are reversed.\n\
There are no guarantees for tables which are non-monotonic or are not\n\
strictly monotonic.\n\
\n\
To get an index value which lies within an interval of the table,\n\
use:\n\
\n\
@example\n\
idx = lookup (table(2:length(table)-1), y) + 1\n\
@end example\n\
\n\
@noindent\n\
This expression puts values before the table into the first\n\
interval, and values after the table into the last interval.\n\
\n\
If complex values are supplied, their magnitudes are used.\n\
@end deftypefn") 
{
  int nargin = args.length ();

  octave_value table = args (0), y = args (1);
  if (nargin == 2 && nargout <= 1 
      && (table = args(0), table.is_numeric_type ()) 
      && (y = args(1), y.is_numeric_type ()) )
    {
      if (table.ndims() == 2 && (table.rows () == 1 || table.columns () == 1))
        {
          // in the case of a complex array, absolute values will be used (though it's
          // not too meaningful).
          NDArray atable = (table.is_complex_type ()) 
            ? table.complex_array_value ().abs ()
            : table.array_value ();
          NDArray ay = (y.is_complex_type ()) 
            ? y.complex_array_value ().abs ()
            : y.array_value ();

          ret_idx_array idx (y.dims ());

          seq_lookup (atable.data (), atable.numel (), ay.data (),
                      ay.numel (), idx.fortran_vec ());

          return octave_value (idx);
        }
      else
        {
          error("lookup: table must be a vector");
          return octave_value_list ();
        }
    }
  else
    {
      print_usage ();
      return octave_value_list ();
    }
}  

/*
%!assert (real(lookup(1:3, 0.5)), 0)     # value before table
%!assert (real(lookup(1:3, 3.5)), 3)     # value after table error
%!assert (real(lookup(1:3, 1.5)), 1)     # value within table error
%!assert (real(lookup(1:3, [3,2,1])), [3,2,1])
%!assert (real(lookup([1:4]', [1.2, 3.5]')), [1, 3]');
%!assert (real(lookup([1:4], [1.2, 3.5]')), [1, 3]');
%!assert (real(lookup([1:4]', [1.2, 3.5])), [1, 3]);
%!assert (real(lookup([1:4], [1.2, 3.5])), [1, 3]);
%!assert (real(lookup(1:3, [3, 2, 1])), [3, 2, 1]);
%!assert (real(lookup([3:-1:1], [3.5, 3, 1.2, 2.5, 2.5])), [0, 1, 2, 1, 1])
%!assert (isempty(lookup([1:3], [])))
%!assert (isempty(lookup([1:3]', [])))
%!assert (real(lookup(1:3, [1, 2; 3, 0.5])), [1, 2; 3, 0]);
*/