fix accepted value

[qmc.git] / main.cpp

/*
 * Copyright (c) 2009 Mauro Iazzi
 *
 * Permission is hereby granted, free of charge, to any person
 * obtaining a copy of this software and associated documentation
 * files (the "Software"), to deal in the Software without
 * restriction, including without limitation the rights to use,
 * copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the
 * Software is furnished to do so, subject to the following
 * conditions:
 *
 * The above copyright notice and this permission notice shall be
 * included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
 * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
 * OTHER DEALINGS IN THE SOFTWARE.
 */

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cstdlib>
#include <cmath>
#include <ctime>
#include <cstdio>

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


// 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;
                ret.j = (ret.i + 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 sigma (int i) const { return m_spins[i]?1:-1; }
        bool at (int i) const { return m_spins[i]; }
        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 halfRandom (int size = DEFAULT_SIZE) {
            SpinChain ret(size);
            for (int i=0;i<size/2;i++) {
                ret.m_spins[i] = rand()/(RAND_MAX+1.0) < 0.5;
                ret.m_spins[size-1-i] = !ret.m_spins[i];
            }
            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':
                    case 'u':
                    case 'U':
                    case '+':
                    case '^':
                        c = true;
                        break;
                    case '0':
                    case 'd':
                    case 'D':
                    case '-':
                    case 'v':
                        c = false;
                        break;
                }
                ret.m_spins[i] = c;
            }
            return ret;
        }

        const std::string& toString () const {
            const int L = size();
            t.resize(L, '\0');
            for (int i=0;i<L;i++) {
                t[i] = at(i)?'1':'0';
            }
            return t;
        }

        friend std::ostream& operator << (std::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:
        static std::string t;
        std::vector<bool> m_spins;
};

std::string SpinChain::t;

std::ostream& operator << (std::ostream& out, const SpinChain& s) {
    out << s.toString ();
    return out;
}

class VarProb {
    private:
        // distance between the two spins in the lattice
        static double _distance (double i, double j) { return 2.0*std::sin(PI*std::abs(i-j)); }
        // coefficient between the two spins in the lattice
        static double _V_ij (double i, double j) { return 2.0*std::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 SpinChain State;
        typedef struct {
            std::vector<double> V;
            std::vector<double> v_ij;
        } 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.sigma(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 std::vector<double>& V = h.V;
            const std::vector<double>& v_ij = h.v_ij;
            const int L = V.size();
            double d = v_ij[abs(c.i-c.j)];
            // 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 {
            std::vector<double>& V = h.V;
            std::vector<double>& v_ij = h.v_ij;
            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[i-k];
                    double dj = v_ij[j-k];
                    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[i-k];
                    double dj = v_ij[k-j];
                    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[k-i];
                    double dj = v_ij[k-j];
                    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[i-j];
            V[i] -= 2.0*d;
            V[j] -= 2.0*d;
        }
        void computeHelper (const SpinChain& s, Helper& h) {
            const int L = s.size();
            std::vector<double>& V = h.V;
            std::vector<double>& v_ij = h.v_ij;
            if (V.size()<L) V.resize(L);
            if (v_ij.size()<L) v_ij.resize(L);
            for (int i=0;i<L;i++) {
                v_ij[i] = _V_ij(double(i)/L, 0.0);
            }
            for (int i=0;i<L;i++) {
                double x = 0.0;
                for (int j=0;j<i;j++) {
                    double d = v_ij[i-j];
                    x += d*(s.parallel(i, j)?1.0:-1.0);
                }
                for (int j=i+1;j<L;j++) {
                    double d = v_ij[j-i];
                    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.V.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:
        Energy() : J(1.0), D(1.0) {}
        Energy(double j) : J(j), D(1.0) {}
        Energy(double j, double d) : J(j), D(d) {}
        Energy(const Energy &e) : J(e.J), D(e.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;
        }

        double J;
        double D;
    private:
};

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

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

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

        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; }
        double weight () const { return 1.0; }

        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;
};

class Walker {
    private:
        void _recompute () {
            m_prob.computeHelper(m_state, m_helper);
            _fillG();
        }

        void _fillG () {
            const double J = m_energy.J;
            const double D = m_energy.D;
            const int L = m_state.size();
            double e_L = 0.0;
            double self = 0.0;
            int j = 1; // 0 is reserved for no change
            _G.resize(L+1, -1.0); // L spins can flip or no flip at all
            _x.resize(L+1, -1); // L spins can flip or no flip at all
            for (int i=0;i<L;i++) {
                if (m_state.parallel(i, (i+1)%L)) {
                    self += J*D;
                } else {
                    Change c = { i, (i+1)%L };
                    const double g_i = 0.5 * J * m_prob.ratioHelper(m_helper, c); // guiding function
                    e_L += g_i; // e_L is the sum of all the green function column
                    _G[j] = -g_i; // the energy due to the change
                    _x[j++] = i; // the spin that is flipped
                    self -= J*D; // the energy of the state
                }
            }
            for (;j<=L;j++) { // overwrite possible leftover values
                _G[j] = -1.0;
                _x[j] = -1;
            }
            self = 0.25*self;
            _G[0] = Lambda-self;
            _x[0] = -1;
            e_L = self+e_L; // e_L is now correct
            m_localEnergy = e_L;
            //std::cout << e_L-m_energy.compute(m_state, m_helper, m_prob) << std::endl;
            //std::cout << e_L-m_energy.compute(m_state, m_prob) << std::endl;
        }
    public:
        typedef SpinChain State;
        typedef /*typename*/ State::Change Change;
        typedef VarProb Probability;
        typedef /*typename*/ Probability::Helper Helper;

        explicit Walker(const Probability &p, const State &s, Energy e = Energy()) : Lambda(-1.0), m_weight(1.0), m_energy(e), m_state(s), m_prob(p) {
            Lambda = s.size() * 0.25 * e.J * (0.1 + e.D); // surely higher than any value
            _recompute();
        }

        const State& state () const { return m_state; }
        void setState (const State& s) { m_state = s; _recompute(); }
        const Probability& probability () const { return m_prob; }
        void setProbability (const Probability& p) { m_prob = p; _recompute(); }
        const Helper& helper () const { return m_helper; }
        double weight () const { return m_weight; }
        double localEnergy () const { return m_localEnergy; }

        bool step () {
            const int L = m_state.size();
            const double b = (Lambda - m_localEnergy); // will be the new weight
            double r = b*rand()/(RAND_MAX+1.0);
            int j = 0;
            while (r>_G[j] && _G[j]>0) {
                //cout << j << " " << r << " " << G[j] << endl;
                r -= _G[j++];
            }
            //if (_G[j]<0.0) throw std::string("WTF?!");
            // now j contains the index of the spin that must be flipped
            const int pos = _x[j];
            //std::cout << "stuck on " << j << " " << r << " " << _G[j] << " " << pos << std::endl;
            //std::cout << pos << "\n"; // to check for translation invariance
            if (pos == -1) {
                // no change needed?
            } else {
                // let's recompute
                const Change f = { pos, (pos+1)%L };
                m_state.applyChange(f);
                m_prob.applyChangeHelper(m_state, m_helper, f); // XXX: leave after applyChange()
                _fillG();
            }
            m_weight = b;
            return pos!=-1;
        }
    private:
        double Lambda;
        double m_weight;
        double m_localEnergy;
        Energy m_energy;
        State m_state;
        Helper m_helper;
        Probability m_prob;
        std::vector<double> _G; // green function for a flip
        std::vector<int> _x; // which spin is flipped
};

using namespace std;

void print_help () {
}

int main (int argc, char **argv) {
    bool print_state = false;
    ofstream out;
    srand(time(NULL));
    string filename("");
    int nspins = 100;
    int niter = 1e6;
    double alpha = 0.25;
    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 's':
                print_state = 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::neel(nspins);
                break;
            case 'h':
            case 'H':
                print_help();
                exit(0);
                break;
            default:
                print_help();
                exit(1);
                break;
        }
    }

    if (filename!="") out.open(filename.c_str(), ofstream::out);
    //QMC<VarProb> mc(VarProb(alpha), init);
    Walker mc(VarProb(alpha), init);
    Energy en;

    bool accepted = true; // first value is always accepted
    double e = 0.0;
    double w = 0.0;
    if (out.is_open()) {
        out << "#!prepend " << alpha << "\n";
        out << "# alpha = " << alpha << "\n";
        out << "# format = '^([01]+)%s+(.*)$'\n";
    }
    std::string state_str = init.toString();
    for (int i=1;i<=niter;i++) {
        if (!accepted) {
            // e still contains the correct value
        } else {
            const SpinChain& s = mc.state();
            //e = en.compute(s, mc.helper(), mc.probability());
            e = mc.localEnergy();
            if (out.is_open()) state_str = s.toString();
        }
        w = std::log(mc.weight());
        if (out.is_open()) {
            if (print_state) {
                out << state_str << " " << (e/nspins) << " " << w << "\n";
            } else {
                out << "* " << (e/nspins) << " " << w << "\n";
            }
        }
        accepted = mc.step();
    }
    if (out.is_open()) {
    }
    return 0;
}