Challenge24 2014EC: Image Compression (author: Psyho)

#include <bits/stdc++.h>
#include <unistd.h>
#include <sys/time.h>
#include <thread>

#include <CImg.h>
//#include <emmintrin.h>

using namespace std;
using namespace cimg_library;

#define FORE(it,c)  for(__typeof((c).begin()) it = (c).begin(); it != (c).end(); ++it)
#define FOR(i,a,b)  for(int i=(a);i<(b);++i)
#define REP(i,a)    FOR(i,0,a)
#define ZERO(m)    memset(m,0,sizeof(m))
#define ALL(x)      x.begin(),x.end()
#define PB          push_back
#define S          size()
#define LL          long long
#define ULL        unsigned long long
#define LD          long double
#define MP          make_pair
#define X          first
#define Y          second
#define VC          vector
#define PII        pair <int, int>
#define VI          VC < int >
#define VPII          VC < PII >
#define VVI        VC < VI >
#define VD          VC < double >
#define VVD        VC < VD >
#define VS          VC < string >
#define DB(a)      cerr << #a << ": " << (a) << endl;

template<class T> void print(VC < T > v) {cerr << "[";if (v.S) cerr << v[0];FOR(i, 1, v.S) cerr << ", " << v[i];cerr << "]\n";}
template<class T> string i2s(T x) {ostringstream o; o << x; return o.str(); }
VS splt(string s, char c = ' ') {VS all; int p = 0, np; while (np = s.find(c, p), np >= 0) {if (np != p) all.PB(s.substr(p, np - p)); p = np + 1;} if (p < s.S) all.PB(s.substr(p)); return all;}

int SX, SY;
int bmp[1000][1000];

int px[2000];
int py[2000];
int pv[2000];
int pno;

double calc() {
    double rv = 0;
    REP(x, SX) REP(y, SY) {
        double v = 0;
        double w = 0;
        REP(i, pno) {
            int d = (px[i] - x) * (px[i] - x) + (py[i] - y) * (py[i] - y);
            if (d == 0) {
                v = pv[i];
                w = 1;
                break;
            }
            w += 1.0 / d;
            v += 1.0 / d * pv[i];
        }
        v /= w;
        rv += (v - bmp[x][y]) * (v - bmp[x][y]);
    }
    rv /= SX * SY;
    rv = sqrt(rv);
    return rv * 1000000;
}

double tw[800][800];
double tv[800][800];
void tadd(int ox, int oy, int v) {
    REP(x, SX) REP(y, SY) {
        int d = (ox - x) * (ox - x) + (oy - y) * (oy - y);
        double w = d == 0 ? 1e4 : 1.0 / d;
        tw[x][y] += w;
        tv[x][y] += w * v;
    }
}

void trem(int ox, int oy, int v) {
    REP(x, SX) REP(y, SY) {
        int d = (ox - x) * (ox - x) + (oy - y) * (oy - y);
        double w = d == 0 ? 1e4 : 1.0 / d;
        tw[x][y] -= w;
        tv[x][y] -= w * v;
    }
}

double tcalc() {
    double rv = 0;
    REP(x, SX) REP(y, SY) {
        double a = tv[x][y] / tw[x][y];
        rv += (a - bmp[x][y]) * (a - bmp[x][y]);
    }
    return sqrt(rv / SX / SY) * 1e6;

}

void save(int test) {
    char fout[100];
    sprintf(fout, "J%d.out", test);
    FILE *f = fopen(fout, "w");
    fprintf(f, "%d\n", pno);
    REP(i, pno) fprintf(f, "%d %d %d\n", px[i], py[i], pv[i]);
    fclose(f);
}

void hc(int test, int tries) {
    const int dx[] = {-1,1,0,0,-1,-1,1,1,0,0};
    const int dy[] = {0,0,-1,1,-1,1,-1,1,0,0};

    REP(i, pno) tadd(px[i], py[i], pv[i]);
    double bv = tcalc();

    int lastImp = -10000;

    REP(cur, tries) {
        int p = rand() % pno;
        int t = rand() % 2;
        int d = rand() % 10;
        int v = 3 - rand() % 7;

        if (v == 0 && d >= 8) continue;

        int ox = px[p];
        int oy = py[p];
        int ov = pv[p];
        int nx = px[p] + dx[d];
        int ny = py[p] + dy[d];
        if (nx < 0 || nx >= SX || ny < 0 || ny >= SY) continue;
        int ok = 0;
        REP(i, pno) ok += px[i] == nx && py[i] == ny;
        if (ok) continue;
        int nv = pv[p] + v;
        if (nv < 0 || nv >= 256) continue;
        trem(px[p], py[p], pv[p]);
        px[p] = nx;
        py[p] = ny;
        pv[p] = nv;
        tadd(px[p], py[p], pv[p]);

        double av = tcalc();
        if (av < bv) {
            if (cur > lastImp + 5000) {
                cerr << cur << ' ' << av << ' ' << calc() << endl;
                save(test);
                lastImp = cur;
            }

            bv = av;
        } else {
            trem(px[p], py[p], pv[p]);
            px[p] = ox;
            py[p] = oy;
            pv[p] = ov;
            tadd(px[p], py[p], pv[p]);
        }
    }

    DB(calc());
}

void load(int test) {
    string fimg = i2s(test) + ".bmp";
    CImg<unsigned char> img(fimg.c_str());
    SX = img.width();
    SY = img.height();
    REP(x, img.width()) REP(y, img.height()) bmp[x][y] = img(x, y);
}

void run(int test, int tries = 0) {
    load(test);


    int bx = 0, by = 0;
    FOR(sx, 1, SX) FOR(sy, 1, SY) {
        if (sx * sy <= SX + SY && ((sx * sy >= bx * by - 10 && min(sx, sy) > min(bx, by)) || sx * sy >= bx * by && min(sx, sy) >= min(bx, by))) bx = sx, by = sy;
    }

    pno = 0;
    ZERO(px);
    ZERO(py);
    ZERO(pv);
    REP(x, bx) REP(y, by) {
        px[pno] = (int)((x + 0.5) * SX / bx + 0.5);
        py[pno] = (int)((y + 0.5) * SY / by + 0.5);
        int no = 0;
        FOR(i, -2, 3) FOR(j, -2, 3) pv[pno] += bmp[px[pno]+i][py[pno]+j], no++;
        pv[pno] = (pv[pno] + no / 2) / no;
        pno++;
    }

    while (pno < SX + SY) {
        px[pno] = rand() % SX;
        py[pno] = rand() % SY;
        pv[pno] = bmp[px[pno]][py[pno]];
        pno++;
    }

    DB(calc());

    save(test);

    if (tries) hc(test, tries);

}

void loadSolution(string fn) {
    fstream fs(fn.c_str());
    fs >> pno;
    REP(i, pno) fs >> px[i] >> py[i] >> pv[i];
    fs.close();
}

int main(int argc, char **argv) {
    int t = atoi(argv[1]);
    if (!isdigit(argv[2][0])) {
        load(t);
        DB(1);
        loadSolution(argv[2]);
        DB(2);
        DB(calc());
    } else {
        run(t, atoi(argv[2]));
    }
}