From a735a3eb295bca0e20229db557f9de9eadde2b40 Mon Sep 17 00:00:00 2001
From: jake <jake@keplr.io>
Date: Tue, 19 Sep 2017 14:03:11 -0400
Subject: [PATCH] style: delete trailing spaces

---
 main.py                   | 70 ++++++++++++++---------------
 rbm.py                    | 12 ++---
 tutorial/ed_tfim1d.py     | 60 ++++++++++++-------------
 tutorial/ising2d.py       | 74 +++++++++++++++----------------
 tutorial/main.py          |  4 +-
 tutorial/plot_observer.py | 70 ++++++++++++++---------------
 tutorial/plot_results.py  | 30 ++++++-------
 tutorial/plot_sampler.py  | 52 +++++++++++-----------
 tutorial/rbm.py           | 16 +++----
 tutorial/sampler.py       | 84 +++++++++++++++++------------------
 tutorial/tfim1d.py        | 92 +++++++++++++++++++--------------------
 11 files changed, 282 insertions(+), 282 deletions(-)

diff --git a/main.py b/main.py
index 6e04593..6b0b6c8 100644
--- a/main.py
+++ b/main.py
@@ -9,26 +9,26 @@
 import json
 
 def main()
-    
+
     # Initialize the command line parser
     parser = argparse.ArgumentParser()
-    
+
     # Read command line arguments
-    parser.add_argument('command',type=str,help='command to execute') 
-    parser.add_argument('-nV',type=int,default=4,help='number of visible nodes')                
-    parser.add_argument('-nH',type=int,default=4,help='number of hidden nodes')   
-    parser.add_argument('-steps',type=int,default=1000000,help='training steps')  
-    parser.add_argument('-lr',type=float,default=1e-3,help='learning rate')   
-    parser.add_argument('-bs',type=int,default=100,help='batch size')   
-    parser.add_argument('-CD',type=int,default=10,help='steps of contrastive divergence') 
-    parser.add_argument('-nC',type=float,default=10,help='number of chains in PCD')  
-    
+    parser.add_argument('command',type=str,help='command to execute')
+    parser.add_argument('-nV',type=int,default=4,help='number of visible nodes')
+    parser.add_argument('-nH',type=int,default=4,help='number of hidden nodes')
+    parser.add_argument('-steps',type=int,default=1000000,help='training steps')
+    parser.add_argument('-lr',type=float,default=1e-3,help='learning rate')
+    parser.add_argument('-bs',type=int,default=100,help='batch size')
+    parser.add_argument('-CD',type=int,default=10,help='steps of contrastive divergence')
+    parser.add_argument('-nC',type=float,default=10,help='number of chains in PCD')
+
     # Parse the arguments
     args = parser.parse_args()
-    
+
     if args.command == 'train':
         train(args)
-    
+
     if args.command == 'sample':
         sample(args)
 
@@ -43,7 +43,7 @@ class Ops(object):
     pass
 
 def train(args):
-   
+
     # Simulation parameters
     num_visible = args.nV           # number of visible nodes
     num_hidden = args.nH            # number of hidden nodes
@@ -57,32 +57,32 @@ def train(args):
     hidden_bias=None                # hidden bias
     bcount=0                        # counter
     epochs_done=1                   # epochs counter
- 
+
     # *************************************************************
     # INSERT HERE THE PATH TO THE TRAINING AND TESTING DATASETS
     trainName = '*******************'
     testName  = '*******************'
-    
+
     # Loading the data
     xtrain = np.loadtxt(trainName)
     xtest = np.loadtxt(testName)
-    
+
     ept=np.random.permutation(xtrain)               # random permutation of training data
     epv=np.random.permutation(xtest)                # random permutation of test data
     iterations_per_epoch = xtrain.shape[0] / bsize  # gradient iteration per epoch
 
     # Initialize RBM class
-    rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples) 
-    
+    rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples)
+
     # Initialize operations and placeholders classes
     ops = Ops()
     placeholders = Placeholders()
 
     placeholders.visible_samples = tf.placeholder(tf.float32, shape=(None, num_visible), name='v') # placeholder for training data
 
-    total_iterations = 0 # starts at zero 
+    total_iterations = 0 # starts at zero
     ops.global_step = tf.Variable(total_iterations, name='global_step_count', trainable=False)
-    
+
     # Decaying learning rate
     learning_rate = tf.train.exponential_decay(
         learning_rate_b,
@@ -98,10 +98,10 @@ def train(args):
     ops.lr=learning_rate
     ops.train = optimizer.minimize(cost, global_step=ops.global_step)
     ops.init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
-    
+
     with tf.Session() as sess:
         sess.run(ops.init)
-        
+
         for ii in range(nsteps):
             if bcount*bsize+ bsize>=xtrain.shape[0]:
                bcount=0
@@ -110,7 +110,7 @@ def train(args):
             batch=ept[ bcount*bsize: bcount*bsize+ bsize,:]
             bcount=bcount+1
             feed_dict = {placeholders.visible_samples: batch}
-            
+
             _, num_steps = sess.run([ops.train, ops.global_step], feed_dict=feed_dict)
 
             if num_steps % iterations_per_epoch == 0:
@@ -119,26 +119,26 @@ def train(args):
                 epochs_done += 1
 
 def sample(args):
-       
+
     num_visible = args.nV   # number of visible nodes
     num_hidden = args.nH    # number of hidden nodes
-    
+
     # *************************************************************
     # INSERT HERE THE PATH TO THE PARAMETERS FILE
     path_to_params = '*******************'
-   
-    # Load the RBM parameters 
+
+    # Load the RBM parameters
     params = np.load(path_to_params)
     weights = params['weights']
     visible_bias = params['visible_bias']
     hidden_bias = params['hidden_bias']
     hidden_bias=np.reshape(hidden_bias,(hidden_bias.shape[0],1))
     visible_bias=np.reshape(visible_bias,(visible_bias.shape[0],1))
-  
+
     # Sampling parameters
     num_samples=1000   # how many independent chains will be sampled
     gibb_updates=100   # how many gibbs updates per call to the gibbs sampler
-    nbins=1000         # number of calls to the RBM sampler      
+    nbins=1000         # number of calls to the RBM sampler
 
     # Initialize RBM class
     rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples)
@@ -146,24 +146,24 @@ def sample(args):
 
     # Initialize tensorflow
     init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
-   
+
     with tf.Session() as sess:
         sess.run(init)
-    
+
         for i in range(nbins):
             print ('bin %d\t' %i)
-            
+
             # Gibbs sampling
             _,samples=sess.run([hsamples,vsamples])
 
 
 def save_parameters(sess,rbm,epochs):
     weights, visible_bias, hidden_bias = sess.run([rbm.weights, rbm.visible_bias, rbm.hidden_bias])
-    
+
     # *************************************************************
     # INSERT HERE THE PATH TO THE PARAMETERS FILE
     parameter_file_path = '*******************'
-    
+
     np.savez_compressed(parameter_file_path, weights=weights, visible_bias=visible_bias, hidden_bias=hidden_bias,
                         epochs=epochs)
 
diff --git a/rbm.py b/rbm.py
index 4ed04ac..58c102e 100644
--- a/rbm.py
+++ b/rbm.py
@@ -3,9 +3,9 @@
 import numpy as np
 
 class RBM(object):
-    
+
     ''' Restricted Boltzmann Machine '''
-    
+
     def __init__(self, num_hidden, num_visible, num_samples=128, weights=None, visible_bias=None, hidden_bias=None):
         ''' Constructor '''
         # number of hidden units
@@ -101,7 +101,7 @@ class variables.
         self.hidden_samples = self.hidden_samples.assign(h_samples)
         self.p_of_v = p_of_v
         return self.hidden_samples, v_samples
-    
+
     def energy(self, hidden_samples, visible_samples):
         # type: (tf.Tensor, tf.Tensor) -> tf.Tensor
         """Compute the energy:
@@ -123,8 +123,8 @@ def free_energy(self, visible_samples):
         free_energy = (tf.matmul(visible_samples, self.visible_bias)
                        + tf.reduce_sum(tf.nn.softplus(tf.matmul(visible_samples, self.weights)
                                                       + tf.transpose(self.hidden_bias)), 1, keep_dims=True))
-        return free_energy   
-    
+        return free_energy
+
     def neg_log_likelihood_grad(self, visible_samples, model_samples=None, num_gibbs=2):
         # type: (tf.Tensor, tf.Tensor, int) -> tf.Tensor
 
@@ -142,7 +142,7 @@ def neg_log_likelihood(self, visible_samples, log_Z):
                        + tf.reduce_sum(tf.nn.softplus(tf.matmul(visible_samples, self.weights)
                                                       + tf.transpose(self.hidden_bias)), 1))
         return -tf.reduce_mean(free_energy - log_Z)
-    
+
     def exact_log_partition_function(self):
         ''' Evaluate the partition function by exact enumerations '''
         with tf.name_scope('exact_log_Z'):
diff --git a/tutorial/ed_tfim1d.py b/tutorial/ed_tfim1d.py
index a3bc1cd..600009f 100755
--- a/tutorial/ed_tfim1d.py
+++ b/tutorial/ed_tfim1d.py
@@ -6,7 +6,7 @@
 Sx = np.array([[0,1.],[1,0]])
 Sz = np.array([[1,0.],[0,-1]])
 
-# Pauli X operator 
+# Pauli X operator
 def sigmaX(L,i):
     ''' Return the many-body operator
         I x I x .. x Sx x I x .. x I
@@ -34,11 +34,11 @@ def sigmaZ(L,i):
 
     return reduce(np.kron,OpList)
 
-# Magnetic interaction term 
+# Magnetic interaction term
 def buildMagneticInteraction(i,B,L):
     return B*sigmaX(L,i)
 
-# Ising interaction term 
+# Ising interaction term
 def buildIsingInteraction(i,j,J,L):
     ''' Return the Ising interaction term
         I x .. x Sz x Sz x .. x I
@@ -56,19 +56,19 @@ def buildIsingInteraction(i,j,J,L):
 
 # Build transverse-field Ising model
 def build1dIsingModel(L,J,B,OBC):
-    
+
     D = 1<<L    # dimension of Hilbert space
     ''' Return the full Hamiltonian '''
-    Ham = np.zeros((D,D))	
+    Ham = np.zeros((D,D))
     for i in range(L-1):
 	Ham = Ham - buildIsingInteraction(i,i+1,J,L)
-	Ham = Ham - buildMagneticInteraction(i,B,L)	
+	Ham = Ham - buildMagneticInteraction(i,B,L)
     Ham = Ham - buildMagneticInteraction(L-1,B,L)
-	
+
     # Periodic Boundary Conditions
     if OBC is True:
         Ham = Ham - buildIsingInteraction(0,L-1,J,L)
-    
+
     return Ham
 
 # Generate the training dataset
@@ -76,17 +76,17 @@ def buildDataset(L,B,psi,Nsamples):
 
     D = 1<<L
     config = np.zeros((D,L))
-    psi2 = np.zeros((D))  
-    
+    psi2 = np.zeros((D))
+
     # Build all spin states and psi^2
     for i in range(D):
         state = (bin(i)[2:].zfill(L)).split()
         for j in range(L):
             config[i,j] = int(state[0][j])
-        psi2[i] = psi[i]**2 
-    
+        psi2[i] = psi[i]**2
+
     config_index = range(D)
-    
+
     # Generate the trainset
     index_samples = np.random.choice(config_index,
                                      Nsamples,
@@ -99,7 +99,7 @@ def buildDataset(L,B,psi,Nsamples):
             dataFile.write("%d " % config[index_samples[i]][j])
         dataFile.write("\n")
     dataFile.close()
-    
+
     # Generate the testset
     index_samples = np.random.choice(config_index,
                                      Nsamples/10,
@@ -114,16 +114,16 @@ def buildDataset(L,B,psi,Nsamples):
 
 # MAIN
 def main(pars):
-	
+
     print ('\n--------------------------------\n')
     print ('EXACT DIAGONALIZATION OF THE 1d-TFIM\n')
-     
+
     # Parameters
     L = pars.L;         # number of spins
-    B = pars.B          # magnetic field 
+    B = pars.B          # magnetic field
     J = pars.J          # ising interaction strength
     Nsamples = 10000    # train dataset size
-    
+
     print ('Number of spins    L = %d' % L)
     print ('Ising interaction  J = %.1f' % J)
     print ('\n\n')
@@ -133,38 +133,38 @@ def main(pars):
     Bmax = 2.0
     Bsteps = 10
     deltaB = (Bmax-Bmin)/float(Bsteps)
-    
+
     # Observables file
     obs_Name = '../data/tfim1d/observables/ed_tfim1d_L' + str(L) + '_observables.txt'
     obs_File = open(obs_Name,'w')
     obs_File.write('# B               E         ')  # write header
     obs_File.write('<|Sz|>           <Sx>\n')       # write header
-    
+
     # Spin-spin correlation file
     corr_Name = '../data/tfim1d/observables/ed_tfim1d_L' + str(L) + '_correlations.txt'
     corr_File = open(corr_Name,'w')
-    
+
     # Loop over magnetic field values
     for b in range(1,Bsteps+1):
-        
+
         B = Bmin + b*deltaB
-        print('Magnetic field B = %.2f' % B) 
+        print('Magnetic field B = %.2f' % B)
 
         # Wavefunction file
         psiName  = '../data/tfim1d/wavefunctions/wavefunction_tfim1d_L' + str(L)
         psiName += '_B' + str(B) + '.txt'
-        
+
         # Diagonalize the Hamiltonian
         print('diagonalizing...')
         H = build1dIsingModel(L,J,B,True)
         (e,psi) = np.linalg.eigh(H)
         psi0 = np.abs(psi[:,0])
         e0 = e[0]
-        
+
         # Save energy and wavefunction
         obs_File.write('%.1f   %.10f   ' % (B,e0/float(L)))
         np.savetxt(psiName,psi0)
-        
+
         # Magnetic observables
         print('computing observables...')
         # Compute <|Sz|>
@@ -184,13 +184,13 @@ def main(pars):
             for j in range(L):
                 for k in range(1<<L):
                     config = (bin(k)[2:].zfill(L)).split()
-                    SzSz[i,j] += (psi0[k]**2)*(2*float(config[0][i])-1)*(2*float(config[0][j])-1) 
+                    SzSz[i,j] += (psi0[k]**2)*(2*float(config[0][i])-1)*(2*float(config[0][j])-1)
                 corr_File.write('%.10f   ' % SzSz[i,j])
             corr_File.write('\n')
-        
+
         # Generate the training datasets
         print('generating the dataset...')
-        buildDataset(L,B,psi0,Nsamples) 
+        buildDataset(L,B,psi0,Nsamples)
         print('\n')
 
 
@@ -198,7 +198,7 @@ def main(pars):
 
 
 if __name__ == "__main__":
-   
+
     # Read arguments from command line
     parser = argparse.ArgumentParser()
     parser.add_argument('-L',type=int,default=10)
diff --git a/tutorial/ising2d.py b/tutorial/ising2d.py
index dfc6a92..f2622c5 100644
--- a/tutorial/ising2d.py
+++ b/tutorial/ising2d.py
@@ -15,7 +15,7 @@ def save_parameters(sess, results_dir, rbm, epochs,L,T):
     parameter_file_path =  '../data/ising2d/parameters/parameters_L' + str(L)
     parameter_file_path += '_T' + str(T)
     np.savez_compressed(parameter_file_path, weights=weights, visible_bias=visible_bias, hidden_bias=hidden_bias,
-                        epochs=epochs) 
+                        epochs=epochs)
 class Args(object):
     pass
 
@@ -26,7 +26,7 @@ class Ops(object):
     pass
 
 def train(args):
-   
+
     # Simulation parameters
     T = args.T                      # temperature
     num_visible = args.L*args.L     # number of visible nodes
@@ -41,27 +41,27 @@ def train(args):
     hidden_bias=None                # hidden bias
     bcount=0                        # counter
     epochs_done=1                   # epochs counter
- 
+
     # Loading the data
     train_dir = '../data/ising2d/datasets/'  # Location of training data.
     trainName = '../data/ising2d/datasets/ising2d_L'+str(args.L)+'_T'+str(T)+'_train.txt'
     testName = '../data/ising2d/datasets/ising2d_L'+str(args.L)+'_T'+str(T)+'_test.txt'
     xtrain = np.loadtxt(trainName)
     xtest = np.loadtxt(testName)
-    
+
     ept=np.random.permutation(xtrain) # random permutation of training data
     epv=np.random.permutation(xtest) # random permutation of test data
-    iterations_per_epoch = xtrain.shape[0] / bsize  
+    iterations_per_epoch = xtrain.shape[0] / bsize
 
     # Initialize RBM class
-    rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples) 
-    
+    rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples)
+
     # Initialize operations and placeholders classes
     ops = Ops()
     placeholders = Placeholders()
     placeholders.visible_samples = tf.placeholder(tf.float32, shape=(None, num_visible), name='v') # placeholder for training data
 
-    total_iterations = 0 # starts at zero 
+    total_iterations = 0 # starts at zero
     ops.global_step = tf.Variable(total_iterations, name='global_step_count', trainable=False)
     learning_rate = tf.train.exponential_decay(
         learning_rate_b,
@@ -69,20 +69,20 @@ def train(args):
         100 * xtrain.shape[0]/bsize,
         1.0 # decay rate =1 means no decay
     )
-    
+
     cost = rbm.neg_log_likelihood_grad(placeholders.visible_samples, num_gibbs=num_gibbs)
     optimizer = tf.train.AdamOptimizer(learning_rate, epsilon=1e-2)
     ops.lr=learning_rate
     ops.train = optimizer.minimize(cost, global_step=ops.global_step)
     ops.init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
     logZ = rbm.exact_log_partition_function()
-    placeholders.logZ = tf.placeholder(tf.float32) 
+    placeholders.logZ = tf.placeholder(tf.float32)
     NLL = rbm.neg_log_likelihood(placeholders.visible_samples,placeholders.logZ)
     path_to_distr = '../data/ising2d/boltzmann_distributions/distribution_ising2d_L4_T'+str(T)+'.txt'
     boltz_distr=np.loadtxt(path_to_distr)
     p_x = tf.exp(rbm.free_energy(placeholders.visible_samples))
     all_v_states= np.array(list(it.product([0, 1], repeat=num_visible)), dtype=np.float32)
-    
+
     # Observer file
     observer_file=open('../data/ising2d/observables/training_observer.txt','w',0)
     observer_file.write('#      O')
@@ -93,7 +93,7 @@ def train(args):
     gibb_updates=10
     observer_samples=rbm.observer_sampling(gibb_updates)
     nbins=100
-    
+
     # Build lattice
     path_to_lattice = '../data/ising2d/lattice2d_L'+str(args.L)+'.txt'
     nn=np.loadtxt(path_to_lattice)
@@ -101,10 +101,10 @@ def train(args):
     N = num_visible
     E=0.0
     M=0.0
-    
+
     with tf.Session() as sess:
         sess.run(ops.init)
-        
+
         for ii in range(nsteps):
             if bcount*bsize+ bsize>=xtrain.shape[0]:
                bcount=0
@@ -113,7 +113,7 @@ def train(args):
             batch=ept[ bcount*bsize: bcount*bsize+ bsize,:]
             bcount=bcount+1
             feed_dict = {placeholders.visible_samples: batch}
-            
+
             _, num_steps = sess.run([ops.train, ops.global_step], feed_dict=feed_dict)
 
             if num_steps % iterations_per_epoch == 0:
@@ -121,24 +121,24 @@ def train(args):
                 lz = sess.run(logZ)
                 nll = sess.run(NLL,feed_dict={placeholders.visible_samples: epv, placeholders.logZ: lz})
                 px = sess.run(p_x,feed_dict={placeholders.visible_samples: all_v_states})
-                
+
                 Ov = 0.0
                 E = 0.0
                 M = 0.0
-                
+
                 for i in range(1<<num_visible):
                     Ov += boltz_distr[i]*m.log(boltz_distr[i])
                     Ov += -boltz_distr[i]*(m.log(px[i])-lz)
-                
+
                 for i in range(nbins):
-                    
+
                     # Gibbs sampling
                     samples=sess.run(observer_samples)
                     spins = np.asarray((2*samples-1))
-                    
+
                     # Compute averages of magnetizations
                     m_avg = np.mean(np.absolute(np.sum(spins,axis=1)))
-                    
+
                     # Compute averages of energies
                     e.fill(0.0)
                     for k in range(num_samples):
@@ -147,35 +147,35 @@ def train(args):
                     e_avg = np.mean(e)
                     E += (e_avg-E)/float(i+1)
                     M += (m_avg-M)/float(i+1)
-              
+
                 # Print observer on screen
                 print ('Ov = %.6f     ' % Ov,end='')
                 print ('NLL = %.6f     ' % nll,end='')
                 print ('<E> = %.6f     ' % (E/float(N)),end='')
                 print ('<|M|> = %.6f     ' % (M/float(N)),end='')
-                
+
                 # Save observer on file
                 observer_file.write('%.6f   ' % Ov)
                 observer_file.write('%.6f   ' % nll)
                 observer_file.write('%.6f   ' % (E/float(N)))
                 observer_file.write('%.6f   ' % (M/float(N)))
                 observer_file.write('\n')
- 
+
                 print()
                 #save_parameters(sess, rbm)
                 epochs_done += 1
 
 
 def sample(args):
-       
+
     T = args.T                      # temperature
-    num_visible = args.L*args.L     # number of visible units 
+    num_visible = args.L*args.L     # number of visible units
     num_hidden = args.nH            # number of hidden units
-    
+
     # Build lattice
     path_to_lattice = '../data/ising2d/lattice2d_L'+str(args.L)+'.txt'
     nn=np.loadtxt(path_to_lattice)
-    
+
     # Load the RBM parameters
     path_to_params = '../data/ising2d/parameters/parameters_nH'+str(num_hidden) + '_L'+str(args.L)+'_T'+str(T)+'.npz'
     params = np.load(path_to_params)
@@ -184,11 +184,11 @@ def sample(args):
     hidden_bias = params['hidden_bias']
     hidden_bias=np.reshape(hidden_bias,(hidden_bias.shape[0],1))
     visible_bias=np.reshape(visible_bias,(visible_bias.shape[0],1))
-  
+
     # Sampling parameters
     num_samples=1000   # how many independent chains will be sampled
     gibb_updates=10    # how many gibbs updates per call to the gibbs sampler
-    nbins=1000         # number of calls to the RBM sampler      
+    nbins=1000         # number of calls to the RBM sampler
 
     # Initialize RBM class
     rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples)
@@ -196,25 +196,25 @@ def sample(args):
 
     # Initialize tensorflow
     init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
-   
+
     N = num_visible
     e = np.zeros((num_samples))
     e2= np.zeros((num_samples))
- 
+
     with tf.Session() as sess:
         sess.run(init)
-    
+
         for i in range(nbins):
             print ('bin %d\t\t' %i,end='')
-            
+
             # Gibbs sampling
             _,samples=sess.run([hsamples,vsamples])
             spins = np.asarray((2*samples-1))
-            
+
             # Compute averages of magnetizations
             m_avg = np.mean(np.absolute(np.sum(spins,axis=1)))
             m2_avg = np.mean(np.multiply(np.sum(spins,axis=1),np.sum(spins,axis=1)))
-            
+
             # Compute averages of energies
             e.fill(0.0)
             e2.fill(0.0)
@@ -224,7 +224,7 @@ def sample(args):
                 e2[k] = e[k]*e[k]
             e_avg = np.mean(e)
             e2_avg= np.mean(e2)
-            
+
             # Compute specific heat and susceptibility
             c = (e2_avg-e_avg*e_avg)
             s = (m2_avg-m_avg*m_avg)
diff --git a/tutorial/main.py b/tutorial/main.py
index 8011497..738247d 100644
--- a/tutorial/main.py
+++ b/tutorial/main.py
@@ -1,4 +1,4 @@
-import ising2d 
+import ising2d
 import tfim1d
 import sampler
 import argparse
@@ -33,7 +33,7 @@
         ising2d.sample(args)
     if args.model == 'tfim1d':
         tfim1d.sample(args)
-    
+
 if args.command == 'sample_all':
     if args.model == 'ising2d':
         sampler.ising2d()
diff --git a/tutorial/plot_observer.py b/tutorial/plot_observer.py
index af01ca9..00fadcf 100644
--- a/tutorial/plot_observer.py
+++ b/tutorial/plot_observer.py
@@ -8,72 +8,72 @@ def observe_ising2d(args):
            2.269,2.524,2.778,3.032,3.286,3.540]
 
     plt.ion()
-    
+
     L=4
     T=args.T
     nH = args.nH
-    
+
     for i in range(len(temps)):
         if temps[i] == T:
             t_index=i
             break
-    
+
     plt.figure(figsize=(12,9), facecolor='w', edgecolor='k')
     lw = 1
     lw_exact = 2
     mS = 6
- 
+
     nameMC  = '../data/ising2d/observables/MC_ising2d_L'
     nameMC += str(L)
     nameMC += '_Observables.txt'
     fileMC = open(nameMC,'r')
     header = fileMC.readline().lstrip('#').split()
     dataMC = np.loadtxt(fileMC)
- 
+
     while(True):
-         
+
         plt.clf()
         # Open MC data file
         exact_energy = dataMC[t_index,header.index('E')]
         exact_magnetization = dataMC[t_index,header.index('M')]
-        
+
         observer = np.loadtxt('../data/ising2d/observables/training_observer.txt')
         x = [i for i in range(observer.shape[0])]
-     
+
         plt.subplot(221)
         plt.ylabel('$KL$',fontsize=25)
         plt.plot(x,observer[:,0],color='red',marker='o',markersize=mS,linewidth=lw)
-        
+
         plt.subplot(222)
         plt.ylabel('$<NLL>$',fontsize=25)
         plt.plot(x,observer[:,1],color='red',marker='o',markersize=mS,linewidth=lw)
-     
+
         plt.subplot(223)
         plt.ylabel('$<E>$',fontsize=25)
         plt.plot(x,observer[:,2],color='red',marker='o',markersize=mS,linewidth=lw)
-        plt.axhline(y=exact_energy, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact') 
-        
+        plt.axhline(y=exact_energy, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact')
+
         plt.subplot(224)
         plt.ylabel('$<M>$',fontsize=25)
         plt.plot(x,observer[:,3],color='red',marker='o',markersize=mS,linewidth=lw)
-        plt.axhline(y=exact_magnetization, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact') 
-        
+        plt.axhline(y=exact_magnetization, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact')
+
         plt.tight_layout()
         plt.pause(0.05)
-    
+
     #plt.show()
 
 def observe_tfim1d(args):
-    
+
     fields = [0.2,0.4,0.6,0.8,1.0,
               1.2,1.4,1.6,1.8,2.0]
 
     plt.ion()
-    
+
     L=10
     B=args.B
     nH = 10
-    
+
     for i in range(len(fields)):
         if fields[i] == B:
             b_index=i
@@ -87,7 +87,7 @@ def observe_tfim1d(args):
     dataMC = np.loadtxt(fileMC)
     exact_Sz = dataMC[b_index,header.index('<|Sz|>')]
     exact_Sx = dataMC[b_index,header.index('<Sx>')]
-    exact_E  = dataMC[b_index,header.index('E')] 
+    exact_E  = dataMC[b_index,header.index('E')]
     #nameCorr = '../data/tfim1d/observables/ed_tfim1d_L'+str(L)
     #nameCorr += '_correlations.txt'
     #fileCorr = open(nameCorr,'r')
@@ -97,7 +97,7 @@ def observe_tfim1d(args):
     #for i in range(L):
     #    for j in range(L):
     #        ed_corr_crit[i][j] = ed_corr_full[b_index*L+i][j]
- 
+
     plt.figure(figsize=(12,9), facecolor='w', edgecolor='k')
     lw = 1
     lw_exact = 2
@@ -106,13 +106,13 @@ def observe_tfim1d(args):
 
 
     while(True):
-         
+
         plt.clf()
         observer = np.loadtxt('../data/tfim1d/observables/training_observer.txt')
-        #correlations= np.loadtxt('../data/tfim1d/observables/training_observer_corr.txt') 
- 
+        #correlations= np.loadtxt('../data/tfim1d/observables/training_observer_corr.txt')
+
         x = [i for i in range(observer.shape[0])]
-    
+
         #plt.subplot(221)
         #xg = np.asarray(range(10))
         #yg = np.asarray(range(10))
@@ -126,7 +126,7 @@ def observe_tfim1d(args):
         #plt.subplot(222)
         #plt.contourf(X,Y,correlations,20, cmap=plt.cm.rainbow)
         #plt.contour(X,Y,correlations,levels,colors=lineCol,linewidths=1.5)
- 
+
         plt.subplot(221)
         plt.ylabel('O',fontsize=25)
         plt.plot(x,observer[:,0],color='red', marker='o',markersize=mS,linewidth=lw)
@@ -135,27 +135,27 @@ def observe_tfim1d(args):
         plt.subplot(222)
         plt.ylabel('$<H>$',fontsize=25)
         plt.plot(x,observer[:,2],color='red',marker='o',markersize=mS,linewidth=lw)
-        plt.axhline(y=exact_E, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact') 
- 
+        plt.axhline(y=exact_E, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact')
+
         plt.subplot(223)
         plt.ylabel('$<|S_z|>$',fontsize=25)
         plt.plot(x,observer[:,3],color='red',marker='o',markersize=mS,linewidth=lw)
-        plt.axhline(y=exact_Sz, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact') 
-        plt.yticks([0,0.25,0.5,0.75,1.0]) 
+        plt.axhline(y=exact_Sz, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact')
+        plt.yticks([0,0.25,0.5,0.75,1.0])
         plt.subplot(224)
         plt.ylabel('$<S_x>$',fontsize=25)
         plt.plot(x,observer[:,4],color='red',marker='o',markersize=mS,linewidth=lw)
-        plt.axhline(y=exact_Sx, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact') 
-        plt.yticks([0,0.25,0.5,0.75,1.0]) 
- 
+        plt.axhline(y=exact_Sx, xmin=0, xmax=x[-1], linewidth=2, color = 'k',label='Exact')
+        plt.yticks([0,0.25,0.5,0.75,1.0])
+
         plt.tight_layout()
         plt.pause(0.1)
-    
+
     plt.show()
 
 
 if __name__ == "__main__":
-    
+
     """ Read command line arguments """
     parser = argparse.ArgumentParser()
 
@@ -165,7 +165,7 @@ def observe_tfim1d(args):
     parser.add_argument('-L',type=int,default=4)
     parser.add_argument('-T',type=float)
     parser.add_argument('-B',type=float)
-     
+
     args = parser.parse_args()
 
     if args.model == 'ising2d':
diff --git a/tutorial/plot_results.py b/tutorial/plot_results.py
index c635f55..dd02c73 100644
--- a/tutorial/plot_results.py
+++ b/tutorial/plot_results.py
@@ -14,10 +14,10 @@ def plot_ising2d_observables(L):
     lineW = 2.0
     lineS = '-'
     markerSize = 12
-   
-    # RBMs used in the plot (sorted with number of hidden nodes) 
+
+    # RBMs used in the plot (sorted with number of hidden nodes)
     hidden = [4,16,64]
-    
+
     # Open MC data file
     nameMC  = '../data/ising2d/observables/MC_ising2d_L'
     nameMC += str(L)
@@ -40,15 +40,15 @@ def plot_ising2d_observables(L):
 
         error = 'd' + obs
         data = dataMC[:,header.index(obs)]
-        plt.plot(x,data,color="k", 
+        plt.plot(x,data,color="k",
             marker='o',
-            linewidth=lineW, 
+            linewidth=lineW,
             linestyle=lineS,
             markersize = markerSize)
-        
+
         # Plot from RBMs with different number of hidden nodes
         for i, n_h in enumerate(hidden):
-        
+
             nameRBM = '../data/ising2d/observables/RBM_nH'
             nameRBM += str(n_h) + '_ising2d_L'
             nameRBM += str(L) + '_Observables.txt'
@@ -57,9 +57,9 @@ def plot_ising2d_observables(L):
             dataRBM = np.loadtxt(fileRBM)
             data = dataRBM[:,header.index(obs)]
             err  = dataRBM[:,header.index(error)]
-            lab = '$n_h$=' 
+            lab = '$n_h$='
             lab += str(n_h)
-            
+
             plt.errorbar(x,data,yerr=err,
                     color = colors[i],
                     linewidth = lineW,
@@ -67,7 +67,7 @@ def plot_ising2d_observables(L):
                     marker = mark,
                     markersize = markerSize,
                     label =lab)
-    
+
         # Set Label and Tickes
         plt.xlabel('$T$', fontsize=20)
         plt.xticks([0,5,10],[1.0,2.269,3.54],fontsize=15)
@@ -75,21 +75,21 @@ def plot_ising2d_observables(L):
         ylab = '$' + obs + '$'
         plt.ylabel(ylab,fontsize=20)
 
-    plt.legend(loc='upper left')    
+    plt.legend(loc='upper left')
     plt.tight_layout()
     plt.show()
     #savefig('observables.pdf', format='pdf', dpi=1000)
 
 
 def plot_tfim1d_observables(L):
-    
+
     fig = plt.figure(figsize=(12,9), facecolor='w', edgecolor='k')
-    
+
     ''' WRITE HERE THE PLOT FUNCTION FOR THE TFIM '''
 
 
 if __name__ == "__main__":
-    
+
     """ Read command line arguments """
     parser = argparse.ArgumentParser()
 
@@ -97,7 +97,7 @@ def plot_tfim1d_observables(L):
     parser.add_argument('-L',type=int)
     parser.add_argument('-T',type=float)
     parser.add_argument('-B',type=float)
-     
+
     par = parser.parse_args()
 
     if par.model == 'ising2d':
diff --git a/tutorial/plot_sampler.py b/tutorial/plot_sampler.py
index 69578b5..0cb7f21 100644
--- a/tutorial/plot_sampler.py
+++ b/tutorial/plot_sampler.py
@@ -14,7 +14,7 @@ def observe_ising2d():
     mS = 9
     mS_exact = 3
     plt.figure(figsize=(12,9), facecolor='w', edgecolor='k')
-    
+
     # Load the MC averages
     nameMC  = '../data/ising2d/observables/MC_ising2d_L'
     nameMC += str(L)
@@ -22,15 +22,15 @@ def observe_ising2d():
     fileMC = open(nameMC,'r')
     header = fileMC.readline().lstrip('#').split()
     dataMC = np.loadtxt(fileMC)
-    xMC = [i for i in range(len(dataMC))] 
- 
+    xMC = [i for i in range(len(dataMC))]
+
     while(True):
-         
+
         plt.clf()
-                
+
         # Load the rbm averages
         observables = np.loadtxt('../data/ising2d/observables/sampler_observer.txt')
-        
+
         # Plot the energy
         plt.subplot(221)
         x = [i for i in range(observables.shape[0])]
@@ -41,7 +41,7 @@ def observe_ising2d():
         plt.xlim(-0.1,10.1)
         plt.ylabel('$<E>$',fontsize=25)
         plt.xlabel('$T$',fontsize=25)
- 
+
         # Plot the magnetization
         plt.subplot(222)
         data = dataMC[:,header.index('M')]
@@ -51,17 +51,17 @@ def observe_ising2d():
         plt.xlim(-0.1,10.1)
         plt.ylabel('$<|M|>$',fontsize=25)
         plt.xlabel('$T$',fontsize=25)
- 
+
         # Plot the specific heat
         plt.subplot(223)
         data = dataMC[:,header.index('C')]
         plt.plot(x,observables[:,2],color='red',marker='o',markersize=mS,linewidth=lw)
         plt.plot(xMC,data,color='k',marker='o',linewidth=lw_exact,markersize=mS_exact)
         plt.ylim(-0.05,1.2)
-        plt.xlim(-0.1,10.1) 
+        plt.xlim(-0.1,10.1)
         plt.ylabel('$<C_V>$',fontsize=25)
         plt.xlabel('$T$',fontsize=25)
- 
+
         # Plot the susceptibility
         plt.subplot(224)
         data = dataMC[:,header.index('S')]
@@ -72,13 +72,13 @@ def observe_ising2d():
         plt.ylabel('$<\chi>$',fontsize=25)
         plt.xlabel('$T$',fontsize=25)
 
-        plt.tight_layout() 
+        plt.tight_layout()
         plt.pause(0.2)
-    
+
 # Plot the average observables for the 1d transverse-field ising model
 def observe_tfim1d():
     plt.ion()
-    
+
     # Plot properties
     L=10
     plt.figure(figsize=(12,9), facecolor='w', edgecolor='k')
@@ -89,7 +89,7 @@ def observe_tfim1d():
     levels = [0.000,0.0005,0.001,0.005,0.01,0.02,0.03,0.04,0.06,0.08,0.1,0.2]
     lineCol =['#1a1a1a','#1a1a1a','#1a1a1a','#1a1a1a','#1a1a1a','#1a1a1a',
               '#1a1a1a','#1a1a1a','#262626','#333333','#404040','#4d4d4d','#595959']
-   
+
     # Load the ed data
     nameMC  = '../data/tfim1d/observables/ed_tfim1d_L'
     nameMC += str(L)
@@ -105,25 +105,25 @@ def observe_tfim1d():
     for i in range(L):
         for j in range(L):
             ed_corr_crit[i][j] = ed_corr_full[4*L+i][j]
- 
+
     while(True):
-         
+
         plt.clf()
-                
+
         # Load the rbm data
         observables = np.loadtxt('../data/tfim1d/observables/sampler_observer.txt')
-        correlations= np.loadtxt('../data/tfim1d/observables/sampler_observer_corr.txt') 
+        correlations= np.loadtxt('../data/tfim1d/observables/sampler_observer_corr.txt')
         plt.subplot(221)
-        
+
         # Plot magnetization Z
         data = dataMC[:,header.index('<|Sz|>')]
-        xMC = [i for i in range(len(data))] 
+        xMC = [i for i in range(len(data))]
         x = [i for i in range(observables.shape[0])]
         plt.plot(xMC,data,color='k',marker='o',linewidth=lw_exact,markersize=mS_exact)
         plt.plot(x,observables[:,0],color='red',marker='o',markersize=mS,linewidth=lw)
         plt.ylabel('$<|S_z|>$',fontsize=25)
         plt.xlabel('$B$',fontsize=25)
-        plt.ylim([-0.05,1.05]) 
+        plt.ylim([-0.05,1.05])
         # Plot magnetization X
         plt.subplot(222)
         data = dataMC[:,header.index('<Sx>')]
@@ -138,21 +138,21 @@ def observe_tfim1d():
         xg = np.asarray(range(10))
         yg = np.asarray(range(10))
         X,Y = np.meshgrid(xg,yg)
-        
+
         plt.contourf(X,Y,ed_corr_crit,50, cmap=plt.cm.rainbow)
         plt.contour(X,Y,ed_corr_crit,levels,colors=lineCol,linewidths=1.5)
-        
+
         # Plot rbm critical ZZ correlations
         plt.subplot(224)
         plt.contourf(X,Y,correlations,50, cmap=plt.cm.rainbow)
         plt.contour(X,Y,correlations,levels,colors=lineCol,linewidths=1.5)
- 
+
         plt.tight_layout()
         plt.pause(0.1)
-    
+
 
 if __name__ == "__main__":
-    
+
     """ Read command line arguments """
     parser = argparse.ArgumentParser()
     parser.add_argument('model',type=str)
diff --git a/tutorial/rbm.py b/tutorial/rbm.py
index b993f32..66e619e 100644
--- a/tutorial/rbm.py
+++ b/tutorial/rbm.py
@@ -3,9 +3,9 @@
 import numpy as np
 
 class RBM(object):
-    
+
     ''' Restricted Boltzmann Machine '''
-    
+
     def __init__(self, num_hidden, num_visible, num_samples=128, weights=None, visible_bias=None, hidden_bias=None):
         ''' Constructor '''
         # number of hidden units
@@ -59,7 +59,7 @@ def _create_parameter_variable(initial_value=None, default=None):
             initial_value = default
         return tf.Variable(initial_value)
 
-     
+
     def p_of_h_given(self, v):
         ''' Conditional probability of hidden layer given visible state '''
         # type: (tf.Tensor) -> tf.Tensor
@@ -105,7 +105,7 @@ class variables.
         self.hidden_samples = self.hidden_samples.assign(h_samples)
         self.p_of_v = p_of_v
         return self.hidden_samples, v_samples
-    
+
     def observer_sampling(self, num_iterations):
         # type: (int) -> (tf.Tensor, tf.Tensor, tf.Tensor)
         """
@@ -121,7 +121,7 @@ class variables.
             h_samples, _ = self.sample_h_given(v_samples)
         self.obs_hid_samples = self.obs_hid_samples.assign(h_samples)
         return v_samples
- 
+
     def energy(self, hidden_samples, visible_samples):
         # type: (tf.Tensor, tf.Tensor) -> tf.Tensor
         """Compute the energy:
@@ -143,8 +143,8 @@ def free_energy(self, visible_samples):
         free_energy = (tf.matmul(visible_samples, self.visible_bias)
                        + tf.reduce_sum(tf.nn.softplus(tf.matmul(visible_samples, self.weights)
                                                       + tf.transpose(self.hidden_bias)), 1, keep_dims=True))
-        return free_energy   
-    
+        return free_energy
+
     def neg_log_likelihood_grad(self, visible_samples, model_samples=None, num_gibbs=2):
         # type: (tf.Tensor, tf.Tensor, int) -> tf.Tensor
 
@@ -162,7 +162,7 @@ def neg_log_likelihood(self, visible_samples, log_Z):
                        + tf.reduce_sum(tf.nn.softplus(tf.matmul(visible_samples, self.weights)
                                                       + tf.transpose(self.hidden_bias)), 1))
         return -tf.reduce_mean(free_energy - log_Z)
-    
+
     def exact_log_partition_function(self):
         ''' Evaluate the partition function by exact enumerations '''
         with tf.name_scope('exact_log_Z'):
diff --git a/tutorial/sampler.py b/tutorial/sampler.py
index 6078305..d59ed44 100644
--- a/tutorial/sampler.py
+++ b/tutorial/sampler.py
@@ -9,18 +9,18 @@
 import argparse
 
 #def main():
-#    
+#
 #    # Initialize the command line parser
 #    parser = argparse.ArgumentParser()
 #    # Read command line arguments
-#    parser.add_argument('model',type=str,help='model to sample') 
+#    parser.add_argument('model',type=str,help='model to sample')
 #    # Parse the arguments
-#    
+#
 #    args = parser.parse_args()
 #
 #    if args.model == 'ising2d':
 #        sample_ising2d()
-#    
+#
 #    if args.model == 'tfim1d':
 #        sample_tfim1d()
 
@@ -32,16 +32,16 @@ def ising2d():
            2.269,2.524,2.778,3.032,3.286,3.540]
     num_visible = L*L     # number of visible nodes
     num_hidden = 4       # number of hidden nodes
-    
+
     # Build lattice
     path_to_lattice = '../data/ising2d/lattice2d_L'+str(L)+'.txt'
     nn=np.loadtxt(path_to_lattice)
-    
+
     # Sampling parameters
     num_samples=500 # how many independent chains will be sampled
     gibb_updates=2  # how many gibbs updates per call to the gibbs sampler
-    nbins=1000      # number of calls to the RBM sampler      
-    
+    nbins=1000      # number of calls to the RBM sampler
+
     rbms = []
     rbm_samples = []
     for i in range(len(temps)):
@@ -53,68 +53,68 @@ def ising2d():
         hidden_bias = params['hidden_bias']
         hidden_bias=np.reshape(hidden_bias,(hidden_bias.shape[0],1))
         visible_bias=np.reshape(visible_bias,(visible_bias.shape[0],1))
-    
+
         # Initialize RBM class
         rbms.append(RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples))
         rbm_samples.append(rbms[i].stochastic_maximum_likelihood(gibb_updates))
-    
+
     # Initialize tensorflow
     init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
-    
+
     # Thermodynamic observables
     N = num_visible
-    
+
     with tf.Session() as sess:
         sess.run(init)
-    
+
         for i in range(nbins):
             print ('bin %d\t' %i,end='')
-            print() 
+            print()
             fout = open('../data/ising2d/observables/sampler_observer.txt','w')
-     
+
             for t in range(len(temps)):
                 _,samples=sess.run(rbm_samples[t])
                 spins = np.asarray((2*samples-1))
-    
+
                 m_avg = np.mean(np.absolute(np.sum(spins,axis=1)))
                 e = np.zeros((num_samples))
                 e2= np.zeros((num_samples))
-    
+
                 for k in range(num_samples):
                     for i in range(N):
                         e[k] += -spins[k,i]*(spins[k,int(nn[i,0])]+spins[k,int(nn[i,1])])
                     e2[k] = e[k]*e[k]
                 e_avg = np.mean(e)
-                e2_avg= np.mean(e2) 
+                e2_avg= np.mean(e2)
                 m2_avg = np.mean(np.multiply(np.sum(spins,axis=1),np.sum(spins,axis=1)))
                 c = (e2_avg-e_avg*e_avg)
                 s = (m2_avg-m_avg*m_avg)
-                
+
                 fout.write('%.6f  ' % (e_avg/float(N)))
                 fout.write('%.6f  ' % (m_avg/float(N)))
                 fout.write('%.6f  ' % (c/float(N*temps[t]**2)))
                 fout.write('%.6f\n' % (s/float(N*temps[t])))
 
 def tfim1d():
-    
+
     L=10
     fields = [0.2,0.4,0.6,0.8,1.0,
               1.2,1.4,1.6,1.8,2.0]
-    
+
     num_visible = L     # number of visible nodes
     num_hidden = 10     # number of hidden nodes
-    
+
     # Sampling parameters
     num_samples=500 # how many independent chains will be sampled
     gibb_updates=2 # how many gibbs updates per call to the gibbs sampler
-    nbins=1000      # number of calls to the RBM sampler      
-    
+    nbins=1000      # number of calls to the RBM sampler
+
     class Placeholders(object):
         pass
-    
+
     placeholders = Placeholders()
     placeholders.visible_samples = tf.placeholder(tf.float32, shape=(None, num_visible), name='v') # placeholder for training data
-    
+
     rbms = []
     rbm_samples = []
     psi_list = []
@@ -129,40 +129,40 @@ class Placeholders(object):
         hidden_bias = params['hidden_bias']
         hidden_bias=np.reshape(hidden_bias,(hidden_bias.shape[0],1))
         visible_bias=np.reshape(visible_bias,(visible_bias.shape[0],1))
-    
+
         # Initialize RBM class
         rbms.append(RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples))
         rbm_samples.append(rbms[i].stochastic_maximum_likelihood(gibb_updates))
         logZ.append(rbms[i].exact_log_partition_function())
         psi_list.append(tf.exp(0.5*rbms[i].free_energy(placeholders.visible_samples)))
-    
+
     # Initialize tensorflow
     init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
-    
+
     psix = []
     lz = []
-    
+
     with tf.Session() as sess:
         sess.run(init)
         for i in range(len(fields)):
             lz.append(sess.run(logZ[i]))
-            psix.append(sess.run(psi_list[i],feed_dict={placeholders.visible_samples: all_v_states})) 
-            for j in range(1<<num_visible): 
+            psix.append(sess.run(psi_list[i],feed_dict={placeholders.visible_samples: all_v_states}))
+            for j in range(1<<num_visible):
                 psix[i][j] /= m.exp(0.5*lz[i])
-    
+
         for i in range(nbins):
             print ('bin %d\t' %i,end='')
-            print() 
+            print()
             fout = open('../data/tfim1d/observables/sampler_observer.txt','w')
-            fcorr = open('../data/tfim1d/observables/sampler_observer_corr.txt','w') 
-    
+            fcorr = open('../data/tfim1d/observables/sampler_observer_corr.txt','w')
+
             for b in range(len(fields)):
                 _,samples=sess.run(rbm_samples[b])
                 spins = np.asarray((2*samples-1))
-                
+
                 # Compute average of longitudinal magnetizations
                 sZ_avg = np.mean(np.absolute(np.sum(spins,axis=1)))
-                
+
                 # Compute averages of energies
                 sX= np.zeros((num_samples))
                 for k in range(num_samples):
@@ -173,16 +173,16 @@ class Placeholders(object):
                         state_flip = int(samples[k].dot(1 << np.arange(samples[k].size)[::-1]))
                         sX[k] += float(psix[b][state_flip])/float(psix[b][state])
                         samples[k,i] = 1 - samples[k,i]
-                    
-                sX_avg = np.mean(sX) 
-                
+
+                sX_avg = np.mean(sX)
+
                 if (b==4):
                     sZsZ = np.zeros((L,L))
                     for i in range(L):
                         for j in range(L):
                             for k in range(num_samples):
                                 sZsZ[i,j] += spins[k,i]*spins[k,j]/float(num_samples)
-                        
+
                             fcorr.write('%.10f   ' % sZsZ[i,j])
                         fcorr.write('\n')
                 fout.write('%.6f  ' % (sZ_avg/float(L)))
diff --git a/tutorial/tfim1d.py b/tutorial/tfim1d.py
index 4adcbed..f9cbee6 100644
--- a/tutorial/tfim1d.py
+++ b/tutorial/tfim1d.py
@@ -15,7 +15,7 @@ def save_parameters(sess, results_dir, rbm, epochs,L,B):
     parameter_file_path =  '../data/tfim1d/parameters/parameters_L' + str(L)
     parameter_file_path += '_B' + str(B)
     np.savez_compressed(parameter_file_path, weights=weights, visible_bias=visible_bias, hidden_bias=hidden_bias,
-                        epochs=epochs) 
+                        epochs=epochs)
 class Args(object):
     pass
 
@@ -26,7 +26,7 @@ class Ops(object):
     pass
 
 def train(args):
-   
+
     # Simulation parameters
     B = args.B                      # magnetic field
     num_visible = args.L            # number of visible nodes
@@ -41,28 +41,28 @@ def train(args):
     hidden_bias=None                # hidden bias
     bcount=0                        # counter
     epochs_done=1                   # epochs counter
- 
+
     # Loading the data
     train_dir = '../data/tfim1d/datasets/'  # Location of training data.
     trainName = '../data/tfim1d/datasets/tfim1d_L'+str(args.L)+'_B'+str(B)+'_train.txt'
     testName = '../data/tfim1d/datasets/tfim1d_L'+str(args.L)+'_B'+str(B)+'_test.txt'
     xtrain = np.loadtxt(trainName)
     xtest = np.loadtxt(testName)
-    
+
     ept=np.random.permutation(xtrain) # random permutation of training data
     epv=np.random.permutation(xtest) # random permutation of test data
-    iterations_per_epoch = xtrain.shape[0] / bsize  
+    iterations_per_epoch = xtrain.shape[0] / bsize
 
     # Initialize RBM class
-    rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples) 
-    
+    rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples)
+
     # Initialize operations and placeholders classes
     ops = Ops()
     placeholders = Placeholders()
 
     placeholders.visible_samples = tf.placeholder(tf.float32, shape=(None, num_visible), name='v') # placeholder for training data
 
-    total_iterations = 0 # starts at zero 
+    total_iterations = 0 # starts at zero
     ops.global_step = tf.Variable(total_iterations, name='global_step_count', trainable=False)
 
     learning_rate = tf.train.exponential_decay(
@@ -79,15 +79,15 @@ def train(args):
     ops.lr=learning_rate
     ops.train = optimizer.minimize(cost, global_step=ops.global_step)
     ops.init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
-    
+
     logZ = rbm.exact_log_partition_function()
-    placeholders.logZ = tf.placeholder(tf.float32) 
+    placeholders.logZ = tf.placeholder(tf.float32)
     NLL = rbm.neg_log_likelihood(placeholders.visible_samples,placeholders.logZ)
     path_to_wf='../data/tfim1d/wavefunctions/wavefunction_tfim1d_L'+str(args.L)+'_B'+str(B)+'.txt'
     wf=np.loadtxt(path_to_wf)
     psi_x = tf.exp(0.5*rbm.free_energy(placeholders.visible_samples))
     all_v_states= np.array(list(it.product([0, 1], repeat=num_visible)), dtype=np.float32)
-    
+
     observer_file=open('../data/tfim1d/observables/training_observer.txt','w',0)
     observer_file.write('#      O')
     observer_file.write('        NLL')
@@ -98,14 +98,14 @@ def train(args):
     gibb_updates=10
     observer_samples=rbm.observer_sampling(gibb_updates)
     nbins=100
-    
+
     e = np.zeros((num_samples))
     sX= np.zeros((num_samples))
     L = num_visible
-    
+
     with tf.Session() as sess:
         sess.run(ops.init)
-        
+
         for ii in range(nsteps):
             if bcount*bsize+ bsize>=xtrain.shape[0]:
                bcount=0
@@ -114,48 +114,48 @@ def train(args):
             batch=ept[ bcount*bsize: bcount*bsize+ bsize,:]
             bcount=bcount+1
             feed_dict = {placeholders.visible_samples: batch}
-            
+
             _, num_steps = sess.run([ops.train, ops.global_step], feed_dict=feed_dict)
 
             if num_steps % iterations_per_epoch == 0:
                 print ('Epoch = %d     ' % epochs_done,end='')
                 lz = sess.run(logZ)
                 nll = sess.run(NLL,feed_dict={placeholders.visible_samples: epv, placeholders.logZ: lz})
-                psix = sess.run(psi_x,feed_dict={placeholders.visible_samples: all_v_states}) 
- 
+                psix = sess.run(psi_x,feed_dict={placeholders.visible_samples: all_v_states})
+
                 Ov = 0.0
                 E = 0.0
                 Sz = 0.0
                 Sx = 0.0
                 SzSz = np.zeros((L,L))
 
-                for i in range(1<<num_visible): 
+                for i in range(1<<num_visible):
                     psix[i] /= m.exp(0.5*lz)
                     Ov += wf[i]*psix[i]
-                
+
                 for i in range(nbins):
-                    
+
                     # Gibbs sampling
                     samples=sess.run(observer_samples)
                     spins = np.asarray((2*samples-1))
-                    
+
                     # Compute average of longitudinal magnetizations
                     sZ_avg = np.mean(np.absolute(np.sum(spins,axis=1)))
-                    
+
                     # Compute averages of energies
                     e.fill(0.0)
                     sX.fill(0.0)
 
                     for k in range(num_samples):
                         state = int(samples[k].dot(1 << np.arange(samples[k].size)[::-1]))
-                        
+
                         # Compute the average of transverse magnetization
                         for i in range(L):
                             samples[k,i] = 1 - samples[k,i]
                             state_flip = int(samples[k].dot(1 << np.arange(samples[k].size)[::-1]))
                             sX[k] += float(psix[state_flip])/float(psix[state])
                             samples[k,i] = 1 - samples[k,i]
-                        
+
                         # Compute the correlations ZZ
                         for i in range(L):
                             for j in range(L):
@@ -173,14 +173,14 @@ def train(args):
                         state_flip = int(samples[k].dot(1 << np.arange(samples[k].size)[::-1]))
                         e[k] += -B*psix[state_flip]/psix[state]
                         samples[k,L-1] = 1 - samples[k,L-1]
-                    sX_avg = np.mean(sX) 
+                    sX_avg = np.mean(sX)
                     e_avg = np.mean(e)
-                    
+
                     E += (e_avg-E)/float(i+1)
                     Sz += (sZ_avg-Sz)/float(i+1)
                     Sx += (sX_avg-Sx)/float(i+1)
- 
-                               
+
+
                 # Print observer on screen
                 print ('Ov = %.6f     ' % Ov,end='')
                 print ('NLL = %.6f     ' % nll,end='')
@@ -188,7 +188,7 @@ def train(args):
                 print ('<|Sz|> = %.6f     ' % (Sz/float(L)),end='')
                 print ('<Sx> = %.6f     ' % (Sx/float(L)),end='')
                 print()
- 
+
                 # Save observer on file
                 observer_file.write('%.6f   ' % Ov)
                 observer_file.write('%.6f   ' % nll)
@@ -196,7 +196,7 @@ def train(args):
                 observer_file.write('%.6f   ' % (Sz/float(L)))
                 observer_file.write('%.6f   ' % (Sx/float(L)))
                 observer_file.write('\n')
-                
+
                 #observer_corr_file=open('../data/tfim1d/observables/training_observer_corr.txt','w')
                 #for i in range(L):
                 #    for j in range(L):
@@ -207,12 +207,12 @@ def train(args):
 
 
 def sample(args):
-       
+
     # Architecture
     B = args.B                      # magnetic field
     num_visible = args.L            # number of visible nodes
     num_hidden = args.nH            # number of hidden nodes
-    
+
     # Load the RBM parameters
     path_to_params = '../data/tfim1d/parameters/parameters_nH'+str(num_hidden) + '_L'+str(args.L)+'_B'+str(B)+'.npz'
     params = np.load(path_to_params)
@@ -221,11 +221,11 @@ def sample(args):
     hidden_bias = params['hidden_bias']
     hidden_bias=np.reshape(hidden_bias,(hidden_bias.shape[0],1))
     visible_bias=np.reshape(visible_bias,(visible_bias.shape[0],1))
-    
+
     # Sampling parameters
     num_samples=1000   # how many independent chains will be sampled
     gibb_updates=10    # how many gibbs updates per call to the gibbs sampler
-    nbins=1000         # number of calls to the RBM sampler      
+    nbins=1000         # number of calls to the RBM sampler
 
     # Initialize RBM class
     rbm = RBM(num_hidden=num_hidden, num_visible=num_visible, weights=weights, visible_bias=visible_bias,hidden_bias=hidden_bias, num_samples=num_samples)
@@ -236,32 +236,32 @@ def sample(args):
     # Initialize tensorflow
     init = tf.group(tf.initialize_all_variables(), tf.initialize_local_variables())
     logZ = rbm.exact_log_partition_function()
-    placeholders.logZ = tf.placeholder(tf.float32) 
+    placeholders.logZ = tf.placeholder(tf.float32)
     psi_x = tf.exp(0.5*rbm.free_energy(placeholders.visible_samples))
     all_v_states= np.array(list(it.product([0, 1], repeat=num_visible)), dtype=np.float32)
 
     sX= np.zeros((num_samples))
     L = num_visible
-    e = np.zeros((num_samples)) 
+    e = np.zeros((num_samples))
     with tf.Session() as sess:
         sess.run(init)
         lz = sess.run(logZ)
-        psix = sess.run(psi_x,feed_dict={placeholders.visible_samples: all_v_states}) 
+        psix = sess.run(psi_x,feed_dict={placeholders.visible_samples: all_v_states})
         #Ov=0.0
-        for i in range(1<<num_visible): 
+        for i in range(1<<num_visible):
             psix[i] /= m.exp(0.5*lz)
         #    Ov += wf[i]*psix[i]
- 
+
         for i in range(nbins):
             print ('bin %d\t\t' %i,end='')
-            
+
             # Gibbs sampling
             _,samples=sess.run([hsamples,vsamples])
             spins = np.asarray((2*samples-1))
-            
+
             # Compute average of longitudinal magnetizations
             sZ_avg = np.mean(np.absolute(np.sum(spins,axis=1)))
-            
+
             # Compute averages of energies
             e.fill(0.0)
             sX.fill(0.0)
@@ -273,7 +273,7 @@ def sample(args):
                     state_flip = int(samples[k].dot(1 << np.arange(samples[k].size)[::-1]))
                     sX[k] += float(psix[state_flip])/float(psix[state])
                     samples[k,i] = 1 - samples[k,i]
-            
+
                 # Compute the Energy
                 for i in range(L-1):
                     e[k] += -spins[k,i]*spins[k,i+1]
@@ -286,10 +286,10 @@ def sample(args):
                 state_flip = int(samples[k].dot(1 << np.arange(samples[k].size)[::-1]))
                 e[k] += -B*psix[state_flip]/psix[state]
                 samples[k,L-1] = 1 - samples[k,L-1]
-            
+
             e_avg = np.mean(e)
-            sX_avg = np.mean(sX) 
- 
+            sX_avg = np.mean(sX)
+
             # Print observer on screen
             print ('<H> = %.6f     ' % (e_avg/float(L)),end='')
             print ('<|Sz|> = %.6f     ' % (sZ_avg/float(L)),end='')