functions to create half-domain and neel state

[qmc.git] / main.cpp

#include <iostream>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <ctime>

#define PI 3.141592653589792
#define DEFAULT_SIZE 10

using namespace std;

class State;
class Zero;
class Probability;
class Energy;

// End of forward decl

class State {
    private:
    public:
        struct Change {
            int i;
            int j;
        };

        State () : m_spins(DEFAULT_SIZE) {}
        State (int size) : m_spins(size) {}
        State (const State& other) : m_spins(other.m_spins) {}

        Change trialChange () {
            Change ret;
            ret.i = (rand()/(RAND_MAX+1.0))*size();
            ret.j = (rand()/(RAND_MAX+1.0))*(size()-1);
            if (ret.j>=ret.i) ret.j++;
            return ret;
        }

        State applyChange (const Change& change) {
            State ret(*this);
            ret.m_spins[change.i] = !ret.m_spins[change.i];
            ret.m_spins[change.j] = !ret.m_spins[change.j];
            return ret;
        }

        int size () const { return m_spins.size(); }
        int at (int i) const { return m_spins[i]?1:-1; }
        bool parallel (int i, int j) const { return m_spins[i] == m_spins[j]; }

        static State up (int size = DEFAULT_SIZE) {
            State ret(size);
            for (int i=0;i<size;i++) {
                ret.m_spins[i] = true;
            }
            return ret;
        }

        static State down (int size = DEFAULT_SIZE) {
            State ret(size);
            for (int i=0;i<size;i++) {
                ret.m_spins[i] = false;
            }
            return ret;
        }

        static State neel (int size = DEFAULT_SIZE) {
            State ret(size);
            for (int i=0;i<size;i++) {
                ret.m_spins[i] = i%2;
            }
            return ret;
        }

        static State half (int size = DEFAULT_SIZE) {
            State ret(size);
            for (int i=0;i<size;i++) {
                ret.m_spins[i] = 2*i < size;
            }
            return ret;
        }

        static State random (int size = DEFAULT_SIZE) {
            State ret(size);
            for (int i=0;i<size;i++) {
                ret.m_spins[i] = rand()/(RAND_MAX+1.0) < 0.5;
            }
            return ret;
        }

        friend ostream& operator << (ostream&, const State&);
        const State& operator = (const State& other) { m_spins = other.m_spins; return *this; }
        State& swap (State& other) { m_spins.swap(other.m_spins); return *this; }

    private:
        std::vector<bool> m_spins;
};

ostream& operator << (ostream& out, const State& s) {
    for (std::vector<bool>::const_iterator iter = s.m_spins.begin();iter != s.m_spins.end();iter++) {
        if (iter!=s.m_spins.begin() && iter!=s.m_spins.end()) out << ", ";
        out << *iter;
    }
    return out;
}

class VarProb {
    private:
        double distance (double i, double j) { return 2*sin(PI*abs(i-j)); }
    public:
        VarProb() : alpha(1.0) {}
        VarProb(double a) : alpha(a) {}

        double compute (const State& state) {
            double x = 0.0;
            int L = state.size();
            int ups = 0;
            for (int i=0;i<L;i++) {
                ups += state.at(i);
                for (int j=i+1;j<L;j++) {
                    double d = distance(double(i)/L, double(j)/L);
                    x += log(d*d)*(state.parallel(i, j)?1.0:-1.0);
                }
            }
            if (ups!=0) return 0.0;
            return exp(alpha*x);
        }
    private:
        double alpha;
};

class Observable {
    private:
    public:
};

template <typename Probability>
class QMC {
    private:
    public:

        QMC (const Probability& p, const Observable &o, const State& s) : m_probability(p), m_obs(o), m_state(s) {
            current_probability = m_probability.compute(m_state);
        }

        State state () const { return m_state; }
        void setState (const State& s) { m_state = s; }
        Probability probability () const { return m_probability; }
        void setProbability (const Probability& p) { m_probability = p; }

        bool step() {
            bool ret = false;
            State::Change trial = m_state.trialChange();
            cout << "original state: " << m_state << endl;
            State t_state = m_state.applyChange(trial);
            cout << " --> proposed:  " << t_state << endl;
            double prob = m_probability.compute(t_state);
            cout << " --> with probability: " << prob << endl;
            double a = prob/current_probability;
            cout << " --> and acceptance ratio: " << a << endl;
            if (a>(rand()/(RAND_MAX+1.0))) {
                // accept
                m_state = t_state;
                current_probability = prob;
                cout << " --> accepted" << endl;
                ret = true;
            } else {
                // reject
                cout << " --> rejected" << endl;
                ret = false;
            }
            return ret;
        }
    private:
        State m_state;
        Probability m_probability;
        double current_probability;
        Observable m_obs;
};

int main (int argc, char **argv) {
    srand(time(NULL));
    QMC<VarProb> qmc(VarProb(+0.26), Observable(), State::half(100));
    for (int i=0;i<1e3;i++) qmc.step();
    return 0;
}