diff --git a/exercises/ex10_solution/Makefile b/exercises/ex10_solution/Makefile
new file mode 100755
index 0000000000000000000000000000000000000000..7bd55db609488df0034ed1d355941ab3ba3589e6
--- /dev/null
+++ b/exercises/ex10_solution/Makefile
@@ -0,0 +1,4 @@
+include $(ALPS_ROOT)/share/alps/include.mk
+
+ising: ising.cpp
+ $(CXX) -Wall $(CPPFLAGS) $(CXXFLAGS) $(LDFLAGS) $(LIBS) -o ising ising.cpp
diff --git a/exercises/ex10_solution/classical_ising-mag.pdf b/exercises/ex10_solution/classical_ising-mag.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..7964e9e202347ff6610101e715a9ac8f1c82fed7
Binary files /dev/null and b/exercises/ex10_solution/classical_ising-mag.pdf differ
diff --git a/exercises/ex10_solution/classical_ising-mag2.pdf b/exercises/ex10_solution/classical_ising-mag2.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..30c7457e9fa9d905afd3453a15b8f7adfb64ca20
Binary files /dev/null and b/exercises/ex10_solution/classical_ising-mag2.pdf differ
diff --git a/exercises/ex10_solution/classical_ising-mag_error.pdf b/exercises/ex10_solution/classical_ising-mag_error.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..da2ef9ce5c21acba312108676b8e10a01be5583d
Binary files /dev/null and b/exercises/ex10_solution/classical_ising-mag_error.pdf differ
diff --git a/exercises/ex10_solution/classical_ising-mag_tau.pdf b/exercises/ex10_solution/classical_ising-mag_tau.pdf
new file mode 100755
index 0000000000000000000000000000000000000000..59a4905418d4dc10e82e47d3c20b81b4f88fef50
Binary files /dev/null and b/exercises/ex10_solution/classical_ising-mag_tau.pdf differ
diff --git a/exercises/ex10_solution/data/run_classical.sh b/exercises/ex10_solution/data/run_classical.sh
new file mode 100755
index 0000000000000000000000000000000000000000..55e055b953dbd156eb7c1b9986a9d7b8d9afe77f
--- /dev/null
+++ b/exercises/ex10_solution/data/run_classical.sh
@@ -0,0 +1,11 @@
+#!/bin/sh
+
+L=16
+N=1000000
+
+for J in 0.1 0.2 0.3 0.4 0.44 0.5 0.6 0.7 0.8 0.9
+do
+ time ../ising $L $L $J $J $N &
+ time ../ising $L $L $J $J $N local &
+ wait
+done
diff --git a/exercises/ex10_solution/ising.cpp b/exercises/ex10_solution/ising.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..0b301f58b26f478536bcb0c33b79bcc02d37d4cc
--- /dev/null
+++ b/exercises/ex10_solution/ising.cpp
@@ -0,0 +1,344 @@
+/*****************************************************************************
+ *
+ * Copyright (C) 2010 - 2012 Jan Gukelberger
+ * Copyright (C) 2010 Brigitte Surer
+ *
+ * This software is based on the ALPS Code-02 tutorial, published under the ALPS
+ * Library License; you can use, redistribute it and/or modify it under
+ * the terms of the license, either version 1 or (at your option) any later
+ * version.
+ *
+ * You should have received a copy of the ALPS Library License along with
+ * the ALPS Libraries; see the file LICENSE.txt. If not, the license is also
+ * available from http://alps.comp-phys.org/.
+ *
+ * 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, TITLE AND NON-INFRINGEMENT. IN NO EVENT
+ * SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
+ * FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
+ * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *
+ *****************************************************************************/
+
+#include <alps/scheduler/montecarlo.h>
+#include <alps/alea.h>
+#include <boost/random.hpp>
+#include <boost/multi_array.hpp>
+#include <iostream>
+#include <sstream>
+#include <vector>
+#include <map>
+#include <cmath>
+using namespace std;
+
+
+class Simulation
+{
+public:
+ Simulation(double Jy, double Jx, size_t M, size_t L, std::string output_file, bool local)
+ : eng_(42)
+ , rng_(eng_, dist_)
+ , localUpdates_(local)
+ , spins_(boost::extents[M][L])
+ , energy_("E")
+ , magnetization_("m")
+ , abs_magnetization_("|m|")
+ , m2_("m^2")
+ , m4_("m^4")
+ , filename_(output_file)
+ {
+ L_[0] = M ; L_[1] = L ;
+ J_[0] = Jy; J_[1] = Jx;
+ exp2j_[0] = exp(-2*Jy);
+ exp2j_[1] = exp(-2*Jx);
+
+ // Init random spin configuration
+ for(unsigned i = 0; i < L_[0]; ++i)
+ {
+ for(unsigned j = 0; j < L_[1]; ++j)
+ spin(make_pair(i,j)) = 2 * randint(2) - 1;
+ }
+ }
+
+ void run(size_t ntherm,size_t n)
+ {
+ thermalization_ = ntherm;
+ sweeps_=n;
+ // Thermalize for ntherm steps
+ while(ntherm--)
+ step();
+
+ // Run n steps
+ while(n--)
+ {
+ step();
+ measure();
+ }
+
+ //save the observables
+ save(filename_);
+
+ // Print observables
+ // std::cout << abs_magnetization_;
+ std::cout << abs_magnetization_.name() << ":\t" << abs_magnetization_.mean()
+ << " +- " << abs_magnetization_.error() << ";\ttau = " << abs_magnetization_.tau()
+ << ";\tconverged: " << alps::convergence_to_text(abs_magnetization_.converged_errors())
+ << std::endl;
+ std::cout << energy_.name() << ":\t" << energy_.mean()
+ << " +- " << energy_.error() << ";\ttau = " << energy_.tau()
+ << ";\tconverged: " << alps::convergence_to_text(energy_.converged_errors())
+ << std::endl;
+ std::cout << magnetization_.name() << ":\t" << magnetization_.mean()
+ << " +- " << magnetization_.error() << ";\ttau = " << magnetization_.tau()
+ << ";\tconverged: " << alps::convergence_to_text(magnetization_.converged_errors())
+ << std::endl;
+ }
+
+ void step()
+ {
+ switch( localUpdates_ )
+ {
+ case false: cluster_step();
+ case true: local_sweep();
+ }
+ }
+
+ void cluster_step()
+ {
+ // printSpins();
+ vector<TodoItem> todo;
+ set<Site> flipped;
+
+ // Pick random site k=(i,j)
+ Site x = randsite();
+ Spin s = spin(x);
+ flip(x);
+ flipped.insert(x);
+ addNeighbors(x,todo);
+
+ while( !todo.empty() )
+ {
+ TodoItem t = todo.back();
+ todo.pop_back();
+
+ if( spin(t.x) != s || flipped.count(t.x) ) continue; // do not do anything if the spins where not parallel or if the potential one (t) is already added to the cluster
+
+ if( rng_() > exp2j_[t.dir] )
+ {
+ flip(t.x);
+ flipped.insert(t.x);
+ addNeighbors(t.x,todo);
+ }
+ }
+ }
+
+ void local_sweep()
+ {
+ for( size_t i = 0; i < nsites(); ++i )
+ {
+ // Pick random site k=(i,j)
+ Site x = randsite();
+ Spin s = spin(x);
+
+ // Measure local energy e = -s_k * sum_{l nn k} s_l
+ Site l(x.first,(x.second+L_[1]-1)%L_[1]), d((x.first+L_[0]-1)%L_[0],x.second),
+ r(x.first,(x.second +1)%L_[1]), u((x.first +1)%L_[0],x.second);
+ int ey = s * (spin(u) + spin(d)); // up + down
+ int ex = s * (spin(l) + spin(r)); // left + right
+ double p = 1;
+ switch( ey )
+ {
+ case -2: p /= (exp2j_[0]*exp2j_[0]); break;
+ case +2: p *= (exp2j_[0]*exp2j_[0]); break;
+ default: break;
+ }
+ switch( ex )
+ {
+ case -2: p /= (exp2j_[1]*exp2j_[1]); break;
+ case +2: p *= (exp2j_[1]*exp2j_[1]); break;
+ default: break;
+ }
+
+ // Flip s_k with probability exp(2 beta e)
+ if(p > 1 || rng_() < p)
+ flip(x);
+ }
+ }
+
+ void measure()
+ {
+ int E = 0; // energy
+ int M = 0; // magnetization
+ for(size_t i = 0; i < L_[0]; ++i)
+ {
+ for(size_t j = 0; j < L_[1]; ++j)
+ {
+ E -= J_[0]*spins_[i][j]*spins_[(i+1)%L_[0]][j] + J_[1]*spins_[i][j]*spins_[i][(j+1)%L_[1]];
+ M += spins_[i][j];
+ }
+ }
+
+ // Add sample to observables
+ energy_ << E/double(nsites());
+ double m = M/double(nsites());
+ magnetization_ << m;
+ abs_magnetization_ << fabs(m);
+ m2_ << m*m;
+ m4_ << m*m*m*m;
+ }
+
+ void save(std::string const & filename){
+ alps::hdf5::archive ar(filename,alps::hdf5::archive::REPLACE);
+ ar << alps::make_pvp("/parameters/L", L_[1]);
+ ar << alps::make_pvp("/parameters/M", L_[0]);
+ ar << alps::make_pvp("/parameters/Jx", J_[1]);
+ ar << alps::make_pvp("/parameters/Jy", J_[0]);
+ ar << alps::make_pvp("/parameters/SWEEPS", sweeps_);
+ ar << alps::make_pvp("/parameters/THERMALIZATION", thermalization_);
+
+ ar << alps::make_pvp("/simulation/results/"+energy_.representation(), energy_);
+ ar << alps::make_pvp("/simulation/results/"+magnetization_.representation(), magnetization_);
+ ar << alps::make_pvp("/simulation/results/"+abs_magnetization_.representation(), abs_magnetization_);
+ ar << alps::make_pvp("/simulation/results/"+m2_.representation(), m2_);
+ ar << alps::make_pvp("/simulation/results/"+m4_.representation(), m4_);
+
+ alps::RealObsevaluator m = abs_magnetization_;
+ alps::RealObsevaluator m2 = m2_;
+ alps::RealObsevaluator m4 = m4_;
+
+ alps::RealObsevaluator u2("Binder Cumulant U2");
+ u2 = m2/(m*m);
+ ar << alps::make_pvp("/simulation/results/"+u2.name(), u2);
+
+ alps::RealObsevaluator chi("Connected Susceptibility");
+ chi=nsites()*(m2-m*m);
+ ar << alps::make_pvp("/simulation/results/"+chi.name(), chi);
+
+ alps::RealObsevaluator u4("Binder Cumulant");
+ u4 = m4/(m2*m2);
+ ar << alps::make_pvp("/simulation/results/"+u4.name(), u4);
+ }
+
+protected:
+ typedef signed char Spin;
+ typedef std::pair<unsigned,unsigned> Site;
+ typedef unsigned Direction;
+ struct TodoItem { Site x; Direction dir; TodoItem(unsigned i=0, unsigned j=0, Direction d=0) : x(i,j), dir(d) {} };
+
+ // Random int from the interval [0,max)
+ int randint(int max) const
+ {
+ return static_cast<int>(max * rng_());
+ }
+ Site randsite() const
+ {
+ return Site(randint(L_[0]),randint(L_[1]));
+ }
+ size_t nsites() const
+ {
+ return L_[0]*L_[1];
+ }
+ Spin& spin(Site x)
+ {
+ return spins_[x.first][x.second];
+ }
+ const Spin& spin(Site x) const
+ {
+ return spins_[x.first][x.second];
+ }
+ void flip(Site x)
+ {
+ spin(x) *= -1;
+ }
+ void addNeighbors(Site x, vector<TodoItem>& n) const
+ {
+ n.push_back(TodoItem(x.first,(x.second+L_[1]-1)%L_[1],1)); // left
+ n.push_back(TodoItem(x.first,(x.second +1)%L_[1],1)); // right
+ n.push_back(TodoItem((x.first+L_[0]-1)%L_[0],x.second,0)); // down
+ n.push_back(TodoItem((x.first +1)%L_[0],x.second,0)); // up
+ }
+
+ void printSpins() const
+ {
+ for( unsigned i = 0; i < L_[0]; ++i )
+ {
+ for( unsigned j = 0; j < L_[1]; ++j )
+ cout << (int(spins_[i][j]) > 0 ? 'x' : '.');
+ cout << "\n";
+ }
+ cout << endl;
+ }
+
+private:
+ typedef boost::mt19937 engine_type;
+ typedef boost::uniform_real<> distribution_type;
+ typedef boost::variate_generator<engine_type&, distribution_type> rng_type;
+ engine_type eng_;
+ distribution_type dist_;
+ mutable rng_type rng_;
+
+ static const unsigned DIMENSIONS = 2;
+ size_t L_[DIMENSIONS];
+ double J_[DIMENSIONS];
+ size_t sweeps_;
+ size_t thermalization_;
+ bool localUpdates_;
+
+ boost::multi_array<Spin,DIMENSIONS> spins_;
+ double exp2j_[DIMENSIONS];
+
+ alps::RealObservable energy_;
+ alps::RealObservable magnetization_;
+ alps::RealObservable abs_magnetization_;
+ alps::RealObservable m2_;
+ alps::RealObservable m4_;
+
+ std::string filename_;
+};
+
+string usage(const std::string& prog)
+{
+ return "usage: " + prog + " M L Jy Jx N";
+}
+int main(int argc, const char** argv)
+{
+ if( argc != 6 && argc != 7 )
+ {
+ cerr << usage(argv[0]) << endl;
+ return -1;
+ }
+
+ size_t M, L, N;
+ double Jy, Jx;
+ bool local = false;
+ try
+ {
+ using boost::lexical_cast;
+ M = lexical_cast<size_t>(argv[1]); // Linear lattice sizes
+ L = lexical_cast<size_t>(argv[2]); //
+ Jy = lexical_cast<double>(argv[3]); // coupling constants
+ Jx = lexical_cast<double>(argv[4]); //
+ N = lexical_cast<size_t>(argv[5]); // # of simulation steps
+ if( argc > 6 && argv[6] == string("local") )
+ local = true;
+ }
+ catch(...)
+ {
+ cerr << usage(argv[0]) << endl;
+ return -1;
+ }
+
+
+ std::cout << "# M: " << M << " L: " << L << " Jy: " << Jy << " Jx: " << Jx << " N: " << N << " local: " << local << std::endl;
+
+ std::stringstream output_name;
+ output_name << "ising_M" << M << "_L" << L << "_Jy" << Jy << "_Jx" << Jx;
+ if( local ) output_name << "_local";
+ output_name <<".h5";
+ Simulation sim(Jy, Jx, M, L, output_name.str(), local);
+ sim.run(N/5,N);
+
+ return 0;
+}
diff --git a/exercises/ex10_solution/plot_classical.py b/exercises/ex10_solution/plot_classical.py
new file mode 100755
index 0000000000000000000000000000000000000000..1c3f9695ece4d55922f6f4c75abebd31307cd2ef
--- /dev/null
+++ b/exercises/ex10_solution/plot_classical.py
@@ -0,0 +1,53 @@
+import sys, os, re
+
+import numpy as np
+import matplotlib.pyplot as plt
+import pyalps
+from pyalps.plot import plot
+
+files = pyalps.getResultFiles(dirname='data')
+data = pyalps.loadMeasurements(files , ['|m|','m^2'])
+
+for d in pyalps.flatten(data):
+ if re.search(r'_local\.h5$',d.props['filename']):
+ d.props['update'] = 'local'
+ else:
+ d.props['update'] = 'cluster'
+
+
+m = pyalps.collectXY(data, 'Jx', '|m|', foreach=['update'])
+m2 = pyalps.collectXY(data, 'Jx', 'm^2', foreach=['update'])
+
+plt.figure()
+plot(m)
+plt.xlabel('$\\beta J$')
+plt.ylabel('|magnetization|')
+plt.legend(loc='best', frameon=False)
+plt.savefig('classical_ising-mag.pdf',transparent=True)
+
+
+plt.figure()
+plot(m2)
+plt.xlabel('$\\beta J$')
+plt.ylabel('|magnetization|^2')
+plt.legend(loc='best', frameon=False)
+plt.savefig('classical_ising-mag2.pdf',transparent=True)
+
+plt.figure()
+for d in m:
+ plt.plot(d.x,[y.error for y in d.y],label=d.props['label'])
+plt.xlabel('$\\beta J$')
+plt.ylabel('error')
+plt.legend(loc='best', frameon=False)
+plt.savefig('classical_ising-mag_error.pdf',transparent=True)
+
+plt.figure()
+for d in m:
+ plt.plot(d.x,[y.tau for y in d.y],label=d.props['label'])
+plt.xlabel('$\\beta J$')
+plt.ylabel('autocorrelation time')
+plt.legend(loc='best', frameon=False)
+plt.savefig('classical_ising-mag_tau.pdf',transparent=True)
+
+
+plt.show()