print the acceptance

[qmc.git] / main.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <ctime>
extern "C" {
#include <lua.h>
#include <lualib.h>
#include <lauxlib.h>
}

#define PI 3.141592653589792
#define DEFAULT_SIZE 10
#define MAX(a, b) ((a)>(b)?(a):(b))
#define MIN(a, b) ((a)<(b)?(a):(b))

using namespace std;

// End of forward decl

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

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

        Change trialChange () {
            const int L = size();
            Change ret;
            do {
                ret.i = (rand()/(RAND_MAX+1.0))*L;
                ret.j = (rand()/(RAND_MAX+1.0))*(L-1);
                if (ret.j>=ret.i) ret.j = (ret.j+1)%L;
            } while (parallel(ret.i, ret.j));
            return ret;
        }

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

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

        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]; }
        int magnetization () const {
            int m = 0;
            int L = m_spins.size();
            for (int i=0;i<L;i++) {
                m += m_spins[i]?1:-1;
            }
            return m;
        }

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

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

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

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

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

        static SpinChain fromString (const char *s, int size) {
            SpinChain ret(size);
            for (int i=0;i<size;i++) {
                bool c;
                switch (s[i]) {
                    case '1':
                        c = true;
                        break;
                    case '0':
                        c = false;
                        break;
                }
                ret.m_spins[i] = c;
            }
            return ret;
        }

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

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

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

class VarProb {
    private:
        // distance between the two spins in the lattice
        static double _distance (double i, double j) { return 2.0*sin(PI*abs(i-j)); }
        // coefficient between the two spins in the lattice
        static double _V_ij (double i, double j) { return 2.0*log(_distance(i, j)); }
        // (-1)^(# of spin ups in even indices)
        static int _sign (const SpinChain& s) {
            int ret = 1;
            const int L = s.size();
            for (int i=0;i<L;i+=2) {
                ret *= s.at(i)?-1:1;
            }
            return ret;
        }
        // ratio between the signs of the two states
        static int _signRatio (const SpinChain::Change& c) {
            return (c.i%2)==(c.j%2)?1:-1; // each flip on even site bears a sign change
        }
    public:
        typedef std::vector<double> Helper;

        VarProb() : alpha(1.0) {}
        VarProb(double a) : alpha(a) {}

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

        double diff (const SpinChain& state, const SpinChain::Change& c) const {
            double ret = 0.0;
            const int L = state.size();
            const int i = MAX(c.i, c.j);
            const int j = MIN(c.i, c.j);
            if (state.parallel(i, j)) return 0.0;
            for (int k=0;k<j;k++) {
                    double di = _V_ij(double(k)/L, double(i)/L);
                    double dj = _V_ij(double(k)/L, double(j)/L);
                    ret += di*(state.parallel(k, i)?-1.0:1.0);
                    ret += dj*(state.parallel(k, j)?-1.0:1.0);
            }
            for (int k=j+1;k<i;k++) {
                    double di = _V_ij(double(k)/L, double(i)/L);
                    double dj = _V_ij(double(k)/L, double(j)/L);
                    ret += di*(state.parallel(k, i)?-1.0:1.0);
                    ret += dj*(state.parallel(k, j)?-1.0:1.0);
            }
            for (int k=i+1;k<L;k++) {
                    double di = _V_ij(double(k)/L, double(i)/L);
                    double dj = _V_ij(double(k)/L, double(j)/L);
                    ret += di*(state.parallel(k, i)?-1.0:1.0);
                    ret += dj*(state.parallel(k, j)?-1.0:1.0);
            }
            return 2.0 * alpha * ret;
        }

        double ratio (const SpinChain& s, const SpinChain::Change& c) const {
            const double exponent = diff(s, c);
            return _signRatio(c)*exp(exponent);
        }

        double diffHelper (const Helper& h, const SpinChain::Change& c) const {
            const vector<double>& v = h;
            const int L = v.size();
            double d = _V_ij(double(c.i)/L, double(c.j)/L);
            // compensate double counting: i, j are surely antiparallel, so had a -1 sign
            // they are counted once in v[i] and once in v[j]
            // with the wrong sign
            double x = 2.0 * d;
            return -2.0* alpha * (v[c.i]+v[c.j]+x);
        }
        double ratioHelper (const Helper& h, const SpinChain::Change& c) const {
            const double exponent = diffHelper(h, c);
            return _signRatio(c)*exp(exponent);
        }
        void applyChangeHelper (const SpinChain& s, Helper& h, const SpinChain::Change& c) const {
            vector<double>& v = h;
            const int L = v.size();
            const int i = MAX(c.i, c.j);
            const int j = MIN(c.i, c.j);
            v[i] = -v[i];
            v[j] = -v[j];
            for (int k=0;k<j;k++) {
                    double di = _V_ij(double(k)/L, double(i)/L);
                    double dj = _V_ij(double(k)/L, double(j)/L);
                    v[k] -= 2.0*di*(s.parallel(k, i)?-1.0:1.0);
                    v[k] -= 2.0*dj*(s.parallel(k, j)?-1.0:1.0);
            }
            for (int k=j+1;k<i;k++) {
                    double di = _V_ij(double(k)/L, double(i)/L);
                    double dj = _V_ij(double(k)/L, double(j)/L);
                    v[k] -= 2.0*di*(s.parallel(k, i)?-1.0:1.0);
                    v[k] -= 2.0*dj*(s.parallel(k, j)?-1.0:1.0);
            }
            for (int k=i+1;k<L;k++) {
                    double di = _V_ij(double(k)/L, double(i)/L);
                    double dj = _V_ij(double(k)/L, double(j)/L);
                    v[k] -= 2.0*di*(s.parallel(k, i)?-1.0:1.0);
                    v[k] -= 2.0*dj*(s.parallel(k, j)?-1.0:1.0);
            }
            double d = _V_ij(double(c.i)/L, double(c.j)/L);
            v[i] -= 2.0*d;
            v[j] -= 2.0*d;
        }
        void computeHelper (const SpinChain& s, Helper& h) {
            const int L = s.size();
            vector<double>& v = h;
            if (v.size()<L) v.resize(L);
            for (int i=0;i<L;i++) {
                double x = 0.0;
                for (int j=0;j<i;j++) {
                    double d = _V_ij(double(i)/L, double(j)/L);
                    x += d*(s.parallel(i, j)?1.0:-1.0);
                }
                for (int j=i+1;j<L;j++) {
                    double d = _V_ij(double(i)/L, double(j)/L);
                    x += d*(s.parallel(i, j)?1.0:-1.0);
                }
                v[i] = x;
            }
        }

    private:
        double alpha;
};

class Energy {
    private:
        double _H_i (const SpinChain& s, const VarProb& p, int i) const {
            double self = 0.0;
            double other = 0.0;
            const int L = s.size();
            if (s.parallel(i, (i+1)%L)) {
                self += D*J;
            } else {
                SpinChain::Change flip = { i, (i+1)%L };
                self -= D*J;
                other += 0.5*J*(p.ratio(s, flip));
            }
            return (0.25*self)+other;
        }

        double _H_i (const SpinChain& s, const VarProb::Helper& h, const VarProb& p, int i) const {
            double self = 0.0;
            double other = 0.0;
            const int L = h.size();
            if (s.parallel(i, (i+1)%L)) {
                self += D*J;
            } else {
                SpinChain::Change flip = { i, (i+1)%L };
                self -= D*J;
                other += 0.5*J*(p.ratioHelper(h, flip));
            }
            return (0.25*self)+other;
        }
    public:
        typedef double ValueType;

        Energy() : J(1.0), D(1.0) {}
        Energy(double j) : J(j), D(1.0) {}
        Energy(double j, double d) : J(j), D(d) {}

        double compute (const SpinChain& s, const VarProb& p) const {
            double ret = 0.0;
            const int L = s.size();
            for (int i=0;i<L;i++) {
                ret += _H_i (s, p, i);
            }
            return ret;
        }

        double compute (const SpinChain& s, const VarProb::Helper& h, const VarProb& p) const {
            double ret = 0.0;
            const int L = s.size();
            for (int i=0;i<L;i++) {
                ret += _H_i (s, h, p, i);
            }
            return ret;
        }

    private:
        double J;
        double D;
};

template <typename Probability, typename State>
class QMC {
    private:
        void _recompute () {
            m_probability.computeHelper(m_state, m_helper);
        }
    public:

        QMC (const Probability& p, const State& s) : m_probability(p), m_state(s) {
            _recompute();
        }

        typedef typename State::Change Change;
        typedef typename Probability::Helper Helper;

        const State& state () const { return m_state; }
        void setState (const State& s) { m_state = s; _recompute(); }
        const Probability& probability () const { return m_probability; }
        void setProbability (const Probability& p) { m_probability = p; _recompute(); }
        const Change& lastChange () const { return last_change; }
        const Helper& helper () const { return m_helper; }

        bool step() {
            bool ret = false; // whether it's accepted
            Change trial = m_state.trialChange(); // description of a trial step
            double pchange_h = m_probability.ratioHelper(m_helper, trial);
            double a = pchange_h*pchange_h; // pchange is in fact related to an amplitude
            if (a>(rand()/(RAND_MAX+1.0))) {
                // accept
                m_state.applyChange(trial);
                m_probability.applyChangeHelper(m_state, m_helper, trial); // XXX: leave after applyChange()
                last_change = trial;
                ret = true;
            } else {
                // reject
                ret = false;
            }
            return ret;
        }
    private:
        State m_state;
        Probability m_probability;
        Change last_change;
        Helper m_helper;
};

void print_help () {
}

void print_state (const SpinChain& s) { cout << s << endl; }

void load_file (int& niter, int& nspins, double& alpha, SpinChain& init, const char* fn) {
    lua_State *L = lua_open();
    luaL_openlibs(L);
    lua_pushnil(L);
    lua_setglobal(L, "debug");
    lua_pushnil(L);
    lua_setglobal(L, "io");
    luaL_dofile(L, fn);
    lua_getglobal(L, "qmc_result");
    if (!lua_isnil(L, -1)) {
        lua_getfield(L, -1, "iterations");
        niter = lua_tointeger(L, -1);
        lua_pop(L, 1);
        lua_getfield(L, -1, "alpha");
        alpha = lua_tonumber(L, -1);
        lua_pop(L, 1);
        lua_getfield(L, -1, "final_state");
        init = SpinChain::fromString(lua_tostring(L, -1), lua_objlen(L, -1));
        nspins = init.size();
        lua_pop(L, 1);
    }
    lua_pop(L, 1);
    lua_close(L);
}

int main (int argc, char **argv) {
    bool append = false;
    ofstream out;
    srand(time(NULL));
    string filename("");
    string summaryname("");
    int nspins = 20;
    int niter = 1e3;
    double alpha = 0.3;
    SpinChain init = SpinChain::neel(nspins);
    for (int i=1;i<argc;i++) {
        if (argv[i][0]!='-') {
            print_help();
            exit(1);
        }
        switch (argv[i][1]) {
            case 'O':
                append = true;
            case 'o':
                filename = argv[++i];
                break;
            case 'r':
                srand(atoi(argv[++i]));
                break;
            case 'i':
                niter = atoi(argv[++i]);
                break;
            case 'a':
                alpha = atof(argv[++i]);
                break;
            case 'n':
                nspins = atoi(argv[++i]);
                init = SpinChain::half(nspins);
                break;
            case 'c':
                summaryname = argv[++i];
                load_file(niter, nspins, alpha, init, argv[i]);
                break;
            case 'h':
            case 'H':
                print_help();
                exit(0);
                break;
            default:
                print_help();
                exit(1);
                break;
        }
    }

    if (filename!="") out.open(filename.c_str(), append?ofstream::app:ofstream::out);
    QMC<VarProb, SpinChain> qmc(VarProb(alpha), init);
    Energy en;

    int nacc = 0;
    bool accepted = true; // first value is always accepted
    double esum = 0.0;
    double e = 0.0;
    vector<double> csum(nspins/2, 0.0);
    vector<double> c(nspins/2, 0.0);
    for (int i=1;i<=niter;i++) {
        if (!accepted) {
            // e still contains the correct value
        } else {
            const SpinChain& s = qmc.state();
            e = en.compute(s, qmc.helper(), qmc.probability());
            c[0] = s.at(0);
            for (int j=1;j<nspins/2;j++) {
                c[j] = s.parallel(0, j)?1.0:-1.0;
            }
        }
        esum += e;
        for (int j=0;j<nspins/2;j++) {
            csum[j] += c[j];
        }
        accepted = qmc.step();
        nacc += accepted?1:0;
        if (out.is_open()) out << i << " " << double(nacc)/(i) << " " << e/nspins << " " << esum/(i*nspins) << "\n";
    }
    ofstream sfile;
    if (summaryname!="") sfile.open(summaryname.c_str(), ofstream::app);
    ostream &summary = (sfile.is_open())?sfile:cout;
    summary << "qmc_result = {\n";
    summary << "\tinitial_state = \"" << init << "\",\n";
    summary << "\tfinal_state   = \"" << qmc.state() << "\",\n";
    summary << "\talpha = " << alpha << ",\n";
    summary << "\tenergy_sum = " << esum << ", -- average per spin: " << (esum/(niter*nspins)) << "\n";
    summary << "\titerations = " << niter << ",\n";
    summary << "\tacceptances = " << nacc << ",\n";
    summary << "\tmagnetization_sum = " << csum[0] << ",\n";
    summary << "\tcorrelation_sum = { " << (niter-csum[0]*csum[0]/niter) << ", ";
    for (int j=1;j<nspins/2;j++) {
        summary << csum[j] << ", ";
    }
    summary << "},\n";
    summary << "}\n";
    summary << "qmcdata = qmcdata or {}\n";
    summary << "table.insert(qmcdata, qmc_result)\n" << endl;
    return 0;
}