ShortCorr.h

Go to the documentation of this file.
00001 #ifndef SHORTCORR_H_
00002 #define SHORTCORR_H_
00003 #include "FatalException.h"
00004 #include "statutils.h"
00005 #include <algorithm>
00006 #include <numeric>
00007 
00008 namespace gplib
00009   {
00010 
00011     /** \addtogroup tstools Time series analysis methods */
00012     /* @{ */
00013 
00014     //! Calculate the correlation between a short time series and a master time series
00015     /*! The short time series will be shifted along the master and the corresponding correlation will be stored in
00016      * the structure referenced by outpebegin, make sure that structure can hold enough data*/
00017     template<typename _InputIterator, typename _OutputIterator>
00018     void ShortCorr(_InputIterator masterbegin, _InputIterator masterend,
00019         _InputIterator shortbegin, _InputIterator shortend,
00020         _OutputIterator outbegin)
00021       {
00022         //check the length of both time series, the "short" time series has to be shorter than the master
00023         if (distance(shortbegin, shortend) > distance(masterbegin, masterend))
00024           throw FatalException("Short time series longer than master !");
00025 
00026         typedef std::vector<
00027             typename std::iterator_traits<_InputIterator>::value_type>
00028             tinvector;
00029         _InputIterator currstart = masterbegin;
00030         _OutputIterator currout = outbegin;
00031         const int shortsize = distance(shortbegin, shortend);
00032 
00033         // create a local copy of the short time series, so we can do some manipulation
00034         tinvector shortts(shortbegin, shortend);
00035         SubMean(shortts.begin(), shortts.end());
00036         //calculate the zero lag autocorrelation of the short time series
00037         tinvector shortauto(shortts.begin(), shortts.end());
00038         std::transform(shortauto.begin(), shortauto.end(), shortauto.begin(),
00039             shortauto.begin(), std::multiplies<double>());
00040         const double shortpower = std::accumulate(shortauto.begin(),
00041             shortauto.end(),
00042             typename std::iterator_traits<_OutputIterator>::value_type());
00043         // go through the time series and calculate the correlation
00044         while (currstart + shortsize < masterend)
00045           {
00046             tinvector currmul(currstart, currstart + shortsize);
00047 
00048             SubMean(currmul.begin(), currmul.end());
00049             tinvector currauto(currmul.begin(), currmul.end());
00050             std::transform(currauto.begin(), currauto.end(), currauto.begin(),
00051                 currauto.begin(), std::multiplies<double>());
00052             double currpower = std::accumulate(currauto.begin(),
00053                 currauto.end(),
00054                 typename std::iterator_traits<_OutputIterator>::value_type());
00055 
00056             std::transform(currmul.begin(), currmul.end(), shortts.begin(),
00057                 currmul.begin(), std::multiplies<double>());
00058             *currout = std::accumulate(currmul.begin(), currmul.end(),
00059                 typename std::iterator_traits<_OutputIterator>::value_type());
00060             *currout /= sqrt(shortpower * currpower);
00061             ++currout;
00062             ++currstart;
00063           }
00064       }
00065   /* @} */
00066   }
00067 #endif /*SHORTCORR_H_*/