arma.cc

#include <algorithm>
#include <chrono>
#include <cstring>
#include <sstream>
#include <sstream>
#include <stdexcept>
#include <string>

#include <vtestbed/config/real_type.hh>
#include <vtestbed/core/arma.hh>
#include <vtestbed/core/convolution.hh>
#include <vtestbed/core/debug.hh>
#include <vtestbed/base/for_loop.hh>
#include <vtestbed/core/polynomial.hh>
#include <vtestbed/profile/profile.hh>

//#define VTB_DEBUG_ARMA

namespace {

    template <class T, int N>
    inline blitz::Array<T,N>
    exponential_decay(
        const vtb::core::Grid<T,N>& grid,
        const blitz::TinyVector<T,N>& alpha
    ) {
        const auto& delta = grid.delta();
        return generate(
            grid,
            [&grid,&delta,alpha] (const blitz::TinyVector<int,N>& idx) {
                using blitz::abs;
                using blitz::sum;
                using std::abs;
                using std::exp;
                return exp(-(sum(abs(grid(idx, delta) * alpha))));
            }
        );
    }

    template <class T, int N>
    inline void
    exponential_decay(
        blitz::Array<T,N>& surface,
        const vtb::core::Grid<T,N>& grid,
        const blitz::TinyVector<T,N>& alpha
    ) {
        using blitz::variance;
        using std::isfinite;
        using std::sqrt;
        const T old_variance = variance(surface);
        surface *= exponential_decay(grid, alpha);
        const T new_variance = variance(surface);
        const T a = sqrt(old_variance / new_variance);
        if (!isfinite(a)) {
            throw std::runtime_error{"too large exponential decay"};
        }
        surface *= a;
    }

}

template <class T, int N>
void
vtb::core::validate_process(blitz::Array<T,N> phi) {
    VTB_PROFILE_BLOCK(__func__);
    const int n = phi.numElements();
    Polynomial<T> characteristic_polynomial(n);
    std::copy_n(phi.data(), n, characteristic_polynomial.begin());
    characteristic_polynomial.coefficients() =
        -characteristic_polynomial.coefficients();
    characteristic_polynomial.coefficients(0) = 1;
    int num_bad_roots = num_roots_inside_unit_disk(characteristic_polynomial);
    if (num_bad_roots > 0) {
        throw Invalid_ARMA_process(num_bad_roots);
    }
}

template <class T, int N>
blitz::Array<T,N>
vtb::core::auto_covariance(const blitz::Array<T,N>& rhs) {
    VTB_PROFILE_BLOCK(__func__);
    typedef std::complex<T> C;
    using blitz::abs;
    using blitz::pow2;
    using blitz::real;
    Chirp_Z_transform<C,N> fft(rhs.shape());
    blitz::Array<C,N> tmp(rhs.shape());
    tmp = rhs;
    fft.forward(tmp);
    tmp = pow2(abs(tmp));
    fft.forward(tmp);
    const int n = rhs.numElements();
    return blitz::Array<T,N>(real(tmp) / n / n);
}

template <class T, int N>
blitz::Array<T,N>
vtb::core::covariance(
    blitz::Array<T,N> lhs,
    blitz::Array<T,N> rhs,
    Chirp_Z_transform<std::complex<T>,N>& fft
) {
    VTB_PROFILE_BLOCK(__func__);
    typedef std::complex<T> C;
    using namespace blitz;
    if (!all(lhs.shape() == rhs.shape())) {
        throw std::invalid_argument{"covariance: bad shape"};
    }
    blitz::Array<C,N> f1(lhs.shape());
    f1 = lhs;
    fft.forward(f1);
    blitz::Array<C,N> f2(rhs.shape());
    f2 = rhs;
    fft.forward(f2);
    f1 = conj(f1) * f2;
    fft.backward(f1);
    const int n = rhs.numElements();
    return blitz::Array<T,N>(real(f1) / n / n);
}


template <class T, int N>
void
vtb::core::ARMA_wavy_surface_generator_base<T,N>::seed() {
    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);
}

template <class T, int N>
void
vtb::core::ARMA_wavy_surface_generator_base<T,N>::generate_white_noise(
    array_type result
) {
    VTB_PROFILE_BLOCK(__func__);
    const T stdev = std::sqrt(this->white_noise_variance());
    std::normal_distribution<T> normal(T{0}, stdev);
    for (auto& x : result) { x = normal(this->_prng); }
}

template <class T, int N>
void
vtb::core::ARMA_wavy_surface_generator_base<T,N>::calculate_coefficients() {
    if (this->has_acf_generator()) {
        grid_type grid = this->_acf_generator->acf_grid();
        array_type surface(grid.shape());
        this->_acf_generator->generate(grid, surface);
        #if defined(VTB_DEBUG_ARMA)
        dbg::write_animated_surface("base_surface", surface);
        #endif
        exponential_decay(surface, grid, this->decay());
        #if defined(VTB_DEBUG_ARMA)
        dbg::write_animated_surface("base_surface_exp", surface);
        #endif
        auto tmp = auto_covariance(surface);
        tmp.resizeAndPreserve(tmp.shape()/2);
        #if defined(VTB_DEBUG_ARMA)
        dbg::write_animated_surface("acf", tmp);
        #endif
        this->_acf.reference(tmp);
    }
    this->_coefficients.reference(this->_solver->solve(this->acf()));
    #if defined(VTB_DEBUG_ARMA)
    dbg::write_animated_surface("coefficients", this->_coefficients);
    #endif
    validate_process(this->_coefficients);
    this->_whitenoisevariance = this->_solver->white_noise_variance();
    if (this->_whitenoisevariance <= T{0}) {
        std::cerr << "white noise variance = " << _whitenoisevariance << std::endl;
        throw std::out_of_range("white noise variance <= 0");
    }
}

template <class T, int N>
void
vtb::core::AR_wavy_surface_generator<T,N>::generate(
    const grid_type& grid_txy,
    array_type& result
) {
    typedef blitz::TinyVector<int,N> index;
    const auto& phi = this->coefficients();
    auto& buffer = this->_buffer;
    auto& offset = this->_offset;
    if (buffer.extent(1) != grid_txy.extent(1) ||
        buffer.extent(2) != grid_txy.extent(2) ||
        buffer.extent(0) < phi.extent(0) + grid_txy.extent(0)) {
        buffer.resize(
            std::max(phi.extent(0) + 1, phi.extent(0) + grid_txy.extent(0)),
            grid_txy.extent(1),
            grid_txy.extent(2)
        );
        buffer = 0;
        offset = 0;
    }
    const auto nt = grid_txy.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;
    }
    index r0{offset,0,0};
    index r1{offset+grid_txy.extent(0), buffer.extent(1), buffer.extent(2)};
    blitz::RectDomain<N> r{r0,r1-1};
    this->generate_white_noise(buffer(r));
    for_loop<N>(
        r,
        [&] (const index& idx) {
            T sum{0};
            index shape{blitz::min(idx+1, phi.shape())};
            for_loop<N>(
                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;
}

template <class T, int N>
void
vtb::core::MA_wavy_surface_generator<T,N>::generate(
    const grid_type& grid,
    array_type& result
) {
    this->generate_white_noise(result);
    typedef std::complex<T> C;
    typedef blitz::Array<C,N> carray_type;
    using blitz::real;
    const auto& theta = this->coefficients();
    carray_type signal(result.shape()), kernel(theta.shape());
    signal = result;
    kernel = theta;
    kernel(0,0,0) = -1;
    kernel = -kernel;
    Convolution<C,N> convolve(result.shape(), theta.shape());
    result = real(convolve(signal, kernel));
}

template void
vtb::core::validate_process(blitz::Array<VTB_REAL_TYPE,3> phi);

template blitz::Array<VTB_REAL_TYPE,3>
vtb::core::auto_covariance(const blitz::Array<VTB_REAL_TYPE,3>& rhs);

template blitz::Array<VTB_REAL_TYPE,1>
vtb::core::covariance(
    blitz::Array<VTB_REAL_TYPE,1> lhs,
    blitz::Array<VTB_REAL_TYPE,1> rhs,
    Chirp_Z_transform<std::complex<VTB_REAL_TYPE>,1>& fft
);

template blitz::Array<VTB_REAL_TYPE,3>
vtb::core::covariance(
    blitz::Array<VTB_REAL_TYPE,3> lhs,
    blitz::Array<VTB_REAL_TYPE,3> rhs,
    Chirp_Z_transform<std::complex<VTB_REAL_TYPE>,3>& fft
);

template class vtb::core::ARMA_wavy_surface_generator_base<VTB_REAL_TYPE,3>;
template class vtb::core::AR_wavy_surface_generator<VTB_REAL_TYPE,3>;
template class vtb::core::MA_wavy_surface_generator<VTB_REAL_TYPE,3>;

#if defined(VTB_REAL_TYPE_FLOAT)
template class vtb::core::ARMA_wavy_surface_generator_base<double,3>;
template class vtb::core::AR_wavy_surface_generator<double,3>;
#endif