wl_num_bits.cpp

// ref: https://users.wfu.edu/natalie/s14phy770/projects/

#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <ctime>

#include <vector>
#include <random>

#include <boost/dynamic_bitset.hpp>

int main() {

  /* Factor to multiply g(E) after move is accepted */
  double f = 2.7;
  double ln_f = log(f);

  /* minimum f */
  double const min_f = 1 + 0.01;

  /* Exponent of power law to decrease f after each run */
  double dec_pow = 0.9;

  /* Number of moves to skip between checks of histogram */
  int const skip = 5;

  /* Minimum number of moves for each f */
  int min_steps = 5;

  /* Threshold for flat histogram */
  double const flat_thres = 0.90;

  /* number of bits (system parameter) */
  int const num_bits = 1024;

  /* random number generator seed */
  unsigned int seed = std::random_device()();
  srand(seed);
  fprintf(stderr, "seed=%u\n", seed);

  /* random number generator */
  std::mt19937 mt(seed);
  std::uniform_real_distribution<> rand_r(0.0, 1.0);

  /* begin system specific */
  int const num_bits = 1024;
  double const inv_num_bits_p = 1.0 / (num_bits + 1);

  std::uniform_int_distribution<> rand_bit_num(0, num_bits - 1);

  boost::dynamic_bitset<> bit_array(num_bits);

  for(int i=0; i<num_bits; i++) { bit_array[i] = 1; }
  int energy_cur = num_bits;
  int energy_next = energy_cur;

  std::vector<double> log_rho(num_bits+1);
  std::vector<uint64_t> histogram(num_bits+1);

  /* end system specific */

  /* Repeat loop until f reaches minimum value */
  while(f > min_f) {

    /* whether to continue MC sweep */
    bool continue_mcmc = 1;

    /* Start with flat histogram */
    histogram.clear();
    histogram.resize(num_bits+1);

    /* Reset counter for number of steps */
    int num_total_sweep_step = 0;
    int skip_counter = skip + 1;
    int mc_steps = 0;

    do {

      /* Perform one Monte Carlo sweep */
      for(int sweep_step = 0; sweep_step < num_bits; sweep_step++) {

        /* Choose random site for flipping */
        int const bit_num = rand_bit_num(mt);
        int const bit_value = bit_array[bit_num];

        /* Number of new energy level */
        if(bit_value) { energy_next = energy_cur - 1; }
        else { energy_next = energy_cur + 1; }

        /* Calculate probability for acceptance */
        double const prob = exp(log_rho[energy_cur] - log_rho[energy_next]);
        
        /* Accept if proposed move has lower g (prob>=1.0) */
        bool accept = prob > 1;

        /* or if random number [0.0,1.0] is less than prob */
        if(!accept) {
          double const acc_rand = rand_r(mt);
          accept = prob > acc_rand;
        }

        if(accept) {
          energy_cur = energy_next;
          bit_array[bit_num] = !bit_array[bit_num];
        }

        /* Multiply g[b] by f: add logarithms */
        log_rho[energy_cur] += ln_f;

        /* Update counters */
        histogram[energy_cur] += 1;
        num_total_sweep_step++;

      }

      /* Update counters for Monte Carlo steps */
      mc_steps++;
      skip_counter++;

      /* Check for flat histogram after skipping skip steps */
      if ((mc_steps >= min_steps) && (skip_counter >= skip)) {

        /* Reset skip counter */
        skip_counter = 0;

        /* Check for flat histogram */
        /* Look into each energy level and continue if fraction is larger than given threshold */
        bool is_flat = true;
        double histo_thres = num_total_sweep_step * flat_thres * inv_num_bits_p;
        for (int energy = 0; energy<num_bits+1; energy++) {
          uint64_t const histo_e = histogram[energy];
          bool const is_flat_e = histo_e > histo_thres;
          if(!is_flat_e) {
            is_flat = false;
            break;
          }
        }

        continue_mcmc = !is_flat;

      } else {
        continue_mcmc = 1;
      }

    } while(continue_mcmc);

    /* Normalize (logarithmic) g values, so g(0)=1 */
    for (int energy = 1; energy < num_bits+1; energy++) {
      log_rho[energy] -= log_rho[0];
    }
    log_rho[0] = 0.0;

    /* Give status message with current f */
    fprintf(stderr, "f-1=%.3e\tlog(f)=%g\tmc_steps=%d\n", f - 1., ln_f, mc_steps);

    /* Decrease factor f by a power law */
    f = pow(f, dec_pow);
    ln_f = log(f);

  } 

  /* Print g distribution and final histogram */
  /* Output of g is normalized to the number of ground states */
  /* We have two ground states in the Ising model, so g(0)=2 */
  // fputs("b\tE(b)\tlog(g(E))\tH(E)\n", out_file);
  printf("energy,log(rho(E)),H(E)\n");
  // fputs("-----------------------------------------\n", out_file);
  for (int energy = 0; energy < num_bits + 1; energy++) {
    printf("%d,%g,%d\n", energy, log_rho[energy], histogram[energy]);
  }

}