arma.hh

#ifndef VTESTBED_CORE_ARMA_HH
#define VTESTBED_CORE_ARMA_HH

#include <complex>
#include <random>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/core/fourier_transform.hh>
#include <vtestbed/core/wavy_surface_generator.hh>
#include <vtestbed/core/yule_walker_solver.hh>

namespace vtb {

    namespace core {

        template <class T>
        inline std::unique_ptr<T>
        copy(const std::unique_ptr<T>& ptr) {
            return std::unique_ptr<T>(ptr ? ptr->copy() : nullptr);
        }

        class Invalid_ARMA_process: public std::exception {

        private:
            int _nbadroots = 0;

        public:

            inline explicit
            Invalid_ARMA_process(int nbadroots):
            _nbadroots(nbadroots) {}

            inline int
            num_bad_roots() const noexcept {
                return this->_nbadroots;
            }

            const char*
            what() const noexcept override {
                return "AR/MA process is not stationary/invertible: "
                    "some roots lie inside unit disk.";
            }

        };

        template <class T, int N>
        void
        validate_process(blitz::Array<T,N> phi);

        template <class T, int N>
        blitz::Array<T,N>
        auto_covariance(const blitz::Array<T,N>& rhs);

        template <class T, int N>
        blitz::Array<T,N>
        covariance(
            blitz::Array<T,N> lhs,
            blitz::Array<T,N> rhs,
            Chirp_Z_transform<std::complex<T>,N>& fft
        );

        template <class T, int N>
        inline blitz::Array<T,N>
        covariance(blitz::Array<T,N> lhs, blitz::Array<T,N> rhs) {
            Chirp_Z_transform<std::complex<T>,N> fft(rhs.shape());
            return covariance<T,N>(lhs, rhs, fft);
        }

        template <class T, int N>
        inline blitz::Array<T,N>
        correlation(
            blitz::Array<T,N> lhs,
            blitz::Array<T,N> rhs,
            Chirp_Z_transform<std::complex<T>,N>& fft
        ) {
            const auto& lhs_stat = statistics(lhs);
            const auto& rhs_stat = statistics(rhs);
            blitz::Array<T,N> result(covariance<T,N>(
                blitz::Array<T,N>(lhs - lhs_stat.mean()),
                blitz::Array<T,N>(rhs - rhs_stat.mean()),
                fft
            ));
            const T denominator = lhs_stat.variance() * rhs_stat.variance();
            if (denominator > 0) {
                result /= std::sqrt(denominator);
            }
            return result;
        }

        template <class T, int N>
        inline blitz::Array<T,N>
        correlation(blitz::Array<T,N> lhs, blitz::Array<T,N> rhs) {
            Chirp_Z_transform<std::complex<T>,N> fft(rhs.shape());
            return correlation<T,N>(lhs, rhs, fft);
        }

        template <class T, int N>
        class ARMA_wavy_surface_generator_base:
            public Wavy_surface_generator<T,N> {

        private:
            typedef Wavy_surface_generator<T,N> base_type;

        public:
            using typename base_type::array_type;
            using typename base_type::grid_type;

        public:
            typedef std::mt19937 prng_type;
            typedef Yule_walker_solver<T,N> solver_type;
            typedef std::unique_ptr<solver_type> solver_ptr;
            typedef Wavy_surface_generator<T,N> generator_type;
            typedef std::unique_ptr<generator_type> generator_ptr;
            typedef blitz::TinyVector<T,N> vec3;
            typedef blitz::RectDomain<N> rect3;

        private:
            array_type _acf;
            generator_ptr _acf_generator;
            array_type _coefficients;
            T _whitenoisevariance{0};
            prng_type _prng;
            solver_ptr _solver;
            vec3 _alpha{0,0,0.01};

        public:

            ARMA_wavy_surface_generator_base() = default;

            virtual
            ~ARMA_wavy_surface_generator_base() = default;

            inline explicit
            ARMA_wavy_surface_generator_base(array_type acf_in) {
                this->acf(acf_in);
            }

            inline bool
            has_coefficients() const {
                return this->_coefficients.size() > 0;
            }

            inline bool
            has_acf() const {
                return this->_acf.size() > 0;
            }

            inline bool
            has_acf_generator() const {
                return static_cast<bool>(this->_acf_generator);
            }

            inline void
            acf_generator(generator_ptr ptr) {
                this->_acf_generator = std::move(ptr);
                this->_acf.free();
            }

            inline const generator_ptr&
            acf_generator() const {
                return this->_acf_generator;
            }

            inline void
            acf(array_type acf_in) {
                this->_acf.reference(acf_in);
                this->_acf_generator = nullptr;
            }

            inline const array_type&
            acf() const noexcept {
                return this->_acf;
            }

            inline const array_type&
            coefficients() const noexcept {
                return this->_coefficients;
            }

            inline T
            white_noise_variance() const noexcept {
                return this->_whitenoisevariance;
            }

            inline T
            variance() const noexcept {
                return *this->_acf.data();
            }

            void
            seed();

            template <class SeedSequence>
            inline void
            seed(SeedSequence& seq) {
                this->_prng.seed(seq);
            }

            inline void
            coefficient_solver(solver_ptr&& ptr) {
                this->_solver = std::move(ptr);
            }

            inline const solver_type*
            coefficient_solver() const noexcept {
                return this->_solver.get();
            }

            inline solver_type*
            coefficient_solver() noexcept {
                return this->_solver.get();
            }

            inline void
            decay(const vec3& rhs) {
                this->_alpha = rhs;
            }

            inline const vec3&
            decay() const noexcept {
                return this->_alpha;
            }

            void
            calculate_coefficients();

        protected:

            void
            generate_white_noise(array_type result);

        };

        template <class T, int N>
        class AR_wavy_surface_generator:
            public ARMA_wavy_surface_generator_base<T,N> {

        private:
            typedef ARMA_wavy_surface_generator_base<T,N> base_type;

        public:
            using typename base_type::array_type;
            using typename base_type::grid_type;
            using typename base_type::solver_ptr;

        private:
            array_type _buffer;
            int _offset = 0;

        public:

            inline
            AR_wavy_surface_generator() = default;

            virtual
            ~AR_wavy_surface_generator() = default;

            void
            generate(const grid_type& grid, array_type& result) override;

        };

        template <class T, int N>
        class MA_wavy_surface_generator:
            public ARMA_wavy_surface_generator_base<T,N> {

        private:
            typedef ARMA_wavy_surface_generator_base<T,N> base_type;

        public:
            typedef Yule_walker_solver<T,N> solver_type;
            typedef std::unique_ptr<solver_type> solver_ptr;
            typedef Fixed_point_iteration_yule_walker_solver<T,N>
                default_coefficient_solver;

        public:
            using typename base_type::array_type;
            using typename base_type::grid_type;

        public:

            inline
            MA_wavy_surface_generator() {
                this->coefficient_solver(
                    solver_ptr{new default_coefficient_solver}
                );
            }

            MA_wavy_surface_generator(const MA_wavy_surface_generator&) = default;

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

            void
            generate(const grid_type& grid, array_type& result) override;

        };

    }

}

#endif // vim:filetype=cpp