/*!
 * \file
 * \brief Class for numerically efficient log-likelihood algebra
 * \author Erik G. Larsson and Martin Senst
 *
 * -------------------------------------------------------------------------
 *
 * IT++ - C++ library of mathematical, signal processing, speech processing,
 *        and communications classes and functions
 *
 * Copyright (C) 1995-2008  (see AUTHORS file for a list of contributors)
 *
 * This program 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 2 of the License, or
 * (at your option) any later version.
 *
 * This program 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 this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 *
 * -------------------------------------------------------------------------
 */

#ifndef LLR_H
#define LLR_H

#include <limits>
#include <itpp/base/vec.h>
#include <itpp/base/mat.h>
#include <itpp/base/specmat.h>
#include <itpp/base/matfunc.h>
#include <limits>

namespace itpp {

  /*! \relates LLR_calc_unit
    The quantized log-likelihood ratio (QLLR) representation, scalar form. See \c LLR_calc_unit.
  */
  typedef signed int QLLR;

  /*!  \relates LLR_calc_unit
    The quantized log-likelihood ratio (QLLR) representation, vector form. See \c LLR_calc_unit.
  */
  typedef Vec<QLLR> QLLRvec;

  /*!  \relates LLR_calc_unit
    The quantized log-likelihood ratio (QLLR) representation, matrix form. See \c LLR_calc_unit.
  */
  typedef Mat<QLLR> QLLRmat;

  /*!  \relates LLR_calc_unit
    The largest possible QLLR value
  */
  const QLLR QLLR_MAX = (std::numeric_limits<QLLR>::max() >> 4);
  // added some margin to make sure the sum of two LLR is still permissible

  /*!
    \brief Log-likelihood algebra calculation unit.

    This class contains functions for algebra with log-likelihood
    ratios (LLRs). The (sole) purpose of this class is to provide
    numerically efficient functions for turbo and LDPC codes, which
    rely on certain nonlinear operations on LLRs.

    An LLR for an information bit b is defined according to \f[ L =
    \log \frac{P(b=0)}{P(b=1)} \f] and it is in general a real number.

    LLR values are represented via the special type, "quantized
    LLR" (QLLR).  The relation between the quantized representation
    and the real (floating-point) LLR value is
    \f[ \mbox{QLLR} = \mbox{round} \left(2^{\mbox{Dint1}}\cdot
    \mbox{LLR}\right) \f]
    The user parameter Dint1 determines the
    granularity of the quantization, and it can be set arbitrarily.
    The functions to_double() and to_qllr() can be used to perform
    conversions between the two representations (QLLR to
    floating-point, and vice versa).

    The class provides functions for the computation of the Jacobian
    logarithm and Hagenauer's "boxplus" operator.  These functions are
    based on a table-lookup.  The resolution of the table is
    determined by the parameters Dint2 and Dint3. See the class
    constructor for more detail.  When an object of LLR_calc_unit is
    created, corresponding lookup-tables are also generated. The
    resolution of these tables can be adjusted by providing parameters
    to the constructor.

    The variable table resolution allows one to study complexity
    versus accuracy (i.e., how different table resolutions would
    degrade performance) to some extent. Yet the main purpose of the
    QLLR representation is to provide a tool for writing efficient
    simulation code, rather than to provide for bit-level
    (fixed-point) simulations. For bit-level simulations, a true fixed
    point representation of LLRs would be preferred/required.  With
    the default setting of the table parameters, using the QLLR type
    is practically as accurate (but much faster) as using "double" to
    represent LLRs.  Decoder implementations may then provide
    functions using QLLR, fixed-point, or double (for compatibility
    reasons) representations of LLR values.

    Note: the QLLR type does not check that the correct quantization
    level is used. I.e., in theory it would be possible to add two
    QLLR types with different quantization (Dint) parameters.  This is
    intentionally implemented this way to achieve maximum runtime
    efficiency.

  */
  class LLR_calc_unit {
  public:
    //! Constructor, using the default table resolution
    LLR_calc_unit();

    /*!
     * \brief Constructor, using a specific table resolution.
     *
     * See init_llr_tables() for more details on the parameters.
     */
    LLR_calc_unit(short int Dint1, short int Dint2, short int Dint3);

    /*! \brief Set the quantization and table parameters

      \param Dint1 Determines the relation between LLR represented as
      real number and as integer.  The relation is
      \f[ \mbox{QLLR} = \mbox{round} \left(2^{\mbox{Dint1}}\cdot
      \mbox{LLR}\right) \f]

      \param Dint2 Number of entries in the table. If this is zero,
      then logmap becomes logmax.

      \param Dint3 Determines the table resolution. The spacing between each
      entry is \f[ 2^{-(Dint1-Dint3)} \f]

      The default parameter values are chosen to give a performance
      practically indistinguishable from that of using floating point
      calculations.

      Example: (recommended settings with "exact" computation via high
      resolution lookup table)
      \code
      LLR_calc_unit lcalc(12, 300, 7);
      \endcode

      Example: (recommended settings with logmax, i.e. no table lookup)
      \code
      LLR_calc_unit lcalc(12, 0, 7);
      \endcode
    */
    void init_llr_tables(short int Dint1 = 12, short int Dint2 = 300,
                         short int Dint3 = 7);

    //! Convert a "real" LLR value to an LLR type
    QLLR to_qllr(double l) const;

    //! Convert a vector of "real" LLR values to an LLR type
    QLLRvec to_qllr(const vec &l) const;

    //! Convert a matrix of "real" LLR values to an LLR type
    QLLRmat to_qllr(const mat &l) const;

    //! Convert an LLR type to a "real" LLR
    double to_double(QLLR l) const;

    //! Convert a vector of LLR types to a "real" LLR
    vec to_double(const QLLRvec &l) const;

    //! Convert a matrix of LLR types to a "real" LLR
    mat to_double(const QLLRmat &l) const;

    /*!
     * \brief Jacobian logarithm.
     *
     * This function computes \f[ \log(\exp(a)+\exp(b)) \f]
     */
    inline QLLR jaclog(QLLR a, QLLR b) const;
    // Note: a version of this function taking "double" values as input
    // is deliberately omitted, because this is rather slow.

    /*!
     * \brief Hagenauer's "Boxplus" operator.
     *
     * This function computes:
     * \f[ \mbox{sign}(a) * \mbox{sign}(b) * \mbox{min}(|a|,|b|)
     * + f(|a+b|) - f(|a-b|) \f]
     * where \f[ f(x) = \log(1+\exp(-x))  \f]
     */
    QLLR Boxplus(QLLR a, QLLR b) const;

    /*!
     * \brief Logexp operator.
     *
     * This function computes \f[ f(x) = \log(1+\exp(-x)) \f]
     */
    inline QLLR logexp(QLLR x) const;

    //! Retrieve the table resolution values
    ivec get_Dint();

    //! Print some properties of the LLR calculation unit in plain text
    friend std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &l);

  private:
    //! Compute the table for \f[ f(x) = \log(1+\exp(-x)) \f]
    ivec construct_logexp_table();

    //! The lookup tables for the decoder
    ivec logexp_table;

    //! Decoder (lookup-table) parameters
    short int Dint1, Dint2, Dint3;
  };

  /*!
    \relatesalso LLR_calc_unit
    \brief Print some properties of the LLR calculation unit in plain text.
  */
  std::ostream &operator<<(std::ostream &os, const LLR_calc_unit &lcu);


  // ----------------------------------------------------------------------
  // implementation of some inline functions
  // ----------------------------------------------------------------------

  inline double LLR_calc_unit::to_double(QLLR l) const
  {
    return static_cast<double>(l) / (1<<Dint1);
  }

  inline QLLR LLR_calc_unit::to_qllr(double l) const
  {
    double QLLR_MAX_double = to_double(QLLR_MAX);
    // Don't abort when overflow occurs, just saturate the QLLR
    if (l > QLLR_MAX_double) {
      it_info_debug("LLR_calc_unit::to_qllr(): LLR overflow");
      return QLLR_MAX;
    }
    if (l < -QLLR_MAX_double) {
      it_info_debug("LLR_calc_unit::to_qllr(): LLR overflow");
      return -QLLR_MAX;
    }
    return static_cast<QLLR>(std::floor(0.5 + (1<<Dint1) * l));
  }


  inline QLLR LLR_calc_unit::logexp(QLLR x) const
  {
    it_assert_debug(x >= 0,"LLR_calc_unit::logexp(): Wrong LLR value");
    int ind = x >> Dint3;
    if (ind >= Dint2) // outside table
      return 0;

    it_assert_debug(ind >= 0,"LLR_calc_unit::logexp(): Internal error");
    it_assert_debug(ind < Dint2,"LLR_calc_unit::logexp(): internal error");

    // With interpolation
    // int delta=x-(ind<<Dint3);
    // return ((delta*logexp_table(ind+1) + ((1<<Dint3)-delta)*logexp_table(ind)) >> Dint3);

    // Without interpolation
    return logexp_table(ind);
  }


  inline QLLR LLR_calc_unit::jaclog(QLLR a, QLLR b) const
  {
    QLLR x, maxab;

    if (a > b) {
      maxab = a;
      x = a - b;
    }
    else {
      maxab = b;
      x = b - a;
    }

    if (maxab >= QLLR_MAX)
      return QLLR_MAX;
    else
      return (maxab + logexp(x));
  }

}

#endif