ian@0: //  Copyright 2011, Ian Hobson.
ian@0: //
ian@0: //  This file is part of gpsynth.
ian@0: //
ian@0: //  gpsynth is free software: you can redistribute it and/or modify
ian@0: //  it under the terms of the GNU General Public License as published by
ian@0: //  the Free Software Foundation, either version 3 of the License, or
ian@0: //  (at your option) any later version.
ian@0: //
ian@0: //  gpsynth is distributed in the hope that it will be useful,
ian@0: //  but WITHOUT ANY WARRANTY; without even the implied warranty of
ian@0: //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
ian@0: //  GNU General Public License for more details.
ian@0: //
ian@0: //  You should have received a copy of the GNU General Public License
ian@0: //  along with gpsynth in the file COPYING. 
ian@0: //  If not, see http://www.gnu.org/licenses/.
ian@0: 
ian@0: // Some useful stats functions
ian@0: 
ian@0: #pragma once
ian@0: 
ian@0: #include "range.hpp"
ian@0: #include "std_ex.hpp"
ian@0: 
ian@0: #include <algorithm>
ian@0: #include <iterator>
ian@0: #include <numeric>
ian@0: #include <sstream>
ian@0: #include <stdexcept>
ian@0: 
ian@0: namespace stats {
ian@0: 
ian@0: // Sum
ian@0: template<typename Iterator>
ian@0: typename std::iterator_traits<Iterator>::value_type Sum(Iterator start,
ian@0:                                                         Iterator end) {
ian@0:   typename std::iterator_traits<Iterator>::value_type initializer(0);
ian@0:   return std::accumulate(start, end, initializer);
ian@0: }
ian@0: 
ian@0: // Sum - container
ian@0: template<typename Container>
ian@0: typename Container::value_type Sum(const Container& container) {
ian@0:   return Sum(container.begin(), container.end());
ian@0: }
ian@0: 
ian@0: // Mean
ian@0: template<typename Iterator>
ian@0: typename std::iterator_traits<Iterator>::value_type Mean(Iterator start,
ian@0:                                                          Iterator end) {
ian@0:   return Sum(start, end) / std::distance(start, end);
ian@0: }
ian@0: 
ian@0: // Mean - container
ian@0: template<typename Container>
ian@0: typename Container::value_type Mean(const Container& container) {
ian@0:   return Mean(container.begin(), container.end());
ian@0: }
ian@0: 
ian@0: // measures the MSE between a range and a target
ian@0: template<typename Iterator1, typename Iterator2>
ian@0: typename std::iterator_traits<Iterator1>::value_type 
ian@0: MeanSquaredError(Iterator1 start1, Iterator1 end1, Iterator2 start2) {
ian@0:   typedef typename std::iterator_traits<Iterator1>::value_type Value;
ian@0:   Value n = std::distance(start1, end1);
ian@0:   Value error(0);
ian@0:   while (start1 != end1) {
ian@0:     error += std::pow(*start1 - *start2, Value(2));
ian@0:     ++start1;
ian@0:     ++start2;
ian@0:   }
ian@0:   return error / n;
ian@0: }
ian@0: 
ian@0: // MeanSquaredError - container adaptor
ian@0: template<typename Container>
ian@0: typename Container::value_type MeanSquaredError(const Container& x,
ian@0:                                                 const Container& y) {
ian@0:   return MeanSquaredError(x.begin(), x.end(), y.begin());
ian@0: }
ian@0: 
ian@0: // root mean square error
ian@0: template<typename Iterator1, typename Iterator2>
ian@0: double RMSE(Iterator1 start, Iterator1 end, Iterator2 target) {
ian@0:   return std::sqrt(MeanSquaredError(start, end, target));
ian@0: }
ian@0: 
ian@0: // normalized root mean square error, using precomputed minima and maxima
ian@0: template<typename Iterator1, typename Iterator2, typename T>
ian@0: T NRMSE(Iterator1 start, Iterator1 end, Iterator2 target,
ian@0:         const stdx::Range<T>& range1, const stdx::Range<T>& range2) {
ian@0:   T range = std::max(range1.Maximum(), range2.Maximum())
ian@0:             - std::min(range1.Minimum(), range2.Minimum());
ian@0:   return RMSE(start, end, target) / range;
ian@0: }
ian@0: 
ian@0: // Takes a set of numbers with sum <= 1 which define a distribution.
ian@0: // operator() returns index chosen randomly with distribution defined by the
ian@0: // provided probabilites.
ian@0: // if the probabilites have sum < 1 then the difference is taken to imply a 
ian@0: // single additional index
ian@0: // e.g. given:
ian@0: // ProbabilitySelector selector(boost::assign::list_of(0.1)(0.8));
ian@0: // selector() will yield 
ian@0: // '0' 10%,
ian@0: // '1' 80%,
ian@0: // '2' 10%
ian@0: class ProbabilitySelector {
ian@0:   std::vector<double> boundaries_;
ian@0:   
ian@0: public:
ian@0:   ProbabilitySelector() {}
ian@0:   ProbabilitySelector(const std::vector<double>& probabilities)
ian@0:   {
ian@0:     SetProbabilities(probabilities);
ian@0:   }
ian@0:   
ian@0:   void SetProbabilities(const std::vector<double>& probabilities) {
ian@0:     if (std::count_if(probabilities.begin(), probabilities.end(),
ian@0:                       stdx::LessThan<double>(0.0))) {
ian@0:       throw std::runtime_error("ProbabilitySelector: negative value found");
ian@0:     }
ian@0:     // assign range boundaries
ian@0:     std::partial_sum(probabilities.begin(), probabilities.end(), 
ian@0:                      std::back_inserter(boundaries_));
ian@0:   }
ian@0:   
ian@0:   std::size_t operator()() const {
ian@0:     double random = rand() / static_cast<double>(RAND_MAX);
ian@0:     std::vector<double>::const_iterator index;
ian@0:     index = std::upper_bound(boundaries_.begin(), boundaries_.end(), random);
ian@0:     return std::distance(boundaries_.begin(), index);
ian@0:   }
ian@0:   
ian@0:   bool Initialized() const { return !boundaries_.empty(); }
ian@0: };
ian@0: 
ian@0: } // stats namespace