#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <math.h>
using namespace std;

int main()
{
    // Part 1: Declare variables and read file

    //// Part 1.1: Declare variables
    struct Zone {                  // zone struct
        int a;                     // y axis position
        int b;                     // x axis position
        double W,                  // demand of this position
               C,                  // build cost
               Wbar,               // total demand of R region
               Q;                  // assessed value = Wbar/C
        bool isQ,                  // is qualified to build new facility
             isF,                  // is selected to build new facility
             isB;                  // is assigned to a specific facility
    };
    int M,                         // y axis range (1~M)
        N,                         // x axis range (1~N)
        U = 0,                     // amount of qualified zones
        L,                         // amount of facilities to build
        R;                         // the range of facilities, using Manhattan Distance
    double Qstar;                  // specified accessed value
    Zone *V;                       // start pointer of zone array
    int *VF;                       // start pointer of facility index array
    bool **E;                      // start pointer of two dimension responsible array, candidate point <-> Zone position

    char filename[] = "data1.txt"; // file name and default filename
    fstream file;                  // file stream
    int remained,                  // remained step (remained Manhattan Distance for x-direction)
        tm,                        // temp y-position
        cm,                        // current y-position
        cn,                        // current x-position
        zi,                        // zone index (m*N+n)
        u_max,                     // facility index of max demand facility
        i_max;                     // zone index of max-demand facility
    double w_bar,                  // temp w_bar
           q,                      // temp q
           w_bar_max,              // w_bar max
           w,                      // temp w
           total_satisfied_w = 0;  // total satisfied w
    bool debug = false;            // debug flag

    //// Part 1.2: Get input file name
    if (!debug) { // show input and get filename when not in debug mode
        cout << "Enter filename: ";
        cin >> filename;
    }

    //// Part 1.3: Read and parse file
    file.open(filename, ios::in);
    file >> M;
    file >> N;
    file >> Qstar;
    file >> L;
    file >> R;
    V = new Zone[M*N];
    VF = new int[M*N];
    for(int i = 0; i < M*N; i++){
        file >> V[i].W;
        V[i].a = floor(i/N);
        V[i].b = i % N;
        V[i].isQ = false;
        V[i].isF = false;
        V[i].isB = false;
    }
    for(int i = 0; i < M*N; i++) file >> V[i].C;
    file.close();

    //// Part 1.4: Print result
    ////// Part 1.4.1: Print variables
    printf("M=%d; N=%d; Q*=%.1lf; L=%d; R=%d;\n", M, N, Qstar, L, R);

    ////// Part 1.4.2: Print W matrix
    printf("%6c", 'W');
    for (int n = 0; n < N; n++) printf("%6d", n+1);
    printf("\n");
    for (int m = 0; m < M; m++) {
        printf("%6d", m+1);
        for (int n = 0; n < N; n++) printf("%6.1lf", V[m*N+n].W);
        printf("\n");
    }
    printf("\n");

    ////// Part 1.4.3: Print C matrix
    printf("%6c", 'C');
    for (int n = 0; n < N; n++) printf("%6d", n+1);
    printf("\n");
    for (int m = 0; m < M; m++) {
        printf("%6d", m+1);
        for (int n = 0; n < N; n++) printf("%6.1lf", V[m*N+n].C);
        printf("\n");
    }
    printf("\n");

    // Part 2: Calculate Wbar, Qi, VF, U, E:
    //// Part 2.1: Calculate Wbar, Qi, VF and U:
    for (int m = 0; m < M; m++) {
        for (int n = 0; n < N; n++) {
            w_bar = 0;
            for(int dm = -R; dm <= R; dm++) {
                tm = m + dm;
                if(tm < 0 || tm >= M) continue;
                remained = R - abs(dm);
                for(int tn = n-remained; tn <= n+remained; tn++) {
                    if (tn < 0 || tn >= N) continue;
                    w_bar += V[tm*N+tn].W;
                }
            }
            zi = m*N+n;
            q = w_bar / V[zi].C;
            if (q > Qstar) {
                VF[U++] = zi;
                V[zi].isQ = true;
            }
            V[zi].Wbar = w_bar;
            V[zi].Q = q;
        }
    }

    //// Part 2.2: Calculate E
    E = new bool*[U];
    for (int u = 0; u < U; u++) {
        E[u] = new bool[M*N];
        zi = VF[u];
        cm = V[zi].a;
        cn = V[zi].b;
        for(int dm = -R; dm <= R; dm++) {
            tm = cm + dm;
            if(tm < 0 || tm >= M) continue;
            remained = R - abs(dm);
            for(int tn = cn-remained; tn <= cn+remained; tn++) {
                if (tn < 0 || tn >= N) continue;
                E[u][tm*N+tn] = true;
            }
        }
    }

    //// Part 2.3: Print result
    printf("U=%d, ", U);
    printf("VF:");
    for (int u = 0; u < U; u++) printf(" %d", VF[u]);
    printf("\n");
    printf("\n");
    printf("%6s", "W_bar");
    for (int n = 0; n < N; n++) printf("%8d", n+1);
    printf("\n");
    for (int m = 0; m < M; m++) {
        printf("%6d", m+1);
        for (int n = 0; n < N; n++) printf("%8.1lf", V[m*N+n].Wbar);
        printf("\n");
    }
    printf("\n");
    printf("%6s", "Q");
    for (int n = 0; n < N; n++) printf("%6d", n+1);
    printf("\n");
    for (int m = 0; m < M; m++) {
        printf("%6d", m+1);
        for (int n = 0; n < N; n++) printf("%6.2lf", V[m*N+n].Q);
        printf("\n");
    }
    printf("\n");
    printf("\n");
    printf("%6s", "E");
    for (int i = 0; i < M*N; i++) printf("%3d", i);
    printf("\n");
    for (int u = 0; u < U; u++) {
        printf("%6d", u);
        for (int i = 0; i < M*N; i++) printf("%3d", E[u][i]);
        printf("\n");
    }
    printf("\n");

    // Part 3: Greedy method
    for (int cl = 1; cl <= L; cl++) {
        // Part 3.1: Select facility with max Wbar in unselected candidate
        w_bar_max = 0;
        u_max = -1;
        for (int u = 0; u < U; u++) {
            zi = VF[u];
            if(V[zi].isF) continue;
            w_bar = V[zi].Wbar;
            if(w_bar <= w_bar_max) continue;
            w_bar_max = w_bar;
            u_max = u;
        }
        i_max = VF[u_max];
        V[i_max].isF = true;
        total_satisfied_w += w_bar_max;

        // Part 3.2: Update isB, Wbar and E
        for (int i = 0; i < M*N; i++) { // for each zone
            if(!E[u_max][i]) continue; // skip if this zone is not near selected facility
            if(V[i].isB) continue; // skip if this zone is assigned
            V[i].isB = true; // this zone is assigned to selected faciity
            w = V[i].W;
            for (int u = 0; u < U; u++) { // for each facility
                if(V[VF[u]].isF) continue; // skip already selected facilities
                if(!E[u][i]) continue; // skip if assigned zone is not near this facility
                V[VF[u]].Wbar -= w; // remove this zone's demand from facility's Wbar
                E[u][i] = false; // drop the relationship between the zone and the facility
            }
        }

        // Part 3.3: print iteration result
        printf("Iteration %3d: select Facility %3d, satisfied %8.1lf demands for zone ", cl, i_max, w_bar_max);
        for (int i = 0; i < M*N; i++) if(E[u_max][i]) printf("%d,", i);
        printf("\n");
    }
    // Part 3.4: Print total result
    printf("Total satisfied demands: %6.1lf\n", total_satisfied_w);
}