mathops.hpp

Go to the documentation of this file.

/*
 * Copyright 2008,2009,2010,2011,2012 Stephan Ossowski, Korbinian Schneeberger,
 * Felix Ott, Joerg Hagmann, Alf Scotland, Sebastian Bender
 * 
 * This file is part of SHORE.
 * 
 * SHORE 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 3 of the License, or
 * (at your option) any later version.
 * 
 * SHORE 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 SHORE.  If not, see <http://www.gnu.org/licenses/>.
 */


#ifndef SHORE_MISC_MATHOPS_HPP__
#define SHORE_MISC_MATHOPS_HPP__

#include <deque>
#include <iterator>
#include <numeric>
#include <string>

#include <boost/random.hpp>

namespace shore
{

template<class U,class T>
U round(T floatval)
{
        return static_cast<U>(floatval+T(0.5));
}

double log2(const double val);

template<class Iter>
double mean(const Iter begin,const Iter end)
{
        if(begin==end) return 0.0;
        const double sum=std::accumulate(begin,end,0.0);
        const double size=end-begin;
        return sum/size;
}

class incremental_mean_var
{
 private:

        size_t m_size;
        double m_mean;
        double m_rec;

 public:

        incremental_mean_var();

        void append(const double d);

        size_t size() const;
        double mean() const;
        double var() const;
        double sd() const;
};

template<class Iter>
double pairs_mean(Iter begin,const Iter end)
{
        if(begin==end) return 0.0;
        double sum=0.0;
        size_t nelem=0;

        while(begin!=end)
        {
                sum+=(static_cast<double>(begin->first)*static_cast<double>(begin->second));
                nelem+=static_cast<size_t>(begin->second);
                ++begin;
        }

        return sum/static_cast<double>(nelem);
}

template<class Iter>
double medianize(const Iter begin,const Iter end)
{
        double ret=0.0;

        if(begin==end)
                return ret;

        const size_t size=end-begin;

        const Iter median_i=begin+(size>>1);

        std::nth_element(begin,median_i,end);

        if((size&1)>0)
        {
                ret=(*median_i);
        }
        else
        {
                const Iter median_i_lower=median_i-1;

                std::iter_swap(median_i_lower,std::max_element(begin,median_i));

                const double med_low= *median_i_lower;
                const double med_high= *median_i;
                ret=0.5*(med_low+med_high);
        }

        return ret;
}

template<class Iter>
double geometric_mean(const Iter begin,const Iter end,const bool skipZ=false)
{
        double logsmean=0.0;
        size_t s=0;

        for(Iter i=begin;i!=end;++i)
        {
                if((*i)==0.0)
                {
                        if(skipZ) continue;
                        return 0.0;
                }

                logsmean+=log2(*i);
                ++s;
        }

        if(s==0) return 0.0;

        logsmean/=double(s);

        return std::pow(2.0,logsmean);
}

template<class Iter>
double harmonic_mean(const Iter begin,const Iter end)
{
        const double amean1=mean(begin,end);
        const double gmean1=geometric_mean(begin,end);
        const double gmean2=gmean1*gmean1;

        return gmean2/amean1;
}

template<class Iter,class Oter>
void skewness_kurtosis(Iter begin,Iter end,Oter dest,const bool var_n)
{
        // use precision of the destination value type
        const double m=mean(begin,end);

        double var=0.0;
        double var3=0.0;
        double var4=0.0;

        const size_t size=end-begin;

        while(begin!=end)
        {
                const double diff=(*begin)-m;
                const double diff2=diff*diff;
                var+=diff2;
                var3+=diff*diff2;
                var4+=diff2*diff2;
                ++begin;
        }

        var/=static_cast<double>(var_n?(size):(size-1));

        (*dest++)=var3/(var*sqrt(var)*static_cast<double>(size));
        (*dest)=var4/(var*var*static_cast<double>(size))-3.0;
}

template<class Iter>
double var(const Iter begin,const Iter end,
  const double mu,const bool var_n=false)
{
        double var=0.0;

        Iter s=begin;

        while(s!=end)
        {
                const double diff=(*s)-mu;
                const double diff2=diff*diff;
                var+=diff2;
                ++s;
        }

        const size_t size=end-begin;

        if(var_n)
        {
                var/=static_cast<double>(size);
        }
        else
        {
                var=(size>1)?(var/static_cast<double>(size-1)):0.0;
        }

        return var;
}

template<class Iter>
double pairs_var(const Iter begin,const Iter end,
  const double mu,const bool var_n=false)
{
        double var=0.0;
        size_t size=0;

        Iter s=begin;

        while(s!=end)
        {
                const double diff=s->first-mu;
                const double diff2=diff*diff;
                var+=(diff2*static_cast<double>(s->second));
                size+=static_cast<size_t>(s->second);
                ++s;
        }

        if(var_n)
        {
                var/=static_cast<double>(size);
        }
        else
        {
                var=(size>1)?(var/(static_cast<double>(size-1))):0.0;
        }

        return var;
}

template<class Iter>
double var(const Iter begin,const Iter end,const bool var_n=false)
{
        const double mu=mean(begin,end);
        return var(begin,end,mu,var_n);
}

template<class Iter>
double stddev(const Iter begin,const Iter end,const bool var_n=false)
{
        return std::sqrt(var(begin,end,var_n));
}

template<class Iter>
double stddev(const Iter begin,const Iter end,
  const double mu,const bool var_n=false)
{
        return std::sqrt(var(begin,end,mu,var_n));
}

template<class Iter>
double pairs_stddev(const Iter begin,const Iter end,
  const double mu,const bool var_n=false)
{
        return std::sqrt(pairs_var(begin,end,mu,var_n));
}

template<class Xter,class Yter>
std::pair<double,double> polycog(
  Xter x_beg,Xter x_end,
  Yter y_beg,double*const ret_area2)
{
        std::pair<double,double> ret(0.0,0.0);
        double area2=0.0;
        Xter xi=x_beg;
        Yter yi=y_beg;

        while(xi!=x_end)
        {
                const double x=double(*xi);
                const double y=double(*yi);

                ++xi;
                ++yi;

                double xnext;
                double ynext;
                if(xi==x_end)
                {
                        xnext=double(*x_beg);
                        ynext=double(*y_beg);
                }
                else
                {
                        xnext=double(*xi);
                        ynext=double(*yi);
                }
                const double xynext=x*ynext;
                const double xnexty=xnext*y;

                area2+=(xynext-xnexty);
                ret.first+=((x+xnext)*(xynext-xnexty));
                ret.second+=((y+ynext)*(xynext-xnexty));
        }

        if(area2!=0.0)
        {
                const double sc=(1.0/(area2+area2+area2));
                ret.first*=sc;
                ret.second*=sc;
        }
        else
        {
                // to-do: correct this? (currently assumes it's a straight line)
                ret.first=0.5*((*(--xi))+(*x_beg));
                ret.second=0.5*((*(--yi))+(*y_beg));
        }

        if(ret_area2!=0) (*ret_area2)=area2;

        return ret;
}

template<class Xter,class Yter>
std::pair<double,double> polycog(
  Xter x_beg,Xter x_end,Yter y_beg)
{
        return polycog(x_beg,x_end,y_beg,0);
}

template<class Iter>
void mean_filter(Iter beg,Iter end,const size_t width)
{
        if(width==1) return; // nothing to do

        if((width&1)==0)
                throw std::runtime_error("mean filter width must be odd");

        Iter w_beg=beg;
        Iter w_cen=beg;
        Iter w_end=beg;

        std::deque<double> w_mod;

        double w_div=width;

        double w_total=0.0;
        for(size_t i=0;i<width;++i)
        {
                if(w_end==end)
                {
                        w_div=i;
                        break;
                }
                w_total+=(*w_end);
                ++w_end;
        }

        // set first half of the first window
        double w_mean=w_total/w_div;
        for(size_t i=0;i<=(width>>1);++i)
        {
                if(w_cen==end)
                {
                        break;
                }
                w_mod.push_back(*w_cen);
                (*w_cen)=w_mean;
                ++w_cen;
        }

        while(w_end!=end)
        {
                w_total+=(double(*w_end)-w_mod.front());

                w_mod.pop_front();
                w_mod.push_back(*w_cen);

                w_mean=w_total/w_div;
                (*w_cen)=w_mean;

                ++w_beg;
                ++w_cen;
                ++w_end;
        }

        while(w_cen!=w_end)
        {
                (*w_cen)=w_mean;
                ++w_cen;
        }
}

template<class T,class Iter>
void seq(T lb,const T ub,const size_t len,Iter dest)
{
        if(len<2)
        {
                if(len==1) (*dest)=lb;
                return;
        }

        const T step=(ub-lb)/(len-1);

        for(size_t i=0;i<len;++i)
        {
                (*dest++)=lb;
                lb+=step;
        }
}

template<class IIter,class OIter,class Rng>
void random_copy(IIter b,IIter e,OIter o,Rng& rng)
{
        if(e<=b) return;

        const size_t size=e-b;

        std::vector<size_t> vi(size);
        seq(size_t(0),size-1,size,vi.begin());

        while(!vi.empty())
        {
                boost::uniform_int<size_t> range(0,vi.size()-1);

                boost::variate_generator<Rng&,boost::uniform_int<size_t> >
                        die(rng,range);

                const size_t id=die();

                (*o)=(*(b+vi[id]));
                vi[id]=vi.back();
                vi.pop_back();
                ++o;
        }
}

template<typename T>
struct index_generator
{
        T m_value;

        index_generator()
        :m_value(T())
        {}

        T operator()()
        {
                return m_value++;
        }
};

template<typename IndexIter>
void index_fill(IndexIter ib,IndexIter ie)
{
        typedef typename std::iterator_traits<IndexIter>::value_type value_type;

        index_generator<value_type> gen;

        std::generate(ib,ie,gen);
}

template<typename DataIter,typename Comparator=void>
struct index_comparator
{
        const DataIter m_data;
        Comparator m_cmp;

        index_comparator(const DataIter data,Comparator cmp)
        :m_data(data),m_cmp(cmp)
        {}

        bool operator()(const size_t i,const size_t j)
        {
                return m_cmp(m_data[i],m_data[j]);
        }
};

template<typename DataIter>
struct index_comparator<DataIter,void>
{
        const DataIter m_data;

        index_comparator(const DataIter data)
        :m_data(data)
        {}

        bool operator()(const size_t i,const size_t j) const
        {
                return m_data[i]<m_data[j];
        }
};

template<typename DataIter,typename IndexIter>
void index_sort(IndexIter ib,IndexIter ie,DataIter db)
{
        const index_comparator<DataIter> cmp(db);

        std::sort(ib,ie,cmp);
}

template<typename DataIter,typename IndexIter,typename Comparator>
void index_sort(IndexIter ib,IndexIter ie,DataIter db,Comparator data_cmp)
{
        const index_comparator<DataIter,Comparator> cmp(db,data_cmp);

        std::sort(ib,ie,cmp);
}

template<typename DataIter,typename OrderingIter,typename RankIter>
void ranks(RankIter first,RankIter last,OrderingIter ob,DataIter db)
{
        const size_t size=last-first;
        if(size>1)
        {
                first[ob[size-1]]=size-1;
                for(size_t i=size-2;i>0;--i)
                {
                        if(db[ob[i]]<db[ob[i+1]]) first[ob[i]]=i;
                        else first[ob[i]]=first[ob[i+1]];
                }
                if(db[ob[0]]<db[ob[1]]) first[ob[0]]=0;
                else first[ob[0]]=first[ob[1]];
        }
}

template<class Abs_t>
struct relabs
{
        Abs_t abs;
        double rel;

        relabs():abs(0),rel(0.0) {}

        double calc(const Abs_t& refval) const
        {
                if(rel!=0.0) return rel*double(refval);
                return abs;
        }

        Abs_t round(const Abs_t& refval) const
        {
                return long(calc(refval)+0.5);
        }

        Abs_t floor(const Abs_t& refval) const
        {
                return std::floor(calc(refval));
        }

        Abs_t ceil(const Abs_t& refval) const
        {
                return std::ceil(calc(refval));
        }

        bool is_absolute() const
        {
                return rel==0.0;
        }

        static relabs<Abs_t> mkabs(const Abs_t v)
        {
                relabs<Abs_t> ret;
                ret.abs=v;
                return ret;
        }

        static relabs<Abs_t> mkrel(const double v)
        {
                relabs<Abs_t> ret;
                ret.rel=v;
                return ret;
        }
};

template<class Abs_t>
std::istream& operator>>(std::istream& is,relabs<Abs_t>& dest)
{
        dest=relabs<Abs_t>();

        double d;
        is>>d;

        if((!is.eof())&&(is.peek()=='%'))
        {
                is.ignore(1);
                dest.rel=0.01*d;
        }
        else dest.abs=d;

        return is;
}

template<class Abs_t>
std::ostream& operator<<(std::ostream& os,const relabs<Abs_t>& val)
{
        if(val.rel!=0.0) os<<(100.0*val.rel)<<'%';
        else os<<val.abs;
        return os;
}

} // namespace shore

#endif // SHORE_MISC_MATHOPS_HPP__