/* boost random/piecewise_linear_distribution.hpp header file
 *
 * Copyright Steven Watanabe 2011
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 *
 * See http://www.boost.org for most recent version including documentation.
 *
 * $Id$
 */

#ifndef BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED
#define BOOST_RANDOM_PIECEWISE_LINEAR_DISTRIBUTION_HPP_INCLUDED

#include <vector>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <boost/assert.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/discrete_distribution.hpp>
#include <boost/random/detail/config.hpp>
#include <boost/random/detail/operators.hpp>
#include <boost/random/detail/vector_io.hpp>

#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
#include <initializer_list>
#endif

#include <boost/range/begin.hpp>
#include <boost/range/end.hpp>

namespace lslboost {
namespace random {

/**
 * The class @c piecewise_linear_distribution models a \random_distribution.
 */
template<class RealType = double>
class piecewise_linear_distribution {
public:
    typedef std::size_t input_type;
    typedef RealType result_type;

    class param_type {
    public:

        typedef piecewise_linear_distribution distribution_type;

        /**
         * Constructs a @c param_type object, representing a distribution
         * that produces values uniformly distributed in the range [0, 1).
         */
        param_type()
        {
            _weights.push_back(RealType(1));
            _weights.push_back(RealType(1));
            _intervals.push_back(RealType(0));
            _intervals.push_back(RealType(1));
        }
        /**
         * Constructs a @c param_type object from two iterator ranges
         * containing the interval boundaries and weights at the boundaries.
         * If there are fewer than two boundaries, then this is equivalent to
         * the default constructor and the distribution will produce values
         * uniformly distributed in the range [0, 1).
         *
         * The values of the interval boundaries must be strictly
         * increasing, and the number of weights must be the same as
         * the number of interval boundaries.  If there are extra
         * weights, they are ignored.
         */
        template<class IntervalIter, class WeightIter>
        param_type(IntervalIter intervals_first, IntervalIter intervals_last,
                   WeightIter weight_first)
          : _intervals(intervals_first, intervals_last)
        {
            if(_intervals.size() < 2) {
                _intervals.clear();
                _weights.push_back(RealType(1));
                _weights.push_back(RealType(1));
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
            } else {
                _weights.reserve(_intervals.size());
                for(std::size_t i = 0; i < _intervals.size(); ++i) {
                    _weights.push_back(*weight_first++);
                }
            }
        }
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
        /**
         * Constructs a @c param_type object from an initializer_list
         * containing the interval boundaries and a unary function
         * specifying the weights at the boundaries.  Each weight is
         * determined by calling the function at the corresponding point.
         *
         * If the initializer_list contains fewer than two elements,
         * this is equivalent to the default constructor and the
         * distribution will produce values uniformly distributed
         * in the range [0, 1).
         */
        template<class T, class F>
        param_type(const std::initializer_list<T>& il, F f)
          : _intervals(il.begin(), il.end())
        {
            if(_intervals.size() < 2) {
                _intervals.clear();
                _weights.push_back(RealType(1));
                _weights.push_back(RealType(1));
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
            } else {
                _weights.reserve(_intervals.size());
                for(typename std::vector<RealType>::const_iterator
                    iter = _intervals.begin(), end = _intervals.end();
                    iter != end; ++iter)
                {
                    _weights.push_back(f(*iter));
                }
            }
        }
#endif
        /**
         * Constructs a @c param_type object from Boost.Range ranges holding
         * the interval boundaries and the weights at the boundaries.  If
         * there are fewer than two interval boundaries, this is equivalent
         * to the default constructor and the distribution will produce
         * values uniformly distributed in the range [0, 1).  The
         * number of weights must be equal to the number of
         * interval boundaries.
         */
        template<class IntervalRange, class WeightRange>
        param_type(const IntervalRange& intervals_arg,
                   const WeightRange& weights_arg)
          : _intervals(lslboost::begin(intervals_arg), lslboost::end(intervals_arg)),
            _weights(lslboost::begin(weights_arg), lslboost::end(weights_arg))
        {
            if(_intervals.size() < 2) {
                _weights.clear();
                _weights.push_back(RealType(1));
                _weights.push_back(RealType(1));
                _intervals.clear();
                _intervals.push_back(RealType(0));
                _intervals.push_back(RealType(1));
            }
        }

        /**
         * Constructs the parameters for a distribution that approximates a
         * function.  The range of the distribution is [xmin, xmax).  This
         * range is divided into nw equally sized intervals and the weights
         * are found by calling the unary function f on the boundaries of the
         * intervals.
         */
        template<class F>
        param_type(std::size_t nw, RealType xmin, RealType xmax, F f)
        {
            std::size_t n = (nw == 0) ? 1 : nw;
            double delta = (xmax - xmin) / n;
            BOOST_ASSERT(delta > 0);
            for(std::size_t k = 0; k < n; ++k) {
                _weights.push_back(f(xmin + k*delta));
                _intervals.push_back(xmin + k*delta);
            }
            _weights.push_back(f(xmax));
            _intervals.push_back(xmax);
        }

        /**  Returns a vector containing the interval boundaries. */
        std::vector<RealType> intervals() const { return _intervals; }

        /**
         * Returns a vector containing the probability densities
         * at all the interval boundaries.
         */
        std::vector<RealType> densities() const
        {
            RealType sum = static_cast<RealType>(0);
            for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
                RealType width = _intervals[i + 1] - _intervals[i];
                sum += (_weights[i] + _weights[i + 1]) * width / 2;
            }
            std::vector<RealType> result;
            result.reserve(_weights.size());
            for(typename std::vector<RealType>::const_iterator
                iter = _weights.begin(), end = _weights.end();
                iter != end; ++iter)
            {
                result.push_back(*iter / sum);
            }
            return result;
        }

        /** Writes the parameters to a @c std::ostream. */
        BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(os, param_type, parm)
        {
            detail::print_vector(os, parm._intervals);
            detail::print_vector(os, parm._weights);
            return os;
        }
        
        /** Reads the parameters from a @c std::istream. */
        BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(is, param_type, parm)
        {
            std::vector<RealType> new_intervals;
            std::vector<RealType> new_weights;
            detail::read_vector(is, new_intervals);
            detail::read_vector(is, new_weights);
            if(is) {
                parm._intervals.swap(new_intervals);
                parm._weights.swap(new_weights);
            }
            return is;
        }

        /** Returns true if the two sets of parameters are the same. */
        BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(param_type, lhs, rhs)
        {
            return lhs._intervals == rhs._intervals
                && lhs._weights == rhs._weights;
        }
        /** Returns true if the two sets of parameters are different. */
        BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(param_type)

    private:
        friend class piecewise_linear_distribution;

        std::vector<RealType> _intervals;
        std::vector<RealType> _weights;
    };

    /**
     * Creates a new @c piecewise_linear_distribution that
     * produces values uniformly distributed in the range [0, 1).
     */
    piecewise_linear_distribution()
    {
        default_init();
    }
    /**
     * Constructs a piecewise_linear_distribution from two iterator ranges
     * containing the interval boundaries and the weights at the boundaries.
     * If there are fewer than two boundaries, then this is equivalent to
     * the default constructor and creates a distribution that
     * produces values uniformly distributed in the range [0, 1).
     *
     * The values of the interval boundaries must be strictly
     * increasing, and the number of weights must be equal to
     * the number of interval boundaries.  If there are extra
     * weights, they are ignored.
     *
     * For example,
     *
     * @code
     * double intervals[] = { 0.0, 1.0, 2.0 };
     * double weights[] = { 0.0, 1.0, 0.0 };
     * piecewise_constant_distribution<> dist(
     *     &intervals[0], &intervals[0] + 3, &weights[0]);
     * @endcode
     *
     * produces a triangle distribution.
     */
    template<class IntervalIter, class WeightIter>
    piecewise_linear_distribution(IntervalIter first_interval,
                                  IntervalIter last_interval,
                                  WeightIter first_weight)
      : _intervals(first_interval, last_interval)
    {
        if(_intervals.size() < 2) {
            default_init();
        } else {
            _weights.reserve(_intervals.size());
            for(std::size_t i = 0; i < _intervals.size(); ++i) {
                _weights.push_back(*first_weight++);
            }
            init();
        }
    }
#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    /**
     * Constructs a piecewise_linear_distribution from an
     * initializer_list containing the interval boundaries
     * and a unary function specifying the weights.  Each
     * weight is determined by calling the function at the
     * corresponding interval boundary.
     *
     * If the initializer_list contains fewer than two elements,
     * this is equivalent to the default constructor and the
     * distribution will produce values uniformly distributed
     * in the range [0, 1).
     */
    template<class T, class F>
    piecewise_linear_distribution(std::initializer_list<T> il, F f)
      : _intervals(il.begin(), il.end())
    {
        if(_intervals.size() < 2) {
            default_init();
        } else {
            _weights.reserve(_intervals.size());
            for(typename std::vector<RealType>::const_iterator
                iter = _intervals.begin(), end = _intervals.end();
                iter != end; ++iter)
            {
                _weights.push_back(f(*iter));
            }
            init();
        }
    }
#endif
    /**
     * Constructs a piecewise_linear_distribution from Boost.Range
     * ranges holding the interval boundaries and the weights.  If
     * there are fewer than two interval boundaries, this is equivalent
     * to the default constructor and the distribution will produce
     * values uniformly distributed in the range [0, 1).  The
     * number of weights must be equal to the number of
     * interval boundaries.
     */
    template<class IntervalsRange, class WeightsRange>
    piecewise_linear_distribution(const IntervalsRange& intervals_arg,
                                  const WeightsRange& weights_arg)
      : _intervals(lslboost::begin(intervals_arg), lslboost::end(intervals_arg)),
        _weights(lslboost::begin(weights_arg), lslboost::end(weights_arg))
    {
        if(_intervals.size() < 2) {
            default_init();
        } else {
            init();
        }
    }
    /**
     * Constructs a piecewise_linear_distribution that approximates a
     * function.  The range of the distribution is [xmin, xmax).  This
     * range is divided into nw equally sized intervals and the weights
     * are found by calling the unary function f on the interval boundaries.
     */
    template<class F>
    piecewise_linear_distribution(std::size_t nw,
                                  RealType xmin,
                                  RealType xmax,
                                  F f)
    {
        if(nw == 0) { nw = 1; }
        RealType delta = (xmax - xmin) / nw;
        _intervals.reserve(nw + 1);
        for(std::size_t i = 0; i < nw; ++i) {
            RealType x = xmin + i * delta;
            _intervals.push_back(x);
            _weights.push_back(f(x));
        }
        _intervals.push_back(xmax);
        _weights.push_back(f(xmax));
        init();
    }
    /**
     * Constructs a piecewise_linear_distribution from its parameters.
     */
    explicit piecewise_linear_distribution(const param_type& parm)
      : _intervals(parm._intervals),
        _weights(parm._weights)
    {
        init();
    }

    /**
     * Returns a value distributed according to the parameters of the
     * piecewise_linear_distribution.
     */
    template<class URNG>
    RealType operator()(URNG& urng) const
    {
        std::size_t i = _bins(urng);
        bool is_in_rectangle = (i % 2 == 0);
        i = i / 2;
        uniform_real<RealType> dist(_intervals[i], _intervals[i+1]);
        if(is_in_rectangle) {
            return dist(urng);
        } else if(_weights[i] < _weights[i+1]) {
            return (std::max)(dist(urng), dist(urng));
        } else {
            return (std::min)(dist(urng), dist(urng));
        }
    }
    
    /**
     * Returns a value distributed according to the parameters
     * specified by param.
     */
    template<class URNG>
    RealType operator()(URNG& urng, const param_type& parm) const
    {
        return piecewise_linear_distribution(parm)(urng);
    }
    
    /** Returns the smallest value that the distribution can produce. */
    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return _intervals.front(); }
    /** Returns the largest value that the distribution can produce. */
    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const
    { return _intervals.back(); }

    /**
     * Returns a vector containing the probability densities
     * at the interval boundaries.
     */
    std::vector<RealType> densities() const
    {
        RealType sum = static_cast<RealType>(0);
        for(std::size_t i = 0; i < _intervals.size() - 1; ++i) {
            RealType width = _intervals[i + 1] - _intervals[i];
            sum += (_weights[i] + _weights[i + 1]) * width / 2;
        }
        std::vector<RealType> result;
        result.reserve(_weights.size());
        for(typename std::vector<RealType>::const_iterator
            iter = _weights.begin(), end = _weights.end();
            iter != end; ++iter)
        {
            result.push_back(*iter / sum);
        }
        return result;
    }
    /**  Returns a vector containing the interval boundaries. */
    std::vector<RealType> intervals() const { return _intervals; }

    /** Returns the parameters of the distribution. */
    param_type param() const
    {
        return param_type(_intervals, _weights);
    }
    /** Sets the parameters of the distribution. */
    void param(const param_type& parm)
    {
        std::vector<RealType> new_intervals(parm._intervals);
        std::vector<RealType> new_weights(parm._weights);
        init(new_intervals, new_weights);
        _intervals.swap(new_intervals);
        _weights.swap(new_weights);
    }
    
    /**
     * Effects: Subsequent uses of the distribution do not depend
     * on values produced by any engine prior to invoking reset.
     */
    void reset() { _bins.reset(); }

    /** Writes a distribution to a @c std::ostream. */
    BOOST_RANDOM_DETAIL_OSTREAM_OPERATOR(
        os, piecewise_linear_distribution, pld)
    {
        os << pld.param();
        return os;
    }

    /** Reads a distribution from a @c std::istream */
    BOOST_RANDOM_DETAIL_ISTREAM_OPERATOR(
        is, piecewise_linear_distribution, pld)
    {
        param_type parm;
        if(is >> parm) {
            pld.param(parm);
        }
        return is;
    }

    /**
     * Returns true if the two distributions will return the
     * same sequence of values, when passed equal generators.
     */
    BOOST_RANDOM_DETAIL_EQUALITY_OPERATOR(
        piecewise_linear_distribution, lhs,  rhs)
    {
        return lhs._intervals == rhs._intervals && lhs._weights == rhs._weights;
    }
    /**
     * Returns true if the two distributions may return different
     * sequences of values, when passed equal generators.
     */
    BOOST_RANDOM_DETAIL_INEQUALITY_OPERATOR(piecewise_linear_distribution)

private:

    /// @cond \show_private

    void init(const std::vector<RealType>& intervals_arg,
              const std::vector<RealType>& weights_arg)
    {
        using std::abs;
        std::vector<RealType> bin_weights;
        bin_weights.reserve((intervals_arg.size() - 1) * 2);
        for(std::size_t i = 0; i < intervals_arg.size() - 1; ++i) {
            RealType width = intervals_arg[i + 1] - intervals_arg[i];
            RealType w1 = weights_arg[i];
            RealType w2 = weights_arg[i + 1];
            bin_weights.push_back((std::min)(w1, w2) * width);
            bin_weights.push_back(abs(w1 - w2) * width / 2);
        }
        typedef discrete_distribution<std::size_t, RealType> bins_type;
        typename bins_type::param_type bins_param(bin_weights);
        _bins.param(bins_param);
    }

    void init()
    {
        init(_intervals, _weights);
    }

    void default_init()
    {
        _intervals.clear();
        _intervals.push_back(RealType(0));
        _intervals.push_back(RealType(1));
        _weights.clear();
        _weights.push_back(RealType(1));
        _weights.push_back(RealType(1));
        init();
    }

    discrete_distribution<std::size_t, RealType> _bins;
    std::vector<RealType> _intervals;
    std::vector<RealType> _weights;

    /// @endcond
};

}
}

#endif