wind_generator.hh

#ifndef VTESTBED_CORE_WIND_GENERATOR_HH
#define VTESTBED_CORE_WIND_GENERATOR_HH

#include <memory>
#include <random>
#include <chrono>
#include <cmath>

#include <vtestbed/base/blitz.hh>
#include <vtestbed/core/debug.hh>
#include <vtestbed/core/grid.hh>
#include <vtestbed/core/types.hh>
#include <vtestbed/core/yule_walker_solver.hh>
#include <vtestbed/linalg/linear_algebra.hh>

//#define VTB_DEBUG_WIND

namespace vtb {

    namespace core {

        template <class T>
        class Wind_param {
            public:
            typedef Grid<T,4> grid_type;
            typedef blitz::TinyVector<int,4> shape_type;

            public:
            T _c1, _c2;
            T _sigma;
            grid_type _grid;
            shape_type _order;

            public:
            Wind_param() {
                _order = shape_type{4,4,4,4};
                _c1 = 0.01;
                _c2 = 1;
                _grid = Grid<T,4>{
                        {0,-50,-50,0},
                        {2,50,50,50},
                        {10,10,10,10}};
                _sigma = 1;
            }
            ~Wind_param() = default;

            Wind_param& operator=(const Wind_param& rhs) {
                this->_c1 = rhs._c1;
                this->_c2 = rhs._c2;
                this->_sigma = rhs._sigma;
                this->_grid = rhs._grid;
                this->_order = rhs._order;
                return *this;
            }
        };

        //template <class T>
        //class Wind_generator {
        //private:
        //public:

        //  Wind_generator() = default;
        //  Wind_generator(const Wind_generator&) = default;
        //  virtual ~Wind_generator() = default;
        //  Wind_generator& operator=(const Wind_generator&) = default;

        //}

        template <class T>
        class Wind_generator {
        public:
            typedef Grid<T,4> grid_type;
            typedef Array<T,4> array_type;
            typedef std::mt19937 prng_type;
            typedef blitz::TinyVector<int,4> shape_type;
        private:
            array_type _buffer;
            int _offset = 0;
            array_type _coefficients;
            T _whitenoisevariance;
            prng_type _prng;
            shape_type _order_acf;
            Wind_param<T> _param;
        public:

            Wind_generator() {
                createSeed();
            }

            ~Wind_generator() = default;
            Wind_generator(const Wind_generator&) = default;
            Wind_generator& operator=(const Wind_generator&) = default;

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

            blitz::Array<T,4>
            calculate_acf(blitz::TinyVector<int,4> shape) {
                blitz::Array<T,4> result(shape);
                T c1 = this->_param._c1;
                //T c2 = this->_param._c2;
                T sigma = this->_param._sigma;
                grid_type grid = this->_param._grid;
                //T U = 5;
                typedef blitz::TinyVector<int,4> index;
                for_loop<4>(
                    shape,
                    [&] (const index& idx) {
                    T a = idx(0) * grid.length(0) / (shape(0) - 1);
                    T b = idx(1) * grid.length(1) / (shape(1) - 1);
                    T c = idx(2) * grid.length(2) / (shape(2) - 1);
                    T d = idx(3) * grid.length(3) / (shape(3) - 1);
                    result(idx) =
                            sigma*sigma*exp(-c1*(a + b + c + d)) /*
                            (2*c2*(U)*sqrt((b*b + c*c + d*d))) /
                            (pow(c2*sqrt(b*b + c*c + d*d), 2) +
                            a*a*(U)*(U))*/;
                    }
                );
                result(0,0,0,0) = sigma;
                return result;
            }

            void setParam(Wind_param<T> param) {
                this->_param = param;
            }

            void calculateCoefficients() {
                using namespace vtb::core;
                using blitz::abs;
                using blitz::all;
                using blitz::max;
                using blitz::shape;
                shape_type order = this->_param._order;
                T variance = 1.0;
                blitz::Array<T,4> acf(variance * calculate_acf(order + 1));
                #if defined(VTB_DEBUG_WIND)
                const auto _ = blitz::Range::all();
                dbg::write_animated_surface("wind.acf", acf(0,_,_,_));
                #endif
                Gauss_yule_walker_solver<T,4> gauss;
                gauss.chop(true);
                _coefficients.reference(gauss.solve(acf));
                _whitenoisevariance = gauss.white_noise_variance();
                #if defined(VTB_DEBUG_WIND)
                std::clog << "_whitenoisevariance=" << _whitenoisevariance << std::endl;
                #endif
            }

            void createSeed() {
                typedef std::chrono::high_resolution_clock clock_type;
                std::random_device dev;
                std::seed_seq seq{{
                    clock_type::now().time_since_epoch().count(),
                    clock_type::rep(dev())
                }};
                this->_prng.seed(seq);
            }

            array_type
            generate() {
                grid_type grid = this->_param._grid;
                array_type result(grid.shape());
                const auto& phi = this->_coefficients;
                auto& buffer = this->_buffer;
                auto& offset = this->_offset;

                if (buffer.extent(1) != grid.extent(1) ||
                    buffer.extent(2) != grid.extent(2) ||
                    buffer.extent(3) != grid.extent(3) ||
                    buffer.extent(0) < phi.extent(0) + grid.extent(0)) {
                    buffer.resize(
                        std::max(
                            phi.extent(0) + 1,
                            phi.extent(0) + grid.extent(0)),
                        grid.extent(1),
                        grid.extent(2),
                        grid.extent(3)
                    );
                    buffer = 0;
                    offset = 0;
                }

                const auto nt = grid.extent(0);
                if (offset + nt >= buffer.extent(0)) {
                    int shift = offset + nt - buffer.extent(0);
                    int slice_size = buffer.extent(1)*buffer.extent(2);
                    std::memmove(
                        buffer.data(),
                        buffer.data() + shift*slice_size,
                        (offset-shift)*slice_size*sizeof(T)
                    );
                    offset -= shift;
                }
                typedef blitz::TinyVector<int,4> index;
                index r0{offset,0,0,0};
                index r1{offset+grid.extent(0), buffer.extent(1),
                    buffer.extent(2), buffer.extent(3)};
                blitz::RectDomain<4> r{r0,r1-1};
                this->generate_white_noise(buffer(r));

                for_loop<4>(
                    r,
                    [&] (const index& idx) {
                        T sum{0};
                        index shape{blitz::min(idx+1, phi.shape())};
                        for_loop<4>(
                            shape,
                            [&] (const index& idx2) {
                                sum += phi(idx2) * buffer(index{idx-idx2});
                            }
                        );
                        buffer(idx) += sum;
                    }
                );

                result = buffer(r);
                offset += nt;
                //auto m = blitz::minmax(result);
                //std::clog << "m.min=" << m.min << std::endl;
                //std::clog << "m.max=" << m.max << std::endl;
                return result;
            }

            void
            generate_white_noise(
                array_type result
            ) {
                const T stdev = std::sqrt(this->_whitenoisevariance);
                //std::clog << "stdev" << std::endl;
                //std::clog << stdev << std::endl;
                std::normal_distribution<T> normal(T{0}, stdev);
                for (auto& x : result) { x = normal(this->_prng); }
            }


        };

    }

}

#endif // vim:filetype=cpp