text:
#pragma once
/*
 * (C) Copyright 2020-2021 UCAR
 *
 * This software is licensed under the terms of the Apache Licence Version 2.0
 * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
 */
/*! \defgroup ioda_cxx_variable Variables, Data Access, and Selections
 * \brief The main data storage methods and objects in IODA.
 * \ingroup ioda_cxx_api
 *
 * @{
 * \file Variable.h
 * \brief @link ioda_cxx_variable Interfaces @endlink for ioda::Variable and related classes.
 */
#include <cstring>
#include <gsl/gsl-lite.hpp>
#include <list>
#include <map>
#include <memory>
#include <string>
#include <tuple>
#include <utility>
#include <vector>

#include "ioda/Attributes/Has_Attributes.h"
#include "ioda/Exception.h"
#include "ioda/Misc/Eigen_Compat.h"
#include "ioda/Python/Var_ext.h"
#include "ioda/Types/Marshalling.h"
#include "ioda/Types/Type.h"
#include "ioda/Types/Type_Provider.h"
#include "ioda/Variables/Fill.h"
#include "ioda/Variables/Selection.h"
#include "ioda/defs.h"

namespace ioda {
class Attribute;
class Variable;
struct VariableCreationParameters;

struct Named_Variable;

namespace detail {
class Attribute_Backend;
class Variable_Backend;

/// \brief Exists to prevent constructor conflicts when passing a backend into
///   a frontend object.
/// \ingroup ioda_cxx_variable
/// \bug Need to add a virtual resize function.
template <class Variable_Implementation = Variable>
class Variable_Base {
protected:
  /// Using an opaque object to implement the backend.
  std::shared_ptr<Variable_Backend> backend_;

  /// @name General Functions
  /// @{

  Variable_Base(std::shared_ptr<Variable_Backend>);

public:
  virtual ~Variable_Base();

  /// @}
  /// @name Metadata manipulation
  /// @{

  /// Attributes
  Has_Attributes atts;

  /// @}

  /// @name General Functions
  /// @{

  /// Gets a handle to the underlying object that implements the backend functionality.
  std::shared_ptr<Variable_Backend> get() const;

  /// @}
  /// @name Type-querying Functions
  /// @{

  /// Get type
  virtual Type getType() const;
  /// Get type
  inline Type type() const { return getType(); }

  /// Query the backend and get the type provider.
  virtual detail::Type_Provider* getTypeProvider() const;

  /// \brief Convenience function to check a Variable's storage type.
  /// \param DataType is the type of the data. I.e. float, int, int32_t, uint16_t, std::string, etc.
  /// \returns True if the type matches
  /// \returns False (0) if the type does not match
  /// \throws ioda::Exception if an error occurred.
  template <class DataType>
  bool isA() const {
    Type templateType = Types::GetType_Wrapper<DataType>::GetType(getTypeProvider());

    return isA(templateType);
  }
  /// Hand-off to the backend to check equivalence
  virtual bool isA(Type lhs) const;

  /// Python compatability function
  inline bool isA(BasicTypes dataType) const { return isA(Type(dataType, getTypeProvider())); }
  /// \internal pybind11
  inline bool _py_isA2(BasicTypes dataType) { return isA(dataType); }
  /// Convenience function to query type
  BasicTypes getBasicType() const;

  /// @}
  /// @name Querying Functions for Fill Values, Chunking and Compression
  /// @{

  /// @brief Convenience function to get fill value, attributes,
  ///   chunk sizes, and compression in a collective call.
  /// @details This function has better performance on some engines
  ///   for bulk operations than calling separately.
  /// @param doAtts includes attributes in the creation parameters.
  /// @param doDims includes dimension scales in the creation parameters.
  /// Although dimensions are attributes on some backends, we treat them
  /// separately in this function.
  /// @returns the filled-in VariableCreationParameters object
  virtual VariableCreationParameters getCreationParameters(bool doAtts = true,
                                                           bool doDims = true) const;

  /// \brief Check if a variable has a fill value set
  /// \returns true if a fill value is set, false otherwise.
  virtual bool hasFillValue() const;

  /// Remap fill value storage type into this class.
  typedef detail::FillValueData_t FillValueData_t;

  /// \brief Retrieve the fill value
  /// \returns an object of type FillValueData_t that stores
  ///   the fill value. If there is no fill value, then
  ///   FillValueData_t::set_ will equal false. The fill
  ///   value data will be stored in
  ///   FillValueData_t::fillValue_ (for simple types), or in
  ///   FillValueData_t::stringFillValue_ (for strings only).
  /// \note Recommend querying isA to make sure that you
  ///   are reading the fill value as the correct type.
  virtual FillValueData_t getFillValue() const;

  /// \brief Retrieve the chunking options for the Variable.
  /// \note Not all backends support chunking, but they should
  ///   all store the desired chunk size information in case the
  ///   Variable is copied to a new backend.
  /// \returns a vector containing the chunk sizes.
  /// \returns an empty vector if chunking is not used.
  virtual std::vector<Dimensions_t> getChunkSizes() const;

  /// \brief Retrieve the GZIP compression options for the Variable.
  /// \note Not all backends support compression (and those that
  ///   do also require chunking support). They should store this
  ///   information anyways, in case the Variable is copied to
  ///   a new backend.
  /// \returns a pair indicating 1) whether GZIP is requested and
  ///   2) the compression level that is desired.
  virtual std::pair<bool, int> getGZIPCompression() const;

  /// \brief Retrieve the SZIP compression options for the Variable.
  /// \note Not all backends support compression (and those that
  ///   do also require chunking support). They should store this
  ///   information anyways, in case the Variable is copied to
  ///   a new backend.
  /// \returns a tuple indicating 1) whether SZIP is requested,
  ///   2) PixelsPerBlock, and 3) general SZIP filter option flags.
  virtual std::tuple<bool, unsigned, unsigned> getSZIPCompression() const;

  /// @}
  /// @name Data Space-Querying Functions
  /// @{

  // Get dataspace
  // JEDI_NODISCARD DataSpace getSpace() const;

  /// Get current and maximum dimensions, and number of total points.
  /// \note In Python, see the dims property.
  virtual Dimensions getDimensions() const;

  /// \brief Resize the variable.
  /// \note Not all variables are resizable. This depends
  ///   on backend support. For HDF5, the variable must be chunked
  ///   and must not exceed max_dims.
  /// \note Bad things may happen if a variable's dimension scales
  ///   have different lengths than its dimensions. Resize them
  ///   together, preferably using the ObsSpace resize function.
  /// \see ObsSpace::resize
  /// \param newDims are the new dimensions.
  virtual Variable resize(const std::vector<Dimensions_t>& newDims);

  /// Attach a dimension scale to this Variable.
  virtual Variable attachDimensionScale(unsigned int DimensionNumber, const Variable& scale);
  /// Detach a dimension scale
  virtual Variable detachDimensionScale(unsigned int DimensionNumber, const Variable& scale);
  /// Set dimensions (convenience function to several invocations of attachDimensionScale).
  Variable setDimScale(const std::vector<Variable>& dims);
  /// Set dimensions (convenience function to several invocations of attachDimensionScale).
  Variable setDimScale(const std::vector<Named_Variable>& dims);
  /// Set dimensions (convenience function to several invocations of attachDimensionScale).
  Variable setDimScale(const Variable& dims);
  /// Set dimensions (convenience function to several invocations of attachDimensionScale).
  Variable setDimScale(const Variable& dim1, const Variable& dim2);
  /// Set dimensions (convenience function to several invocations of attachDimensionScale).
  Variable setDimScale(const Variable& dim1, const Variable& dim2, const Variable& dim3);

  /// Is this Variable used as a dimension scale?
  virtual bool isDimensionScale() const;

  /// Designate this table as a dimension scale
  virtual Variable setIsDimensionScale(const std::string& dimensionScaleName);
  /// Get the name of this Variable's defined dimension scale
  inline std::string getDimensionScaleName() const {
    std::string r;
    getDimensionScaleName(r);
    return r;
  }
  virtual Variable getDimensionScaleName(std::string& res) const;

  /// Is a dimension scale attached to this Variable in a certain position?
  virtual bool isDimensionScaleAttached(unsigned int DimensionNumber, const Variable& scale) const;

  /// \brief Which dimensions are attached at which positions? This function may offer improved
  /// performance on some backends compared to serial isDimensionScaleAttached calls.
  /// \param scalesToQueryAgainst is a vector containing the scales. You can pass in
  ///   "tagged" strings that map the Variable to a name.
  ///   If you do not pass in a scale to check against, then this scale will not be checked and
  ///   will not be present in the function output.
  /// \param firstOnly is specified when only one dimension can be attached to each axis (the
  /// default).
  /// \returns a vector with the same length as the variable's dimensionality.
  ///   Each variable dimension in the vector can have one or more attached scales. These scales are
  ///   returned as their own, inner, vector of pair<string, Variable>.
  virtual std::vector<std::vector<Named_Variable>> getDimensionScaleMappings(
    const std::list<Named_Variable>& scalesToQueryAgainst,
    bool firstOnly = true) const;

  /// @}
  /// @name Writing Data
  /// @{

  /// \brief The fundamental write function. Backends overload this function to implement
  ///   all write operations.
  ///
  /// \details This function writes a span of bytes (characters) to the backend attribute
  ///   storage. No type conversions take place here (see the templated conversion function,
  ///   below).
  ///
  /// \param data is a span of data.
  /// \param in_memory_datatype is an opaque (backend-level) object that describes the
  ///   placement of the data in memory. Usually ignorable - needed for complex data structures.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \throws ioda::xError if data has the wrong size.
  /// \returns The variable (for chaining).
  virtual Variable write(gsl::span<const char> data, const Type& in_memory_dataType,
                         const Selection& mem_selection  = Selection::all,
                         const Selection& file_selection = Selection::all);

  /// \brief Write the Variable
  /// \note Ensure that the correct dimension ordering is preserved.
  /// \note With default parameters, the entire Variable is written.
  /// \tparam DataType is the type of the data to be written.
  /// \tparam Marshaller is a class that serializes / deserializes data.
  /// \tparam TypeWrapper translates DataType into a form that the backend understands.
  /// \param data is a span of data.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \throws ioda::xError if data has the wrong size.
  /// \returns The variable (for chaining).
  template <class DataType, class Marshaller = Object_Accessor<DataType>,
            class TypeWrapper = Types::GetType_Wrapper<DataType>>
  Variable_Implementation write(const gsl::span<DataType> data,
                                const Selection& mem_selection  = Selection::all,
                                const Selection& file_selection = Selection::all) {
    try {
      Marshaller m;
      auto d = m.serialize(data);
      return write(gsl::make_span<const char>(
                      reinterpret_cast<const char*>(d->DataPointers.data()),
                     d->DataPointers.size() * Marshaller::bytesPerElement_),
                   TypeWrapper::GetType(getTypeProvider()), mem_selection, file_selection);
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  /// \brief Write the Variable
  /// \note Ensure that the correct dimension ordering is preserved.
  /// \note With default parameters, the entire Variable is written.
  /// \tparam DataType is the type of the data to be written.
  /// \tparam Marshaller is a class that serializes / deserializes data.
  /// \tparam TypeWrapper translates DataType into a form that the backend understands.
  /// \param data is a span of data.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \throws ioda::xError if data has the wrong size.
  /// \returns The variable (for chaining).
  template <class DataType, class Marshaller = Object_Accessor<DataType>,
            class TypeWrapper = Types::GetType_Wrapper<DataType>>
  Variable_Implementation write(const gsl::span<const DataType> data,
                                const Selection& mem_selection  = Selection::all,
                                const Selection& file_selection = Selection::all) {
    try {
      Marshaller m;
      auto d = m.serialize(data);
      return write(gsl::make_span<const char>(
                      reinterpret_cast<const char*>(d->DataPointers.data()),
                     d->DataPointers.size() * Marshaller::bytesPerElement_),
                   TypeWrapper::GetType(getTypeProvider()), mem_selection, file_selection);
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  /// \brief Write the variable
  /// \tparam DataType is the type of the data to be written.
  /// \tparam Marshaller is a class that serializes / deserializes data.
  /// \tparam TypeWrapper translates DataType into a form that the backend understands.
  /// \param data is a span of data.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \throws ioda::xError if data has the wrong size.
  /// \returns The variable (for chaining).

  template <class DataType, class Marshaller = Object_Accessor<DataType>,
            class TypeWrapper = Types::GetType_Wrapper<DataType>>
  Variable_Implementation write(const std::vector<DataType>& data,
                                const Selection& mem_selection  = Selection::all,
                                const Selection& file_selection = Selection::all) {
    try {
      return this->write<DataType, Marshaller, TypeWrapper>(gsl::make_span(data), mem_selection,
                                                            file_selection);
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  /// \brief Write an Eigen object (a Matrix, an Array, a Block, a Map).
  /// \tparam EigenClass is the type of the Eigen object being written.
  /// \param d is the data to be written.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \throws ioda::xError on a dimension mismatch.
  /// \returns the variable
  template <class EigenClass>
  Variable_Implementation writeWithEigenRegular(const EigenClass& d,
                                                const Selection& mem_selection  = Selection::all,
                                                const Selection& file_selection = Selection::all) {
#if 1  //__has_include("Eigen/Dense")
    try {
      typedef typename EigenClass::Scalar ScalarType;
      // If d is already in Row Major form, then this is optimized out.
      Eigen::Array<ScalarType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> dout;
      dout.resize(d.rows(), d.cols());
      dout               = d;
      const auto& dconst = dout;  // To make some compilers happy.
      auto sp            = gsl::make_span(dconst.data(), static_cast<int>(d.rows() * d.cols()));

      return write<ScalarType>(sp, mem_selection, file_selection);
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
#else
    static_assert(false, "The Eigen headers cannot be found, so this function cannot be used.");
#endif
  }

  /// \brief Write an Eigen Tensor-like object
  /// \tparam EigenClass is the type of the Eigen object being written.
  /// \param d is the data to be written.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \throws ioda::xError on a dimension mismatch.
  /// \returns the variable
  template <class EigenClass>
  Variable_Implementation writeWithEigenTensor(const EigenClass& d,
                                               const Selection& mem_selection  = Selection::all,
                                               const Selection& file_selection = Selection::all) {
#if 1  //__has_include("unsupported/Eigen/CXX11/Tensor")
    try {
      ioda::Dimensions dims = detail::EigenCompat::getTensorDimensions(d);

      auto sp  = (gsl::make_span(d.data(), dims.numElements));
      auto res = write(sp, mem_selection, file_selection);
      return res;
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
#else
    static_assert(
      false, "The Eigen unsupported/ headers cannot be found, so this function cannot be used.");
#endif
  }

  /// @}
  /// @name Reading Data
  /// @{

  /// \brief Read the Variable - as char array. Ordering is row-major.
  /// \details This is the fundamental read function that has to be implemented.
  /// \param data is a byte-array that will hold the read data.
  /// \param in_memory_dataType describes how ioda should arrange the read data in memory.
  ///   As floats? As doubles? Strings?
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \note Ensure that the correct dimension ordering is preserved
  /// \note With default parameters, the entire Variable is read
  virtual Variable read(gsl::span<char> data, const Type& in_memory_dataType,
                        const Selection& mem_selection  = Selection::all,
                        const Selection& file_selection = Selection::all) const;

  /// \brief Read the variable into a span (range) or memory. Ordering is row-major.
  /// \tparam DataType is the type of the data to be written.
  /// \tparam Marshaller is a class that serializes / deserializes data.
  /// \tparam TypeWrapper translates DataType into a form that the backend understands.
  /// \param data is a byte-array that will hold the read data.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \todo Add in the dataspaces!
  template <class DataType, class Marshaller = ioda::Object_Accessor<DataType>,
            class TypeWrapper = Types::GetType_Wrapper<DataType>>
  Variable_Implementation read(gsl::span<DataType> data,
                               const Selection& mem_selection  = Selection::all,
                               const Selection& file_selection = Selection::all) const {
    try {
      const size_t numObjects = data.size();

      detail::PointerOwner pointerOwner = getTypeProvider()->getReturnedPointerOwner();
      Marshaller m(pointerOwner);
      auto p = m.prep_deserialize(numObjects);
      read(gsl::make_span<char>(
             reinterpret_cast<char*>(p->DataPointers.data()),
             // Logic note: sizeof mutable data type. If we are
             // reading in a string, then mutable data type is char*,
             // which works because address pointers have the same size.
             p->DataPointers.size() * Marshaller::bytesPerElement_),
           TypeWrapper::GetType(getTypeProvider()), mem_selection, file_selection);
      m.deserialize(p, data);

      return Variable_Implementation{backend_};
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
  }

  /// \brief Read the variable into a vector. Resize if needed. For a non-resizing
  ///   version, use a gsl::span.
  /// \details Ordering is row-major.
  /// \tparam DataType is the type of the data to be written.
  /// \tparam Marshaller is a class that serializes / deserializes data.
  /// \tparam TypeWrapper translates DataType into a form that the backend understands.
  /// \param data is a byte-array that will hold the read data.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \bug Resize only if needed, and resize to the proper extent depending on
  ///   mem_selection and file_selection.
  template <class DataType, class Marshaller = ioda::Object_Accessor<DataType>,
            class TypeWrapper = Types::GetType_Wrapper<DataType>>
  Variable_Implementation read(std::vector<DataType>& data,
                               const Selection& mem_selection  = Selection::all,
                               const Selection& file_selection = Selection::all) const {
    data.resize(getDimensions().numElements); // TODO(Ryan): remove
    return read<DataType, Marshaller, TypeWrapper>(gsl::make_span(data.data(), data.size()),
                                                   mem_selection, file_selection);
  }

  /// \brief Read the variable into a new vector. Python convenience function.
  /// \bug Get correct size based on selection operands.
  /// \tparam DataType is the type of the data to be written.
  /// \tparam Marshaller is a class that serializes / deserializes data.
  /// \tparam TypeWrapper translates DataType into a form that the backend understands.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  template <class DataType, class Marshaller = ioda::Object_Accessor<DataType>,
            class TypeWrapper = Types::GetType_Wrapper<DataType>>
  std::vector<DataType> readAsVector(const Selection& mem_selection  = Selection::all,
                                     const Selection& file_selection = Selection::all) const {
    std::vector<DataType> data(getDimensions().numElements);
    read<DataType, Marshaller, TypeWrapper>(gsl::make_span(data.data(), data.size()), mem_selection,
                                            file_selection);
    return data;
  }

  /// \brief Valarray read convenience function. Resize if needed.
  ///   For a non-resizing version, use a gsl::span.
  /// \tparam DataType is the type of the data to be written.
  /// \tparam Marshaller is a class that serializes / deserializes data.
  /// \tparam TypeWrapper translates DataType into a form that the backend understands.
  /// \param data is a valarray acting as a data buffer that is filled with the
  ///   metadata's contents. It gets resized as needed.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \returns Another instance of this Attribute. Used for operation chaining.
  /// \note data will be stored in row-major order.
  template <class DataType, class Marshaller = ioda::Object_Accessor<DataType>,
            class TypeWrapper = Types::GetType_Wrapper<DataType>>
  Variable_Implementation read(std::valarray<DataType>& data,
                               const Selection& mem_selection  = Selection::all,
                               const Selection& file_selection = Selection::all) const {
    /// \bug Resize only if needed, and resize to the proper extent depending on
    /// mem_selection and file_selection.
    data.resize(getDimensions().numElements); // TODO(Ryan): remove
    return read<DataType, Marshaller, TypeWrapper>(gsl::make_span(std::begin(data), std::end(data)),
                                                   mem_selection, file_selection);
  }

  /// \brief Read data into an Eigen::Array, Eigen::Matrix, Eigen::Map, etc.
  /// \tparam EigenClass is a template pointing to the Eigen object.
  ///   This template must provide the EigenClass::Scalar typedef.
  /// \tparam Resize indicates whether the Eigen object should be resized
  ///   if there is a dimension mismatch. Not all Eigen objects can be resized.
  /// \param res is the Eigen object.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \returns Another instance of this Variable. Used for operation chaining.
  /// \throws ioda::xError if the variable's dimensionality is
  ///   too high.
  /// \throws ioda::xError if resize = false and there is a dimension mismatch.
  /// \note When reading in a 1-D object, the data are read as a column vector.
  template <class EigenClass, bool Resize = detail::EigenCompat::CanResize<EigenClass>::value>
  Variable_Implementation readWithEigenRegular(EigenClass& res,
                                               const Selection& mem_selection = Selection::all,
                                               const Selection& file_selection
                                               = Selection::all) const {
    /// \bug Resize only if needed, and resize to the proper extent depending on
    /// mem_selection and file_selection.
#if 1  //__has_include("Eigen/Dense")
    try {
      typedef typename EigenClass::Scalar ScalarType;

      static_assert(
        !(Resize && !detail::EigenCompat::CanResize<EigenClass>::value),
        "This object cannot be resized, but you have specified that a resize is required.");

      // Check that the dimensionality is 1 or 2.
      const auto dims = getDimensions();
      if (dims.dimensionality > 2)
        throw Exception("Dimensionality too high for a regular Eigen read. Use "
          "Eigen::Tensor reads instead.", ioda_Here());

      int nDims[2] = {1, 1};
      if (dims.dimsCur.size() >= 1) nDims[0] = gsl::narrow<int>(dims.dimsCur[0]);
      if (dims.dimsCur.size() >= 2) nDims[1] = gsl::narrow<int>(dims.dimsCur[1]);

      // Resize if needed.
      if (Resize)
        detail::EigenCompat::DoEigenResize(res, nDims[0],
                                           nDims[1]);  // nullop if the size is already correct.
      else if (dims.numElements != (size_t)(res.rows() * res.cols()))
        throw Exception("Size mismatch", ioda_Here());

      // Array copy to preserve row vs column major format.
      // Should be optimized away by the compiler if unneeded.
      // Note to the reader: We are reading in the data to a temporary object.
      // We can size _this_ temporary object however we want.
      // The temporary is used to swap row / column indices if needed.
      // It should be optimized away if not needed... making sure this happens is a todo.
      Eigen::Array<ScalarType, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> data_in(res.rows(),
                                                                                res.cols());

      auto ret = read<ScalarType>(gsl::span<ScalarType>(data_in.data(), dims.numElements),
                                  mem_selection, file_selection);
      res      = data_in;
      return ret;
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
#else
    static_assert(false, "The Eigen headers cannot be found, so this function cannot be used.");
#endif
  }

  /// \brief Read data into an Eigen::Array, Eigen::Matrix, Eigen::Map, etc.
  /// \tparam EigenClass is a template pointing to the Eigen object.
  ///   This template must provide the EigenClass::Scalar typedef.
  /// \param res is the Eigen object.
  /// \param mem_selection is the user's memory layout representing the location where
  ///   the data is read from.
  /// \param file_selection is the backend's memory layout representing the
  ///   location where the data are written to.
  /// \returns Another instance of this Variable. Used for operation chaining.
  /// \throws ioda::Exception if there is a size mismatch.
  /// \note When reading in a 1-D object, the data are read as a column vector.
  template <class EigenClass>
  Variable_Implementation readWithEigenTensor(EigenClass& res,
                                              const Selection& mem_selection = Selection::all,
                                              const Selection& file_selection
                                              = Selection::all) const {
#if 1  //__has_include("unsupported/Eigen/CXX11/Tensor")
    try {
      // Check dimensionality of source and destination
      const auto ioda_dims  = getDimensions();
      const auto eigen_dims = ioda::detail::EigenCompat::getTensorDimensions(res);
      if (ioda_dims.numElements != eigen_dims.numElements)
        throw Exception("Size mismatch for Eigen Tensor-like read.", ioda_Here());

      auto sp = (gsl::make_span(res.data(), eigen_dims.numElements));
      return read(sp, mem_selection, file_selection);
    } catch (...) {
      std::throw_with_nested(Exception(ioda_Here()));
    }
#else
    static_assert(
      false, "The Eigen unsupported/ headers cannot be found, so this function cannot be used.");
#endif
  }

  /// \internal Python binding function
  template <class EigenClass>
  EigenClass _readWithEigenRegular_python(const Selection& mem_selection  = Selection::all,
                                          const Selection& file_selection = Selection::all) const {
    EigenClass data;
    readWithEigenRegular(data, mem_selection, file_selection);
    return data;
  }

  /// @brief Convert a selection into its backend representation
  /// @param sel is the frontend selection.
  /// @return The cached backend selection.
  virtual Selections::SelectionBackend_t instantiateSelection(const Selection& sel) const;

  /// @}

private:
  /// \brief get the fill value from the netcdf specification (_FillValue attribute)
  FillValueData_t getNcFillValue() const;

  /// \brief check if fill data objects match, print warning if they don't match
  /// \param hdfFill fill value obtained from the hdf fill value property
  /// \param ncFill fill value obtained from the netcdf fill value attribute
  void checkWarnFillValue(FillValueData_t & hdfFill, FillValueData_t & ncFill) const;

  /// \brief run an action according to the current variable data type
  /// \details this function will call the function passed to it through the action
  /// parameter, giving this action a typed entity that can be identified by decltype.
  /// At this point, all of the types supported by Variable_Base::getBasicType() except
  /// for bool are handled in this function.
  /// \param action Function object callable with a single argument which is identifiable
  /// by decltype
  // TODO(srh) Additional development is required to properly handle the bool data type.
  template <typename Action>
  auto runForVarType(const Action & action) const {
    if (isA<int>()) {
      int typeMe;
      return action(typeMe);
    } else if (isA<unsigned int>()) {
      unsigned int typeMe;
      return action(typeMe);
    } else if (isA<float>()) {
      float typeMe;
      return action(typeMe);
    } else if (isA<double>()) {
      double typeMe;
      return action(typeMe);
    } else if (isA<std::string>()) {
      std::string typeMe;
      return action(typeMe);
    } else if (isA<long>()) {
      long typeMe;
      return action(typeMe);
    } else if (isA<unsigned long>()) {
      unsigned long typeMe;
      return action(typeMe);
    } else if (isA<short>()) {
      short typeMe;
      return action(typeMe);
    } else if (isA<unsigned short>()) {
      unsigned short typeMe;
      return action(typeMe);
    } else if (isA<long long>()) {
      long long typeMe;
      return action(typeMe);
    } else if (isA<unsigned long long>()) {
      unsigned long long typeMe;
      return action(typeMe);
    } else if (isA<int32_t>()) {
      int32_t typeMe;
      return action(typeMe);
    } else if (isA<uint32_t>()) {
      uint32_t typeMe;
      return action(typeMe);
    } else if (isA<int16_t>()) {
      int16_t typeMe;
      return action(typeMe);
    } else if (isA<uint16_t>()) {
      uint16_t typeMe;
      return action(typeMe);
    } else if (isA<int64_t>()) {
      int64_t typeMe;
      return action(typeMe);
    } else if (isA<uint64_t>()) {
      uint64_t typeMe;
      return action(typeMe);
    } else if (isA<long double>()) {
      long double typeMe;
      return action(typeMe);
    } else if (isA<char>()) {
      char typeMe;
      return action(typeMe);
    } else if (isA<unsigned char>()) {
      unsigned char typeMe;
      return action(typeMe);
    } else {
      std::throw_with_nested(Exception("Unsupported variable data type", ioda_Here()));
    }
  }

};
// extern template class Variable_Base<Variable>;
}  // namespace detail

/** \brief Variables store data!
 * \ingroup ioda_cxx_variable
 *
 * A variable represents a single field of data. It can be multi-dimensional and
 * usually has one or more attached **dimension scales**.
 *
 * Variables have Metadata, which describe the variable (i.e. valid_range, long_name, units).
 * Variables can have different data types (i.e. int16_t, float, double, string, datetime).
 * Variables can be resized. Depending on the backend, the data in a variable can be stored
 * using chunks, and may also be compressed.
 *
 * The backend manages how variables are stored in memory or on disk. The functions in the
 * Variable class provide methods to query and set data. The goal is to have data transfers
 * involve as few copies as possible.
 *
 * Variable objects themselves are lightweight handles that may be passed across different
 * parts of a program. Variables are always stored somewhere in a Group (or ObsSpace), so you
 * can always re-open a handle.
 *
 * \note Thread and MPI safety depend on the specific backends used to implement a variable.
 * \note A variable may be linked to multiple groups and listed under multiple names, so long as
 * the storage backends are all the same.
 **/
class IODA_DL Variable : public detail::Variable_Base<Variable> {
public:
  /// @name General Functions
  /// @{

  Variable();
  Variable(std::shared_ptr<detail::Variable_Backend> b);
  Variable(const Variable&);
  Variable& operator=(const Variable&);
  virtual ~Variable();

  /// @}
  /// @name Python compatability objects
  /// @{

  detail::python_bindings::VariableIsA<Variable> _py_isA;

  detail::python_bindings::VariableReadVector<Variable> _py_readVector;
  detail::python_bindings::VariableReadNPArray<Variable> _py_readNPArray;

  detail::python_bindings::VariableWriteVector<Variable> _py_writeVector;
  detail::python_bindings::VariableWriteNPArray<Variable> _py_writeNPArray;

  detail::python_bindings::VariableScales<Variable> _py_scales;

  /// @}
};

namespace detail {
/// \brief Variable backends inherit from this.
class IODA_DL Variable_Backend : public Variable_Base<Variable> {
public:
  virtual ~Variable_Backend();

  /// Default, trivial implementation. Customizable by backends for performance.
  std::vector<std::vector<Named_Variable>> getDimensionScaleMappings(
    const std::list<Named_Variable>& scalesToQueryAgainst,
    bool firstOnly = true) const override;

  /// Default implementation. Customizable by backends for performance.
  VariableCreationParameters getCreationParameters(bool doAtts = true,
                                                   bool doDims = true) const override;

protected:
  Variable_Backend();

  /// @brief This function de-encapsulates an Attribute's backend storage object.
  ///   This function is used by Variable_Backend's derivatives when accessing a
  ///   Variable's Attributes. IODA-internal use only.
  /// @tparam Attribute_Implementation is a dummy parameter used by Attributes.
  /// @param base is the Attribute whose backend you want.
  /// @return The encapsulated backend object.
  template <class Attribute_Implementation = Attribute>
  static std::shared_ptr<Attribute_Backend> _getAttributeBackend(
    const Attribute_Implementation& att) {
    return att.backend_;
  }

  /// @brief This function de-encapsulates a Has_Attributes backend storage object.
  ///   IODA-internal use only.
  /// @tparam Has_Attributes_Implementation is a dummy parameter for template resolution.
  /// @param hatts is the Has_Attributes whose backend you want.
  /// @return The encapsulated backend object.
  template <class Has_Attributes_Implementation = Has_Attributes>
  static std::shared_ptr<Has_Attributes_Backend> _getHasAttributesBackend(
    const Has_Attributes_Implementation& hatts) {
    return hatts.backend_;
  }
};
}  // namespace detail

/// @brief A named pair of (variable_name, ioda::Variable).
struct Named_Variable {
  std::string name;
  ioda::Variable var;
  bool operator<(const Named_Variable& rhs) const { return name < rhs.name; }

  Named_Variable() = default;
  Named_Variable(const std::string& name, const ioda::Variable& var) : name(name), var(var) {}
};

}  // namespace ioda

/// @}

max_stars_repo_name: NOAA-EMC/ioda

max_issues_repo_name: NOAA-EMC/ioda

max_forks_repo_name: NOAA-EMC/ioda


text:
#include "emulator_struct.h"
#include "libEmu/emulator.h"
#include "libEmu/regression.h"
#include "libEmu/emulate-fns.h"
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

/*********************************************************************
Allocates emulator_struct.
*
* ccs, should be thread safe now
*********************************************************************/
emulator_struct * alloc_emulator_struct(modelstruct * model)
{
	emulator_struct * e = (emulator_struct *)malloc(sizeof(emulator_struct));
	double determinant_c	= 0.0;
	e->nparams = model->options->nparams;
	e->nmodel_points = model->options->nmodel_points;
	e->nregression_fns = model->options->nregression_fns;
	e->nthetas = model->options->nthetas;
	e->model = model;
	gsl_matrix * c_matrix = gsl_matrix_alloc(e->nmodel_points, e->nmodel_points);
	e->cinverse = gsl_matrix_alloc(e->nmodel_points, e->nmodel_points);
	e->beta_vector = gsl_vector_alloc(e->nregression_fns);
	e->h_matrix = gsl_matrix_alloc(e->nmodel_points, e->nregression_fns);
	gsl_matrix * temp_matrix = gsl_matrix_alloc(e->nmodel_points, e->nmodel_points);

	makeCovMatrix_es(c_matrix, e);
	gsl_matrix_memcpy(temp_matrix, c_matrix);
	chol_inverse_cov_matrix(model->options, temp_matrix, e->cinverse, &determinant_c);
	makeHMatrix_es(e->h_matrix, e);
	estimateBeta_es(e->beta_vector, e);

	gsl_matrix_free(c_matrix);
	gsl_matrix_free(temp_matrix);
	return e;
}

/*********************************************************************
Free emulator_struct.
*********************************************************************/
void free_emulator_struct(emulator_struct * e)
{
	gsl_matrix_free(e->cinverse);
	gsl_matrix_free(e->h_matrix);
	gsl_vector_free(e->beta_vector);
	free((void *)e);
}

/**
 * the following fns are wrapped versions of commonly used fns in emulator.h etc, this
 * should make manipulations a little cleaner without breaking libRbind
 * they also use the emulator_struct->model->fnptrs instead
 */

/**
 * make the matrix of regression fns HMatrix, (thread safe?)
 * uses modelstruct fnptr: (void)(makeHvector)(gsl_vector*, gsl_vector*, int);
 *
 * although emulator_struct contains a h_matrix, we might want this for somethign else, so it is given as the 1st arg
 * the h_matrix in e is *not set* by this function
 */
void makeHMatrix_es(gsl_matrix *h_matrix, emulator_struct *e)
{
	makeHMatrix_fnptr(h_matrix,
							e->model->xmodel,
							e->nmodel_points,
							e->nparams,
							e->nregression_fns,
							e->model->makeHVector);
}

/**
 * make a covariance matrix: makeCovMatrix_es, (thread safe?)
 * uses the modelstruct fnptr (double)(*covariance_fn)(gsl_vector*, gsl_vector*, gsl_vector*, int, int)
 */
void makeCovMatrix_es(gsl_matrix *cov_matrix, emulator_struct *e)
{
	makeCovMatrix_fnptr(cov_matrix,
											e->model->xmodel,
											e->model->thetas,
											e->nmodel_points,
											e->nthetas,
											e->nparams,
											e->model->covariance_fn);
}

/**
 * compute a kvector at point k(point, emulator_struct) := { c(point, x_1), ... , c(point, x_n) }^T
 *
 * threadsafe
 *
 * requires: e->model to be init
 * e->model->thetas needs to contain sensible values
 */
void makeKVector_es(gsl_vector *kvector, gsl_vector *point, emulator_struct *e)
{
	makeKVector_fnptr(kvector, e->model->xmodel, point,
										e->model->thetas, e->nmodel_points, e->nthetas, e->nparams,
										e->model->covariance_fn);
}

/**
 * estimate betas from an emulator_struct, no need to modify anything to make this threadsafe
 *
 * betas are set into the first vector, not in the emulator_struct
 *
 * @require e->cinverse and e->h_matrix to be init
 */
void estimateBeta_es(gsl_vector *beta_vector, emulator_struct *e)
{
	estimateBeta(beta_vector, e->h_matrix, e->cinverse, e->model->training_vector,
							 e->nmodel_points, e->nregression_fns);
}

/*********************************************************************
return mean and variance at a point.
thread-safe, uses fnptrs defined in e->modelstruct
*********************************************************************/
void emulate_point(emulator_struct* e, gsl_vector * point,  double * mean, double * variance)
{
	gsl_vector * kplus = gsl_vector_alloc(e->nmodel_points); //allocate variable storage
	gsl_vector * h_vector = gsl_vector_alloc(e->nregression_fns); //allocate variable storage

	makeKVector_es(kplus, point, e);
	/* use the fnptr in modelstruct */
	e->model->makeHVector(h_vector, point, e->nparams);

	(*mean) = makeEmulatedMean(e->cinverse, e->model->training_vector,
		kplus, h_vector, e->h_matrix, e->beta_vector, e->nmodel_points);
	double kappa = e->model->covariance_fn(point, point, e->model->thetas, e->nthetas, e->nparams);
	(*variance) = makeEmulatedVariance(e->cinverse, kplus, h_vector,
		e->h_matrix, kappa, e->nmodel_points, e->nregression_fns);
	gsl_vector_free(kplus);
	gsl_vector_free(h_vector);

	// can we really return on void?
	return;
}

ext: c

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
/* stable/stable_integration.h
 *
 * Functions to perform numerical integration used when calculating
 * the PDF and CDF of alpha-stable distributions. Based on GSL
 * numerical quadrature methods.
 *
 * Copyright (C) 2013. Javier Royuela del Val
 *                     Federico Simmross Wattenberg
 *
 * 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; version 3 of the License.
 *
 * 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, see <http://www.gnu.org/licenses/>.
 *
 *
 *  Javier Royuela del Val.
 *  E.T.S.I. Telecomunicación
 *  Universidad de Valladolid
 *  Paseo de Belén 15, 47002 Valladolid, Spain.
 *  jroyval@lpi.tel.uva.es
 *
 */
#ifndef _STABLE_INTEGRATION_H_
#define _STABLE_INTEGRATION_H_

#include "stable.h"

#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_math.h>

int stable_integration_METHODNAME(unsigned short method, char *s);

void
stable_integration(StableDist *dist,double(function)(double, void*),
                   double a, double b,
                   double epsabs, double epsrel, unsigned short limit,
                   double *result, double *abserr, unsigned short method);

#endif

ext: h

max_stars_repo_name: TantraLabs/gonum

max_issues_repo_name: TantraLabs/gonum

max_forks_repo_name: TantraLabs/gonum


text:
/* rng/ran0.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 *
 * 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 Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_rng.h>

/* This is an implementation of the algorithm used in Numerical
   Recipe's ran0 generator. It is the same as MINSTD with an XOR mask
   of 123459876 on the seed.

   The period of this generator is 2^31.

   Note, if you choose a seed of 123459876 it would give a degenerate
   series 0,0,0,0, ...  I've made that into an error. */

static inline unsigned long int ran0_get (void *vstate);
static double ran0_get_double (void *vstate);
static void ran0_set (void *state, unsigned long int s);

static const long int m = 2147483647, a = 16807, q = 127773, r = 2836;
static const unsigned long int mask = 123459876;

typedef struct
  {
    unsigned long int x;
  }
ran0_state_t;

static inline unsigned long int
ran0_get (void *vstate)
{
  ran0_state_t *state = (ran0_state_t *) vstate;

  const unsigned long int x = state->x;

  const long int h = x / q;
  const long int t = a * (x - h * q) - h * r;

  if (t < 0)
    {
      state->x = t + m;
    }
  else
    {
      state->x = t;
    }

  return state->x;
}

static double
ran0_get_double (void *vstate)
{
  return ran0_get (vstate) / 2147483647.0 ;
}

static void
ran0_set (void *vstate, unsigned long int s)
{
  ran0_state_t *state = (ran0_state_t *) vstate;

  if (s == mask)
    {
      GSL_ERROR_VOID ("ran0 should not use seed == mask",
                                GSL_EINVAL);
    }

  state->x = s ^ mask;

  return;
}

static const gsl_rng_type ran0_type =
{"ran0",                        /* name */
 2147483646,                    /* RAND_MAX */
 1,                             /* RAND_MIN */
 sizeof (ran0_state_t),
 &ran0_set,
 &ran0_get,
 &ran0_get_double};

const gsl_rng_type *gsl_rng_ran0 = &ran0_type;

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: manggoguy/parsec-modified

max_forks_repo_name: manggoguy/parsec-modified


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_math.h> // gsl_finite

#ifdef _OPENMP
#include <omp.h>
#endif

const size_t N=10000;
const size_t T=1000;
const size_t K=20;
const size_t L=2;
const double alpha_k = 0.5;

/*
 * 0: thompson sampling
 * 1: uniform subset
 *
 */
const unsigned int strategy = 0;

#include "helpers.h"
#include "model.h"
#include "set_counter.h"

#define to_string(arr, k) \
    for (size_t i = 0; i < k-1; i++) { \
        printf("%.16lf,", arr[i]); \
    } \
    printf("%.16lf", arr[k-1]);

unsigned int
move_gibbs(double *restrict random_numbers, double logth[K], size_t ngames,
           const size_t games[ngames][L], const size_t game_counts[ngames],
           const size_t win_counts[K])
{
    double ll, ll_p;
    double alpha;
    double logth_comp_old;
    double logth_Km1_old;
    unsigned int accepted = 0;

    for (size_t comp=0; comp<K-1; comp++) {
        assert(gsl_fcmp(log_sum_exp(logth, K), 1.0, 1e-15) == 0);
        ll = fullcond(comp, logth, ngames, games, game_counts, win_counts);

        /* sample a suitable value for the current component */
        logth_comp_old = logth[comp];
        logth_Km1_old = logth[K-1];

        if (logth_comp_old > logth[K-1]) {
            logth[comp] = log(*random_numbers++) + logth_comp_old + log1p(exp(logth[K-1] - logth_comp_old));
            logth[K-1] = logth_comp_old + log1p(exp(logth_Km1_old - logth_comp_old) - exp(logth[comp] - logth_comp_old));
        } else {
            logth[comp] = log(*random_numbers++) + logth_Km1_old + log1p(exp(logth_comp_old - logth_Km1_old));
            logth[K-1] = logth_Km1_old + log1p(exp(logth_comp_old - logth_Km1_old) - exp(logth[comp] - logth_Km1_old));
        }

        /* compute full conditional density at th_p */
        ll_p = fullcond(comp, logth, ngames, games, game_counts, win_counts);

        alpha = *random_numbers++;
        if (log(alpha) < ll_p - ll) {
            /* accept */
            accepted = 1;
        }
        else {
            /* reject */
            /* reset the proposed component back to its original value */
            /* th[K-1] += th[comp] - th_comp_old; */
            /* th[comp] = th_comp_old; */
            /* assert(th[K-1] >= 0); */
            logth[comp] = logth_comp_old;
            logth[K-1] = logth_Km1_old;
        }
    }

    return accepted;
}

void
resample_move(const gsl_rng *r, double logtheta[N][K], const double w[N],
              const struct set_counter *games_counter, const size_t wins[K])
{
    unsigned int cnt[N];
    size_t accepted = 0;

    gsl_ran_multinomial(r, N, N, w, cnt);

    /* populate particles */
    double (*logtheta_new)[K] = malloc(N * sizeof *logtheta_new);
    size_t n_new = 0;
    for (size_t n=0; n<N; n++) {
        for (size_t i=0; i < cnt[n]; i++) {
            memcpy(logtheta_new[n_new++], logtheta[n], sizeof *logtheta_new);
        }
    }

    /* pre-generate random numbers to avoid thread synchronization */
    double (*random_numbers)[2*(K-1)] = malloc(N * sizeof *random_numbers);
    for (size_t n=0; n<N; n++) {
        for (size_t k=0; k<2*(K-1); k++) {
            random_numbers[n][k] = gsl_rng_uniform_pos(r);
        }
    }

    /* extract game counts from the counter */
    size_t ngames = games_counter->size;
    size_t (*games)[L] = malloc(ngames * sizeof *games);
    size_t *game_counts = malloc(ngames * sizeof *game_counts);
    set_counter_keys(games_counter, games);
    set_counter_values(games_counter, game_counts);

#pragma omp parallel for reduction(+:accepted)
    for (size_t n=0; n<N; n++) {
        accepted += move_gibbs(random_numbers[n], logtheta_new[n], ngames, games, game_counts, wins);
    }

    printf("# to_move = %zu\n", N);
    printf("# accepted = %zu\n", accepted);
    printf("# acceptance ratio = %lf\n", (double) accepted / N);

    memcpy(logtheta, logtheta_new, N*K*sizeof(double));
    free(logtheta_new);
    free(random_numbers);
    free(games);
    free(game_counts);
}

void
sample_theta_star(const gsl_rng *r, double theta_star[K])
{
    double a[K];
    for (size_t k=0; k<K; k++) {
        a[k] = alpha_k;
    }

    gsl_ran_dirichlet(r, K, a, theta_star);
}

void
read_theta_star(const char *file_name, double theta_star[K])
{
    char buf[80];
    FILE *ts = fopen(file_name, "r");

    if (!ts) {
        fprintf(stderr, "error reading %s\n", file_name);
        exit(EXIT_FAILURE);
    }

    for (size_t k = 0; k<K; k++) {
        if (!fgets(buf, 80, ts)) {
            fprintf(stderr, "error reading %s\n", file_name);
            exit(EXIT_FAILURE);
        }
        theta_star[k] = atof(buf);
    };
    fclose(ts);
}

void
sim(const gsl_rng *r, const double theta_star[K])
{
    double (*logtheta)[K] = malloc(N * sizeof *logtheta);
    double *w = malloc(N * sizeof(double));
    double *logw = malloc(N * sizeof(double));

    ones(w, N);
    zeros(logw, N);

    size_t *wins = calloc(K, sizeof *wins);
    struct set_counter *games_counter = set_counter_alloc();

    /* general info */
    printf("# generator type: %s\n", gsl_rng_name(r));
    printf("# seed = %lu\n", gsl_rng_default_seed);
    printf("\n");

    /* sample N particles from the `uniform` prior */
    {
        double alpha[K];
        ones(alpha, K);
        double theta[K];
        for (size_t n = 0; n < N; n++) {
            gsl_ran_dirichlet(r, K, alpha, theta);
#pragma omp simd
            for (size_t k=0; k<K; k++)
                logtheta[n][k] = log(theta[k]);
        }
    }

    for(size_t t = 0; t < T; t++) {
        fprintf(stderr, "t = %zu\r", t); /* for progress monitoring */
        printf("# iteration = %zu\n", t);

        size_t players[L];
        if (strategy == 0) {
            /* presentation strategy: thompson sampling */
            printf("# strategy: thompson\n");
            /* sample a theta from the current posterior */
            gsl_ran_discrete_t *g = gsl_ran_discrete_preproc(N, w);
            size_t theta_sample_idx = gsl_ran_discrete(r, g);
            gsl_ran_discrete_free(g);

            printf("# sampled theta: ");
            to_string(logtheta[theta_sample_idx], K);
            printf("\n");

            /* pick L elements from current sample */
            gsl_sort_largest_index(players, L, logtheta[theta_sample_idx], 1, K);
        } else if (strategy == 1) {
            /* presentation strategy: uniform subset */
            printf("# strategy: uniform subset\n");
            size_t idx[K];
            for (size_t k=0; k<K; k++)
                idx[k] = k;

            gsl_ran_choose(r, players, L, idx, K, sizeof(size_t));
        }

        set_counter_add(games_counter, players);
        printf("# number of unique subsets so far: %zu\n", games_counter->size);

        double player_w[L];

        printf("# game: ");
        for (size_t l=0; l<L-1; l++) {
            printf("%zu,", players[l]);
            player_w[l] = theta_star[players[l]];
        }
        printf("%zu\n", players[L-1]);
        player_w[L-1] = theta_star[players[L-1]];

        printf("# player weights = ");
        to_string(player_w, L);
        printf("\n");

        /* determine outcome using theta_star */
        size_t winner;
        {
            gsl_ran_discrete_t *g = gsl_ran_discrete_preproc(L, player_w);
            size_t wn = gsl_ran_discrete(r, g);
            winner = players[wn];
            gsl_ran_discrete_free(g);
        }
        printf("# winner: %zu\n", winner);

        wins[winner]++;

        /* update weights */
#pragma omp parallel for
        for(size_t n = 0; n < N; n++) {
            double logtheta_winner = logtheta[n][winner];
            double lth_game[L];
            for (size_t l=0; l<L; l++) {
                lth_game[l] = logtheta[n][players[l]];
            }
            logw[n] += logtheta_winner - log_sum_exp(lth_game, L);
        }

        /* compute w from logw */
        {
            double lse = log_sum_exp(logw, N);
            for (size_t n=0; n<N; n++) {
                w[n] = exp(logw[n] - lse);
                assert(gsl_finite(w[n]) == 1);
            }
        }

        /* compute ess and perform resampling if necessary */
        {
            double two_logw[N];
            for (size_t n=0; n<N; n++)
                two_logw[n] = 2*logw[n];

            double ess = exp(2*log_sum_exp(logw, N) - log_sum_exp(two_logw, N));

            printf("# ess = %lf\n", ess);

            if (ess < .5*N) {
                printf("# resampling at iteration %zu\n", t);
                resample_move(r, logtheta, w, games_counter, wins);
                ones(w, N);
                zeros(logw, N);
            }
        }

        printf("\n");
    }
    fprintf(stderr, "\n");

    /* resample at the end  */
    printf("# resampling at iteration %zu\n", T);
    resample_move(r, logtheta, w, games_counter, wins);
    /* no need to reset the weights at this point but just to be safe... */
    ones(w, N);
    zeros(logw, N);

    for(size_t n = 0; n < N; n++) {
        to_string(logtheta[n], K);
        printf("\n");
    }

    /* cleanup */
    free(logtheta);
    free(w);
    free(logw);
    free(wins);
    set_counter_free(games_counter);
}

int
main(int argc, char *argv[])
{
    const gsl_rng_type *t;
    gsl_rng *r;

    gsl_rng_env_setup();

    t = gsl_rng_default;
    r = gsl_rng_alloc(t);

    gsl_set_error_handler_off();

    double *theta_star = malloc(K * sizeof(double));
    /* read theta_star from a file */
    if (argc==2)
        read_theta_star(argv[1], theta_star);
    else
        sample_theta_star(r, theta_star);

    printf("# K = %zu\n", K);
    printf("# N = %zu\n", N);
    printf("# T = %zu\n", T);
    printf("# L = %zu\n", L);
    printf("# theta_star = ");
    to_string(theta_star, K);
    printf("\n");

    // perform simulation
    sim(r, theta_star);

    // cleanup
    free(theta_star);
    gsl_rng_free(r);

    exit(EXIT_SUCCESS);
}

max_stars_repo_name: ilkerg/preference-sampler

max_issues_repo_name: ilkerg/preference-sampler

max_forks_repo_name: ilkerg/preference-sampler


text:
#ifndef LDA_DOC_TOPIC_PROB_SELECT_H
#define LDA_DOC_TOPIC_PROB_SELECT_H

#include "Lambda.h"
#include "LambdaCreationFunctions.h"
#include "SelectionComp.h"
#include "PDBVector.h"
#include "IntIntVectorPair.h"
#include "IntDoubleVectorPair.h"
#include "DocAssignment.h"
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <random>
#include <gsl/gsl_vector.h>

/* This class samples the topic probability for the document */
using namespace pdb;
class LDADocTopicProbSelection : public SelectionComp<IntDoubleVectorPair, DocAssignment> {

private:
    Vector<double> prior;
    Handle<Vector<char>> myMem;
    int topicNum;

public:
    ENABLE_DEEP_COPY

    LDADocTopicProbSelection() {}

    LDADocTopicProbSelection(Vector<double>& fromPrior) {

        this->prior = fromPrior;
        this->topicNum = fromPrior.size();

	/* Set up the random number generator */
        gsl_rng* src = gsl_rng_alloc(gsl_rng_mt19937);
        std::random_device rd;
        std::mt19937 gen(rd());
        gsl_rng_set(src, gen());

        int spaceNeeded = sizeof(gsl_rng) + src->type->size;
        myMem = makeObject<Vector<char>>(spaceNeeded, spaceNeeded);

        memcpy(myMem->c_ptr(), src, sizeof(gsl_rng));
        memcpy(myMem->c_ptr() + sizeof(gsl_rng), src->state, src->type->size);

        gsl_rng_free(src);
    }

    Lambda<bool> getSelection(Handle<DocAssignment> checkMe) override {
        return makeLambda(checkMe, [](Handle<DocAssignment>& checkMe) { return true; });
    }

    Lambda<Handle<IntDoubleVectorPair>> getProjection(Handle<DocAssignment> checkMe) override {
        return makeLambda(checkMe, [&](Handle<DocAssignment>& checkMe) {

            gsl_rng* rng = getRng();
            Vector<unsigned>& topicCount = checkMe->getVector();
            Handle<Vector<double>> mySamples = makeObject<Vector<double>>(topicNum, topicNum);
            double* totalProb = new double[topicNum];
            for (int i = 0; i < topicNum; i++) {
                totalProb[i] = (this->prior)[i] + topicCount[i];
            }

	    /* Sample the topic probability */
            gsl_ran_dirichlet(rng, topicNum, totalProb, mySamples->c_ptr());

            Handle<IntDoubleVectorPair> result =
                makeObject<IntDoubleVectorPair>(checkMe->getKey(), mySamples);

            delete[] totalProb;

            return result;
        });
    }

    gsl_rng* getRng() {
        gsl_rng* dst = (gsl_rng*)myMem->c_ptr();
        dst->state = (void*)(myMem->c_ptr() + sizeof(gsl_rng));
        dst->type = gsl_rng_mt19937;
        return dst;
    }
};

#endif


text:
/**
 *
 * @file qwrapper_dlaswp.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mathieu Faverge
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:44:58 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlaswp(Quark *quark, Quark_Task_Flags *task_flags,
                       int n, double *A, int lda,
                       int i1,  int i2, const int *ipiv, int inc)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_dlaswp_quark, task_flags,
        sizeof(int),                      &n,    VALUE,
        sizeof(double)*lda*n,  A,        INOUT | LOCALITY,
        sizeof(int),                      &lda,  VALUE,
        sizeof(int),                      &i1,   VALUE,
        sizeof(int),                      &i2,   VALUE,
        sizeof(int)*n,                     ipiv,     INPUT,
        sizeof(int),                      &inc,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlaswp_quark = PCORE_dlaswp_quark
#define CORE_dlaswp_quark PCORE_dlaswp_quark
#endif
void CORE_dlaswp_quark(Quark *quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    double *A;

    quark_unpack_args_7(quark, n, A, lda, i1, i2, ipiv, inc);
    LAPACKE_dlaswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlaswp_f2(Quark *quark, Quark_Task_Flags *task_flags,
                          int n, double *A, int lda,
                          int i1,  int i2, const int *ipiv, int inc,
                          double *fake1, int szefake1, int flag1,
                          double *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_dlaswp_f2_quark, task_flags,
        sizeof(int),                        &n,     VALUE,
        sizeof(double)*lda*n,    A,         INOUT | LOCALITY,
        sizeof(int),                        &lda,   VALUE,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*n,                       ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(double)*szefake1, fake1,     flag1,
        sizeof(double)*szefake2, fake2,     flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlaswp_f2_quark = PCORE_dlaswp_f2_quark
#define CORE_dlaswp_f2_quark PCORE_dlaswp_f2_quark
#endif
void CORE_dlaswp_f2_quark(Quark* quark)
{
    int n, lda, i1, i2, inc;
    int *ipiv;
    double *A;
    void *fake1, *fake2;

    quark_unpack_args_9(quark, n, A, lda, i1, i2, ipiv, inc, fake1, fake2);
    LAPACKE_dlaswp_work(LAPACK_COL_MAJOR, n, A, lda, i1, i2, ipiv, inc );
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlaswp_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, double *Aij,
                              int i1,  int i2, const int *ipiv, int inc, double *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_dlaswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(double)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(double)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_dlaswp_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(double)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(double)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlaswp_ontile_quark = PCORE_dlaswp_ontile_quark
#define CORE_dlaswp_ontile_quark PCORE_dlaswp_ontile_quark
#endif
void CORE_dlaswp_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    double *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_dlaswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlaswp_ontile_f2(Quark *quark, Quark_Task_Flags *task_flags,
                                 PLASMA_desc descA, double *Aij,
                                 int i1,  int i2, const int *ipiv, int inc,
                                 double *fake1, int szefake1, int flag1,
                                 double *fake2, int szefake2, int flag2)
{
    DAG_CORE_LASWP;
    QUARK_Insert_Task(
        quark, CORE_dlaswp_ontile_f2_quark, task_flags,
        sizeof(PLASMA_desc),                &descA, VALUE,
        sizeof(double)*1,        Aij,       INOUT | LOCALITY,
        sizeof(int),                        &i1,    VALUE,
        sizeof(int),                        &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),      ipiv,      INPUT,
        sizeof(int),                        &inc,   VALUE,
        sizeof(double)*szefake1, fake1, flag1,
        sizeof(double)*szefake2, fake2, flag2,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlaswp_ontile_f2_quark = PCORE_dlaswp_ontile_f2_quark
#define CORE_dlaswp_ontile_f2_quark PCORE_dlaswp_ontile_f2_quark
#endif
void CORE_dlaswp_ontile_f2_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    double *A;
    PLASMA_desc descA;
    void *fake1, *fake2;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, fake1, fake2);
    CORE_dlaswp_ontile(descA, i1, i2, ipiv, inc);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dswptr_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, double *Aij,
                              int i1,  int i2, const int *ipiv, int inc,
                              const double *Akk, int ldak)
{
    DAG_CORE_TRSM;
    QUARK_Insert_Task(
        quark, CORE_dswptr_ontile_quark, task_flags,
        sizeof(PLASMA_desc),              &descA, VALUE,
        sizeof(double)*1,      Aij,       INOUT | LOCALITY,
        sizeof(int),                      &i1,    VALUE,
        sizeof(int),                      &i2,    VALUE,
        sizeof(int)*(i2-i1+1)*abs(inc),    ipiv,      INPUT,
        sizeof(int),                      &inc,   VALUE,
        sizeof(double)*ldak,   Akk,       INPUT,
        sizeof(int),                      &ldak,  VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dswptr_ontile_quark = PCORE_dswptr_ontile_quark
#define CORE_dswptr_ontile_quark PCORE_dswptr_ontile_quark
#endif
void CORE_dswptr_ontile_quark(Quark *quark)
{
    int i1, i2, inc, ldak;
    int *ipiv;
    double *A, *Akk;
    PLASMA_desc descA;

    quark_unpack_args_8(quark, descA, A, i1, i2, ipiv, inc, Akk, ldak);
    CORE_dswptr_ontile(descA, i1, i2, ipiv, inc, Akk, ldak);
}

/***************************************************************************//**
 *
 **/
void QUARK_CORE_dlaswpc_ontile(Quark *quark, Quark_Task_Flags *task_flags,
                              PLASMA_desc descA, double *Aij,
                              int i1,  int i2, const int *ipiv, int inc, double *fakepanel)
{
    DAG_CORE_LASWP;
    if (fakepanel == Aij) {
        QUARK_Insert_Task(
            quark, CORE_dlaswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(double)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(double)*1,      fakepanel,     SCRATCH,
            0);
    } else {
        QUARK_Insert_Task(
            quark, CORE_dlaswpc_ontile_quark, task_flags,
            sizeof(PLASMA_desc),              &descA,     VALUE,
            sizeof(double)*1,      Aij,           INOUT | LOCALITY,
            sizeof(int),                      &i1,        VALUE,
            sizeof(int),                      &i2,        VALUE,
            sizeof(int)*(i2-i1+1)*abs(inc),   ipiv,           INPUT,
            sizeof(int),                      &inc,       VALUE,
            sizeof(double)*1,      fakepanel,     INOUT,
            0);
    }
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_dlaswpc_ontile_quark = PCORE_dlaswpc_ontile_quark
#define CORE_dlaswpc_ontile_quark PCORE_dlaswpc_ontile_quark
#endif
void CORE_dlaswpc_ontile_quark(Quark *quark)
{
    int i1, i2, inc;
    int *ipiv;
    double *A, *fake;
    PLASMA_desc descA;

    quark_unpack_args_7(quark, descA, A, i1, i2, ipiv, inc, fake);
    CORE_dlaswpc_ontile(descA, i1, i2, ipiv, inc);
}


text:
/*	$Id$ */
/*
 * Copyright (c) 2014 Kristaps Dzonsons <kristaps@kcons.eu>
 *
 * Permission to use, copy, modify, and distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */
#include <assert.h>
#include <ctype.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <gtk/gtk.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_histogram.h>
#include <kplot.h>

#include "extern.h"

#define	BUFSZ 1024

/*
 * All possible input tokens, including some which are "virtual" tokens
 * in that they don't map to an exact character representation.
 */
enum	token {
	TOKEN_PAREN_OPEN,
	TOKEN_PAREN_CLOSE,
	TOKEN_ADD,
	TOKEN_SUB,
	TOKEN_MUL,
	TOKEN_DIV,
	TOKEN_EXP,
	TOKEN_END,
	TOKEN_NUMBER,
	TOKEN_SQRT,
	TOKEN_EXPF,
	TOKEN_P1,
	TOKEN_PN,
	TOKEN_N,
	TOKEN_ERROR,
	TOKEN_SKIP,
	TOKEN_POSITIVE,
	TOKEN_NEGATIVE,
	TOK__MAX
};

enum	arity {
	ARITY_NONE = 0,
	ARITY_ONE,
	ARITY_TWO
};

enum	assoc {
	ASSOC_NONE = 0,
	ASSOC_L,
	ASSOC_R,
};

struct	tok {
	char		 key; /* input identifier */
	int		 prec; /* precedence (if applicable) */
	int		 oper; /* is operator? */
	int		 func; /* is function? */
	enum assoc	 assoc; /* associativity */
	enum arity	 arity; /* n-ary-ess */
	enum htype	 map; /* map to htype (if applicable) */
};

static	const struct tok toks[TOK__MAX] = {
	{ '(',  -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE__MAX }, /* TOKEN_PAREN_OPEN */
	{ ')',  -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE__MAX }, /* TOKEN_PAREN_CLOSE */
	{ '+',  2, 1, 0, ASSOC_L, ARITY_TWO, HNODE_ADD }, /* TOKEN_ADD */
	{ '-',  2, 1, 0, ASSOC_L, ARITY_TWO, HNODE_SUB }, /* TOKEN_SUB */
	{ '*',  3, 1, 0, ASSOC_L, ARITY_TWO, HNODE_MUL }, /* TOKEN_MUL */
	{ '/',  3, 1, 0, ASSOC_L, ARITY_TWO, HNODE_DIV }, /* TOKEN_DIV */
	{ '^',  5, 1, 0, ASSOC_R, ARITY_TWO, HNODE_EXP }, /* TOKEN_EXP */
	{ '\0', -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE__MAX }, /* TOKEN_END */
	{ '\0', -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE_NUMBER }, /* TOKEN_NUMBER */
	{ '\0', -1, 0, 1, ASSOC_NONE, ARITY_NONE, HNODE_SQRT }, /* TOKEN_SQRT */
	{ '\0', -1, 0, 1, ASSOC_NONE, ARITY_NONE, HNODE_EXPF }, /* TOKEN_EXPF */
	{ 'x',  -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE_P1 }, /* TOKEN_P1 */
	{ 'X',  -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE_PN }, /* TOKEN_PN */
	{ 'n',  -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE_N }, /* TOKEN_N */
	{ '\0', -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE__MAX }, /* TOKEN_ERROR */
	{ ' ',  -1, 0, 0, ASSOC_NONE, ARITY_NONE, HNODE__MAX }, /* TOKEN_SKIP */
	{ '+',  4, 1, 0, ASSOC_R, ARITY_ONE, HNODE_POSITIVE }, /* TOKEN_POSITIVE */
	{ '-',  4, 1, 0, ASSOC_R, ARITY_ONE, HNODE_NEGATIVE }, /* TOKEN_NEGATIVE */
};

#if 0
static void
hnode_print(struct hnode **p)
{
	struct hnode	**pp;

	if (NULL == p)
		return;

	for (pp = p; NULL != *pp; pp++) {
		switch ((*pp)->type) {
		case (HNODE_VAR):
			putchar('v');
			break;
		case (HNODE_P1):
			putchar('x');
			break;
		case (HNODE_P2):
			putchar('y');
			break;
		case (HNODE_NUMBER):
			printf("%g", (*pp)->real);
			break;
		case (HNODE_SQRT):
			printf("sqrt");
			break;
		case (HNODE_ADD):
			putchar('+');
			break;
		case (HNODE_EXP):
			putchar('^');
			break;
		case (HNODE_SUB):
			putchar('-');
			break;
		case (HNODE_MUL):
			putchar('*');
			break;
		case (HNODE_DIV):
			putchar('/');
			break;
		case (HNODE_POSITIVE):
			printf("+'");
			break;
		case (HNODE_NEGATIVE):
			printf("-'");
			break;
		default:
			abort();
		}
		putchar(' ');
	}
	putchar('\n');
}
#endif

/*
 * Check if the current token (which has ARITY_ONE, we assume) is in
 * fact unary.
 * We do this by checking the previous token: if it's a binary operator
 * or an open parenthesis, then we're a unary operator.
 */
static int
check_unary(enum token lasttok)
{
	switch (lasttok) {
	case (TOK__MAX):
		/* FALLTHROUGH */
	case (TOKEN_PAREN_OPEN):
		/* FALLTHROUGH */
	case (TOKEN_ADD):
		/* FALLTHROUGH */
	case (TOKEN_SUB):
		/* FALLTHROUGH */
	case (TOKEN_MUL):
		/* FALLTHROUGH */
	case (TOKEN_DIV):
		/* FALLTHROUGH */
	case (TOKEN_EXP):
		return(1);
	default:
		break;
	}
	return(0);
}

/*
 * Check if the current token (which has ARITY_TWO, we assume) is in
 * fact binary.
 * We do this by checking the previous token: if it's an operand or a
 * close-paren (i.e., an operand), then we're a binary operator.
 */
static int
check_binary(enum token lasttok)
{
	switch (lasttok) {
	case (TOKEN_PAREN_CLOSE):
		/* FALLTHROUGH */
	case (TOKEN_NUMBER):
		/* FALLTHROUGH */
	case (TOKEN_P1):
		/* FALLTHROUGH */
	case (TOKEN_PN):
		/* FALLTHROUGH */
	case (TOKEN_N):
		return(1);
	default:
		break;
	}
	return(0);
}

/*
 * Extra a token from input.
 * If TOKEN_ERROR, something bad has happened in the attempt to do so.
 * Otherwise, this returns a valid token (possibly end of input).
 */
static enum token
tokenise(const char **v, char *buf, enum token lasttok)
{
	size_t		i, sz;

	/* Short-circuit this case. */
	if ('\0' == **v)
		return(TOKEN_END);

	/* Look for token in our predefined inputs. */
	for (i = 0; i < TOK__MAX; i++)
		if ('\0' != toks[i].key && **v == toks[i].key) {
			switch (toks[i].arity) {
			case (ARITY_NONE):
				(*v)++;
				return((enum token)i);
			case (ARITY_ONE):
				if ( ! check_unary(lasttok))
					continue;
				break;
			case (ARITY_TWO):
				if ( ! check_binary(lasttok))
					continue;
				break;
			}
			(*v)++;
			return((enum token)i);
		}

	/* See if we're a real number or identifier. */
	if (isdigit((int)**v) || '.' == **v) {
		sz = 0;
		for ( ; isdigit((int)**v) || '.' == **v; (*v)++) {
			assert(sz < BUFSZ);
			buf[sz++] = **v;
		}
		buf[sz] = '\0';
		return(TOKEN_NUMBER);
	} else if (isalpha((int)**v)) {
		sz = 0;
		for ( ; isalpha((int)**v); (*v)++) {
			assert(sz < BUFSZ);
			buf[sz++] = **v;
		}
		buf[sz] = '\0';
		if (0 == strcmp(buf, "sqrt"))
			return(TOKEN_SQRT);
		else if (0 == strcmp(buf, "exp"))
			return(TOKEN_EXPF);
	}

	/* Eh... */
	return(TOKEN_ERROR);
}

/*
 * Queue something onto the output RPN queue.
 */
static void
enqueue(struct hnode ***q, size_t *qsz, enum token tok)
{

	*q = realloc(*q, ++(*qsz) * sizeof(struct hnode *));
	(*q)[*qsz - 1] = calloc(1, sizeof(struct hnode));
	(*q)[*qsz - 1]->type = toks[tok].map;
	assert(HNODE__MAX != toks[tok].map);
}

/*
 * Make sure that all operators and functions are matched with
 * arguments.
 * The Shunting-Yard algorithm itself doesn't do this, so we need to do
 * it now.
 */
static int
check(struct hnode **p)
{
	struct hnode	**pp;
	size_t		  ssz;

	for (ssz = 0, pp = p; NULL != *pp; pp++)
		switch ((*pp)->type) {
		case (HNODE_P1):
			/* FALLTHROUGH */
		case (HNODE_PN):
			/* FALLTHROUGH */
		case (HNODE_N):
			/* FALLTHROUGH */
		case (HNODE_NUMBER):
			ssz++;
			break;
		case (HNODE_POSITIVE):
			/* FALLTHROUGH */
		case (HNODE_NEGATIVE):
			/* FALLTHROUGH */
		case (HNODE_SQRT):
			/* FALLTHROUGH */
		case (HNODE_EXPF):
			if (0 == ssz)
				return(0);
			break;
		case (HNODE_ADD):
			/* FALLTHROUGH */
		case (HNODE_SUB):
			/* FALLTHROUGH */
		case (HNODE_MUL):
			/* FALLTHROUGH */
		case (HNODE_DIV):
			/* FALLTHROUGH */
		case (HNODE_EXP):
			if (ssz < 2)
				return(0);
			ssz--;
			break;
		default:
			abort();
		}

	return(1 == ssz);
}

/*
 * Dijkstra's Shunting-Yard algorithm for converting an infix-order
 * expression to prefix order.
 * This returns a NULL-terminated list of expressions to evaluate in
 * post-fix order.
 */
struct hnode **
hnode_parse(const char **v)
{
	enum token	  tok, lasttok;
	enum token	  stack[STACKSZ];
	char		  buf[BUFSZ + 1];
	struct hnode	**q;
	int		  found;
	size_t		  i, qsz, ssz;

	q = NULL;
	qsz = ssz = 0;
	lasttok = TOK__MAX;

	while (TOKEN_END != (tok = tokenise(v, buf, lasttok))) {
		switch (tok) {
		case (TOKEN_ERROR):
			goto err;
		case (TOKEN_SKIP):
			break;
		case (TOKEN_P1):
			/* FALLTHROUGH */
		case (TOKEN_PN):
			/* FALLTHROUGH */
		case (TOKEN_N):
			enqueue(&q, &qsz, tok);
			break;
		case (TOKEN_NUMBER):
			enqueue(&q, &qsz, tok);
			q[qsz - 1]->real = atof(buf);
			break;
		case (TOKEN_SQRT):
			/* FALLTHROUGH */
		case (TOKEN_EXPF):
			/* FALLTHROUGH */
		case (TOKEN_PAREN_OPEN):
			assert(ssz < STACKSZ);
			stack[ssz++] = tok;
			break;
		case (TOKEN_PAREN_CLOSE):
			assert(ssz > 0);
			found = 0;
			do {
				if (TOKEN_PAREN_OPEN == stack[--ssz])
					found = 1;
				else
					enqueue(&q, &qsz, stack[ssz]);
			} while ( ! found && ssz > 0);
			assert(found);
			if (ssz > 0 && toks[stack[ssz - 1]].func)
				enqueue(&q, &qsz, stack[--ssz]);
			break;
		case (TOKEN_POSITIVE):
			/* FALLTHROUGH */
		case (TOKEN_NEGATIVE):
			/* FALLTHROUGH */
		case (TOKEN_SUB):
			/* FALLTHROUGH */
		case (TOKEN_MUL):
			/* FALLTHROUGH */
		case (TOKEN_DIV):
			/* FALLTHROUGH */
		case (TOKEN_EXP):
			/* FALLTHROUGH */
		case (TOKEN_ADD):
			assert(toks[tok].prec >= 0);
			assert(ASSOC_NONE != toks[tok].assoc);
			while (ssz > 0 && toks[stack[ssz - 1]].oper) {
				assert(toks[stack[ssz - 1]].prec >= 0);
				assert(ASSOC_NONE !=
					toks[stack[ssz - 1]].assoc);
				if ((ASSOC_L == toks[tok].assoc
					 && toks[tok].prec ==
					 toks[stack[ssz - 1]].prec) ||
					 (toks[tok].prec <
					  toks[stack[ssz - 1]].prec))
					enqueue(&q, &qsz, stack[--ssz]);
				else
					break;
			}
			assert(ssz < STACKSZ);
			stack[ssz++] = tok;
			break;
		default:
			abort();
		}
		if (TOKEN_SKIP != tok)
			lasttok = tok;
	}

	while (ssz > 0) {
		--ssz;
		if (TOKEN_PAREN_OPEN == stack[ssz])
			goto err;
		enqueue(&q, &qsz, stack[ssz]);
	}

	q = realloc(q, ++qsz * sizeof(struct hnode *));
	q[qsz - 1] = NULL;

	if ( ! check(q))
		goto err;

#if 0
	hnode_print(q);
#endif
	return(q);
err:
	for (i = 0; i < qsz; i++)
		free(q[i]);
	free(q);
	return(NULL);
}

struct hnode **
hnode_copy(struct hnode **p)
{
	struct hnode	**pp;
	size_t		  sz;

	for (sz = 0, pp = p; NULL != *pp; pp++)
		sz++;

	pp = calloc(sz + 1, sizeof(struct hnode *));
	for (sz = 0; NULL != *p; p++, sz++) {
		pp[sz] = calloc(1, sizeof(struct hnode));
		pp[sz]->type = (*p)->type;
		pp[sz]->real = (*p)->real;
	}

	return(pp);
}

void
hnode_free(struct hnode **p)
{
	struct hnode	**pp;

	if (NULL == p)
		return;

	for (pp = p; NULL != *pp; pp++)
		free(*pp);

	free(p);
}

/*
 * Execute a function in prefix order given variables x (given player's
 * strategy), y (other player's strategy), and var (given player's
 * morality index).
 */
double
hnode_exec(const struct hnode *const *p, double x, double X, size_t n)
{
	double		  stack[STACKSZ];
	const struct hnode *const *pp;
	size_t		  ssz;
	double		  val;

	for (ssz = 0, pp = p; NULL != *pp; pp++)
		switch ((*pp)->type) {
		case (HNODE_P1):
			assert(ssz < STACKSZ);
			stack[ssz++] = x;
			break;
		case (HNODE_PN):
			assert(ssz < STACKSZ);
			stack[ssz++] = X;
			break;
		case (HNODE_N):
			assert(ssz < STACKSZ);
			stack[ssz++] = (double)n;
			break;
		case (HNODE_NUMBER):
			assert(ssz < STACKSZ);
			stack[ssz++] = (*pp)->real;
			break;
		case (HNODE_SQRT):
			assert(ssz > 0);
			val = sqrt(stack[--ssz]);
			stack[ssz++] = val;
			break;
		case (HNODE_EXPF):
			assert(ssz > 0);
			val = exp(stack[--ssz]);
			stack[ssz++] = val;
			break;
		case (HNODE_ADD):
			assert(ssz > 1);
			val = stack[ssz - 2] + stack[ssz - 1];
			ssz -= 2;
			stack[ssz++] = val;
			break;
		case (HNODE_SUB):
			assert(ssz > 1);
			val = stack[ssz - 2] - stack[ssz - 1];
			ssz -= 2;
			stack[ssz++] = val;
			break;
		case (HNODE_MUL):
			assert(ssz > 1);
			val = stack[ssz - 2] * stack[ssz - 1];
			ssz -= 2;
			stack[ssz++] = val;
			break;
		case (HNODE_DIV):
			assert(ssz > 1);
			val = stack[ssz - 2] / stack[ssz - 1];
			ssz -= 2;
			stack[ssz++] = val;
			break;
		case (HNODE_EXP):
			assert(ssz > 1);
			val = pow(stack[ssz - 2], stack[ssz - 1]);
			ssz -= 2;
			stack[ssz++] = val;
			break;
		case (HNODE_POSITIVE):
			assert(ssz > 0);
			break;
		case (HNODE_NEGATIVE):
			assert(ssz > 0);
			val = -(stack[--ssz]);
			stack[ssz++] = val;
			break;
		default:
			abort();
		}

	assert(1 == ssz);
	return(stack[0]);
}

static void
hnode_test_expect(double x, double X,
	double n, const char *expf, double vexp)
{
	struct hnode	 **exp;
	double		   v;
	const char	  *expfp;

	expfp = expf;
	exp = hnode_parse((const char **)&expfp);
	assert(NULL != exp);
	v = hnode_exec
		((const struct hnode *const *)exp, x, X, n);
	g_debug("pi(x=%g, X=%g, n=%g) = %s = %g (want %g)",
		x, X, n, expf, v, vexp);
	hnode_free(exp);
}

void
hnode_test(void)
{
	double	x, X, n;

	x = 10.0;
	X = 20.0;
	n = 2.0;
	hnode_test_expect(x, X, n,
		"(1 - exp(-X)) - x",
		(1.0 - exp(-X)) - x);
	hnode_test_expect(x, X, n,
		"sqrt(1 / n * X) - 0.5 * x^2",
		sqrt(1.0 / n * X) - 0.5 * pow(x, 2.0));
	hnode_test_expect(x, X, n,
		"x - (X - x) * x - x^2",
		x - (X - x) * x - pow(x, 2.0));
	hnode_test_expect(x, X, n,
		"x * (1 / X) - x",
		x * (1.0 / X) - x);
}

max_stars_repo_name: kristapsdz/bmigrate

max_issues_repo_name: kristapsdz/bmigrate

max_forks_repo_name: kristapsdz/bmigrate


text:
/* specfunc/dilog.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include "gsl_sf_clausen.h"
#include "gsl_sf_trig.h"
#include "gsl_sf_log.h"
#include "gsl_sf_dilog.h"

/* Evaluate series for real dilog(x)
 * Sum[ x^k / k^2, {k,1,Infinity}]
 *
 * Converges rapidly for |x| < 1/2.
 */
static
int
dilog_series(const double x, gsl_sf_result * result)
{
  const int kmax = 1000;
  double sum  = x;
  double term = x;
  int k;
  for(k=2; k<kmax; k++) {
    double rk = (k-1.0)/k;
    term *= x;
    term *= rk*rk;
    sum += term;
    if(fabs(term/sum) < GSL_DBL_EPSILON) break;
  }

  result->val  = sum;
  result->err  = 2.0 * fabs(term);
  result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);

  if(k == kmax)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
}

/* Assumes x >= 0.0
 */
static
int
dilog_xge0(const double x, gsl_sf_result * result)
{
  if(x > 2.0) {
    const double log_x = log(x);
    gsl_sf_result ser;
    int stat_ser = dilog_series(1.0/x, &ser);
    double t1 = M_PI*M_PI/3.0;
    double t2 = ser.val;
    double t3 = 0.5*log_x*log_x;
    result->val  = t1 - t2 - t3;
    result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
    result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
    result->val += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_ser;
  }
  else if(x > 1.01) {
    const double log_x    = log(x);
    const double log_term = log_x * (log(1.0-1.0/x) + 0.5*log_x);
    gsl_sf_result ser;
    int stat_ser = dilog_series(1.0 - 1.0/x, &ser);
    double t1 = M_PI*M_PI/6.0;
    double t2 = ser.val;
    double t3 = log_term;
    result->val  = t1 + t2 - t3;
    result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
    result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_ser;
  }
  else if(x > 1.0) {
    /* series around x = 1.0 */
    const double eps = x - 1.0;
    const double lne = log(eps);
    const double c0 = M_PI*M_PI/6.0;
    const double c1 =   1.0 - lne;
    const double c2 = -(1.0 - 2.0*lne)/4.0;
    const double c3 =  (1.0 - 3.0*lne)/9.0;
    const double c4 = -(1.0 - 4.0*lne)/16.0;
    const double c5 =  (1.0 - 5.0*lne)/25.0;
    const double c6 = -(1.0 - 6.0*lne)/36.0;
    const double c7 =  (1.0 - 7.0*lne)/49.0;
    const double c8 = -(1.0 - 8.0*lne)/64.0;
    result->val = c0+eps*(c1+eps*(c2+eps*(c3+eps*(c4+eps*(c5+eps*(c6+eps*(c7+eps*c8)))))));
    result->err = 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_SUCCESS;
  }
  else if(x == 1.0) {
    result->val = M_PI*M_PI/6.0;
    result->err = 2.0 * GSL_DBL_EPSILON * M_PI*M_PI/6.0;
    return GSL_SUCCESS;
  }
  else if(x > 0.5) {
    const double log_x = log(x);
    gsl_sf_result ser;
    int stat_ser = dilog_series(1.0-x, &ser);
    double t1 = M_PI*M_PI/6.0;
    double t2 = ser.val;
    double t3 = log_x*log(1.0-x);
    result->val  = t1 - t2 - t3;
    result->err  = GSL_DBL_EPSILON * fabs(log_x) + ser.err;
    result->err += GSL_DBL_EPSILON * (fabs(t1) + fabs(t2) + fabs(t3));
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return stat_ser;
  }
  else if(x > 0.0) {
    return dilog_series(x, result);
  }
  else {
    /* x == 0.0 */
    result->val = 0.0;
    result->err = 0.0;
    return GSL_SUCCESS;
  }
}

/* Evaluate the series representation for Li2(z):
 *
 *   Li2(z) = Sum[ |z|^k / k^2 Exp[i k arg(z)], {k,1,Infinity}]
 *   |z|    = r
 *   arg(z) = theta
 *
 * Assumes 0 < r < 1.
 */
static
int
dilogc_series_1(double r, double cos_theta, double sin_theta,
                gsl_sf_result * real_result, gsl_sf_result * imag_result)
{
  double alpha = 1.0 - cos_theta;
  double beta  = sin_theta;
  double ck = cos_theta;
  double sk = sin_theta;
  double rk = r;
  double real_sum = r*ck;
  double imag_sum = r*sk;
  int kmax = 50 + (int)(22.0/(-log(r))); /* tuned for double-precision */
  int k;
  for(k=2; k<kmax; k++) {
    double ck_tmp = ck;
    ck = ck - (alpha*ck + beta*sk);
    sk = sk - (alpha*sk - beta*ck_tmp);
    rk *= r;
    real_sum += rk/((double)k*k) * ck;
    imag_sum += rk/((double)k*k) * sk;
  }

  real_result->val = real_sum;
  real_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(real_sum);
  imag_result->val = imag_sum;
  imag_result->err = 2.0 * kmax * GSL_DBL_EPSILON * fabs(imag_sum);

  return GSL_SUCCESS;
}

/* Evaluate a series for Li_2(z) when |z| is near 1.
 * This is uniformly good away from z=1.
 *
 *   Li_2(z) = Sum[ a^n/n! H_n(theta), {n, 0, Infinity}]
 *
 * where
 *   H_n(theta) = Sum[ e^(i m theta) m^n / m^2, {m, 1, Infinity}]
 *   a = ln(r)
 *
 *  H_0(t) = Gl_2(t) + i Cl_2(t)
 *  H_1(t) = 1/2 ln(2(1-c)) + I atan2(-s, 1-c)
 *  H_2(t) = -1/2 + I/2 s/(1-c)
 *  H_3(t) = -1/2 /(1-c)
 *  H_4(t) = -I/2 s/(1-c)^2
 *  H_5(t) = 1/2 (2 + c)/(1-c)^2
 *  H_6(t) = I/2 s/(1-c)^5 (8(1-c) - s^2 (3 + c))
 *
 *  assumes: 0 <= theta <= 2Pi
 */
static
int
dilogc_series_2(double r, double theta, double cos_theta, double sin_theta,
                gsl_sf_result * real_result, gsl_sf_result * imag_result)
{
  double a = log(r);
  double omc = 1.0 - cos_theta;
  double H_re[7];
  double H_im[7];
  double an, nfact;
  double sum_re, sum_im;
  gsl_sf_result Him0;
  int n;

  H_re[0] = M_PI*M_PI/6.0 + 0.25*(theta*theta - 2.0*M_PI*fabs(theta));
  gsl_sf_clausen_e(theta, &Him0);
  H_im[0] = Him0.val;

  H_re[1] = -0.5*log(2.0*omc);
  H_im[1] = -atan2(-sin_theta, omc);

  H_re[2] = -0.5;
  H_im[2] = 0.5 * sin_theta/omc;

  H_re[3] = -0.5/omc;
  H_im[3] = 0.0;

  H_re[4] = 0.0;
  H_im[4] = -0.5*sin_theta/(omc*omc);

  H_re[5] = 0.5 * (2.0 + cos_theta)/(omc*omc);
  H_im[5] = 0.0;

  H_re[6] = 0.0;
  H_im[6] = 0.5 * sin_theta/(omc*omc*omc*omc*omc)
            * (8*omc - sin_theta*sin_theta*(3 + cos_theta));

  sum_re = H_re[0];
  sum_im = H_im[0];
  an = 1.0;
  nfact = 1.0;
  for(n=1; n<=6; n++) {
    double t;
    an *= a;
    nfact *= n;
    t = an/nfact;
    sum_re += t * H_re[n];
    sum_im += t * H_im[n];
  }

  real_result->val = sum_re;
  real_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_re) + fabs(an/nfact);
  imag_result->val = sum_im;
  imag_result->err = 2.0 * 6.0 * GSL_DBL_EPSILON * fabs(sum_im) + Him0.err + fabs(an/nfact);

  return GSL_SUCCESS;
}

/* complex dilogarithm in the unit disk
 * assumes:  r < 1  and  0 <= theta <= 2Pi
 */
static
int
dilogc_unitdisk(double r, double theta, gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
{
  const double zeta2 = M_PI*M_PI/6.0;
  int stat_dilog;
  gsl_sf_result cos_theta;
  gsl_sf_result sin_theta;
  int stat_cos = gsl_sf_cos_e(theta, &cos_theta);
  int stat_sin = gsl_sf_sin_e(theta, &sin_theta);
  gsl_sf_result x;
  gsl_sf_result y;
  gsl_sf_result x_tmp, y_tmp, r_tmp;
  gsl_sf_result result_re_tmp, result_im_tmp;
  double cos_theta_tmp;
  double sin_theta_tmp;
  x.val = r * cos_theta.val;
  x.err = r * cos_theta.err;
  y.val = r * sin_theta.val;
  y.err = r * sin_theta.err;

  /* Reflect away from z = 1 if
   * we are too close.
   */
  if(x.val > 0.5) {
    x_tmp.val = 1.0 - x.val;
    x_tmp.err = GSL_DBL_EPSILON * (1.0 + fabs(x.val)) + x.err;
    y_tmp.val = -y.val;
    y_tmp.err = y.err;
    r_tmp.val = sqrt(x_tmp.val*x_tmp.val + y_tmp.val*y_tmp.val);
    r_tmp.err = (x_tmp.err*fabs(x_tmp.val) + y_tmp.err*fabs(y_tmp.val))/fabs(r_tmp.val);
  }
  else {
    x_tmp.val = x.val;
    x_tmp.err = x.err;
    y_tmp.val = y.val;
    y_tmp.err = y.err;
    r_tmp.val = r;
    r_tmp.err = r * GSL_DBL_EPSILON;
  }

  cos_theta_tmp = x_tmp.val / r_tmp.val;
  sin_theta_tmp = y_tmp.val / r_tmp.val;

  /* Calculate dilog of the transformed variable.
   */
  if(r_tmp.val < 0.98) {
    stat_dilog = dilogc_series_1(r_tmp.val, cos_theta_tmp, sin_theta_tmp,
				 &result_re_tmp, &result_im_tmp
				 );
  }
  else {
    double theta_tmp = atan2(y_tmp.val, x_tmp.val);
    stat_dilog = dilogc_series_2(r_tmp.val, theta_tmp, cos_theta_tmp, sin_theta_tmp,
                                 &result_re_tmp, &result_im_tmp
				 );
  }

  /* Unwind reflection if necessary.
   *
   * Li2(z) = -Li2(1-z) + zeta(2) - ln(z) ln(1-z)
   */
  if(x.val > 0.5) {
    double lnz    =  log(r);                         /*  log(|z|)   */
    double lnomz  =  log(r_tmp.val);                 /*  log(|1-z|) */
    double argz   =  theta;                          /*  arg(z)     */
    double argomz =  atan2(y_tmp.val, x_tmp.val);    /*  arg(1-z)   */
    real_dl->val  = -result_re_tmp.val + zeta2 - lnz*lnomz + argz*argomz;
    real_dl->err  =  result_re_tmp.err;
    real_dl->err +=  GSL_DBL_EPSILON * (zeta2 + fabs(lnz*lnomz) + fabs(argz*argomz));
    real_dl->err +=  2.0 * GSL_DBL_EPSILON * fabs(real_dl->val);
    imag_dl->val  = -result_im_tmp.val - argz*lnomz - argomz*lnz;
    imag_dl->err  =  result_im_tmp.err;
    imag_dl->err +=  GSL_DBL_EPSILON * (fabs(argz*lnomz) + fabs(argomz*lnz));
    imag_dl->err +=  2.0 * GSL_DBL_EPSILON * fabs(imag_dl->val);
  }
  else {
    real_dl->val = result_re_tmp.val;
    real_dl->err = result_re_tmp.err;
    imag_dl->val = result_im_tmp.val;
    imag_dl->err = result_im_tmp.err;
  }

  return GSL_ERROR_SELECT_3(stat_dilog, stat_sin, stat_cos);
}

/*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/

int
gsl_sf_dilog_e(const double x, gsl_sf_result * result)
{
  /* CHECK_POINTER(result) */

  if(x >= 0.0) {
    return dilog_xge0(x, result);
  }
  else {
    gsl_sf_result d1, d2;
    int stat_d1 = dilog_xge0( -x, &d1);
    int stat_d2 = dilog_xge0(x*x, &d2);
    result->val  = -d1.val + 0.5 * d2.val;
    result->err  =  d1.err + 0.5 * d2.err;
    result->err += 2.0 * GSL_DBL_EPSILON * fabs(result->val);
    return GSL_ERROR_SELECT_2(stat_d1, stat_d2);
  }
}

int
gsl_sf_complex_dilog_e(const double r, double theta,
                          gsl_sf_result * real_dl, gsl_sf_result * imag_dl)
{
  /* CHECK_POINTER(real_dl) */
  /* CHECK_POINTER(imag_dl) */

  if(r == 0.0) {
    real_dl->val = 0.0;
    real_dl->err = 0.0;
    imag_dl->val = 0.0;
    imag_dl->err = 0.0;
    return GSL_SUCCESS;
  }

/*
  if(theta < 0.0 || theta > 2.0*M_PI) {
    gsl_sf_angle_restrict_pos_e(&theta);
  }
*/

  /* Trap cases of real-valued argument.
   */
  if(theta == 0.0) {
    int stat_d;
    imag_dl->val = ( r > 1.0 ? -M_PI * log(r) : 0.0 );
    imag_dl->err = 2.0 * GSL_DBL_EPSILON * fabs(imag_dl->val);
    stat_d = gsl_sf_dilog_e(r, real_dl);
    return stat_d;
  }
  if(theta == M_PI) {
    int stat_d;
    imag_dl->val = 0.0;
    imag_dl->err = 0.0;
    stat_d = gsl_sf_dilog_e(-r, real_dl);
    return stat_d;
  }

  /* Trap unit circle case.
   */
  if(r == 1.0) {
    gsl_sf_result theta_restrict;
    int stat_r = gsl_sf_angle_restrict_pos_err_e(theta, &theta_restrict);
    int stat_c;
    const double term1 = theta_restrict.val*theta_restrict.val;
    const double term2 = 2.0*M_PI*fabs(theta_restrict.val);
    const double term1_err = 2.0 * fabs(theta_restrict.val * theta_restrict.err);
    const double term2_err = 2.0*M_PI*fabs(theta_restrict.err);
    real_dl->val  = M_PI*M_PI/6.0 + 0.25*(term1 - term2);
    real_dl->err  = 2.0 * GSL_DBL_EPSILON * (M_PI*M_PI/6.0 + 0.25 * (fabs(term1) + fabs(term2)));
    real_dl->err += 0.25 * (term1_err + term2_err);
    real_dl->err += 2.0 * GSL_DBL_EPSILON * fabs(real_dl->val);
    stat_c = gsl_sf_clausen_e(theta, imag_dl);
    stat_r = 0;  /* discard restrict status */
    return stat_c;
  }

  /* Generic case.
   */
  {
    int stat_dilog;
    double r_tmp, theta_tmp;
    gsl_sf_result result_re_tmp, result_im_tmp;

    /* Reduce argument to unit disk.
     */
    if(r > 1.0) {
      r_tmp     = 1.0 / r;
      theta_tmp = /* 2.0*M_PI */ - theta;
    }
    else {
      r_tmp     = r;
      theta_tmp = theta;
    }

    /* Calculate in the unit disk.
     */
    stat_dilog = dilogc_unitdisk(r_tmp, theta_tmp,
                                 &result_re_tmp, &result_im_tmp
				 );

    /* Unwind the inversion if necessary. We calculate
     * the imaginary part explicitly if using the inversion
     * because there is no simple relationship between
     * arg(1-z) and arg(1 - 1/z), which is the "omega"
     * term in [Lewin A.2.5 (1)].
     */
    if(r > 1.0) {
      const double zeta2 = M_PI*M_PI/6.0;
      double x = r * cos(theta);
      double y = r * sin(theta);
      double omega = atan2(y, 1.0-x);
      double lnr = log(r);
      double pmt = M_PI - theta;
      gsl_sf_result Cl_a, Cl_b, Cl_c;
      double r1, r2, r3, r4, r5;
      int stat_c1 = gsl_sf_clausen_e(2.0*omega, &Cl_a);
      int stat_c2 = gsl_sf_clausen_e(2.0*theta, &Cl_b);
      int stat_c3 = gsl_sf_clausen_e(2.0*(omega+theta), &Cl_c);
      int stat_c  = GSL_ERROR_SELECT_3(stat_c1, stat_c2, stat_c3);
      r1 = -result_re_tmp.val;
      r2 = -0.5*lnr*lnr;
      r3 =  0.5*pmt*pmt;
      r4 = -zeta2;
      r5 =  omega*lnr;
      real_dl->val  = r1 + r2 + r3 + r4;
      real_dl->err  = result_re_tmp.err;
      real_dl->err += GSL_DBL_EPSILON * (fabs(r1) + fabs(r2) + fabs(r3) + fabs(r4));
      real_dl->err += 2.0 * GSL_DBL_EPSILON * fabs(real_dl->val);
      imag_dl->val  = r5 + 0.5*(Cl_a.val + Cl_b.val - Cl_c.val);
      imag_dl->err  = GSL_DBL_EPSILON * fabs(r5);
      imag_dl->err += GSL_DBL_EPSILON * 0.5*(fabs(Cl_a.val) + fabs(Cl_b.val) + fabs(Cl_c.val));
      imag_dl->err += 0.5*(Cl_a.err + Cl_b.err + Cl_c.err);
      imag_dl->err += 2.0*GSL_DBL_EPSILON * fabs(imag_dl->val);
      return GSL_ERROR_SELECT_2(stat_dilog, stat_c);
    }
    else {
      real_dl->val = result_re_tmp.val;
      real_dl->err = result_re_tmp.err;
      imag_dl->val = result_im_tmp.val;
      imag_dl->err = result_im_tmp.err;
      return stat_dilog;
    }
  }
}

/*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/

#include "eval.h"

double gsl_sf_dilog(const double x)
{
  EVAL_RESULT(gsl_sf_dilog_e(x, &result));
}