main.cpp

   1 /*
   2  * Copyright (c) 2009 Mauro Iazzi
   3  *
   4  * Permission is hereby granted, free of charge, to any person
   5  * obtaining a copy of this software and associated documentation
   6  * files (the "Software"), to deal in the Software without
   7  * restriction, including without limitation the rights to use,
   8  * copy, modify, merge, publish, distribute, sublicense, and/or sell
   9  * copies of the Software, and to permit persons to whom the
  10  * Software is furnished to do so, subject to the following
  11  * conditions:
  12  *
  13  * The above copyright notice and this permission notice shall be
  14  * included in all copies or substantial portions of the Software.
  15  *
  16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
  18  * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
  19  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
  20  * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
  21  * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  22  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
  23  * OTHER DEALINGS IN THE SOFTWARE.
  24  */
  25
  26 #include <iostream>
  27 #include <fstream>
  28 #include <string>
  29 #include <vector>
  30 #include <cstdlib>
  31 #include <cmath>
  32 #include <ctime>
  33 #include <cstdio>
  34
  35 #define PI 3.141592653589792
  36 #define DEFAULT_SIZE 10
  37 #define MAX(a, b) ((a)>(b)?(a):(b))
  38 #define MIN(a, b) ((a)<(b)?(a):(b))
  39
  40
  41 // End of forward decl
  42
  43 class SpinChain {
  44     private:
  45     public:
  46         struct Change {
  47             int i;
  48             int j;
  49         };
  50
  51         SpinChain () : m_spins(DEFAULT_SIZE) {}
  52         SpinChain (int size) : m_spins(size) {}
  53         SpinChain (const SpinChain& other) : m_spins(other.m_spins) {}
  54
  55         Change trialChange () {
  56             const int L = size();
  57             Change ret;
  58             do {
  59                 ret.i = (rand()/(RAND_MAX+1.0))*L;
  60                 //ret.j = (rand()/(RAND_MAX+1.0))*(L-1);
  61                 //if (ret.j>=ret.i) ret.j = (ret.j+1)%L;
  62                 ret.j = (ret.i + 1) % L;
  63             } while (parallel(ret.i, ret.j));
  64             return ret;
  65         }
  66
  67         SpinChain tryChange (const Change& change) const {
  68             SpinChain ret(*this);
  69             ret.m_spins[change.i] = !ret.m_spins[change.i];
  70             ret.m_spins[change.j] = !ret.m_spins[change.j];
  71             return ret;
  72         }
  73
  74         SpinChain& applyChange (const Change& change) {
  75             this->m_spins[change.i] = !this->m_spins[change.i];
  76             this->m_spins[change.j] = !this->m_spins[change.j];
  77             return *this;
  78         }
  79
  80         int size () const { return m_spins.size(); }
  81         int sigma (int i) const { return m_spins[i]?1:-1; }
  82         bool at (int i) const { return m_spins[i]; }
  83         bool parallel (int i, int j) const { return m_spins[i] == m_spins[j]; }
  84         int magnetization () const {
  85             int m = 0;
  86             int L = m_spins.size();
  87             for (int i=0;i<L;i++) {
  88                 m += m_spins[i]?1:-1;
  89             }
  90             return m;
  91         }
  92
  93         static SpinChain up (int size = DEFAULT_SIZE) {
  94             SpinChain ret(size);
  95             for (int i=0;i<size;i++) {
  96                 ret.m_spins[i] = true;
  97             }
  98             return ret;
  99         }
 100
 101         static SpinChain down (int size = DEFAULT_SIZE) {
 102             SpinChain ret(size);
 103             for (int i=0;i<size;i++) {
 104                 ret.m_spins[i] = false;
 105             }
 106             return ret;
 107         }
 108
 109         static SpinChain neel (int size = DEFAULT_SIZE) {
 110             SpinChain ret(size);
 111             for (int i=0;i<size;i++) {
 112                 ret.m_spins[i] = i%2;
 113             }
 114             return ret;
 115         }
 116
 117         static SpinChain half (int size = DEFAULT_SIZE) {
 118             SpinChain ret(size);
 119             for (int i=0;i<size;i++) {
 120                 ret.m_spins[i] = 2*i < size;
 121             }
 122             return ret;
 123         }
 124
 125         static SpinChain random (int size = DEFAULT_SIZE) {
 126             SpinChain ret(size);
 127             for (int i=0;i<size;i++) {
 128                 ret.m_spins[i] = rand()/(RAND_MAX+1.0) < 0.5;
 129             }
 130             return ret;
 131         }
 132
 133         static SpinChain halfRandom (int size = DEFAULT_SIZE) {
 134             SpinChain ret(size);
 135             for (int i=0;i<size/2;i++) {
 136                 ret.m_spins[i] = rand()/(RAND_MAX+1.0) < 0.5;
 137                 ret.m_spins[size-1-i] = !ret.m_spins[i];
 138             }
 139             return ret;
 140         }
 141
 142         static SpinChain fromString (const char *s, int size) {
 143             SpinChain ret(size);
 144             for (int i=0;i<size;i++) {
 145                 bool c;
 146                 switch (s[i]) {
 147                     case '1':
 148                     case 'u':
 149                     case 'U':
 150                     case '+':
 151                     case '^':
 152                         c = true;
 153                         break;
 154                     case '0':
 155                     case 'd':
 156                     case 'D':
 157                     case '-':
 158                     case 'v':
 159                         c = false;
 160                         break;
 161                 }
 162                 ret.m_spins[i] = c;
 163             }
 164             return ret;
 165         }
 166
 167         const std::string& toString () const {
 168             const int L = size();
 169             t.resize(L, '\0');
 170             for (int i=0;i<L;i++) {
 171                 t[i] = at(i)?'1':'0';
 172             }
 173             return t;
 174         }
 175
 176         friend std::ostream& operator << (std::ostream&, const SpinChain&);
 177         const SpinChain& operator = (const SpinChain& other) { m_spins = other.m_spins; return *this; }
 178         bool operator == (const SpinChain& other) { return m_spins == other.m_spins; }
 179         SpinChain& swap (SpinChain& other) { m_spins.swap(other.m_spins); return *this; }
 180
 181     private:
 182         static std::string t;
 183         std::vector<bool> m_spins;
 184 };
 185
 186 std::string SpinChain::t;
 187
 188 std::ostream& operator << (std::ostream& out, const SpinChain& s) {
 189     out << s.toString ();
 190     return out;
 191 }
 192
 193 class VarProb {
 194     private:
 195         // distance between the two spins in the lattice
 196         static double _distance (double i, double j) { return 2.0*std::sin(PI*std::abs(i-j)); }
 197         // coefficient between the two spins in the lattice
 198         static double _V_ij (double i, double j) { return 2.0*std::log(_distance(i, j)); }
 199         // (-1)^(# of spin ups in even indices)
 200         static int _sign (const SpinChain& s) {
 201             int ret = 1;
 202             const int L = s.size();
 203             for (int i=0;i<L;i+=2) {
 204                 ret *= s.at(i)?-1:1;
 205             }
 206             return ret;
 207         }
 208         // ratio between the signs of the two states
 209         static int _signRatio (const SpinChain::Change& c) {
 210             return (c.i%2)==(c.j%2)?1:-1; // each flip on even site bears a sign change
 211         }
 212     public:
 213         typedef SpinChain State;
 214         typedef struct {
 215             std::vector<double> V;
 216             std::vector<double> v_ij;
 217         } Helper;
 218
 219         VarProb() : alpha(1.0) {}
 220         VarProb(double a) : alpha(a) {}
 221
 222         double compute (const SpinChain& state) const {
 223             const int L = state.size();
 224             double x = 0.0;
 225             int ups = 0;
 226             for (int i=0;i<L;i++) {
 227                 ups += state.sigma(i);
 228                 for (int j=i+1;j<L;j++) {
 229                     double d = _V_ij(double(i)/L, double(j)/L);
 230                     x += d*(state.parallel(i, j)?1.0:-1.0);
 231                 }
 232             }
 233             if (ups!=0) return 0.0;
 234             return alpha*x;
 235         }
 236
 237         double diff (const SpinChain& state, const SpinChain::Change& c) const {
 238             double ret = 0.0;
 239             const int L = state.size();
 240             const int i = MAX(c.i, c.j);
 241             const int j = MIN(c.i, c.j);
 242             if (state.parallel(i, j)) return 0.0;
 243             for (int k=0;k<j;k++) {
 244                     double di = _V_ij(double(k)/L, double(i)/L);
 245                     double dj = _V_ij(double(k)/L, double(j)/L);
 246                     ret += di*(state.parallel(k, i)?-1.0:1.0);
 247                     ret += dj*(state.parallel(k, j)?-1.0:1.0);
 248             }
 249             for (int k=j+1;k<i;k++) {
 250                     double di = _V_ij(double(k)/L, double(i)/L);
 251                     double dj = _V_ij(double(k)/L, double(j)/L);
 252                     ret += di*(state.parallel(k, i)?-1.0:1.0);
 253                     ret += dj*(state.parallel(k, j)?-1.0:1.0);
 254             }
 255             for (int k=i+1;k<L;k++) {
 256                     double di = _V_ij(double(k)/L, double(i)/L);
 257                     double dj = _V_ij(double(k)/L, double(j)/L);
 258                     ret += di*(state.parallel(k, i)?-1.0:1.0);
 259                     ret += dj*(state.parallel(k, j)?-1.0:1.0);
 260             }
 261             return 2.0 * alpha * ret;
 262         }
 263
 264         double ratio (const SpinChain& s, const SpinChain::Change& c) const {
 265             const double exponent = diff(s, c);
 266             return _signRatio(c)*exp(exponent);
 267         }
 268
 269         double diffHelper (const Helper& h, const SpinChain::Change& c) const {
 270             const std::vector<double>& V = h.V;
 271             const std::vector<double>& v_ij = h.v_ij;
 272             const int L = V.size();
 273             double d = v_ij[abs(c.i-c.j)];
 274             // compensate double counting: i, j are surely antiparallel, so had a -1 sign
 275             // they are counted once in v[i] and once in v[j]
 276             // with the wrong sign
 277             double x = 2.0 * d;
 278             return -2.0* alpha * (V[c.i]+V[c.j]+x);
 279         }
 280         double ratioHelper (const Helper& h, const SpinChain::Change& c) const {
 281             const double exponent = diffHelper(h, c);
 282             return _signRatio(c)*exp(exponent);
 283         }
 284         void applyChangeHelper (const SpinChain& s, Helper& h, const SpinChain::Change& c) const {
 285             std::vector<double>& V = h.V;
 286             std::vector<double>& v_ij = h.v_ij;
 287             const int L = V.size();
 288             const int i = MAX(c.i, c.j);
 289             const int j = MIN(c.i, c.j);
 290             V[i] = -V[i];
 291             V[j] = -V[j];
 292             for (int k=0;k<j;k++) {
 293                     double di = v_ij[i-k];
 294                     double dj = v_ij[j-k];
 295                     V[k] -= 2.0*di*(s.parallel(k, i)?-1.0:1.0);
 296                     V[k] -= 2.0*dj*(s.parallel(k, j)?-1.0:1.0);
 297             }
 298             for (int k=j+1;k<i;k++) {
 299                     double di = v_ij[i-k];
 300                     double dj = v_ij[k-j];
 301                     V[k] -= 2.0*di*(s.parallel(k, i)?-1.0:1.0);
 302                     V[k] -= 2.0*dj*(s.parallel(k, j)?-1.0:1.0);
 303             }
 304             for (int k=i+1;k<L;k++) {
 305                     double di = v_ij[k-i];
 306                     double dj = v_ij[k-j];
 307                     V[k] -= 2.0*di*(s.parallel(k, i)?-1.0:1.0);
 308                     V[k] -= 2.0*dj*(s.parallel(k, j)?-1.0:1.0);
 309             }
 310             double d = v_ij[i-j];
 311             V[i] -= 2.0*d;
 312             V[j] -= 2.0*d;
 313         }
 314         void computeHelper (const SpinChain& s, Helper& h) {
 315             const int L = s.size();
 316             std::vector<double>& V = h.V;
 317             std::vector<double>& v_ij = h.v_ij;
 318             if (V.size()<L) V.resize(L);
 319             if (v_ij.size()<L) v_ij.resize(L);
 320             for (int i=0;i<L;i++) {
 321                 v_ij[i] = _V_ij(double(i)/L, 0.0);
 322             }
 323             for (int i=0;i<L;i++) {
 324                 double x = 0.0;
 325                 for (int j=0;j<i;j++) {
 326                     double d = v_ij[i-j];
 327                     x += d*(s.parallel(i, j)?1.0:-1.0);
 328                 }
 329                 for (int j=i+1;j<L;j++) {
 330                     double d = v_ij[j-i];
 331                     x += d*(s.parallel(i, j)?1.0:-1.0);
 332                 }
 333                 V[i] = x;
 334             }
 335         }
 336
 337     private:
 338         double alpha;
 339 };
 340
 341 class Energy {
 342     private:
 343         double _H_i (const SpinChain& s, const VarProb& p, int i) const {
 344             double self = 0.0;
 345             double other = 0.0;
 346             const int L = s.size();
 347             if (s.parallel(i, (i+1)%L)) {
 348                 self += D*J;
 349             } else {
 350                 SpinChain::Change flip = { i, (i+1)%L };
 351                 self -= D*J;
 352                 other += 0.5*J*(p.ratio(s, flip));
 353             }
 354             return (0.25*self)+other;
 355         }
 356
 357         double _H_i (const SpinChain& s, const VarProb::Helper& h, const VarProb& p, int i) const {
 358             double self = 0.0;
 359             double other = 0.0;
 360             const int L = h.V.size();
 361             if (s.parallel(i, (i+1)%L)) {
 362                 self += D*J;
 363             } else {
 364                 SpinChain::Change flip = { i, (i+1)%L };
 365                 self -= D*J;
 366                 other += 0.5*J*(p.ratioHelper(h, flip));
 367             }
 368             return (0.25*self)+other;
 369         }
 370     public:
 371         Energy() : J(1.0), D(1.0) {}
 372         Energy(double j) : J(j), D(1.0) {}
 373         Energy(double j, double d) : J(j), D(d) {}
 374         Energy(const Energy &e) : J(e.J), D(e.D) {}
 375
 376         double compute (const SpinChain& s, const VarProb& p) const {
 377             double ret = 0.0;
 378             const int L = s.size();
 379             for (int i=0;i<L;i++) {
 380                 ret += _H_i (s, p, i);
 381             }
 382             return ret;
 383         }
 384
 385         double compute (const SpinChain& s, const VarProb::Helper& h, const VarProb& p) const {
 386             double ret = 0.0;
 387             const int L = s.size();
 388             for (int i=0;i<L;i++) {
 389                 ret += _H_i (s, h, p, i);
 390             }
 391             return ret;
 392         }
 393
 394         double J;
 395         double D;
 396     private:
 397 };
 398
 399 template <typename Probability>
 400 class QMC {
 401     private:
 402         void _recompute () {
 403             m_probability.computeHelper(m_state, m_helper);
 404         }
 405     public:
 406
 407         typedef typename Probability::State State;
 408         typedef typename State::Change Change;
 409         typedef typename Probability::Helper Helper;
 410
 411         QMC (const Probability& p, const State& s) : m_probability(p), m_state(s) {
 412             _recompute();
 413         }
 414
 415         const State& state () const { return m_state; }
 416         void setState (const State& s) { m_state = s; _recompute(); }
 417         const Probability& probability () const { return m_probability; }
 418         void setProbability (const Probability& p) { m_probability = p; _recompute(); }
 419         const Change& lastChange () const { return last_change; }
 420         const Helper& helper () const { return m_helper; }
 421         double weight () const { return 1.0; }
 422
 423         bool step() {
 424             bool ret = false; // whether it's accepted
 425             Change trial = m_state.trialChange(); // description of a trial step
 426             double pchange_h = m_probability.ratioHelper(m_helper, trial);
 427             double a = pchange_h*pchange_h; // pchange is in fact related to an amplitude
 428             if (a>(rand()/(RAND_MAX+1.0))) {
 429                 // accept
 430                 m_state.applyChange(trial);
 431                 m_probability.applyChangeHelper(m_state, m_helper, trial); // XXX: leave after applyChange()
 432                 last_change = trial;
 433                 ret = true;
 434             } else {
 435                 // reject
 436                 ret = false;
 437             }
 438             return ret;
 439         }
 440     private:
 441         State m_state;
 442         Probability m_probability;
 443         Change last_change;
 444         Helper m_helper;
 445 };
 446
 447 class Walker {
 448     private:
 449         void _recompute () {
 450             m_prob.computeHelper(m_state, m_helper);
 451             _fillG();
 452         }
 453
 454         void _fillG () {
 455             const double J = m_energy.J;
 456             const double D = m_energy.D;
 457             const int L = m_state.size();
 458             double e_L = 0.0;
 459             double self = 0.0;
 460             int j = 0;
 461             _G.resize(L+1, -1.0); // L spins can flip or no flip at all
 462             _x.resize(L+1, -1); // L spins can flip or no flip at all
 463             for (int i=0;i<L;i++) {
 464                 if (m_state.parallel(i, (i+1)%L)) {
 465                     self += J*D;
 466                 } else {
 467                     Change c = { i, (i+1)%L };
 468                     const double g_i = 0.5 * J * m_prob.ratioHelper(m_helper, c); // guiding function
 469                     e_L += g_i; // e_L is the sum of all the green function column
 470                     _G[j] = -g_i; // the energy due to the change
 471                     _x[j++] = i; // the spin that is flipped
 472                     self -= J*D; // the energy of the state
 473                 }
 474             }
 475             self = 0.25*self; // spin-1/2
 476             _G[j] = Lambda-self;
 477             _x[j++] = -1;
 478             for (;j<=L;j++) { // overwrite possible leftover values
 479                 _G[j] = -1.0;
 480                 _x[j] = -1;
 481             }
 482             e_L = self+e_L; // e_L is now correct
 483             m_localEnergy = e_L;
 484             //std::cout << e_L-m_energy.compute(m_state, m_helper, m_prob) << std::endl;
 485             //std::cout << e_L-m_energy.compute(m_state, m_prob) << std::endl;
 486         }
 487     public:
 488         typedef SpinChain State;
 489         typedef /*typename*/ State::Change Change;
 490         typedef VarProb Probability;
 491         typedef /*typename*/ Probability::Helper Helper;
 492
 493         explicit Walker(const Probability &p, const State &s, Energy e = Energy()) : Lambda(0.0), m_weight(1.0), m_energy(e), m_state(s), m_prob(p) {
 494             // Lambda = s.size() * 0.25 * e.J * (0.01 + e.D); // surely higher than any value in H
 495             _recompute();
 496         }
 497
 498         const State& state () const { return m_state; }
 499         void setState (const State& s) { m_state = s; _recompute(); }
 500         const Probability& probability () const { return m_prob; }
 501         void setProbability (const Probability& p) { m_prob = p; _recompute(); }
 502         const Helper& helper () const { return m_helper; }
 503         double weight () const { return m_weight; }
 504         double localEnergy () const { return m_localEnergy; }
 505
 506         bool step () {
 507             const int L = m_state.size();
 508             const double b = (Lambda - m_localEnergy); // will be the new weight
 509             double r = b*rand()/(RAND_MAX+1.0);
 510             int j = 0;
 511             while (r>_G[j] && _G[j]>0) {
 512                 //cout << j << " " << r << " " << G[j] << endl;
 513                 r -= _G[j++];
 514             }
 515             //if (_G[j]<0.0) throw std::string("WTF?!");
 516             // now j contains the index of the spin that must be flipped
 517             const int pos = _x[j];
 518             //std::cout << "stuck on " << j << " " << r << " " << _G[j] << " " << pos << std::endl;
 519             //std::cout << pos << "\n"; // to check for translation invariance
 520             if (pos == -1) {
 521                 // no change needed?
 522             } else {
 523                 // let's recompute
 524                 const Change f = { pos, (pos+1)%L };
 525                 m_state.applyChange(f);
 526                 m_prob.applyChangeHelper(m_state, m_helper, f); // XXX: leave after applyChange()
 527                 _fillG();
 528             }
 529             m_weight = b;
 530             return pos!=-1;
 531         }
 532     private:
 533         double Lambda;
 534         double m_weight;
 535         double m_localEnergy;
 536         Energy m_energy;
 537         State m_state;
 538         Helper m_helper;
 539         Probability m_prob;
 540         std::vector<double> _G; // green function for a flip
 541         std::vector<int> _x; // which spin is flipped
 542 };
 543
 544 using namespace std;
 545
 546 void print_help () {
 547 }
 548
 549 int main (int argc, char **argv) {
 550     bool print_state = false;
 551     ofstream out;
 552     srand(time(NULL));
 553     string filename("");
 554     int nspins = 100;
 555     int niter = 1e6;
 556     double alpha = 0.25;
 557     SpinChain init = SpinChain::neel(nspins);
 558     for (int i=1;i<argc;i++) {
 559         if (argv[i][0]!='-') {
 560             print_help();
 561             exit(1);
 562         }
 563         switch (argv[i][1]) {
 564             case 's':
 565                 print_state = true;
 566             case 'o':
 567                 filename = argv[++i];
 568                 break;
 569             case 'r':
 570                 srand(atoi(argv[++i]));
 571                 break;
 572             case 'i':
 573                 niter = atoi(argv[++i]);
 574                 break;
 575             case 'a':
 576                 alpha = atof(argv[++i]);
 577                 break;
 578             case 'n':
 579                 nspins = atoi(argv[++i]);
 580                 init = SpinChain::neel(nspins);
 581                 break;
 582             case 'h':
 583             case 'H':
 584                 print_help();
 585                 exit(0);
 586                 break;
 587             default:
 588                 print_help();
 589                 exit(1);
 590                 break;
 591         }
 592     }
 593
 594     if (filename!="") out.open(filename.c_str(), ofstream::out);
 595     Energy en;
 596     //QMC<VarProb> mc(VarProb(alpha), init);
 597     Walker mc(VarProb(alpha), init, en);
 598
 599     bool accepted = true; // first value is always accepted
 600     double e = 0.0;
 601     double w = 0.0;
 602     if (out.is_open()) {
 603         out.precision(10);
 604         out << "#!prepend " << alpha << "\n";
 605         out << "# alpha = " << alpha << "\n";
 606         out << "# format = '^([01]+)%s+(.*)$'\n";
 607     }
 608     std::string state_str = init.toString();
 609     for (int i=1;i<=niter;i++) {
 610         if (!accepted) {
 611             // e still contains the correct value
 612         } else {
 613             //e = en.compute(mc.state(), mc.helper(), mc.probability());
 614             e = mc.localEnergy();
 615             if (print_state) state_str = mc.state().toString();
 616         }
 617         w = mc.weight();
 618         if (out.is_open()) {
 619             if (print_state) {
 620                 out << state_str << " " << (e/nspins) << " " << w << "\n";
 621             } else {
 622                 out << "* " << (e/nspins) << " " << w << "\n";
 623             }
 624         }
 625         accepted = mc.step();
 626     }
 627     if (out.is_open()) {
 628         out.close();
 629     }
 630     return 0;
 631 }
 632