yule_walker_solver.hh

#ifndef VTESTBED_CORE_YULE_WALKER_SOLVER_HH
#define VTESTBED_CORE_YULE_WALKER_SOLVER_HH

#include <cmath>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/linalg/types.hh>

namespace vtb {

    namespace core {

        template <class T>
        blitz::TinyVector<int,3>
        chop_right(const blitz::Array<T,3>& rhs, T eps);
        template <class T>
        blitz::TinyVector<int,4>
        chop_right(const blitz::Array<T,4>& rhs, T eps);

        template <class T, int N>
        inline T
        autoregression_white_noise_variance(
            const blitz::Array<T,N>& acf,
            const blitz::Array<T,N>& phi
        ) {
            using blitz::RectDomain;
            using blitz::sum;
            RectDomain<N> rect{{0}, phi.shape() - 1};
            const T& variance = *acf.data();
            return variance - sum(phi*acf(rect));
        }

        template <class T, int N>
        T
        moving_average_white_noise_variance(
            const blitz::Array<T,N>& acf,
            const blitz::Array<T,N>& theta
        ) {
            using blitz::sum;
            using blitz::pow2;
            const T& variance = *acf.data();
            return variance / sum(pow2(theta));
        }

        template <class T, int N>
        inline T
        acf_variance(const blitz::Array<T,N>& acf) {
            return *acf.data();
        }

        template <class T, int N>
        class Yule_walker_solver {

        public:
            typedef T value_type;
            typedef blitz::Array<T,N> array_type;
            typedef std::unique_ptr<Yule_walker_solver<T,N>> pointer;

        private:
            T _varwn = T(0);

        public:

            Yule_walker_solver() = default;

            Yule_walker_solver(const Yule_walker_solver&) = default;

            Yule_walker_solver&
            operator=(const Yule_walker_solver&) = default;

            virtual ~Yule_walker_solver() = default;

            virtual array_type
            solve(array_type acf) = 0;

            virtual pointer
            copy() const = 0;

            inline array_type
            operator()() {
                return this->solve();
            }

            inline value_type
            white_noise_variance() const noexcept {
                return this->_varwn;
            }

        protected:

            inline void
            white_noise_variance(value_type rhs) {
                this->_varwn = rhs;
            }

        };

        template <class T, int N>
        class Choi_yule_walker_solver: public Yule_walker_solver<T,N> {

        private:
            typedef Yule_walker_solver<T,N> base_type;

        public:
            using typename base_type::value_type;
            using typename base_type::array_type;
            using typename base_type::pointer;

        private:
            T _varepsilon = T(1e-5);
            bool _determinetheorder = true;
            T _chopepsilon = T(1e-10);
            bool _chop = true;
            int _maxorder = 0;
            T _minvariance{0};

        public:

            Choi_yule_walker_solver() = default;

            ~Choi_yule_walker_solver() = default;

            Choi_yule_walker_solver(const Choi_yule_walker_solver&) = default;

            Choi_yule_walker_solver&
            operator=(const Choi_yule_walker_solver&) = default;

            array_type
            solve(array_type acf) override;

            inline pointer
            copy() const override {
                return pointer(new Choi_yule_walker_solver<T,N>(*this));
            }

            inline bool
            determine_the_order() const noexcept {
                return this->_determinetheorder;
            }

            inline void
            determine_the_order(bool rhs) noexcept {
                this->_determinetheorder = rhs;
            }

            inline void
            do_not_determine_the_order() noexcept {
                this->_determinetheorder = false;
            }

            inline void
            chop(bool rhs) noexcept {
                this->_chop = rhs;
            }

            inline bool
            chop() const noexcept {
                return this->_chop;
            }

            inline void
            do_not_chop() noexcept {
                this->_chop = false;
            }

            inline value_type
            chop_epsilon() const noexcept {
                return this->_chopepsilon;
            }

            inline void
            chop_epsilon(value_type rhs) noexcept {
                this->_chopepsilon = rhs;
            }

            inline value_type
            var_epsilon() const noexcept {
                return this->_varepsilon;
            }

            inline void
            var_epsilon(value_type rhs) noexcept {
                this->_varepsilon = rhs;
            }

            inline void
            max_order(int rhs) {
                this->_maxorder = rhs;
            }

            inline int
            max_order() const noexcept {
                return this->_maxorder;
            }

            inline void
            min_variance(T rhs) noexcept {
                this->_minvariance = rhs;
            }

            inline T
            min_variance() const noexcept {
                return this->_minvariance;
            }

        private:

            inline bool
            variance_has_not_changed_much(T var0, T var1) noexcept {
                return (this->_determinetheorder &&
                       std::abs(var1-var0) < this->_varepsilon) ||
                       (var1 > 0 && !(var1 > this->_minvariance));
            }

        };

        template <class T, int N>
        class Autocovariance_matrix_generator;

        template <class T>
        class Autocovariance_matrix_generator<T,3> {
        private:
            typedef blitz::Array<T,3> array_type;
            typedef blitz::TinyVector<int,3> shape_type;
            typedef vtb::linalg::Square_matrix<T> matrix_type;

        private:
            const array_type& _acf;
            const shape_type& _arorder;

        public:
            Autocovariance_matrix_generator(
                const array_type& acf,
                const shape_type& ar_order
            ): _acf(acf), _arorder(ar_order) {}

            matrix_type
            AC_matrix_block(int i0, int j0);

            matrix_type
            AC_matrix_block(int i0);

            matrix_type
            operator()();
        };

        template <class T>
        class Autocovariance_matrix_generator<T,4> {

        private:
            typedef blitz::Array<T,4> array_type;
            typedef blitz::TinyVector<int,4> shape_type;
            typedef vtb::linalg::Square_matrix<T> matrix_type;

        private:
            const array_type& _acf;
            const shape_type& _arorder;

        public:
            Autocovariance_matrix_generator(
                const array_type& acf,
                const shape_type& ar_order
            ): _acf(acf), _arorder(ar_order) {}

            matrix_type
            AC_matrix_block(int i0, int j0, int k0);

            matrix_type
            AC_matrix_block(int i0, int j0);

            matrix_type
            AC_matrix_block(int i0);

            matrix_type
            operator()();
        };

        template <class T, int N>
        class Gauss_yule_walker_solver: public Yule_walker_solver<T, N> {

        private:
            typedef Yule_walker_solver<T,N> base_type;
            typedef blitz::TinyVector<int,N> shape_type;

        public:
            using typename base_type::value_type;
            using typename base_type::array_type;
            using typename base_type::pointer;

        private:
            shape_type _order{0};
            T _chopepsilon = T(1e-10);
            bool _chop = true;

        public:

            Gauss_yule_walker_solver() = default;

            Gauss_yule_walker_solver(const Gauss_yule_walker_solver&) = default;

            Gauss_yule_walker_solver&
            operator=(const Gauss_yule_walker_solver&) = default;

            ~Gauss_yule_walker_solver() = default;

            array_type
            solve(array_type acf) override;

            inline pointer
            copy() const override {
                return pointer(new Gauss_yule_walker_solver<T,N>(*this));
            }

            inline void order(const shape_type& rhs) { this->_order = rhs; }
            inline const shape_type& order() const { return this->_order; }
            inline void chop(bool rhs) { this->_chop = rhs; }
            inline value_type chop_epsilon() const { return this->_chopepsilon; }
            inline void chop_epsilon(value_type rhs) { this->_chopepsilon = rhs; }

        };

        template <class T, int N>
        class Fixed_point_iteration_yule_walker_solver:
            public Yule_walker_solver<T,N> {

        private:
            typedef Yule_walker_solver<T,N> base_type;

        public:
            using typename base_type::value_type;
            using typename base_type::array_type;
            using typename base_type::pointer;
            typedef blitz::TinyVector<int,N> shape_type;

        private:
            shape_type _order;
            int _maxiterations = 1000;
            T _maxresidual = T(1e-5);
            T _minvarwn = T(1e-6);
            T _epsvarwn = T(1e-5);

        public:

            array_type
            solve(array_type acf) override;

            inline pointer
            copy() const override {
                return pointer(
                    new Fixed_point_iteration_yule_walker_solver<T,N>(*this)
                );
            }

            inline void
            order(const shape_type& rhs) {
                this->_order = rhs;
            }

            inline const shape_type&
            order() const noexcept {
                return this->_order;
            }

            inline int
            max_iterations() const noexcept {
                return this->_maxiterations;
            }

            inline void
            max_iterations(int rhs) noexcept {
                this->_maxiterations = rhs;
            }

            inline T
            max_residual() const noexcept {
                return this->_maxresidual;
            }

            inline void
            max_residual(T rhs) noexcept {
                this->_maxresidual = rhs;
            }

            inline void
            min_white_noise_variance(T rhs) noexcept {
                this->_minvarwn = rhs;
            }

            inline T
            min_white_noise_variance() const noexcept {
                return this->_minvarwn;
            }

            inline void
            min_white_noise_variance_delta(T rhs) noexcept {
                this->_epsvarwn = rhs;
            }

            inline T
            min_white_noise_variance_delta() const noexcept {
                return this->_epsvarwn;
            }

        };

    }

}

#endif // vim:filetype=cpp