main.cpp

   1 #include <iostream>
   2 #include <vector>
   3 #include <cstdlib>
   4 #include <cmath>
   5 #include <ctime>
   6
   7 #define PI 3.141592653589792
   8 #define DEFAULT_SIZE 10
   9
  10 using namespace std;
  11
  12 class State;
  13 class Zero;
  14 class Probability;
  15 class Energy;
  16
  17 // End of forward decl
  18
  19 class State {
  20     private:
  21     public:
  22         struct Change {
  23             int i;
  24             int j;
  25         };
  26
  27         State () : m_spins(DEFAULT_SIZE) {}
  28         State (int size) : m_spins(size) {}
  29         State (const State& other) : m_spins(other.m_spins) {}
  30
  31         Change trialChange () {
  32             Change ret;
  33             ret.i = (rand()/(RAND_MAX+1.0))*size();
  34             ret.j = (rand()/(RAND_MAX+1.0))*(size()-1);
  35             if (ret.j>=ret.i) ret.j++;
  36             return ret;
  37         }
  38
  39         State applyChange (const Change& change) {
  40             State ret(*this);
  41             ret.m_spins[change.i] = !ret.m_spins[change.i];
  42             ret.m_spins[change.j] = !ret.m_spins[change.j];
  43             return ret;
  44         }
  45
  46         int size () const { return m_spins.size(); }
  47         int at (int i) const { return m_spins[i]?1:-1; }
  48         bool parallel (int i, int j) const { return m_spins[i] == m_spins[j]; }
  49
  50         static State up (int size = DEFAULT_SIZE) {
  51             State ret(size);
  52             for (int i=0;i<size;i++) {
  53                 ret.m_spins[i] = true;
  54             }
  55             return ret;
  56         }
  57
  58         static State down (int size = DEFAULT_SIZE) {
  59             State ret(size);
  60             for (int i=0;i<size;i++) {
  61                 ret.m_spins[i] = false;
  62             }
  63             return ret;
  64         }
  65
  66         static State random (int size = DEFAULT_SIZE) {
  67             State ret(size);
  68             for (int i=0;i<size;i++) {
  69                 ret.m_spins[i] = rand()/(RAND_MAX+1.0) < 0.5;
  70             }
  71             return ret;
  72         }
  73
  74         friend ostream& operator << (ostream&, const State&);
  75         const State& operator = (const State& other) { m_spins = other.m_spins; return *this; }
  76
  77     private:
  78         std::vector<bool> m_spins;
  79 };
  80
  81 ostream& operator << (ostream& out, const State& s) {
  82     for (std::vector<bool>::const_iterator iter = s.m_spins.begin();iter != s.m_spins.end();iter++) {
  83         if (iter!=s.m_spins.begin() && iter!=s.m_spins.end()) out << ", ";
  84         out << *iter;
  85     }
  86     return out;
  87 }
  88
  89 class VarProb {
  90     private:
  91         double distance (double i, double j) { return 2*sin(PI*abs(i-j)); }
  92     public:
  93         VarProb() : alpha(1.0) {}
  94         VarProb(double a) : alpha(a) {}
  95
  96         double compute (const State& state) {
  97             double x = 0.0;
  98             int L = state.size();
  99             int ups = 0;
 100             for (int i=0;i<L;i++) {
 101                 ups += state.at(i);
 102                 for (int j=i+1;j<L;j++) {
 103                     double d = distance(double(i)/L, double(j)/L);
 104                     x += log(d*d)*(state.parallel(i, j)?1.0:-1.0);
 105                 }
 106             }
 107             if (ups!=0) return 0.0;
 108             return exp(alpha*x);
 109         }
 110     private:
 111         double alpha;
 112 };
 113
 114 class Observable {
 115     private:
 116     public:
 117 };
 118
 119 template <typename Probability>
 120 class QMC {
 121     private:
 122     public:
 123
 124         QMC (const Probability& p, const Observable &o, const State& s) : m_probability(p), m_obs(o), m_state(s) {
 125             current_probability = m_probability.compute(m_state);
 126         }
 127
 128         State state () const { return m_state; }
 129         void setState (const State& s) { m_state = s; }
 130         Probability probability () const { return m_probability; }
 131         void setProbability (const Probability& p) { m_probability = p; }
 132
 133         bool step() {
 134             bool ret = false;
 135             State::Change trial = m_state.trialChange();
 136             cout << "original state: " << m_state << endl;
 137             State t_state = m_state.applyChange(trial);
 138             cout << " --> proposed:  " << t_state << endl;
 139             double prob = m_probability.compute(t_state);
 140             cout << " --> with probability: " << prob << endl;
 141             double a = prob/current_probability;
 142             cout << " --> and acceptance ratio: " << a << endl;
 143             if (a>(rand()/(RAND_MAX+1.0))) {
 144                 // accept
 145                 m_state = t_state;
 146                 current_probability = prob;
 147                 cout << " --> accepted" << endl;
 148                 ret = true;
 149             } else {
 150                 // reject
 151                 cout << " --> rejected" << endl;
 152                 ret = false;
 153             }
 154             return ret;
 155         }
 156     private:
 157         State m_state;
 158         Probability m_probability;
 159         double current_probability;
 160         Observable m_obs;
 161 };
 162
 163 int main (int argc, char **argv) {
 164     srand(time(NULL));
 165     QMC<VarProb> qmc(VarProb(+0.26), Observable(), State::random(100));
 166     for (int i=0;i<1e5;i++) qmc.step();
 167     cout << State::random(30) << endl;
 168     return 0;
 169 }
 170