Untitled

mail@pastecode.io avatar
unknown
c_cpp
21 days ago
16 kB
4
Indexable
Never
#include <armadillo>
#include <cassert>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iostream>
#include <random>

#include "graph.hpp"
#include "pcg_random.hpp"
#include <unordered_set>
#include <vector> 
#include <cassert>  


using namespace std;

std::ostream& log_out = std::cout;

void
Graph::load(potts_model* model)
{
  params = model;
}

ostream&
Graph::print_distribution(ostream& os)
{
  vector<size_t> conf(n);

  double norm = 0;
  while (true) {
    double x = 0;
    for (size_t i = 0; i < n; ++i) {
      x += params->h(conf[i], i);
    }
    for (size_t i = 0; i < n; ++i) {
      for (size_t j = i + 1; j < n; ++j) {
        x += params->J(i, j)(conf[i], conf[j]);
      }
    }
    norm += exp(x);
    size_t j = 0;
    while (j < n && ++conf[j] == q) {
      conf[j] = 0;
      j++;
    }
    if (j == n) {
      break;
    }
  }
  while (true) {
    double x = 0;
    for (size_t i = 0; i < n; ++i) {
      x += params->h(conf[i], i);
    }
    for (size_t i = 0; i < n; ++i) {
      for (size_t j = i + 1; j < n; ++j) {
        x += params->J(i, j)(conf[i], conf[j]);
      }
    }
    os << "G2 " << exp(x) / norm << endl;
    size_t j = 0;
    while (j < n && ++conf[j] == q) {
      conf[j] = 0;
      j++;
    }
    if (j == n) {
      break;
    }
  }
  return os;
};

void
Graph::sample_mcmc(arma::Mat<int>* ptr,
                   size_t m,
                   size_t mc_iters0,
                   size_t mc_iters,
                   long int seed,
                   double temperature)
{
  pcg32 rng;
  rng.seed(seed);
  std::uniform_real_distribution<double> uniform(0, 1);

  arma::Col<size_t> conf = arma::Col<size_t>(n);
  std::unordered_set<size_t> sample_positions = {251, 252, 253, 254, 255, 257, 258, 259, 262, 263, 264, 265, 266, 267, 268};
  assert(sample_positions.size() == 15);

  // Hardcoded WT states for non-sampled positions
  std::vector<int> wt_states = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 11, 15, 9, 15, 0, 2, 0, 20, 16, 17, 16, 1, 1, 18, 10, 1, 1, 18, 0, 17, 10, 5, 18, 0, 10, 16, 18, 3, 0, 0, 0, 15, 6, 0, 0, 18, 8, 1, 3, 16, 5, 0, 0, 16, 1, 12, 14, 0, 0, 0, 0, 0, 4, 8, 15, 20, 16, 0, 0, 4, 18, 0, 0, 17, 13, 0, 20, 0, 0, 7, 18, 17, 16, 18, 19, 17, 9, 6, 18, 17, 13, 13, 1, 0, 0, 0, 0, 0, 12, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 17, 0, 14, 6, 0, 0, 4, 3, 0, 0, 0, 0, 0, 0, 0, 0, 18, 0, 5, 7, 1, 13, 20, 18, 1, 12, 14, 0, 0, 0, 0, 0, 6, 0, 19, 20, 3, 8, 17, 9, 17, 0, 5, 0, 12, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 9, 0, 0, 0, 0, 3, 3, 10, 10, 2, 6, 1, 1, 17, 1, 6, 12, 11, 10, 7, 19, 19, 5, 3, 14, 12, 9, 3, 14, 8, 9, 15, 20, 0, 10, 4, 4, 7, 13, 0, 0, 4, 0, 0, 0, 0, 0, 0, 0, 9, 14, 9, 8, 12, 12, 6, 4, 14, 0, 11, 5, 3, 18, 9, 4, 1, 8, 3, 17, 0, 0, 9, 12, 7, 14, 10, 3, 0, 16, 9, 10, 0, 0, 5, 4, 20, 5, 9, 4, 9, 1, 5, 13, 20, 10, 7, 10, 6, 18, 5, 13, 3, 7, 18, 8, 3, 11, 5, 8, 12, 6, 20, 0, 0, 15, 10, 0, 16, 10, 7, 6, 0, 13, 17, 13, 18, 9, 4, 0, 0, 6, 16, 9, 0, 3, 13, 0, 15, 6, 6, 8, 5, 3, 0, 0, 0, 0, 1, 0, 0, 0, 0, 18, 5, 17, 15, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 3, 14, 16, 9, 10, 10, 17, 0, 16, 0, 15, 7, 3, 0, 0, 0, 0, 5, 9, 4, 9, 0, 0, 0, 12, 10, 9, 0, 0, 0, 4, 8, 16, 3, 10, 8, 9, 9, 4, 10, 0, 17, 4, 6, 9, 0, 1, 10, 6, 10, 16, 7, 17, 0, 0, 0, 0, 20, 1, 0, 12, 0, 18, 0, 15, 8, 12, 7, 18, 8, 12, 10, 19, 6, 1, 3, 5, 3, 16, 0, 0, 12, 6, 0, 0, 12, 10, 9, 1, 8, 20, 18, 17, 3, 16, 3, 0, 16, 0, 12, 0, 0, 0, 0, 0, 1, 0, 0, 16, 8, 0, 0, 6, 0, 0, 0, 11, 9, 9, 20, 5, 18, 6, 0, 18, 12, 0, 0, 16, 0, 1, 0, 6, 0, 9, 18, 1, 8, 0, 16, 1, 0, 0, 9, 0, 0, 0, 0, 0, 4, 8, 9, 0, 4, 0, 3, 0, 0, 12, 8, 6, 1, 0, 0, 0, 14, 0, 18, 10, 0, 0, 6, 10, 5, 17, 10, 16, 17, 0, 6, 14, 3, 16, 19, 12, 14, 17, 0, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0};
  assert(wt_states.size() == 546);

  // Initialize configuration
  for (size_t i = 0; i < n; ++i) {
      if (sample_positions.find(i) != sample_positions.end()) {
          conf(i) = size_t(q * uniform(rng));
          assert(conf(i) < q);
      } else {
          conf(i) = wt_states[i];
      }
  }

  std::vector<size_t> sample_positions_vector(sample_positions.begin(), sample_positions.end());
  std::uniform_int_distribution<size_t> uniform_index(0, sample_positions_vector.size() - 1);




  for (size_t k = 0; k < mc_iters0; ++k) {
    size_t i = sample_positions_vector[uniform_index(rng)];
    size_t dq = 1 + size_t((q - 1) * uniform(rng));

    size_t q0 = conf(i);
    size_t q1 = (q0 + dq) % q;

    double e0 = -params->h(q0, i);
    double e1 = -params->h(q1, i);
    for (size_t j = 0; j < n; ++j) {
      if (i > j) {
        e0 -= params->J(j, i)(conf(j), q0);
        e1 -= params->J(j, i)(conf(j), q1);
      } else if (i < j) {
        e0 -= params->J(i, j)(q0, conf(j));
        e1 -= params->J(i, j)(q1, conf(j));
      }
    }
    double de = e1 - e0;
    if ((de < 0) || (uniform(rng) < exp(-de / temperature))) {
      conf(i) = q1;
    }
  }

  for (size_t s = 0; s < m; ++s) {
    for (size_t k = 0; k < mc_iters; ++k) {
      size_t i = sample_positions_vector[uniform_index(rng)];
      size_t dq = 1 + size_t((q - 1) * uniform(rng));

      size_t q0 = conf(i);
      size_t q1 = (q0 + dq) % q;

      double e0 = -params->h(q0, i);
      double e1 = -params->h(q1, i);
      for (size_t j = 0; j < n; ++j) {
        if (i > j) {
          e0 -= params->J(j, i)(conf(j), q0);
          e1 -= params->J(j, i)(conf(j), q1);
        } else if (i < j) {
          e0 -= params->J(i, j)(q0, conf(j));
          e1 -= params->J(i, j)(q1, conf(j));
        }
      }
      double de = e1 - e0;
      if ((de < 0) || (uniform(rng) < exp(-de / temperature))) {
        conf(i) = q1;
      }
    }
    for (size_t i = 0; i < n; ++i) {
      (*ptr)(s, i) = conf(i);
    }
  }
  return;
};

void
Graph::sample_mcmc_init(arma::Mat<int>* ptr,
                        size_t m,
                        size_t mc_iters0,
                        size_t mc_iters,
                        arma::Col<int>* init_ptr,
                        long int seed,
                        double temperature)
{
  pcg32 rng;
  rng.seed(seed);
  std::uniform_real_distribution<double> uniform(0, 1);

  arma::Col<size_t> conf = arma::Col<size_t>(n);
  for (size_t i = 0; i < n; ++i) {
    conf(i) = (*ptr)(i);
    assert(conf(i) < q);
  }

  for (size_t k = 0; k < mc_iters0; ++k) {
    size_t i = size_t(n * uniform(rng));
    size_t dq = 1 + size_t((q - 1) * uniform(rng));

    size_t q0 = conf(i);
    size_t q1 = (q0 + dq) % q;

    double e0 = -params->h(q0, i);
    double e1 = -params->h(q1, i);
    for (size_t j = 0; j < n; ++j) {
      if (i > j) {
        e0 -= params->J(j, i)(conf(j), q0);
        e1 -= params->J(j, i)(conf(j), q1);
      } else if (i < j) {
        e0 -= params->J(i, j)(q0, conf(j));
        e1 -= params->J(i, j)(q1, conf(j));
      }
    }
    double de = e1 - e0;
    if ((de < 0) || (uniform(rng) < exp(-de / temperature))) {
      conf(i) = q1;
    }
  }

  for (size_t s = 0; s < m; ++s) {
    for (size_t k = 0; k < mc_iters; ++k) {
      size_t i = size_t(n * uniform(rng));
      size_t dq = 1 + size_t((q - 1) * uniform(rng));

      size_t q0 = conf(i);
      size_t q1 = (q0 + dq) % q;

      double e0 = -params->h(q0, i);
      double e1 = -params->h(q1, i);
      for (size_t j = 0; j < n; ++j) {
        if (i > j) {
          e0 -= params->J(j, i)(conf(j), q0);
          e1 -= params->J(j, i)(conf(j), q1);
        } else if (i < j) {
          e0 -= params->J(i, j)(q0, conf(j));
          e1 -= params->J(i, j)(q1, conf(j));
        }
      }
      double de = e1 - e0;
      if ((de < 0) || (uniform(rng) < exp(-de / temperature))) {
        conf(i) = q1;
      }
    }
    for (size_t i = 0; i < n; ++i) {
      (*ptr)(s, i) = conf(i);
    }
  }
  return;
};

void
Graph::sample_mcmc_zanella(arma::Mat<int>* ptr,
                           size_t m,
                           size_t mc_iters0,
                           size_t mc_iters,
                           long int seed,
                           std::string mode,
                           double temperature)
{
  pcg32 rng;
  rng.seed(seed);
  std::uniform_real_distribution<double> uniform(0, 1);

  arma::Col<size_t> conf = arma::Col<size_t>(n);
  for (size_t i = 0; i < n; ++i) {
    conf(i) = size_t(q * uniform(rng));
    assert(conf(i) < q);
  }

  arma::Mat<double> de = arma::Mat<double>(n, q, arma::fill::zeros);
  arma::Mat<double> g = arma::Mat<double>(n, q, arma::fill::zeros);
  double lambda = 0.0;

  arma::Mat<double> de_old = arma::Mat<double>(n, q, arma::fill::zeros);
  arma::Mat<double> g_old = arma::Mat<double>(n, q, arma::fill::zeros);
  double lambda_old = 0.0;

  // Compute initial neighborhood.
  for (size_t i = 0; i < n; i++) {
    size_t q0 = conf(i);
    double e0 = -params->h(q0, i);
    for (size_t j = 0; j < n; ++j) {
      if (i > j) {
        e0 -= params->J(j, i)(conf(j), q0);
      } else if (i < j) {
        e0 -= params->J(i, j)(q0, conf(j));
      }
    }

    for (size_t q1 = 0; q1 < q; q1++) {
      if (q0 == q1) {
        de(i, q1) = 0.0;
      } else {
        double e1 = -params->h(q1, i);
        for (size_t j = 0; j < n; ++j) {
          if (i > j) {
            e1 -= params->J(j, i)(conf(j), q1);
          } else if (i < j) {
            e1 -= params->J(i, j)(q1, conf(j));
          }
        }
        de(i, q1) = e1 - e0;
      }
    }
  }

  if (mode == "sqrt") {
    g = arma::exp(de * -0.5 / temperature);
    lambda = arma::accu(g) - n; // n*exp(0) needs to be subtracted.
  } else if (mode == "barker") {
    g = 1.0 / (1.0 + arma::exp(de / temperature));
    lambda = arma::accu(g) - .5 * n;
  }

  de_old = de;
  g_old = g;
  lambda_old = lambda;

  for (size_t k = 0; k < mc_iters0; ++k) {
    double rand = uniform(rng) * lambda;
    double r_sum = 0.0;
    size_t i = 0;
    size_t q0 = 0;
    size_t q1 = 0;
    bool done = false;
    for (i = 0; i < n; i++) {
      for (q1 = 0; q1 < q; q1++) {
        if (conf(i) == q1) {
          continue;
        } else {
          r_sum += g(i, q1);
        }
        if (r_sum > rand) {
          done = true;
          break;
        }
      }
      if (done) {
        break;
      }
    }

    double tmp = de(i, q1);
    q0 = conf(i);
    conf(i) = q1;
    for (size_t pos = 0; pos < n; pos++) {
      for (size_t aa = 0; aa < q; aa++) {
        if (pos < i) {
          de(pos, aa) += params->J(pos, i)(conf(pos), q1) -
                         params->J(pos, i)(conf(pos), q0) -
                         params->J(pos, i)(aa, q1) + params->J(pos, i)(aa, q0);
        } else if (pos > i) {
          de(pos, aa) += params->J(i, pos)(q1, conf(pos)) -
                         params->J(i, pos)(q0, conf(pos)) -
                         params->J(i, pos)(q1, aa) + params->J(i, pos)(q0, aa);
        } else {
          if (q1 == aa) {
            de(pos, aa) = 0;
          } else if (q0 == aa) {
            de(pos, aa) = -tmp;
          } else {
            de(pos, aa) += params->h(q1, pos) - params->h(q0, pos);
            for (size_t pos2 = 0; pos2 < n; pos2++) {
              if (pos2 < i) {
                de(pos, aa) += params->J(pos2, i)(conf(pos2), q1) -
                               params->J(pos2, i)(conf(pos2), q0);
              } else if (pos2 > i) {
                de(pos, aa) += params->J(i, pos2)(q1, conf(pos2)) -
                               params->J(i, pos2)(q0, conf(pos2));
              }
            }
          }
        }
      }
    }

    if (mode == "sqrt") {
      g = arma::exp(de * -0.5 / temperature);
      lambda = arma::accu(g) - n; // n*exp(0) needs to be subtracted.
    } else if (mode == "barker") {
      g = 1.0 / (1.0 + arma::exp(de / temperature));
      lambda = arma::accu(g) - .5 * n;
    }

    if (uniform(rng) < lambda_old / lambda) {
      conf(i) = q1;
      de_old = de;
      g_old = g;
      lambda_old = lambda;
    } else {
      conf(i) = q0;
      g = g_old;
      lambda = lambda_old;
      de = de_old;
    }
  }
  for (size_t s = 0; s < m; ++s) {
    for (size_t k = 0; k < mc_iters; ++k) {
      double rand = uniform(rng) * lambda;
      double r_sum = 0.0;
      size_t i = 0;
      size_t q0 = 0;
      size_t q1 = 0;
      bool done = false;
      for (i = 0; i < n; i++) {
        for (q1 = 0; q1 < q; q1++) {
          if (conf(i) == q1) {
            continue;
          } else {
            r_sum += g(i, q1);
          }
          if (r_sum > rand) {
            done = true;
            break;
          }
        }
        if (done) {
          break;
        }
      }

      double tmp = de(i, q1);
      q0 = conf(i);
      conf(i) = q1;
      for (size_t pos = 0; pos < n; pos++) {
        for (size_t aa = 0; aa < q; aa++) {
          if (pos < i) {
            de(pos, aa) += params->J(pos, i)(conf(pos), q1) -
                           params->J(pos, i)(conf(pos), q0) -
                           params->J(pos, i)(aa, q1) +
                           params->J(pos, i)(aa, q0);
          } else if (pos > i) {
            de(pos, aa) += params->J(i, pos)(q1, conf(pos)) -
                           params->J(i, pos)(q0, conf(pos)) -
                           params->J(i, pos)(q1, aa) +
                           params->J(i, pos)(q0, aa);
          } else {
            if (q1 == aa) {
              de(pos, aa) = 0;
            } else if (q0 == aa) {
              de(pos, aa) = -tmp;
            } else {
              de(pos, aa) += params->h(q1, pos) - params->h(q0, pos);
              for (size_t pos2 = 0; pos2 < n; pos2++) {
                if (pos2 < i) {
                  de(pos, aa) += params->J(pos2, i)(conf(pos2), q1) -
                                 params->J(pos2, i)(conf(pos2), q0);
                } else if (pos2 > i) {
                  de(pos, aa) += params->J(i, pos2)(q1, conf(pos2)) -
                                 params->J(i, pos2)(q0, conf(pos2));
                }
              }
            }
          }
        }
      }

      if (mode == "sqrt") {
        g = arma::exp(de * -0.5 / temperature);
        lambda = arma::accu(g) - n; // n*exp(0) needs to be subtracted.
      } else if (mode == "barker") {
        g = 1.0 / (1.0 + arma::exp(de / temperature));
        lambda = arma::accu(g) - .5 * n;
      }

      if (uniform(rng) < lambda_old / lambda) {
        conf(i) = q1;
        de_old = de;
        g_old = g;
        lambda_old = lambda;
      } else {
        conf(i) = q0;
        de = de_old;
        g = g_old;
        lambda = lambda_old;
      }
    }

    for (size_t i = 0; i < n; ++i) {
      (*ptr)(s, i) = (int)conf(i);
    }
  }
  return;
};

ostream&
Graph::print_parameters(ostream& os)
{
  log_out << "printing parameters J" << endl;
  for (size_t i = 0; i < n; ++i) {
    for (size_t j = i + 1; j < n; ++j) {
      for (size_t yi = 0; yi < q; ++yi) {
        for (size_t yj = 0; yj < q; ++yj) {
          os << "J " << i << " " << j << " " << yi << " " << yj << " "
             << params->J(i, j)(yi, yj) << endl;
        }
      }
    }
  }
  log_out << "printing parameters h" << endl;
  for (size_t i = 0; i < n; ++i) {
    for (size_t yi = 0; yi < q; ++yi) {
      os << "h " << i << " " << yi << " " << params->h(yi, i) << endl;
    }
  }
  log_out << "done" << endl;
  return os;
};

void
Graph::print_parameters(FILE* of)
{
  log_out << "printing parameters J" << endl;
  for (size_t i = 0; i < n; ++i) {
    for (size_t j = i + 1; j < n; ++j) {
      for (size_t yi = 0; yi < q; ++yi) {
        for (size_t yj = 0; yj < q; ++yj) {
          fprintf(of,
                  "J %lu %lu %lu %lu %g\n",
                  i,
                  j,
                  yi,
                  yj,
                  params->J(i, j)(yi, yj));
        }
      }
    }
  }
  log_out << "printing parameters h" << endl;
  for (size_t i = 0; i < n; ++i) {
    for (size_t yi = 0; yi < q; ++yi) {
      fprintf(of, "h %lu %lu %g\n", i, yi, params->h(yi, i));
    }
  }
  log_out << "done" << endl;
};
Leave a Comment