# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     name: python3
# ---

# + [markdown] id="view-in-github" colab_type="text"
# <a href="https://colab.research.google.com/github/tatha04/Machine-Learning-Phase-Transitions/blob/main/Machine_Learning_Phase_Transitions.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>

# + [markdown] id="Lx0MPElRyqyB"
# In this notebook, we will use artificial neural networks and supervised machine learning to classify phases of the nearest-neighbor Ising model with ferromagnetic interactions. The training and test datasets have been generated separately. 

# + [markdown] id="3-cTCMEBgtwU"
# ### Import necessary libraries

# + id="ib9-cuY2gYpC"
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import tensorflow as tf
from tensorflow import keras
from keras.models import Sequential
from keras.layers import Dense
from tensorflow.keras import regularizers
from keras.utils import to_categorical

# + [markdown] id="6F7Of1TMgzmK"
# ### Read the data file for training the model.

# + id="0jbMWAKMgqrh"
ising_data = pd.read_csv('/content/drive/MyDrive/Data Science/Data_Sets/Ising/Ising_data1.csv')

# + colab={"base_uri": "https://localhost:8080/", "height": 269} id="xw5xgLzthism" outputId="13983a14-6f89-4b2c-d9c2-635d7a9008d9"
ising_data.head()

# + [markdown] id="GOvBFQKCjc2_"
# The data file contains lattice configurations obtained using Monte-Carlo simulations of the nearest-neighbour ferromagnetic Ising model. For each initial high-temperature configuration, the temperature is gradually reduced from $T/J = 4.9$ to $T/J = 0.5$ in steps of $\delta T/J = 0.1$. (Number of temperature values *ntemp* = 45). At each temperature value, *n_meas* = 10 lattice configurations are saved. Finally, this process is repeated for *n_init_config* = 10 initial configurations. Thus the data file contains a total of 10 x 10 x 45 = 4500 lattice configurations (100 configurations at each temperature value). The square lattice consisting of 20 x 20 sites are reshaped into an array of size = 400  containing values of $\pm1$, depending on whether the site contains a spin $\uparrow$ or $\downarrow$.

# + colab={"base_uri": "https://localhost:8080/"} id="LzGiLd-YgKAX" outputId="bbbb10c7-dfdd-41dc-b2c4-6ab8e827f28c"
X = ising_data.drop(['Temp','Unnamed: 0'],axis=1)
X.shape

# + [markdown] id="kyNGE9nynUMp"
# Obtain the magnetization and magnetic susceptibility for each lattice configuration.

# + id="b8QxPD95inxg"
L = 20
m = X.sum(axis=1)/L**2; m2 = ((X*X).sum(axis=1))/L**2; 
temp = ising_data['Temp']; chi = (m2 - m**2)/temp

mag = pd.DataFrame({'m': m,
                     'chi' : chi,
                     'Temp' : temp})

# + colab={"base_uri": "https://localhost:8080/", "height": 421} id="Q8AoF5l6ijKv" outputId="b825ebd5-49f0-4cde-d0a7-487d055c330c"
plt.figure(figsize=(12,6))
plt.subplot(121)
plt.title('Magnetic Susceptibility')
plt.xlabel('Temp')
plt.scatter(temp,chi)
plt.subplot(122)
plt.title('Magnetization')
plt.xlabel('Temp')
plt.scatter(temp,m)

# + [markdown] id="Q4qpfZDbpS8c"
# ### Critical Temperature ($T=T_c$)
#
# For each Monte-Carlo simulation, plot the magnetic susceptibility ($\chi$) 
# as a function of temperature. The critical temperature $T_c$ for this particular run can be approximated as the temperature at which the $\chi$ is maximum. $T_c$ is obtained by taking an average over all these values.

# + colab={"base_uri": "https://localhost:8080/", "height": 366} id="NDwS6f1ioh69" outputId="0fa30bec-20ad-4110-f3bb-6f04d840eabb"
n_init_config = 10; n_meas = 10;
thi = 4.9; tlo = 0.5; tstep = 0.1; ntemp = int((thi-tlo)/tstep)+1;
temperatures = np.linspace(thi,tlo,ntemp);
chi_config = np.zeros((ntemp,n_init_config))
m_config   = np.zeros((ntemp,n_init_config))
for j in range(n_init_config):
  for i in range(ntemp):
    chi_config[i,j] = np.average(mag['chi'][j*n_meas*ntemp+i*n_meas:j*n_meas*ntemp+(i+1)*n_meas-1])
    m_config[i,j]   = np.average(mag['m'][j*n_meas*ntemp+i*n_meas:j*n_meas*ntemp+(i+1)*n_meas-1])

plt.figure()
plt.ylabel('Susceptibility')
plt.xlabel('Temperature')

for i in range(n_init_config):
  plt.scatter(temperatures,chi_config[:,i])

plt.show()

Tcs = temperatures[np.argmax(chi_config,axis=0)]
Tc = np.average(Tcs)

print('\n Critical temperatures:')
print(Tcs)
print('\n Average critical temperature <Tc> = ',Tc)

# + [markdown] id="4Hfc2uEyrDpZ"
# ## Plot the average magnetization and susceptibility.

# + id="WcYEpmEZkDfB" colab={"base_uri": "https://localhost:8080/", "height": 297} outputId="c668c87b-821d-4895-b220-2616b5a964a7"
mag.sort_values(by='Temp', ascending=False, inplace = True)
# mag['Temp'].value_counts()

nvalues = n_init_config*n_meas
chi_avg = np.zeros((ntemp))
m_avg = np.zeros((ntemp))
for i in range(ntemp):
  chi_avg[i] = np.average(mag['chi'][i*nvalues:(i+1)*nvalues-1])
  m_avg[i]   = np.average(abs(mag['m'][i*nvalues:(i+1)*nvalues-1]))

fig, ax1 = plt.subplots()

plt.title('Avg. Magnetization and Susceptibility')
plt.scatter(temperatures,m_avg,color='b')
plt.ylabel('<|m|>')
plt.xlabel('Temperature')

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
plt.ylabel('<chi>')
plt.scatter(temperatures,chi_avg,color='r')

fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.show()

# + [markdown] id="7M5y7FNA5djv"
# ### Lattice Configurations
#
# Let us plot sample lattice configurations from the three cases:
# *   High-temperature ($T \gt T_c$): There is no long-range order. The spins $\uparrow$ and $\downarrow$ are randomly distributed.
# *   Critical temperature ($T \approx T_c \approx 2.5 $): The spins begin to demonsrate long-range order. Hence, domains of either spin $\uparrow$ or $\downarrow$ appears.
# *   Low-temperature ($T \lt T_c$): All the spins are aligned either $\uparrow$ or $\downarrow$ exhibitng long range order.

# + colab={"base_uri": "https://localhost:8080/", "height": 257} id="3T9h8Rw05beN" outputId="84878509-193f-4337-980f-22c32d3fa2cd"
# Plot T<Tc, T>Tc, and T=Tc

lattice_T_gt_Tc = ising_data.loc[ising_data['Temp']==4.0].sample().drop(['Temp','Unnamed: 0'],axis=1).values.reshape((20,20))
lattice_T_eq_Tc = ising_data.loc[ising_data['Temp']==2.5].sample().drop(['Temp','Unnamed: 0'],axis=1).values.reshape((20,20))
lattice_T_lt_Tc = ising_data.loc[ising_data['Temp']==1.0].sample().drop(['Temp','Unnamed: 0'],axis=1).values.reshape((20,20))

fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(12,4))

axes[0].matshow(lattice_T_gt_Tc)
axes[0].set_title("T =4.0 > Tc")
axes[1].matshow(lattice_T_eq_Tc)
axes[1].set_title("T = T c = 2.5")
axes[2].matshow(lattice_T_lt_Tc)
axes[2].set_title("T = 1.0 < Tc")
for i in range(3):
  axes[i].set_xticks([])
  axes[i].set_yticks([])
plt.show()

# + [markdown] id="oY8Yf0nvKrg_"
# ### Label the training data
#
# Assign a column 'Phase' and assign it to 1 if $T< T_c$ (ordered phase) and 0 if $T > T_c$ (disordered phase).

# + id="uTu3thkXhkHh" colab={"base_uri": "https://localhost:8080/"} outputId="7f782442-aa8b-4ce5-f325-f940a42c700f"
ising_data['Phase'] = (ising_data['Temp']<Tc).astype(int)
Y = ising_data['Phase']
temp = ising_data['Temp']

print(X.shape)
print(Y.shape)

# + [markdown] id="K5rpFRt_M_gm"
# ### Model Building
#
# Build a Neural Network with 1 input layer (400 input dimension), 1 hidden layer (100 neurons) and 1 output layer (2 output nodes).

# + id="-hPpmEXuqtWU" colab={"base_uri": "https://localhost:8080/"} outputId="fdfc2665-f316-4409-9df8-c9e94bb0faad"
# Build the model
# Neural Network with 3 layers
# Input layer: 400 Input nodes
# 1 hidden layers:  100 neurons + L1 regularization (to reduce overfitting) + ReLu activation function 
# Output layer: 2 neuron and softmax function

model = Sequential()
model.add(Dense(100,activation='relu',kernel_regularizer= regularizers.L1(0.01),input_dim = 400))
model.add(Dense(2,activation='softmax'))

# Complie the model
model.compile(
    optimizer = 'adam',
    loss = 'categorical_crossentropy',
    metrics = ['accuracy']
    )

# Train the model
X_Train_fit = model.fit(
    X,
    to_categorical(Y),
    epochs = 10,
)

# + colab={"base_uri": "https://localhost:8080/", "height": 308} id="Kg-f5UHRrh8b" outputId="5654f7d2-93b8-40bf-886e-049478aca936"
# Plot the training performance
plt.figure(figsize=(12,4))
plt.suptitle('Training Performance')

plt.subplot(121)
plt.plot(X_Train_fit.epoch,X_Train_fit.history['accuracy'])
plt.xlabel('Epochs')
plt.ylabel('Accuracy')

plt.subplot(122)
plt.plot(X_Train_fit.epoch,X_Train_fit.history['loss'])
plt.xlabel('Epochs')
plt.ylabel('Loss')

plt.show()

# + [markdown] id="cC680MUx0zvS"
# ### Testing the model
# For testing the model, we will use a separate set of lattice configurations generated using the same process as the training set.
# At first we will repeat the same process as the training set to measure the average critical teperature $T_c$ of the test data. In order to estimate the accuracy of the model on the test data, we will label the 'Phase' of the configuration as 1 or 0 depending on whether the temperature $T < T_c$ or $T > T_c$.

# + id="7YGZUPdq6N65"
ising_data_test = pd.read_csv('/content/drive/MyDrive/Data Science/Data_Sets/Ising/Ising_data2.csv')

# + colab={"base_uri": "https://localhost:8080/", "height": 434} id="5fbwZzns6XyZ" outputId="a62d2993-1079-4f37-e989-3814d3b4e864"
X_test = ising_data_test.drop(['Temp','Unnamed: 0'],axis=1)

m = X_test.sum(axis=1)/400; m2 = ((X_test*X_test).sum(axis=1))/400; 
temp = ising_data['Temp']; chi = (m2 - m**2)/temp

mag_test = pd.DataFrame({'m': m,
                     'chi' : chi,
                     'Temp' : temp})

n_init_config = 10; n_meas = 10;
thi = 4.9; tlo = 0.5; tstep = 0.1; ntemp = int((thi-tlo)/tstep)+1;
temperatures = np.linspace(thi,tlo,ntemp);
chi_config_test = np.zeros((ntemp,n_init_config))
m_config_test   = np.zeros((ntemp,n_init_config))
for j in range(n_init_config):
  for i in range(ntemp):
    chi_config_test[i,j] = np.average(mag_test['chi'][j*n_meas*ntemp+i*n_meas:j*n_meas*ntemp+(i+1)*n_meas-1])
    m_config_test[i,j]   = np.average(mag_test['m'][j*n_meas*ntemp+i*n_meas:j*n_meas*ntemp+(i+1)*n_meas-1])

plt.figure()
plt.title('Test data')
plt.ylabel('Susceptibility')
plt.xlabel('Temperature')

for i in range(n_init_config):
  plt.scatter(temperatures,chi_config_test[:,i])

plt.show()

Tcs_test = temperatures[np.argmax(chi_config_test,axis=0)]
Tc_test = np.average(Tcs_test)

print('\n Critical temperatures:')
print(Tcs_test)
print('\nCritical temperature <Tc> = ',Tc_test,'\n')

Y_test = (ising_data_test['Temp']<Tc_test).astype(int)
print(X_test.shape)
print(Y_test.shape)

# + colab={"base_uri": "https://localhost:8080/"} id="ZGI-gWhwryfe" outputId="2bd18080-ce08-4009-96e2-a6363c40ff79"
# Evaluate the model
model.evaluate(
    X_test,
    to_categorical(Y_test)
)

# + id="KYvKoTBusp7X"
Y_pred = model.predict(X_test)

# + colab={"base_uri": "https://localhost:8080/"} id="eScWZ_1EtSro" outputId="eedb402e-cda9-48a5-d8a1-3b734564ae0c"
Y_pred

# + colab={"base_uri": "https://localhost:8080/"} id="15BqVzlmtUMa" outputId="cc0e5479-bec2-465f-dff8-155962124155"
Y_pred.shape

# + colab={"base_uri": "https://localhost:8080/", "height": 449} id="YKlhYANJv8Ci" outputId="390a2c10-452c-4315-b4af-bb9461654b8e"
plt.figure(figsize = (16,6))

plt.suptitle('Scatter plot of output weights')

plt.subplot(121)
plt.scatter(temp,Y_pred[:,0],label='W_0')
plt.scatter(temp,Y_pred[:,1],label='W_1')
plt.xlabel('Temp.')
plt.legend()

plt.subplot(122)
plt.scatter(mag_test['m'],Y_pred[:,0],label='W0')
plt.scatter(mag_test['m'],Y_pred[:,1],label='W1')
plt.xlabel('Magnetization')
plt.legend()

# + [markdown] id="WZRHziVv7fRl"
# From the above plot, it is clear that for configurations corresponding to $T < T_c$, the model identifies them (overwhelmingly) belonging to the ordered phase (i.e. W1 = 1 and W0 = 0); whereas for $T > T_c$, the model identifies them belonging to the dirordered phase (i.e. W1 = 0 and W0 = 1). For $T \approx T_c$, there is an overlapping of the weights, however there is a gradual crossover of the weights.
#
# If we plot the weights versus magnetization, it becomes clearer that weights $W1 \approx 1$ if $m = \pm 1$.

# + colab={"base_uri": "https://localhost:8080/", "height": 297} id="IXc6mAjl4WUq" outputId="5582c220-32e8-4c43-d455-e117b9655615"
fig, ax1 = plt.subplots()

ax1.scatter(temp,np.abs(mag_test['m']),color='m',label='|m|')
plt.legend()

ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis
ax2.scatter(temp,Y_pred[:,0],color='b',label='W0')
ax2.scatter(temp,Y_pred[:,1],color='g',label='W1')
fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.legend()
plt.xlabel('Temp.')
plt.show()

# + id="yeq6bMNtG1jr"
# Sort test data w.r.t. temperature
ising_data_test.sort_values(by='Temp', ascending=False,inplace=True)
Y_test = (ising_data_test['Temp']<Tc_test).astype(int)
X_test = ising_data_test.drop(['Temp','Unnamed: 0'],axis=1)
Y_pred = model.predict(X_test)

# + colab={"base_uri": "https://localhost:8080/"} id="BjvFy3NK2YNS" outputId="1b6e83bc-352a-4c91-df8f-8e5f8bdc3a3b"
Y_pred

# + id="VVXEqcUyVuL6"
nvalues = n_init_config*n_meas

Y_pred_comb = np.zeros((ntemp,2))
for i in range(ntemp):
  Y_pred_comb[i,0] = np.average(Y_pred[i*nvalues:(i+1)*nvalues-1,0])
  Y_pred_comb[i,1] = np.average(Y_pred[i*nvalues:(i+1)*nvalues-1,1])


# + colab={"base_uri": "https://localhost:8080/", "height": 449} id="vCa8bfSUV7_I" outputId="87d80f71-9989-4b02-86ea-1d36efb52f6b"
plt.figure(figsize=(18,6))
plt.suptitle('Average Output weights vs Temperature')

plt.subplot(121)
plt.plot(temperatures,Y_pred_comb[:,0],'.-',label='W0')
plt.plot(temperatures,Y_pred_comb[:,1],'.-',label='W1')
plt.xlabel('Temp.')
plt.legend()

plt.subplot(122)
plt.xlim((2.0,3.0))
plt.plot(temperatures,Y_pred_comb[:,0],'.-',label='<W0>')
plt.plot(temperatures,Y_pred_comb[:,1],'.-',label='<W1>')
plt.vlines(x=2.4,ymin=0,ymax=1,label='T=2.4')
plt.xlabel('Temp.')
plt.legend()

# + [markdown] id="25eZsMm28_sj"
# As shown in the above graphs, the average weights $<W0>$ and $<W1>$ obtained from the test dataset of the lattice configurations exhibit a crossover at $T\approx 2.4 \approx T_c$.
#
# The above data suffers severly from finite size effcts. Much better results can be obtained by using higher lattice sizes.
#

# + [markdown] id="qxyj7P2i-Iuc"
# The beauty of neural networks and deep learning lies in the fact that the trained model can also be used to locate the phase transitions of other models (e.g. triangular lattice) without any apriori information about the Hamiltonian.

StarcoderdataJupyterScriptsDedupFiltered