/*
  This code uses C++11 
  Compile with g++-4.9 -std=c++11 
  Jeremy thinks this code is written by Yuri and is a basic spin flip algorithm 
  simulator of a conserved-order-parameter Ising model or a Schelling model.
  */
#include <string>
#include <iostream>
#include <random>
#include <fstream>
int main(int argc, char** argv) {
  std::random_device rd; //truly random, but slow
  std::mt19937 mt(rd());  //Mersenne Twister, seeding with truly random number
  std::uniform_real_distribution<double> urnd(0., 1.); //uniform double on [0, 1)

  const double      p = std::stod(argv[1]); //density

  const int         n = 200, 
                    nsmp = 100,
                    n_sq = n*n,
                    nswp_per_smp = 2,
                    n_mc_per_smp = nswp_per_smp*n_sq,
                    n_ones = p*n_sq;

  int               lat[n][n],
                    shuffle[n_sq];

  const std::string fl_prefix = "frame_",
                    fl_suffix = ".bin";

  for(int i = 0; i < n; i++) {
    for(int j = 0; j < n; j++) {
      lat[i][j] = 0;
    }
  }

  for(int i = 0; i < n_sq; i++) { //Knuth shuffle
    int j = (i + 1)*urnd(mt); 
    shuffle[i] = shuffle[j];
    shuffle[j] = i;
  }

  for(int i = 0; i < n_ones; i++) {
    lat[shuffle[i]/n][shuffle[i] % n] = 1;
  }

  for(int t = 0; t < nsmp; t++) {
    for(int mv_id = 0; mv_id < n_mc_per_smp; mv_id++) {
      int rnd_i = urnd(mt)*n,
          rnd_j = urnd(mt)*n,
          rnd_i2 = urnd(mt)*n,
          rnd_j2 = urnd(mt)*n,
          nbr_count = 0,
          nbr_count2 = 0,
          val_at_site = lat[rnd_i][rnd_j];

      for(int i = -1; i < 2; i++) {
        for(int j = -1; j < 2; j++) {
          int   nbr_i = (rnd_i + i + n) % n, //with periodic boundaries
                nbr_j = (rnd_j + j + n) % n,
                nbr_i2 = (rnd_i2 + i + n) % n,
                nbr_j2 = (rnd_j2 + j + n) % n;

          if(lat[nbr_i][nbr_j] == val_at_site) {
            nbr_count++;
          }
          if(lat[nbr_i2][nbr_j2] == val_at_site) {
            nbr_count2++;
          }
        }
      }

      if(nbr_count < nbr_count2) {
        lat[rnd_i][rnd_j] = lat[rnd_i2][rnd_j2];
        lat[rnd_i2][rnd_j2] = val_at_site;
      } else if(nbr_count == nbr_count2 && urnd(mt) < 0.5) {
        lat[rnd_i][rnd_j] = lat[rnd_i2][rnd_j2];
        lat[rnd_i2][rnd_j2] = val_at_site;
      }
    }
    std::cout << "Got frame " << t << "\n";
    std::string  file_nm = fl_prefix + std::to_string(t) + fl_suffix;
    std::ofstream ofile(file_nm.c_str(), std::ios::binary); //binary file writing
    ofile.write((const char*) &lat[0][0], sizeof(int)*n_sq);
    ofile.close();
  }
  return 0;
}