l_algorithm.h

Go to the documentation of this file.

/* $Id: l_algorithm.h 18301 2014-11-26 19:17:13Z ksimek $ */
/* {{{=========================================================================== *
   |
   |  Copyright (c) 1994-2011 by Kobus Barnard (author)
   |
   |  Personal and educational use of this code is granted, provided that this
   |  header is kept intact, and that the authorship is not misrepresented, that
   |  its use is acknowledged in publications, and relevant papers are cited.
   |
   |  For other use contact the author (kobus AT cs DOT arizona DOT edu).
   |
   |  Please note that the code in this file has not necessarily been adequately
   |  tested. Naturally, there is no guarantee of performance, support, or fitness
   |  for any particular task. Nonetheless, I am interested in hearing about
   |  problems that you encounter.
   |
   |  Author:  Kyle Simek
 * =========================================================================== }}}*/

// vim: tabstop=4 shiftwidth=4 foldmethod=marker

#ifndef KJB_L_CPP_ALGORITHM_H
#define KJB_L_CPP_ALGORITHM_H

#include <l_cpp/l_functors.h>
#include <l_cpp/l_std_parallel.h>
#include <l_cpp/l_exception.h>

#include <map>
#include <iterator>
#include <algorithm>
#include <assert.h>
#include <cfloat>
#include <algorithm>
#include <vector>
#include <utility>

//#include <boost/lambda/lambda.hpp>
#include <l_cpp/l_functors.h>
//#include <boost/bind.hpp>
//#include <boost/function.hpp>
#include <boost/concept_check.hpp>
//#include <boost/iterator/transform_iterator.hpp>
//#include <boost/iterator/indirect_iterator.hpp>

namespace kjb
{

template <class Iterator, class Real>
typename std::iterator_traits<Iterator>::value_type lerp(Iterator begin, Iterator end, Real x)
{
    typedef typename std::iterator_traits<Iterator>::value_type value_type; 
    x = std::max(Real(0.0), x);
    x = std::min(x, Real(1.0));

    int max_index = std::distance(begin, end) - 1;
    int integer = x * max_index;
    Real rational = x * max_index - integer;

    assert(integer >= 0);
    assert(integer <= (int) max_index);
    assert(rational >= 0.0);
    assert(rational < 1.0);

    if(integer == (int) max_index) 
    {
        assert(rational < FLT_EPSILON);
        Iterator it = begin;
        std::advance(it, integer);
        return *it;
    }

    Iterator lower = begin;
    std::advance(lower, integer);

    Iterator upper = lower;
    ++upper;

    value_type interp = (1.0 - rational) * (*lower) + rational * (*upper);

    return interp;
}

template <class IIterator, class OIterator, class QueryType>
typename std::iterator_traits<OIterator>::value_type
lerp(IIterator ibegin, IIterator iend, OIterator obegin, const QueryType& query_point)
{
    size_t n = std::distance(ibegin, iend);
    OIterator oend = obegin + n;

    typedef typename std::iterator_traits<IIterator>::value_type InputType ;
    typedef typename std::iterator_traits<OIterator>::value_type OutputType;

    BOOST_CONCEPT_ASSERT((boost::Convertible<QueryType, InputType>));

    if(n == 0)
        KJB_THROW_2(Illegal_argument, "Function is empty.");

    IIterator iit = ibegin;
    std::advance(iit, n-1);

    const InputType& lower_bound = *ibegin;
    const InputType& upper_bound = *iit;
    // out of range
    if(!(lower_bound <= query_point && query_point <= upper_bound))
        KJB_THROW(Index_out_of_bounds);

    iit = std::upper_bound(ibegin, iend, query_point);
    const InputType& iupper = *iit;
    const InputType& ilower = *(--iit);

    size_t dist = std::distance(ibegin, iit);
    // get output lower/upper values
    OIterator oit = obegin;
    std::advance(oit, dist);
    const OutputType& olower = *oit;

    if(query_point == ilower)
        return olower;

    const OutputType& oupper = *(++oit);

    const InputType offset = query_point - ilower;
    const InputType  range_size = iupper - ilower;
    return olower + offset * (oupper - olower) / range_size;
}

template <class IIterator, class OIterator, class I2Iterator, class O2Iterator>
void ordered_lerp(
        IIterator ibegin,
        IIterator iend,
        OIterator obegin,
        I2Iterator i2begin,
        I2Iterator i2end,
        O2Iterator o2begin)
{
    size_t n = std::distance(ibegin, iend);
    OIterator oend = obegin + n;

    size_t n2 = std::distance(i2begin, i2end);
    O2Iterator o2end = o2begin + n2; 

    typedef typename std::iterator_traits<IIterator>::value_type InputType;
    typedef typename std::iterator_traits<OIterator>::value_type OutputType;
    typedef typename std::iterator_traits<I2Iterator>::value_type Input2Type;
    typedef typename std::iterator_traits<O2Iterator>::value_type Output2Type;

    BOOST_CONCEPT_ASSERT((boost::Convertible<InputType, Input2Type>));
//    BOOST_CONCEPT_ASSERT((boost::Convertible<OutputType, Output2Type>));

    if(n == 0)
        KJB_THROW_2(Illegal_argument, "Function is empty.");
    if(n2 == 0)
        return;

    IIterator iit = ibegin;
    std::advance(iit, n-1);

    IIterator i2it = i2begin;
    std::advance(i2it, n2-1);

    {
        const InputType& lower_bound = *ibegin;
        const InputType& upper_bound = *iit;

        // get range of query inputs
        const Input2Type& query_lower_bound = *i2begin;
        const Input2Type& query_upper_bound = *i2it;

        // out of range
        if(!(lower_bound <= query_lower_bound && query_upper_bound <= upper_bound))
            KJB_THROW(Index_out_of_bounds);

    }

    iit = ibegin;
    OIterator oit = obegin;

    i2it = i2begin;
    O2Iterator o2it = o2begin;

    // handle one or more query points equal to the lower bound.
    while(*i2it == *iit)
    {
        *o2it++ = *oit;
        *i2it++;
    }

    InputType* prev_input;
    OutputType* prev_output;
    Input2Type range_size;

    // fact: *i2it > *iit
    
    for(; i2it != i2end; ++i2it)
    {
        const Input2Type& query = *i2it;

        // make sure iit is the upper bound of the range we're inside
        // will always be true for the first iteration:
        if(query > *iit)
        {
            while(query > *iit)
            {
                prev_input = &*iit;
                prev_output = &*oit;
                ++iit;
                ++oit;
            }
            range_size = *iit - *prev_input;
        }

        // special case: exact hit
        if(query == *iit)
        {
            *o2it++ = *oit;
            continue;
        }

        // we are now between the lower and upper bound
        // 1. find out how far above minimum
        Input2Type distance = query - static_cast<Input2Type>(*prev_input);

        // 2. convert to percent
        distance /= range_size;

        // 3. lerp it
        *o2it++ = *prev_output + distance * (*oit - *prev_output);
    }
}

template <class Key_type, class Value_type, class Query_type>
Value_type lerp(const std::map<Key_type, Value_type>& piecewise_function, const Query_type& query_point)
{
    typedef std::map<Key_type, Value_type> Data_set;
    typedef typename Data_set::const_iterator Iterator;
    const Data_set& data_set = piecewise_function;

    BOOST_CONCEPT_ASSERT((boost::Convertible<Query_type, Key_type>));

    if(data_set.empty())
        KJB_THROW_2(Illegal_argument, "Function is empty.");

    // out of range
    if(!(data_set.begin()->first <= query_point && query_point <= data_set.rbegin()->first))
        KJB_THROW(Index_out_of_bounds);

    Iterator upper = data_set.upper_bound(query_point);
    Iterator lower = upper;
    advance(upper,  -1);

    if(query_point == lower->first)
        return lower->second;

    const Key_type offset = query_point - lower->first;
    const Key_type range_size = upper->first - lower->first;
    return lower->second + offset * (upper->second - lower->second) / range_size;
}

template <class ForwardIterator>
void linspace(double min, double max, size_t n, ForwardIterator begin, bool endpoint = true)
{
#if 0 /* This version accumulates precision error the longer the list */
    double incr;
    if(endpoint)
        incr = (max - min) / (n - 1);
    else
        incr = (max - min) / n;

    ForwardIterator end = begin;
    std::advance(end, n);
    std::generate(begin, end, Increase_by<double>(min, incr));
#else /* This version is a port of matlab's linspace, incurs no drift */
    double last = (endpoint ? n-1 : n);

    ForwardIterator it = begin;
    for(size_t i = 0; i < n-1; ++i)
        *it++ = min + i / last * (max - min);
    *it++ = max;
#endif
}

template <class ForwardIterator>
void logspace(double min, double max, size_t n, ForwardIterator begin)
{
    linspace(min, max, n, begin);

    ForwardIterator end = begin;
    std::advance(end, n);
    std::transform(begin, end, begin, std::bind1st(std::ptr_fun(static_cast<double(*)(double, double)>(pow)), 10));
}

template <class InType, class OutIterator, class UnaryOperator>
void linspace(
        const InType& lower,
        const InType& upper,
        size_t N,
        const OutIterator& out_begin,
        const UnaryOperator& f)
{
    const InType* range[2];
    range[0] = &lower;
    range[1] = &upper;

    // create real-valued inputs
    std::vector<double> x(N);
    linspace(0.0, 1.0, N, x.begin());
//
//    // create a transform function converts real-valued inputs to
//    // InType-valued inputs.
//    boost::function1<InType, double> lerp = boost::bind(
//            static_cast< InType (*)(boost::indirect_iterator<InType const**> , boost::indirect_iterator<InType const**>, double)> (&::kjb::lerp),
//            boost::indirect_iterator<InType const **>(range),
//            boost::indirect_iterator<InType const **>(range+2),
//            _1);
//
//    // evaluate function at each point
//    transform(
//        boost::make_transform_iterator(x.begin(), lerp),
//        boost::make_transform_iterator(x.end(), lerp),
//        out_begin,
//        f);
//
//    OutIterator out_it = out_begin;
//    for(size_t i = 0; i < x.size(); ++i)
//    {
//        *out_it++ = f(::kjb::lerp(lower, upper, x[i]));
//    }

    OutIterator it = out_begin;
    InType diff = upper - lower;
#pragma omp parallel for
    for(size_t i = 0; i < x.size(); ++i)
        *it++ = f(lower + x[i] * diff);
}


template <class InType, class OutIterator>
void linspace(InType min, InType max, size_t n, OutIterator begin)
{
    ::kjb::linspace(
            min,
            max,
            n,
            begin,
            Identity<InType>()); // identity function
}
} // namespace kjb


namespace kjb_parallel
{
template <class InType, class OutIterator, class UnaryOperator>
void linspace(
        const InType& lower,
        const InType& upper,
        size_t N,
        const OutIterator& out_begin,
        const UnaryOperator& f)
{
    const InType* range[2];
    range[0] = &lower;
    range[1] = &upper;

    // create real-valued inputs
    std::vector<double> x(N);
    ::kjb::linspace(0.0, 1.0, N, x.begin());
//
//    // create a transform function converts real-valued inputs to
//    // InType-valued inputs.
//    boost::function1<InType, double> lerp = boost::bind(
//            static_cast< InType (*)(boost::indirect_iterator<InType const**> , boost::indirect_iterator<InType const**>, double)> (&::kjb::lerp),
//            boost::indirect_iterator<InType const **>(range),
//            boost::indirect_iterator<InType const **>(range+2),
//            _1);
//
//    // evaluate function at each point
//    transform(
//        boost::make_transform_iterator(x.begin(), lerp),
//        boost::make_transform_iterator(x.end(), lerp),
//        out_begin,
//        f);
//
//#pragma omp parallel for
//    for(size_t i = 0; i < x.size(); ++i)
//    {
//        out_begin[i] = f(lerp(lower, upper, x[i]));
//    }

    OutIterator it = out_begin;
    InType diff = upper - lower;
#pragma omp parallel for
    for(size_t i = 0; i < x.size(); ++i)
        *it++ = f(lower + x[i] * diff);
}

} // namespace kjb_parallel

#endif