#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <cmath>
using namespace std;

const int N = 32;

struct Node {
    double value_;
    bool fixed_;
    Node()
        : value_(0)
        , fixed_(false)
    {}
    void fix(double value) {
        value_ = value;
        fixed_ = true;
    }
    operator double() const { return value_; }
};

struct Grid {
    int width_, height_;
    vector<Node> nodes_;
    Grid(int width, int height)
        : width_ (width)
        , height_(height)
        , nodes_ (width * height)
    {}
    const Node& operator()(int i, int j) const { return nodes_[width_ * j + i]; }
          Node& operator()(int i, int j)       { return nodes_[width_ * j + i]; }
};

void smooth2d_c1() {
    Grid f(N, N);
    Grid gu(N - 1, N);
    Grid gv(N, N - 1);
    
    ifstream ifs_f("gpu_f.txt");
    for (int j = 0; j < N; ++j) for (int i = 0; i < N; ++i) {
        double d0, d1;
        ifs_f >> d0 >> d1;
        if (d1 != 0.0)
            f(i, j).fix(d0);
    }
    
    ifstream ifs_gu("gpu_gu.txt");
    for (int j = 0; j < N; ++j) for (int i = 0; i < N - 1; ++i) {
        double d0, d1;
        ifs_gu >> d0 >> d1;
        if (d1 != 0.0)
            gu(i, j).fix(d0);
    }
    
    ifstream ifs_gv("gpu_gv.txt");
    for (int j = 0; j < N - 1; ++j) for (int i = 0; i < N; ++i) {
        double d0, d1;
        ifs_gv >> d0 >> d1;
        if (d1 != 0.0)
            gv(i, j).fix(d0);
    }
    
    for (int cnt = 0; cnt < 4; ++cnt) {
        cout << "iteration #: " << cnt << endl;
        for (int k = 0; k < 10000; ++k) {
            Grid f_old = f;
            for (int j = 0; j < N; ++j) for (int i = 0; i < N; ++i) {
                if (f(i, j).fixed_)
                    continue;
                double f_00 = f_old(i, j);
                double f_m0 = i == 0     ? f_00 : (f_old(i - 1, j) + gu(i - 1, j));
                double f_p0 = i == N - 1 ? f_00 : (f_old(i + 1, j) - gu(i    , j));
                double f_0m = j == 0     ? f_00 : (f_old(i, j - 1) + gv(i, j - 1));
                double f_0p = j == N - 1 ? f_00 : (f_old(i, j + 1) - gv(i, j    ));
                f(i, j).value_ =
                    0.5   * f_00 +
                    0.125 * f_m0 +
                    0.125 * f_p0 +
                    0.125 * f_0m +
                    0.125 * f_0p;
            }
            for (int j = 0; j < N; ++j) for (int i = 0; i < N - 1; ++i) {
                if (gu(i, j).fixed_)
                    continue;
                double gu_00 = f_old(i + 1, j) - f_old(i, j);
                double gu_m0 = i == 0     ? gu_00 : gu(i - 1, j).fixed_ ? gu(i - 1, j) : (f_old(i    , j) - f_old(i - 1, j));
                double gu_p0 = i == N - 2 ? gu_00 : gu(i + 1, j).fixed_ ? gu(i + 1, j) : (f_old(i + 2, j) - f_old(i + 1, j));
                double gu_0m = j == 0     ? gu_00 : gu(i, j - 1).fixed_ ? gu(i, j - 1) : (f_old(i + 1, j - 1) - f_old(i, j - 1));
                double gu_0p = j == N - 1 ? gu_00 : gu(i, j + 1).fixed_ ? gu(i, j + 1) : (f_old(i + 1, j + 1) - f_old(i, j + 1));
                gu(i, j).value_ =
                    0.5 * gu_00 +
                    0.125 * gu_m0 +
                    0.125 * gu_p0 +
                    0.125 * gu_0m +
                    0.125 * gu_0p;
            }
            for (int j = 0; j < N - 1; ++j) for (int i = 0; i < N; ++i) {
                if (gv(i, j).fixed_)
                    continue;
                double gv_00 = f_old(i, j + 1) - f_old(i, j);
                double gv_0m = j == 0     ? gv_00 : gv(i, j - 1).fixed_ ? gv(i, j - 1) : (f_old(i, j    ) - f_old(i, j - 1));
                double gv_0p = j == N - 2 ? gv_00 : gv(i, j + 1).fixed_ ? gv(i, j + 1) : (f_old(i, j + 2) - f_old(i, j + 1));
                double gv_m0 = i == 0     ? gv_00 : gv(i - 1, j).fixed_ ? gv(i - 1, j) : (f_old(i - 1, j + 1) - f_old(i - 1, j));
                double gv_p0 = i == N - 1 ? gv_00 : gv(i + 1, j).fixed_ ? gv(i + 1, j) : (f_old(i + 1, j + 1) - f_old(i + 1, j));
                gv(i, j).value_ =
                    0.5 * gv_00 +
                    0.125 * gv_0m +
                    0.125 * gv_0p +
                    0.125 * gv_m0 +
                    0.125 * gv_p0;
            }
        }
        {
            stringstream fname;
            fname << "f" << cnt << ".txt";
            ofstream ofs(fname.str().c_str());
            for (int j = 0; j < N; ++j) {
                for (int i = 0; i < N; ++i)
                    ofs << f(i, j) <<  endl;
                ofs << endl;
            }
        }
        {
            stringstream fname;
            fname << "gu" << cnt << ".txt";
            ofstream ofs(fname.str().c_str());
            for (int j = 0; j < N; ++j) {
                for (int i = 0; i < N - 1; ++i)
                    ofs << gu(i, j) <<  endl;
                ofs << endl;
            }
        }
        {
            stringstream fname;
            fname << "gv" << cnt << ".txt";
            ofstream ofs(fname.str().c_str());
            for (int j = 0; j < N - 1; ++j) {
                for (int i = 0; i < N; ++i)
                    ofs << gv(i, j) <<  endl;
                ofs << endl;
            }
        }
    }
}

void smooth2d_c0() {
    Grid f(N, N);
    
    ifstream ifs("gpu_f.txt");
    for (int j = 0; j < N; ++j) for (int i = 0; i < N; ++i) {
        double d0, d1;
        ifs >> d0 >> d1;
        if (d1 != 0.0)
            f(i, j).fix(d0);
    }
    
    for (int cnt = 0; cnt < 4; ++cnt) {
        cout << "iteration #: " << cnt << endl;
        for (int k = 0; k < 10000; ++k) {
            Grid f_old = f;
            for (int j = 0; j < N; ++j) for (int i = 0; i < N; ++i) {
                if (f(i, j).fixed_)
                    continue;
                double f_00 = f_old(i, j);
                double f_m0 = i == 0     ? f_00 : f_old(i - 1, j);
                double f_p0 = i == N - 1 ? f_00 : f_old(i + 1, j);
                double f_0m = j == 0     ? f_00 : f_old(i, j - 1);
                double f_0p = j == N - 1 ? f_00 : f_old(i, j + 1);
                f(i, j).value_ =
                    0.5   * f_00 +
                    0.125 * f_m0 +
                    0.125 * f_p0 +
                    0.125 * f_0m +
                    0.125 * f_0p;
            }
        }
        stringstream fname;
        fname << "c0_f" << cnt << ".txt";
        ofstream ofs(fname.str().c_str());
        for (int j = 0; j < N; ++j) {
            for (int i = 0; i < N; ++i)
                ofs << f(i, j) <<  endl;
            ofs << endl;
        }
    }
}

void smooth1d_c1() {
    vector<Node> f(N);
    vector<Node> g(N - 1);
    f[10].fix(0);
    g[13].fix(1);
    f[50].fix(0);
    g[53].fix(-1);
    
    for (int cnt = 0; cnt < 4; ++cnt) {
        cout << "iteration #: " << cnt << endl;
        for (int k = 0; k < 100000; ++k) {
            vector<Node> f_old = f;
            for (int i = 0; i < N; ++i) {
                if (f[i].fixed_)
                    continue;
                double f_0 = f_old[i];
                double f_m = i == 0     ? f_0 : (f_old[i - 1].value_ + g[i - 1].value_);
                double f_p = i == N - 1 ? f_0 : (f_old[i + 1].value_ - g[i    ].value_);
                f[i].value_ =
                    //0.5 * f_m + 0.5 * f_p;          // without damping
                    0.25 * f_m +
                    0.5  * f_0 +
                    0.25 * f_p;
            }
            for (int i = 0; i < N - 1; ++i) {
                if (g[i].fixed_)
                    continue;
                double g_0 = f_old[i + 1] - f_old[i];
                double g_m = i == 0     ? g_0 : g[i - 1].fixed_ ? g[i - 1] : (f_old[i    ] - f_old[i - 1]);
                double g_p = i == N - 2 ? g_0 : g[i + 1].fixed_ ? g[i + 1] : (f_old[i + 2] - f_old[i + 1]);
                g[i].value_ =
                    //0.5 * g_m + 0.5 * g_p;        // without damping
                    0.25 * g_m +
                    0.5  * g_0 +
                    0.25 * g_p;
            }
        }
        {
            stringstream fname;
            fname << "1d_f" << cnt << ".txt";
            ofstream ofs(fname.str().c_str());
            for (int i = 0; i < N; ++i)
                ofs << i << " " << f[i].value_ << endl;
        }
        {
            stringstream fname;
            fname << "1d_g" << cnt << ".txt";
            ofstream ofs(fname.str().c_str());
            for (int i = 0; i < N - 1; ++i)
                ofs << i << " " << g[i].value_ << endl;
        }
    }
}
void main() {
    //smooth1d_c1();
    //smooth2d_c0();
    smooth2d_c1();
}
