SimpleFits.h

Go to the documentation of this file.

#ifndef SIMPLEFITS_H
#define SIMPLEFITS_H 1

// C/C++ standard libraries
#include <algorithm> // std::for_each()
#include <array>
#include <cmath>     // std::sqrt()
#include <iterator>  // std::begin(), std::end()
#include <ostream>   // std::endl
#include <stdexcept> // std::range_error
#include <tuple>
#include <type_traits> // std::enable_if<>, std::is_const<>

#include "lardata/Utilities/FastMatrixMathHelper.h" // lar::util::details::FastMatrixOperations
#include "lardataalg/Utilities/StatCollector.h"     // lar::util::identity

namespace lar {
  namespace util {

    namespace details {

      template <typename T, unsigned int D>
      class FitDataCollector {

        template <unsigned int Power, unsigned int N>
        struct SumExtractor {
          static T Sum(WeightTracker<T> const&, DataTracker<Power, T> const& sums)
          {
            return sums.template SumN<N>();
          }
        }; // SumExtractor

        // specialization
        template <unsigned int Power>
        struct SumExtractor<Power, 0U> {
          static T Sum(WeightTracker<T> const& weight, DataTracker<Power, T> const&)
          {
            return weight.Weights();
          }
        }; // SumExtractor<>

      public:
        static constexpr unsigned int Degree = D;
        static constexpr unsigned int NParams = Degree + 1;

        using Data_t = T; 

        using Measurement_t = std::tuple<Data_t, Data_t>;

        using MeasurementAndUncertainty_t = std::tuple<Data_t, Data_t, Data_t>;

        // default constructor, destructor and all the rest


        bool add(Data_t x, Data_t y, Data_t sy = Data_t(1.0));

        bool add(Measurement_t value, Data_t sy = Data_t(1.0))
        {
          return add(std::get<0>(value), std::get<1>(value), sy);
        }

        bool add(MeasurementAndUncertainty_t value)
        {
          return add(std::get<0>(value), std::get<1>(value), std::get<2>(value));
        }

        template <typename Iter>
        void add_without_uncertainty(Iter begin, Iter end)
        {
          add_without_uncertainty(begin, end, identity());
        }

        template <typename Iter, typename Pred>
        void add_without_uncertainty(Iter begin, Iter end, Pred extractor);

        template <typename Cont, typename Pred>
        void add_without_uncertainty(Cont cont, Pred extractor)
        {
          add_without_uncertainty(std::begin(cont), std::end(cont), extractor);
        }

        template <typename Cont>
        void add_without_uncertainty(Cont cont)
        {
          add_without_uncertainty(std::begin(cont), std::end(cont));
        }

        template <typename VIter, typename UIter, typename VPred, typename UPred = identity>
        unsigned int add_with_uncertainty(VIter begin_value,
                                          VIter end_value,
                                          UIter begin_uncertainty,
                                          VPred value_extractor,
                                          UPred uncertainty_extractor = UPred());

        template <typename Iter>
        unsigned int add_with_uncertainty(Iter begin, Iter end);

        template <typename Cont>
        unsigned int add_with_uncertainty(Cont cont)
        {
          return add_with_uncertainty(std::begin(cont), std::end(cont));
        }


        void clear();


        int N() const { return s2.N(); }

        Data_t AverageUncertainty() const { return WeightToUncertainty(s2.AverageWeight()); }

        template <typename V>
        static constexpr V sqr(V const& v)
        {
          return v * v;
        }

        Data_t XN(unsigned int n) const { return (n == 0) ? s2.Weights() : x.Sum(n); }

        Data_t XNY(unsigned int n) const { return (n == 0) ? y.Weights() : xy.Sum(n); }

        template <unsigned int N>
        Data_t XN() const
        {
          return SumExtractor<decltype(x)::Power, N>::Sum(s2, x);
        }

        template <unsigned int N>
        Data_t XNY() const
        {
          return SumExtractor<decltype(xy)::Power, N>::Sum(y, xy);
        }

        Data_t Y2() const { return y2.template SumN<1>(); }


        template <typename Stream>
        void Print(Stream& out) const;

        static Data_t UncertaintyToWeight(Data_t s) { return Data_t(1.) / sqr(s); }

        static Data_t WeightToUncertainty(Data_t w) { return Data_t(1.) / std::sqrt(w); }

      protected:
        WeightTracker<Data_t> s2;          
        DataTracker<Degree * 2, Data_t> x; 
        WeightTracker<Data_t> y;           
        DataTracker<1, Data_t> y2;         
        DataTracker<Degree, Data_t> xy;    

      }; // class FitDataCollector<>

      template <typename T, unsigned int D>
      inline std::ostream& operator<<(std::ostream& out, FitDataCollector<T, D> const& stats)
      {
        stats.Print(out);
        return out;
      }

      template <typename T, unsigned int D>
      class SimplePolyFitterDataBase {
        using Collector_t = FitDataCollector<T, D>; 

      public:
        static constexpr unsigned int Degree = Collector_t::Degree;

        using Data_t = typename Collector_t::Data_t; 

        using Measurement_t = typename Collector_t::Measurement_t;

        using MeasurementAndUncertainty_t = typename Collector_t::MeasurementAndUncertainty_t;

        // default constructor, destructor and all the rest


        bool add(Data_t x, Data_t y, Data_t sy = Data_t(1.0)) { return stats.add(x, y, sy); }

        bool add(Measurement_t value, Data_t sy = Data_t(1.0)) { return stats.add(value, sy); }

        bool add(MeasurementAndUncertainty_t value) { return stats.add(value); }

        template <typename Iter>
        void add_without_uncertainty(Iter begin, Iter end)
        {
          stats.add_without_uncertainty(begin, end);
        }

        template <typename Iter, typename Pred>
        void add_without_uncertainty(Iter begin, Iter end, Pred extractor)
        {
          stats.add_without_uncertainty(begin, end, extractor);
        }

        template <typename Cont, typename Pred>
        void add_without_uncertainty(Cont cont, Pred extractor)
        {
          stats.add_without_uncertainty(cont, extractor);
        }

        template <typename Cont>
        void add_without_uncertainty(Cont cont)
        {
          stats.add_without_uncertainty(cont);
        }

        template <typename VIter, typename UIter, typename VPred, typename UPred = identity>
        unsigned int add_with_uncertainty(VIter begin_value,
                                          VIter end_value,
                                          UIter begin_uncertainty,
                                          VPred value_extractor,
                                          UPred uncertainty_extractor = UPred())
        {
          return stats.add_with_uncertainty(
            begin_value, end_value, begin_uncertainty, value_extractor, uncertainty_extractor);
        }

        template <typename Iter>
        unsigned int add_with_uncertainty(Iter begin, Iter end)
        {
          return stats.add_with_uncertainty(begin, end);
        }

        template <typename Cont>
        unsigned int add_with_uncertainty(Cont cont)
        {
          return stats.add_with_uncertainty(cont);
        }


        void clear() { stats.clear(); }


        int N() const { return stats.N(); }

        Data_t AverageUncertainty() const { return stats.AverageUncertainty(); }


        template <typename Stream>
        void PrintStats(Stream& out) const
        {
          stats.Print(out);
        }

        template <typename V>
        static constexpr V sqr(V const& v)
        {
          return Collector_t::sqr(v);
        }

      protected:
        Collector_t stats; 

        Data_t XN(unsigned int n) const { return stats.XN(n); }

        Data_t XNY(unsigned int n) const { return stats.XNY(n); }

      }; // class SimplePolyFitterDataBase<>

      template <typename T, unsigned int N>
      class SimpleFitterInterface {
      public:
        static constexpr unsigned int NParams = N;

        using Data_t = T; 

        using MatrixOps = FastMatrixOperations<Data_t, NParams>;

        using FitParameters_t = std::array<Data_t, NParams>;

        using FitMatrix_t = typename MatrixOps::Matrix_t;

        virtual ~SimpleFitterInterface() = default;

        virtual bool isValid() const = 0;


        virtual FitParameters_t FitParameters() const = 0;

        virtual FitParameters_t FitParameterErrors() const = 0;

        virtual FitMatrix_t FitParameterCovariance() const = 0;

        virtual Data_t FitParameter(unsigned int n) const { return FitParameters()[n]; }

        virtual Data_t FitParameterError(unsigned int n) const
        {
          return std::sqrt(FitParameterCovariance()[n * (NParams + 1)]);
        }

        virtual Data_t ChiSquare() const = 0;

        virtual int NDF() const = 0;


        virtual bool FillResults(FitParameters_t& params,
                                 FitMatrix_t& Xmat,
                                 Data_t& det,
                                 FitMatrix_t& Smat) const = 0;

        virtual bool FillResults(FitParameters_t& params,
                                 FitParameters_t& paramerrors,
                                 FitMatrix_t& Xmat,
                                 Data_t& det,
                                 FitMatrix_t& Smat) const = 0;

        virtual bool FillResults(FitParameters_t& params, FitParameters_t& paramerrors) const = 0;

        virtual Data_t Evaluate(Data_t x) const = 0;

        Data_t operator()(Data_t x) const { return Evaluate(x); }

        static constexpr Data_t sqr(Data_t v) { return v * v; }

        static constexpr Data_t cube(Data_t v) { return v * v * v; }

      protected:
        virtual Data_t Determinant(FitMatrix_t const& mat) const
        {
          return MatrixOps::Determinant(mat);
        }

        virtual FitMatrix_t InvertMatrix(FitMatrix_t const& mat, Data_t det) const
        {
          return MatrixOps::InvertSymmetricMatrix(mat, det);
        }

        virtual FitMatrix_t InvertMatrix(FitMatrix_t const& mat) const
        {
          return MatrixOps::InvertSymmetricMatrix(mat);
        }

        virtual FitParameters_t MatrixProduct(FitMatrix_t const& mat,
                                              FitParameters_t const& vec) const
        {
          return MatrixOps::MatrixVectorProduct(mat, vec);
        }

      }; // class SimpleFitterInterface<>

      template <typename T, unsigned int D>
      class SimplePolyFitterBase : public SimpleFitterInterface<T, D + 1>,
                                   public SimplePolyFitterDataBase<T, D> {
        using Interface_t = SimpleFitterInterface<T, D + 1>; 
        using Base_t = SimplePolyFitterDataBase<T, D>;       

      public:
        using Base_t::sqr;

        static constexpr unsigned int Degree = Base_t::Degree;

        static constexpr unsigned int NParams = Interface_t::NParams;

        using Data_t = typename Base_t::Data_t; 

        using FitParameters_t = typename Interface_t::FitParameters_t;

        using FitMatrix_t = typename Interface_t::FitMatrix_t;

        virtual bool isValid() const override;


        virtual FitParameters_t FitParameters() const override;

        virtual FitParameters_t FitParameterErrors() const override;

        virtual FitMatrix_t FitParameterCovariance() const override;

        virtual Data_t FitParameter(unsigned int n) const override;

        virtual Data_t FitParameterError(unsigned int n) const override;

        virtual Data_t ChiSquare() const override;

        virtual int NDF() const override { return Base_t::N() - NParams; }


        bool FillResults(FitParameters_t& params,
                         FitMatrix_t& Xmat,
                         Data_t& det,
                         FitMatrix_t& Smat) const override;

        bool FillResults(FitParameters_t& params,
                         FitParameters_t& paramerrors,
                         FitMatrix_t& Xmat,
                         Data_t& det,
                         FitMatrix_t& Smat) const override;

        bool FillResults(FitParameters_t& params, FitParameters_t& paramerrors) const override;

        virtual Data_t Evaluate(Data_t x) const override;

        static FitParameters_t ExtractParameterErrors(FitMatrix_t const& Smat);

      protected:
        virtual FitMatrix_t MakeMatrixX() const;

        virtual FitParameters_t MakeMatrixY() const;

        virtual FitParameters_t FitParameterErrors(FitMatrix_t const& Smat) const
        {
          return ExtractParameterErrors(Smat);
        }

        virtual FitParameters_t FitParameters(FitMatrix_t const& Xmat) const;

        virtual FitParameters_t FitParameters(FitMatrix_t const& Smat, Data_t /* det */) const;

        virtual Data_t Param(unsigned int n, FitMatrix_t const& Xmat) const;
        virtual Data_t Param(unsigned int n, FitMatrix_t const& Xmat, Data_t detXmat) const;

        // import the declaration from the interface
        using Interface_t::Determinant;
        using Interface_t::InvertMatrix;
        using Interface_t::MatrixProduct;

      }; // class SimplePolyFitterBase<>

    } // namespace details

    template <typename T>
    class LinearFit : public details::SimplePolyFitterBase<T, 1U> {
      using Base_t = details::SimplePolyFitterBase<T, 1U>;

    public:
      using Base_t::Degree;
      using Base_t::NParams;
      using typename Base_t::Data_t;

      using typename Base_t::Measurement_t;

      using typename Base_t::MeasurementAndUncertainty_t;

      using Base_t::sqr;

      using FitParameters_t = std::array<Data_t, NParams>;
      using FitMatrix_t = std::array<Data_t, sqr(NParams)>;

      // default constructor, destructor and all the rest

      Data_t Intercept() const { return this->FitParameter(0); }

      Data_t Slope() const { return this->FitParameter(1); }

      Data_t InterceptError() const { return this->FitParameterError(0); }

      Data_t SlopeError() const { return this->FitParameterError(1); }

      Data_t InterceptSlopeCovariance() const
      {
        return this->FitParameterCovariance()[0 * NParams + 1];
      }

      virtual Data_t ChiSquare() const override;

    protected:
      Data_t I() const { return Base_t::stats.template XN<0>(); }
      Data_t X() const { return Base_t::stats.template XN<1>(); }
      Data_t X2() const { return Base_t::stats.template XN<2>(); }
      Data_t Y() const { return Base_t::stats.template XNY<0>(); }
      Data_t XY() const { return Base_t::stats.template XNY<1>(); }
      Data_t Y2() const { return Base_t::stats.Y2(); }

    }; // class LinearFit<>

    template <typename T>
    class QuadraticFit : public details::SimplePolyFitterBase<T, 2U> {
      using Base_t = details::SimplePolyFitterBase<T, 2U>;

    public:
      using Base_t::Degree;
      using Base_t::NParams;
      using typename Base_t::Data_t;

      using typename Base_t::Measurement_t;

      using typename Base_t::MeasurementAndUncertainty_t;

      using Base_t::sqr;

      using FitParameters_t = std::array<Data_t, NParams>;
      using FitMatrix_t = std::array<Data_t, sqr(NParams)>;

      // default constructor, destructor and all the rest

      virtual Data_t ChiSquare() const override;

    protected:
      Data_t I() const { return Base_t::stats.template XN<0>(); }
      Data_t X() const { return Base_t::stats.template XN<1>(); }
      Data_t X2() const { return Base_t::stats.template XN<2>(); }
      Data_t X3() const { return Base_t::stats.template XN<3>(); }
      Data_t X4() const { return Base_t::stats.template XN<4>(); }
      Data_t Y() const { return Base_t::stats.template XNY<0>(); }
      Data_t XY() const { return Base_t::stats.template XNY<1>(); }
      Data_t X2Y() const { return Base_t::stats.template XNY<2>(); }
      Data_t Y2() const { return Base_t::stats.Y2(); }

    }; // class QuadraticFit<>

    template <typename T>
    class GaussianFit : public details::SimpleFitterInterface<T, 3> {
      using Base_t = details::SimpleFitterInterface<T, 3>; 
      using Fitter_t = QuadraticFit<T>;                    

    public:
      static constexpr unsigned int NParams = Base_t::NParams;

      using Data_t = typename Base_t::Data_t; 

      using Measurement_t = typename Fitter_t::Measurement_t;

      using MeasurementAndUncertainty_t = typename Fitter_t::MeasurementAndUncertainty_t;

      using FitParameters_t = typename Fitter_t::FitParameters_t;
      using FitMatrix_t = typename Fitter_t::FitMatrix_t;

      using Base_t::cube;
      using Base_t::sqr;

      // default constructor, destructor and all the rest


      bool add(Data_t x, Data_t y, Data_t sy = Data_t(1.0));

      bool add(Measurement_t value, Data_t sy = Data_t(1.0))
      {
        return add(std::get<0>(value), std::get<1>(value), sy);
      }

      bool add(MeasurementAndUncertainty_t value)
      {
        return add(std::get<0>(value), std::get<1>(value), std::get<2>(value));
      }

      template <typename Iter>
      void add_without_uncertainty(Iter begin, Iter end)
      {
        add_without_uncertainty(begin, end, identity());
      }

      template <typename Iter, typename Pred>
      void add_without_uncertainty(Iter begin, Iter end, Pred extractor);

      template <typename Cont, typename Pred>
      void add_without_uncertainty(Cont cont, Pred extractor)
      {
        add_without_uncertainty(std::begin(cont), std::end(cont), extractor);
      }

      template <typename Cont>
      void add_without_uncertainty(Cont cont)
      {
        add_without_uncertainty(std::begin(cont), std::end(cont));
      }

      template <typename VIter, typename UIter, typename VPred, typename UPred = identity>
      unsigned int add_with_uncertainty(VIter begin_value,
                                        VIter end_value,
                                        UIter begin_uncertainty,
                                        VPred value_extractor,
                                        UPred uncertainty_extractor = UPred());

      template <typename Iter>
      unsigned int add_with_uncertainty(Iter begin, Iter end);

      template <typename Cont>
      unsigned int add_with_uncertainty(Cont cont)
      {
        return add_with_uncertainty(std::begin(cont), std::end(cont));
      }

      void clear() { fitter.clear(); }

      int N() const { return fitter.N(); }

      template <typename Stream>
      void PrintStats(Stream& out) const
      {
        fitter.PrintStats(out);
      }



      virtual bool isValid() const override { return fitter.isValid(); }


      virtual FitParameters_t FitParameters() const override;

      virtual FitParameters_t FitParameterErrors() const override;

      virtual FitMatrix_t FitParameterCovariance() const override;

      virtual Data_t ChiSquare() const override { return fitter.ChiSquare(); }

      virtual int NDF() const override { return fitter.NDF(); }


      virtual bool FillResults(FitParameters_t& params,
                               FitMatrix_t& Xmat,
                               Data_t& det,
                               FitMatrix_t& Smat) const override;
      virtual bool FillResults(FitParameters_t& params,
                               FitParameters_t& paramerrors,
                               FitMatrix_t& Xmat,
                               Data_t& det,
                               FitMatrix_t& Smat) const override;

      virtual bool FillResults(FitParameters_t& params,
                               FitParameters_t& paramerrors) const override;

      virtual Data_t Evaluate(Data_t x) const override
      {
        return Evaluate(x, FitParameters().data());
      }

      virtual Fitter_t const& Fitter() const { return fitter; }

      static Data_t Evaluate(Data_t x, Data_t const* params);

    protected:
      Fitter_t fitter; 


      struct Value_t : public std::tuple<Data_t, Data_t> {
        using Base_t = std::tuple<Data_t, Data_t>;

        Value_t(Data_t v, Data_t e) : Base_t(v, e) {}
        Value_t(MeasurementAndUncertainty_t meas) : Base_t(std::get<1>(meas), std::get<2>(meas)) {}

        constexpr Data_t value() const { return std::get<0>(*this); }
        constexpr Data_t error() const { return std::get<1>(*this); }

        Data_t& value() { return std::get<0>(*this); }
        Data_t& error() { return std::get<1>(*this); }

      }; // Value_t

      static Data_t EncodeValue(Data_t value) { return std::log(value); }

      static Data_t DecodeValue(Data_t value) { return std::exp(value); }

      static Value_t EncodeValue(Data_t value, Data_t error)
      {
        return {std::log(value), error / std::abs(value)};
      }

      static Value_t EncodeValue(Value_t const& value)
      {
        return EncodeValue(value.value(), value.error());
      }

      static Value_t DecodeValue(Value_t const& value)
      {
        const Data_t v = std::exp(value.value());
        return {v, v * value.error()};
      } // DecodeValue()

      static Measurement_t EncodeValue(Measurement_t const& meas)
      {
        return Measurement_t(std::get<0>(meas), EncodeValue(std::get<1>(meas)));
      }

      static MeasurementAndUncertainty_t EncodeValue(MeasurementAndUncertainty_t const& meas)
      {
        Value_t value = EncodeValue(Value_t(meas));
        return {std::get<0>(meas), value.value(), value.error()};
      } // EncodeValue(MeasurementAndUncertainty_t)

      static MeasurementAndUncertainty_t EncodeValue(Measurement_t const& meas, Data_t error)
      {
        Value_t value = EncodeValue(Value_t(std::get<1>(meas), error));
        return {std::get<0>(meas), value.value(), value.error()};
      } // EncodeValue(Measurement_t, Data_t)

      template <typename VPred, typename UPred = void>
      struct EncodeExtractor {
        EncodeExtractor(VPred& vpred, UPred& upred) : value_extractor(vpred), error_extractor(upred)
        {}

        // extractor takes whatever dereferencing Iter gives and returns
        // Measurement_t or MeasurementAndUncertainty_t,
        // the one the extractor returns
        template <typename Elem>
        auto operator()(Elem elem)
        {
          // use explicit casts to make sure we know what we are doing
          return EncodeValue(static_cast<Measurement_t&&>(value_extractor(elem)),
                             static_cast<Data_t&&>(error_extractor(elem)));
        } // operator()

        VPred& value_extractor; 
        UPred& error_extractor; 
      };                        // struct EncodeExtractor

      template <typename Pred>
      struct EncodeExtractor<Pred, void> {
        EncodeExtractor(Pred& pred) : extractor(pred) {}

        // extractor takes whatever dereferencing Iter gives and returns
        // Measurement_t or MeasurementAndUncertainty_t,
        // the one the extractor returns;
        // as usual with enable_if, I am not sure it makes sense
        template <typename Elem, typename = std::enable_if<std::is_const<Pred>::value, void>>
        auto operator()(Elem elem) const
        {
          return EncodeValue(extractor(elem));
        }

        template <typename Elem, typename = std::enable_if<!std::is_const<Pred>::value, void>>
        auto operator()(Elem elem)
        {
          return EncodeValue(extractor(elem));
        }

        Pred& extractor;
      }; // struct EncodeExtractor<Pred>

      template <typename Pred>
      static EncodeExtractor<Pred> Encoder(Pred& pred)
      {
        return {pred};
      }

      template <typename VPred, typename UPred>
      static EncodeExtractor<VPred, UPred> Encoder(VPred& vpred, UPred& upred)
      {
        return {vpred, upred};
      }


      static FitParameters_t ConvertParameters(FitParameters_t const& qpars);

      static void ConvertParametersAndErrors(FitParameters_t const& qpars,
                                             FitMatrix_t const& qparerrmat,
                                             FitParameters_t& params,
                                             FitParameters_t& paramerrors);

      static void ConvertParametersAndVariances(FitParameters_t const& qpars,
                                                FitMatrix_t const& qparerrmat,
                                                FitParameters_t& params,
                                                FitParameters_t& paramvariances);

      static void ConvertParametersAndErrorMatrix(FitParameters_t const& qpars,
                                                  FitMatrix_t const& qparerrmat,
                                                  FitParameters_t& params,
                                                  FitMatrix_t& Smat);

      static bool isValid(FitParameters_t const& params, FitParameters_t const& qpars);

      static void ThrowNotImplemented [[noreturn]] (std::string method)
      {
        throw std::logic_error("Method " + method + "() not implemented");
      }

    }; // class GaussianFit<>

  } // namespace util
} // namespace lar

//==============================================================================
//=== template implementation
//==============================================================================

//******************************************************************************
//***  FitDataCollector<>
//***

template <typename T, unsigned int D>
bool lar::util::details::FitDataCollector<T, D>::add(Data_t x_value,
                                                     Data_t y_value,
                                                     Data_t sy /* = Data_t(1.0) */)
{
  Data_t w = UncertaintyToWeight(sy);
  if (!std::isnormal(w)) return false;
  // the x section has a 1/s^2 weight; we track that weight separately
  s2.add(w);
  x.add(x_value, w);
  // we treat the y section as if it were a x section with a y/s^2 weight;
  // we track that weight separately
  Data_t yw = y_value * w;
  y.add(yw);
  y2.add(sqr(y_value), w); // used only for chi^2
  xy.add(x_value, yw);

  return true; // we did add the value
} // FitDataCollector<>::add()

template <typename T, unsigned int D>
template <typename Iter, typename Pred>
void lar::util::details::FitDataCollector<T, D>::add_without_uncertainty(Iter begin,
                                                                         Iter end,
                                                                         Pred extractor)
{
  std::for_each(begin, end, [this, extractor](auto item) { this->add(extractor(item)); });
} // FitDataCollector<>::add_without_uncertainty(Iter, Pred)

template <typename T, unsigned int D>
template <typename VIter, typename UIter, typename VPred, typename UPred>
unsigned int lar::util::details::FitDataCollector<T, D>::add_with_uncertainty(
  VIter begin_value,
  VIter end_value,
  UIter begin_uncertainty,
  VPred value_extractor,
  UPred uncertainty_extractor /* = UPred() */
)
{
  unsigned int n = 0;
  while (begin_value != end_value) {
    if (add(value_extractor(*begin_value), uncertainty_extractor(*begin_uncertainty))) ++n;
    ++begin_value;
    ++begin_uncertainty;
  } // while
  return n;
} // FitDataCollector<>::add_with_uncertainty(VIter, VIter, UIter, VPred, UPred)

template <typename T, unsigned int D>
template <typename Iter>
unsigned int lar::util::details::FitDataCollector<T, D>::add_with_uncertainty(Iter begin, Iter end)
{
  unsigned int old_n = N();
  std::for_each(begin, end, [this](auto p) { this->add(p); });
  return N() - old_n;
} // FitDataCollector<>::add_with_uncertainty(Iter, Iter)

template <typename T, unsigned int D>
inline void lar::util::details::FitDataCollector<T, D>::clear()
{
  s2.clear();
  x.clear();
  y.clear();
  y2.clear();
  xy.clear();
} // FitDataCollector<>::clear()

template <typename T, unsigned int D>
template <typename Stream>
void lar::util::details::FitDataCollector<T, D>::Print(Stream& out) const
{

  out << "Sums  1/s^2=" << s2.Weights() << "\n      x/s^2=" << x.template SumN<1>();
  for (unsigned int degree = 2; degree <= x.Power; ++degree)
    out << "\n    x^" << degree << "/s^2=" << x.Sum(degree);
  out << "\n      y/s^2=" << y.Weights() << "\n    y^2/s^2=" << y2.Sum();
  if (xy.Power >= 1) out << "\n     xy/s^2=" << xy.template SumN<1>();
  for (unsigned int degree = 2; degree <= xy.Power; ++degree)
    out << "\n   x^" << degree << "y/s^2=" << xy.Sum(degree);
  out << std::endl;
} // FitDataCollector<>::Print()

//******************************************************************************
//***  SimplePolyFitterBase<>
//***
template <typename T, unsigned int D>
inline bool lar::util::details::SimplePolyFitterBase<T, D>::isValid() const
{
  return (Base_t::N() > (int)Degree) && std::isnormal(Determinant(MakeMatrixX()));
} // SimplePolyFitterBase<>::isValid()

template <typename T, unsigned int D>
inline auto lar::util::details::SimplePolyFitterBase<T, D>::FitParameter(unsigned int n) const
  -> Data_t
{
  return Param(n, MakeMatrixX());
} // SimplePolyFitterBase<>::FitParameter(unsigned int)

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::FitParameterError(unsigned int n) const
  -> Data_t
{
  if (n > Degree) return Data_t(0); // no parameter, no error
  return std::sqrt(FitParameterCovariance()[n * (NParams + 1)]);
} // SimplePolyFitterBase<>::FitParameterError()

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::FitParameters() const -> FitParameters_t
{
  FitMatrix_t Xmat = MakeMatrixX();
  FitParameters_t fit_params;
  for (unsigned int iParam = 0; iParam < NParams; ++iParam)
    fit_params[iParam] = Param(iParam, Xmat);
  return fit_params;
} // SimplePolyFitterBase<>::FitParameters()

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::FitParameterErrors() const -> FitParameters_t
{
  return FitParameterErrors(FitParameterCovariance());
} // SimplePolyFitterBase<>::FitParameterErrors()

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::FitParameterCovariance() const -> FitMatrix_t
{
  FitMatrix_t Xmat = MakeMatrixX();
  Data_t det = Determinant(Xmat);
  if (!std::isnormal(det)) {
    throw std::range_error(
      "SimplePolyFitterBase::FitParameterCovariance(): determinant 0 while fitting");
  }
  return InvertMatrix(Xmat, det);
} // SimplePolyFitterBase<>::FitParameterCovariance()

template <typename T, unsigned int D>
bool lar::util::details::SimplePolyFitterBase<T, D>::FillResults(FitParameters_t& params,
                                                                 FitMatrix_t& Xmat,
                                                                 Data_t& det,
                                                                 FitMatrix_t& Smat) const
{

  Xmat = MakeMatrixX();
  det = Determinant(Xmat);
  if (!std::isnormal(det)) {
    Smat.fill(Data_t(0));
    params.fill(Data_t(0));
    return false;
  }
  Smat = InvertMatrix(Xmat, det);
  params = FitParameters(Smat, det);
  return true;
} // SimplePolyFitterBase<>::FillResults(params, matrices, determinant)

template <typename T, unsigned int D>
bool lar::util::details::SimplePolyFitterBase<T, D>::FillResults(FitParameters_t& params,
                                                                 FitParameters_t& paramerrors,
                                                                 FitMatrix_t& Xmat,
                                                                 Data_t& det,
                                                                 FitMatrix_t& Smat) const
{

  if (!this->FillResults(params, Xmat, det, Smat)) return false;
  paramerrors = ExtractParameterErrors(Smat);
  return true;
} // SimplePolyFitterBase<>::FillResults(params, errors, matrices, determinant)

template <typename T, unsigned int D>
bool lar::util::details::SimplePolyFitterBase<T, D>::FillResults(FitParameters_t& params,
                                                                 FitParameters_t& paramerrors) const
{
  // to compute the parameters, we need all the stuff;
  // we just keep it local and discard it in the end. Such a waste.
  FitMatrix_t Xmat, Smat;
  Data_t det;
  return FillResults(params, paramerrors, Xmat, det, Smat);
} // SimplePolyFitterBase<>::FillResults(params, errors)

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::Evaluate(Data_t x) const -> Data_t
{
  FitParameters_t params = FitParameters();
  unsigned int iParam = NParams - 1; // point to last parameter (highest degree)
  Data_t v = params[iParam];
  while (iParam > 0)
    v = v * x + params[--iParam];
  return v;
} // SimplePolyFitterBase<>::Evaluate()

// --- protected methods follow ---
template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::MakeMatrixX() const -> FitMatrix_t
{
  FitMatrix_t Xmat;
  for (unsigned int i = 0; i < NParams; ++i) {   // row
    for (unsigned int j = i; j < NParams; ++j) { // column
      Xmat[j * NParams + i] = Xmat[i * NParams + j] = Base_t::XN(i + j);
    } // for j
  }   // for i
  return Xmat;
} // SimplePolyFitterBase<>::MakeMatrixX()

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::MakeMatrixY() const -> FitParameters_t
{
  FitParameters_t Ymat;
  for (unsigned int i = 0; i < NParams; ++i)
    Ymat[i] = Base_t::XNY(i);
  return Ymat;
} // SimplePolyFitterBase<>::MakeMatrixY()

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::FitParameters(FitMatrix_t const& Xmat) const
  -> FitParameters_t
{
  FitParameters_t fit_params;
  for (unsigned int iParam = 0; iParam < NParams; ++iParam)
    fit_params[iParam] = Param(iParam, Xmat);
  return fit_params;
} // SimplePolyFitterBase<>::FitParameters(FitMatrix_t)

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::FitParameters(FitMatrix_t const& Smat,
                                                                   Data_t /* det */) const
  -> FitParameters_t
{
  return MatrixProduct(Smat, MakeMatrixY());
} // SimplePolyFitterBase<>::FitParameters(FitMatrix_t)

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::Param(unsigned int n,
                                                           FitMatrix_t const& Xmat) const -> Data_t
{
  if (n > Degree) return Data_t(0); // no such a degree, its coefficient is 0

  Data_t detXmat = Determinant(Xmat);
  if (!std::isnormal(detXmat)) {
    throw std::range_error("SimplePolyFitterBase::Param(): Determinant 0 while fitting");
  }
  return Param(n, Xmat, detXmat);
} // SimplePolyFitterBase<>::Param(unsigned int, FitMatrix_t)

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::ExtractParameterErrors(FitMatrix_t const& Smat)
  -> FitParameters_t
{
  FitParameters_t fit_errors;
  for (unsigned int iParam = 0; iParam <= Degree; ++iParam)
    fit_errors[iParam] = std::sqrt(Smat[iParam * (NParams + 1)]);
  return fit_errors;
} // SimplePolyFitterBase<>::FitParameterErrors(FitMatrix_t)

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::Param(unsigned int n,
                                                           FitMatrix_t const& Xmat,
                                                           Data_t detXmat) const -> Data_t
{
  if (n > Degree) return Data_t(0); // no such a degree, its coefficient is 0
  // XYmat is as Xmat...
  FitMatrix_t XYmat(Xmat);
  // ... except that the N-th column is replaced with { sum x^i y }
  for (unsigned int i = 0; i < NParams; ++i)
    XYmat[i * NParams + n] = Base_t::XNY(i);

  return Determinant(XYmat) / detXmat;
} // SimplePolyFitterBase<>::Param(unsigned int, FitMatrix_t, Data_t)

template <typename T, unsigned int D>
auto lar::util::details::SimplePolyFitterBase<T, D>::ChiSquare() const -> Data_t
{
  // the generic implementation of ChiSquare from sums is complex enough that
  // I freaked out
  throw std::logic_error("SimplePolyFitterBase::ChiSquare() not implemented for generic fit");
} // SimplePolyFitterBase<>::ChiSquare()

//******************************************************************************
//***  LinearFit<>
//***

template <typename T>
auto lar::util::LinearFit<T>::ChiSquare() const -> Data_t
{
  FitParameters_t fit_params = this->FitParameters();
  const Data_t b = fit_params[0];
  const Data_t a = fit_params[1];
  return Y2() + sqr(a) * X2() + sqr(b) * I() + Data_t(2) * (a * b * X2() - a * XY() - b * Y());
} // LinearFit<T>::ChiSquare()

//******************************************************************************
//***  QuadraticFit<>
//***

template <typename T>
auto lar::util::QuadraticFit<T>::ChiSquare() const -> Data_t
{
  FitParameters_t a = this->FitParameters();
  return Y2() - Data_t(2) * (a[0] * Y() + a[1] * XY() + a[2] * X2Y()) + sqr(a[0]) * I() +
         Data_t(2) * a[0] * (a[1] * X() + a[2] * X2()) + sqr(a[1]) * X2() +
         Data_t(2) * a[1] * (a[2] * X3()) + sqr(a[2]) * X4();
} // QuadraticFit<T>::ChiSquare()

//******************************************************************************
//***  GaussianFit<>
//***

//
// data interface
//
template <typename T>
bool lar::util::GaussianFit<T>::add(Data_t x, Data_t y, Data_t sy /* = Data_t(1.0) */)
{
  if (y <= Data_t(0)) return false; // ignore the non-positive values
  Value_t value = EncodeValue(Value_t(y, sy));
  return fitter.add(x, value.value(), value.error());
} // GaussianFit<T>::add(Data_t, Data_t, Data_t)

template <typename T>
template <typename Iter, typename Pred>
void lar::util::GaussianFit<T>::add_without_uncertainty(Iter begin, Iter end, Pred extractor)
{
  return fitter.add_without_uncertainty(begin, end, Encoder(extractor));
} // GaussianFit<>::add_without_uncertainty(Iter, Iter, Pred)

template <typename T>
template <typename VIter, typename UIter, typename VPred, typename UPred>
unsigned int lar::util::GaussianFit<T>::add_with_uncertainty(
  VIter begin_value,
  VIter end_value,
  UIter begin_uncertainty,
  VPred value_extractor,
  UPred uncertainty_extractor /* = UPred() */
)
{
  return add_with_uncertainty(
    begin_value, end_value, begin_uncertainty, Encoder(value_extractor, uncertainty_extractor));
} // GaussianFit<T>::add_with_uncertainty()

template <typename T>
template <typename Iter>
unsigned int lar::util::GaussianFit<T>::add_with_uncertainty(Iter begin, Iter end)
{
  unsigned int old_n = N();
  std::for_each(begin, end, [this](auto p) { this->add(p); });
  return N() - old_n;
} // GaussianFit<T>::add_with_uncertainty()

//
// fitting interface
//
template <typename T>
auto lar::util::GaussianFit<T>::FitParameters() const -> FitParameters_t
{
  return ConvertParameters(fitter.FitParameters());
} // GaussianFit<>::FitParameters()

template <typename T>
auto lar::util::GaussianFit<T>::FitParameterErrors() const -> FitParameters_t
{
  FitParameters_t qpars, qparerrors;
  if (!FillResults(qpars, qparerrors)) {
    throw std::runtime_error("GaussianFit::FitParameterErrors() yielded invalid results");
  }
  return qparerrors;
} // GaussianFit<>::FitParameterErrors()

template <typename T>
auto lar::util::GaussianFit<T>::FitParameterCovariance() const -> FitMatrix_t
{
  // we need to go through the whole chain to get the error matrix
  FitParameters_t params;
  FitMatrix_t Xmat;
  Data_t det;
  FitMatrix_t Smat;
  if (!FillResults(params, Xmat, det, Smat)) {
    throw std::runtime_error("GaussianFit::FitParameterCovariance() yielded invalid results");
  }
  return Smat;
} // SimplePolyFitterBase<>::FitParameterCovariance()

template <typename T>
bool lar::util::GaussianFit<T>::FillResults(FitParameters_t& params,
                                            FitParameters_t& paramerrors) const
{
  FitParameters_t qpars;
  FitMatrix_t qparerrmat;
  FitMatrix_t Xmat; // not used
  Data_t det;       // not used
  if (!fitter.FillResults(qpars, Xmat, det, qparerrmat)) return false;
  ConvertParametersAndErrors(qpars, qparerrmat, params, paramerrors);
  return isValid(params, qpars);
} // GaussianFit<>::FillResults()

template <typename T>
bool lar::util::GaussianFit<T>::FillResults(FitParameters_t& params,
                                            FitMatrix_t& Xmat,
                                            Data_t& det,
                                            FitMatrix_t& Smat) const
{
  FitParameters_t qpars;
  FitMatrix_t qparerrmat;
  if (!fitter.FillResults(qpars, Xmat, det, qparerrmat)) return false;
  ConvertParametersAndErrorMatrix(qpars, qparerrmat, params, Smat);
  return isValid(params, qpars);
} // GaussianFit::FillResults()

template <typename T>
bool lar::util::GaussianFit<T>::FillResults(FitParameters_t& params,
                                            FitParameters_t& paramerrors,
                                            FitMatrix_t& Xmat,
                                            Data_t& det,
                                            FitMatrix_t& Smat) const
{
  if (!FillResults(params, Xmat, det, Smat)) return false;
  paramerrors = fitter.ExtractParameterErrors(Smat);
  return true;
} // GaussianFit::FillResults()

template <typename T>
auto lar::util::GaussianFit<T>::Evaluate(Data_t x, Data_t const* params) -> Data_t
{
  Data_t z = (x - params[1]) / params[2];
  return params[0] * std::exp(-0.5 * sqr(z));
} // GaussianFit<>::Evaluate()

template <typename T>
auto lar::util::GaussianFit<T>::ConvertParameters(FitParameters_t const& qpars) -> FitParameters_t
{
  FitParameters_t params;

  Data_t sigma2 = -0.5 / qpars[2]; // sigma^2 = -1 / (2 a2)
  params[2] = std::sqrt(sigma2);   // sigma

  params[1] = sigma2 * qpars[1]; // mean = sigma2 a1

  params[0] = std::exp(qpars[0] - 0.25 * sqr(qpars[1]) / qpars[2]);

  return params;
} // GaussianFit<>::ConvertParameters()

template <typename T>
void lar::util::GaussianFit<T>::ConvertParametersAndVariances(FitParameters_t const& qpars,
                                                              FitMatrix_t const& qparerrmat,
                                                              FitParameters_t& params,
                                                              FitParameters_t& paramvariances)
{
  params = ConvertParameters(qpars);

  FitParameters_t const& a = qpars;
  Data_t const& A = params[0];
  Data_t const& mu = params[1];
  Data_t const& sigma = params[2];

  // error on sigma
  paramvariances[2] = qparerrmat[3 * 2 + 2] / sqr(cube(sigma));

  // error on mu
  paramvariances[1] =
    sqr(mu * (+qparerrmat[3 * 1 + 1] / sqr(a[1]) - 2. * qparerrmat[3 * 2 + 1] / (a[1] * a[2]) +
              qparerrmat[3 * 2 + 2] / sqr(a[2])));

  // error on A
  paramvariances[0] =
    sqr(A * (+qparerrmat[3 * 0 + 0] + 2. * qparerrmat[3 * 0 + 1] * mu +
             (qparerrmat[3 * 1 + 1] + 2. * qparerrmat[3 * 0 + 2]) * sqr(mu) +
             2. * qparerrmat[3 * 1 + 2] * cube(mu) + qparerrmat[3 * 2 + 2] * sqr(sqr(mu))));

} // GaussianFit<>::ConvertParametersAndVariances()

template <typename T>
void lar::util::GaussianFit<T>::ConvertParametersAndErrors(FitParameters_t const& qpars,
                                                           FitMatrix_t const& qparerrmat,
                                                           FitParameters_t& params,
                                                           FitParameters_t& paramerrors)
{
  ConvertParametersAndVariances(qpars, qparerrmat, params, paramerrors);
  // paramerrors actually stores the square of the error; fix it:
  for (Data_t& paramerror : paramerrors)
    paramerror = std::sqrt(paramerror);
} // GaussianFit<>::ConvertParametersAndErrors()

template <typename T>
void lar::util::GaussianFit<T>::ConvertParametersAndErrorMatrix(FitParameters_t const& qpars,
                                                                FitMatrix_t const& qparerrmat,
                                                                FitParameters_t& params,
                                                                FitMatrix_t& Smat)
{
  FitParameters_t paramvariances;
  ConvertParametersAndVariances(qpars, qparerrmat, params, paramvariances);

  // let's call things with their names
  FitParameters_t const& a = qpars;
  Data_t const& A = params[0];
  Data_t const& mu = params[1];
  Data_t const& sigma = params[2];

  // variance on sigma
  Smat[3 * 2 + 2] = paramvariances[2];

  // variance on mu
  Smat[3 * 1 + 1] = paramvariances[1];

  // variance on A
  Smat[3 * 0 + 0] = paramvariances[0];

  // covariance on sigma and mu
  Smat[3 * 1 + 2] = Smat[3 * 2 + 1] =
    (qparerrmat[3 * 1 + 2] + 2 * mu * qparerrmat[3 * 2 + 2]) / sigma;

  // this is the sum of the derivatives of A vs. all a parameters, each one
  // multiplied by the covariance of that parameter with a2
  const Data_t dA_dak_cov_aka2 =
    A * (qparerrmat[3 * 0 + 2] + qparerrmat[3 * 1 + 2] * mu + qparerrmat[3 * 2 + 2] * sqr(mu));
  // covariance on A and sigma
  Smat[3 * 0 + 2] = Smat[3 * 2 + 0] = dA_dak_cov_aka2 / cube(sigma);

  // this other is the same as dA_dak_cov_aka2, but for a1
  const Data_t dA_dak_cov_aka1 =
    A * (qparerrmat[3 * 0 + 1] + qparerrmat[3 * 1 + 1] * mu + qparerrmat[3 * 2 + 1] * sqr(mu));

  // covariance on A and mu
  Smat[3 * 0 + 1] = Smat[3 * 1 + 0] = mu * (dA_dak_cov_aka1 / a[1] - dA_dak_cov_aka2 / a[2]);

} // GaussianFit<>::ConvertParametersAndErrors()

template <typename T>
bool lar::util::GaussianFit<T>::isValid(FitParameters_t const& params, FitParameters_t const& qpars)
{
  return (qpars[2] < Data_t(0)) && (params[0] >= Data_t(0));
} // GaussianFit<>::isValid(FitParameters_t)

//******************************************************************************

#endif // SIMPLEFITS_H