import matplotlib.pyplot as plt
from pylab import *
from numpy import *

N = 100									# Number of spins
M = 100									# Number of systems in ensemble
T = 150									# Number of iterations
set_beta = array([10, 1, 0.1, 0.01])	# Betas

# Compute the global energy of the chain
def global_energy(S, J):
	S_shift = hstack([S[1:], S[0]])
	energy = -1*dot(J*S, S_shift)
	return energy

# Initialize spins: each row a string in ensemble, each column a particular spin
S_initial = arange(N*M)
for i in range(N*M):
	if random.random() > 0.5:
		S_initial[i] = 1
	else:
		S_initial[i] = -1
S_initial.resize(M, N)

# Initialize Gaussian weights
f = open('J.txt', 'r')
J = loadtxt(f)[:N].astype(float)
f.close()

# Compute minimum energy
Emin_bound = -sqrt(J*J).sum()

# Initialize energies and relative probabilities
E = 0*arange(M).astype(float)
p = 0*arange(M).astype(float)

# Initialize output
output = 0*arange(T*set_beta.size)
output.resize(T, set_beta.size)

for b in range(set_beta.size):
	# Set beta
	beta = set_beta[b]
	# Reset spins
	S_ori = copy(S_initial)
	S_ter = copy(S_initial)
	for t in range(T):
		# Calculate energies and probabilities
		for i in range(M):
			E[i] = global_energy(S_ori[i,:], J)
			p[i] = exp(-beta*(E[i] - Emin_bound))
		# Compute minimum energy
		output[t,b] = E.min()
		# Calculate normalized cumulative sum
		p_draw = cumsum(p)
		p_draw /= p_draw[-1]
		for i in range(M):
			# Choose two strings at random
			str1_idx = where(p_draw > random.random())[0][0]
			str1 = S_ori[str1_idx,:]
			str2_idx = where(p_draw > random.random())[0][0]
			str2 = S_ori[str2_idx,:]
			idx = int(N*random.random())
			# Build new string
			str = hstack([str1[:idx],str2[idx:]])
			flip_idx = int(N*random.random())
			# Mutate string
			str[flip_idx] *= -1
			# Write string to output matrix
			S_ter[i,:] = str
		# Replace input matrix with output matrix
		S_ori[:] = S_ter[:]

# Plot output
fig = plt.figure()
ax = fig.gca()
X = arange(T)
ax.plot(X, output[:,0], label = 'beta = 10')
ax.plot(X, output[:,1], label = 'beta = 1')
ax.plot(X, output[:,2], label = 'beta = 0.1')
ax.plot(X, output[:,3], label = 'beta = 0.01')
flr = 0*arange(T) + Emin_bound
ax.plot(X, flr, 'k--')
legend()

# Render plot
plt.show()