#include <vector>
#include <iostream>
#include <iomanip>
#include <cassert>
#include <numeric>

using real = double ;

constexpr real KILL_RATE { 0.062l };
constexpr real FEED_RATE { 0.03l };
constexpr real DT { 1.0l };
constexpr real DIFFUSION_RATE_U { 0.1l };
constexpr real DIFFUSION_RATE_V { 0.05l };

template<typename Collection>
void drop( Collection & data, std::size_t nb_rows, std::size_t nb_cols, real value ) {
    const std::size_t center_row { nb_rows/2ul };
    const std::size_t center_col { nb_cols/2ul };
    for (std::size_t row = (center_row-2ul) ; row < (center_row+2ul); row++) {
        for (std::size_t col = (center_col-3ul); col < (center_col+3ul); col++) {
            data[row*nb_cols+col] = value ;
        }
    }
}

void transform(
    std::vector<real> const & iu, std::vector<real> const & iv,
    std::vector<real> & ou, std::vector<real> & ov,
    std::size_t nb_rows, std::size_t nb_cols ) {
 
    for (std::size_t r = 0 ; r < (nb_rows-2); r++) {
        for (std::size_t c = 0; c < (nb_cols-2); c++) {

            real u = iu[(r+1)*nb_cols+c+1];
            real v = iv[(r+1)*nb_cols+c+1];
            real uvv = u*v*v;

            real full_u = 0.l;
            real full_v = 0.l;
            for(std::size_t dr = 0ul; dr < 3ul; ++dr){
                for(std::size_t dc = 0ul; dc < 3ul; ++dc){
                    full_u += iu[(r+dr)*nb_cols+c+dc] - u;
                    full_v += iv[(r+dr)*nb_cols+c+dc] - v;
                }
            }

            real du = DIFFUSION_RATE_U*full_u - uvv + FEED_RATE*(static_cast<real>(1.0l) - u);
            real dv = DIFFUSION_RATE_V*full_v + uvv - (FEED_RATE + KILL_RATE)*v;

            ou[(r+1)*nb_cols+c+1] = u + du*DT;
            ov[(r+1)*nb_cols+c+1] = v + dv*DT;

        }
    }
}

int main( int argc, char * argv[] ) {

    // runtime parameters
    assert(argc==5) ;
    std::size_t nb_rows {std::stoul(argv[1])};
    std::size_t nb_cols {std::stoul(argv[2])};
    std::size_t nb_images {std::stoul(argv[3])};
    std::size_t nb_iterations {std::stoul(argv[4])};
    std::cout << std::fixed << std::setprecision(6);
    assert(nb_iterations%2==0);
    assert(nb_images<=1000);

    try {

        // Temporary work images
        std::vector<real> u1(nb_rows*nb_cols,1.l);
        std::vector<real> v1(nb_rows*nb_cols,0.l);
        std::vector<real> u2(nb_rows*nb_cols,1.l);
        std::vector<real> v2(nb_rows*nb_cols,0.l);
        drop(u1, nb_rows, nb_cols, 0.l);
        drop(v1, nb_rows, nb_cols, 1.l);

        // Final images sequence
        std::vector<std::vector<real>> images;
        images.push_back(u1);

        // Iterations
        std::cout<<std::endl ;
        for ( std::size_t image = 0 ; image < nb_images ; ++image ) {
            for ( std::size_t iter = 0 ; iter < nb_iterations ; ++iter ) {
                transform( u1, v1, u2, v2, nb_rows, nb_cols );
                std::swap(u1,u2) ; std::swap (v1,v2) ;
            }
            std::cout << '.' << std::flush;
            images.push_back(u1);
        }
        std::cout << std::endl;
        
        // Final product
        std::cout << std::endl;
        real product = std::inner_product(u1.begin(), u1.end(), v1.begin(),static_cast<real>(0.l));
        std::cout
          << "product : " << product
          << std::endl;
        
    }
    catch (std::exception & e) {
      std::cout << e.what() << std::endl;
    }
    catch (const char * e) {
      std::cout << e << std::endl;
    }

    // End
    return 0;
}