<jupyter_start><jupyter_text>Data loading and pre-processingWe can use either entanglement spectrum or the absolute values of wavefunction as the input<jupyter_code># load the data by specifying the data type, interaction strengths and number of samples
def load_data(**args):
    # This is only useful if you run the en_spectrum.py script to generate your own data set
    d_type = args['type']
    L = int(args['len'])
    j1 = float(args['j1'])
    j2 = float(args['j2'])
    runs = int(para['run_1']) # number of samples generated
    loaded = np.load("data_set/j_"+str(j1)+"/1_en_spectrum_L="+str(L)+"_j="+str(j1)+".npz")
    # for ent_specta, the size is determined by the cut position
    # for the wavefunction, it's by the system size
    dim = loaded[d_type].shape
    x = np.zeros((runs, dim[0], dim[1]), dtype='float64')
    
    # load data for the phase 1
    for i in range(runs):
        filename = "data_set/j_"+str(j1)+"/"+str(i+1)+"_en_spectrum_L="+str(L)+"_j="+str(j1)+".npz"
        loaded = np.load(filename)
        if d_type == 'wave':
            wav = loaded[d_type]
            x[i] = np.multiply(wav, np.conj(wav)).real
        else:
            x[i] = loaded[d_type]
            
    p1 = np.concatenate((x[:])) 
    
    # load data for the phase 2
    runs = int(para['run_2'])    
    x = np.zeros((runs, dim[0], dim[1]), dtype='float64')
    for i in range(runs):
        filename = "data_set/j_"+str(j2)+"/"+str(i+1)+"_en_spectrum_L="+str(L)+"_j="+str(j2)+".npz"
        loaded = np.load(filename)
        if d_type == 'wave':
            wav = loaded[d_type]
            x[i] = np.multiply(wav, np.conj(wav)).real
        else:
            x[i] = loaded[d_type]
            
    p2 = np.concatenate((x[:])) 

    # cast the original data to single precision to speed up
    p1 = p1.astype('float32')
    p2 = p2.astype('float32')
    
    return p1, p2
para = {'type':'ent','len':'12', 'run_1':'35', 'run_2':'32', 'j1':'0.01', 'j2':'5.0'}

# load data according to parameters
p1, p2 = load_data(**para)
para = {'type':'ent','len':'12', 'run_1':'20', 'run_2':'21', 'j1':'0.05', 'j2':'1.5'}
# load data according to parameters
pp1, pp2 = load_data(**para)
p1 = np.concatenate((p1, pp1))
p2 = np.concatenate((p2, pp2))
print(p1.shape, p2.shape)
from sklearn.model_selection import train_test_split

num_classes = 2
# processing data 
def data_pipeline(p1, p2, test_size = 0.2, encode = True):
    # create labels for the two phases
    l1 = np.ones(p1.shape[0],dtype=np.int8)
    l2 = np.zeros(p2.shape[0],dtype=np.int8)
    # combine labels and the data
    ph1 = np.column_stack((l1, p1))
    ph2 = np.column_stack((l2, p2))
    # train_test split
    dat = np.concatenate((ph1, ph2))
    X_train, X_test, Y_train, Y_test = train_test_split(dat[:,1:], dat[:,:1], test_size=test_size, random_state=1)
    # convert class vectors to binary class matrices
    if encode == True:
        Y_train = tf.keras.utils.to_categorical(Y_train, num_classes)
        Y_test = tf.keras.utils.to_categorical(Y_test, num_classes)
    
    return X_train, X_test, Y_train, Y_test

# train test split
X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2)
from tensorflow.keras import Sequential
from tensorflow.keras.layers import Dense, Dropout, Flatten, BatchNormalization

def create_model():
    # instantiate model
    model = Sequential()
    # add a normalization layer to improve robustness
    model.add(BatchNormalization(axis=-1, input_shape=(X_test.shape[1],)))
    # add a dense all-to-all relu layer
    model.add(Dense(400, activation='relu'))
    # add a dense all-to-all relu layer
    model.add(Dense(50, activation='relu'))
    # apply dropout with rate 0.4
    model.add(Dropout(0.4))
    # soft-max layer
    model.add(Dense(num_classes, activation='softmax'))
    
    # compile the model with loss function and optimizer
    model.compile(loss=tf.keras.losses.binary_crossentropy,
                  optimizer=tf.keras.optimizers.Adam(),
                  metrics=['accuracy'])
    
    return model


from tensorflow.keras.callbacks import EarlyStopping
#early termination to prevent overfitting
early_stopping = EarlyStopping(monitor='val_loss', patience=3)
# training parameters
batch_size = 200
epochs = 14

# create the deep neural net
model_DNN = create_model()

# train DNN and store training info in history
history = model_DNN.fit(X_train, Y_train,
          batch_size=batch_size,
          epochs=epochs,
          verbose=0,
          validation_data=(X_test, Y_test), callbacks=[early_stopping])
# evaluate model
score = model_DNN.evaluate(X_test, Y_test, verbose=0)

# print performance
print('Test loss:', score[0])
print('Test accuracy:', score[1])
# look into training history
# summarize history for accuracy
plt.plot(history.history['accuracy'])
plt.plot(history.history['val_accuracy'])
plt.ylabel('model accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='best')
plt.show()

# summarize history for loss
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.ylabel('model loss')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>Hyperparameters tuningUse [Keras Tuner](https://blog.tensorflow.org/2020/01/hyperparameter-tuning-with-keras-tuner.html) to search for the ideal parameters, including the discrete parameters such number of layers, number neurons in each layer etc., and continuous parameters like the dropout rate and  learning rate which is perphaps the most important one.<jupyter_code>import kerastuner as kt
from kerastuner.tuners import Hyperband, BayesianOptimization


def create_model(hp):
    # instantiate model
    model = Sequential()
    # add a normalization layer to improve robustness
    model.add(BatchNormalization(axis=-1, input_shape=(X_test.shape[1],)))
    # add a dense all-to-all relu layer
#     for i in range(hp.Int('num_layers', 2, 6)):
#         model.add(layers.Dense(units=hp.Int('units_' + str(i), 32, 512, 32),
#                                activation='relu'))
    model.add(Dense(units=hp.Int('layer_1_units',
                                        min_value=32,
                                        max_value=512,
                                        step=64),
                           activation='relu'))
    model.add(Dense(units=hp.Int('layer_2_units',
                                        min_value=32,
                                        max_value=512,
                                        step=64),
                           activation='relu'))
    model.add(Dropout(
      hp.Float('dropout_rate', 0, 0.5, step=0.1, default=0.5)))

    # soft-max layer
    model.add(Dense(num_classes, activation='softmax'))
    
    # compile the model with loss function and optimizer    
    model.compile(
      optimizer=tf.keras.optimizers.Adam(
        hp.Float('learning_rate', 1e-4, 1e-1, sampling='log')),
      loss='binary_crossentropy',
      metrics=['accuracy'])
    
    return model
tuner = kt.Hyperband(
    create_model,
    objective='val_accuracy',
    max_epochs=3,
    hyperband_iterations=2)
b_tuner = BayesianOptimization(
    create_model,
    objective='val_accuracy',
    max_trials=50)
X_train.shape
# using partial dateset for hp tuning
tuner.search(X_train[:20000], Y_train[:20000],
             epochs=2,
             validation_data=(X_test[:4000], Y_test[:4000]))
hp = tuner.get_best_hyperparameters(1)[0]
print(hp.values)<jupyter_output><empty_output><jupyter_text>Re-train the model with optimized hyper-parametersNow that we've found the (hopefully) best hyper parameters, let's retrain the model to see its performance<jupyter_code>def create_model_hp():
    # instantiate model
    model = Sequential()
    # add a normalization layer to improve robustness
    model.add(BatchNormalization(axis=-1, input_shape=(X_test.shape[1],)))
    # add a dense all-to-all relu layer
    model.add(Dense(288, activation='relu'))
    # add a dense all-to-all relu layer
    model.add(Dense(352, activation='relu'))
    model.add(Dropout(0.1))
    # soft-max layer
    model.add(Dense(num_classes, activation='softmax'))
    
    # compile the model with loss function and optimizer
    model.compile(loss=tf.keras.losses.binary_crossentropy,
                  optimizer=tf.keras.optimizers.Adam(learning_rate = 0.003746593),
                  metrics=['accuracy'])
    
    return model
# training parameters
batch_size = 200
epochs = 14

# create the deep neural net
model_hp = create_model_hp()

# train DNN and store training info in history
history = model_hp.fit(X_train, Y_train,
          batch_size=batch_size,
          epochs=epochs,
          verbose=0,
          validation_data=(X_test, Y_test), callbacks=[early_stopping])
# evaluate model
new_score = model_hp.evaluate(X_test, Y_test, verbose=0)

# print performance
print('Test loss:', new_score[0])
print('Test accuracy:', new_score[1])

# look into training history
# summarize history for accuracy
plt.plot(history.history['accuracy'])
plt.plot(history.history['val_accuracy'])
plt.ylabel('model accuracy')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='best')
plt.show()

# summarize history for loss
plt.plot(history.history['loss'])
plt.plot(history.history['val_loss'])
plt.ylabel('model loss')
plt.xlabel('epoch')
plt.legend(['train', 'test'], loc='best')
plt.show()<jupyter_output><empty_output><jupyter_text>We can see a slight improve of accurary over the old model and a reduction in the gap between the test and training data  Load other data set with different $j$ to see how well the network generalizesAccording to the energy level statistics, the critical point occurs around $j_c\approx 0.23$.So we use different $j$ values as data set to see how our network works.<jupyter_code>para = {'type':'ent','len':'12', 'run_1':'8', 'run_2':'8', 'j1':'0.1', 'j2':'0.5'}

# load data according to parameters
p1, p2 = load_data(**para)

# train test split
X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2, 0.99)
print(X_train.shape, X_test.shape)
# evaluate model
score = model_DNN.evaluate(X_test, Y_test, verbose=0)
score_hp = model_hp.evaluate(X_test, Y_test, verbose=0)
# print performance
print('Test loss:', score[0], score_hp[0])
print('Test accuracy:', score[1], score_hp[1])
# reverse transform the one-hot matrix to 1-D list
_, _, _, Y_test = data_pipeline(p1, p2, 0.99, encode = False)
Y_pred = model_hp.predict_classes(X_test, verbose=0)
from sklearn.preprocessing import OneHotEncoder
enc = OneHotEncoder(handle_unknown='ignore')
enc.fit(Y_test)
Y_test = enc.transform(Y_test).toarray()

Y_true = enc.inverse_transform(Y_test)
Y_true = Y_true.flatten()
from sklearn.metrics import confusion_matrix
cm = confusion_matrix(Y_true, Y_pred)
print(cm)<jupyter_output><empty_output><jupyter_text>According to the convention of sci-kit learn, the $c_{0,0}$ component of the matrix represents the true negatives. So the total positive rate is<jupyter_code>cm[1,1]+cm[1,0]
para = {'type':'ent','len':'12', 'run_1':'8', 'run_2':'8', 'j1':'0.15', 'j2':'0.31'}

# load data according to parameters
p1, p2 = load_data(**para)

# train test split
X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2, 0.99)
print(X_train.shape, X_test.shape)
# evaluate model
score = model_DNN.evaluate(X_test, Y_test, verbose=0)
score_hp = model_hp.evaluate(X_test, Y_test, verbose=0)
# print performance
print('Test loss:', score[0], score_hp[0])
print('Test accuracy:', score[1], score_hp[1])<jupyter_output><empty_output><jupyter_text>Using critical data set for identifying critical pointsAccording to the energy level statistics, the critical point occurs around $j_c\approx 0.23$.To benchmark this result, we use $j=0.21$ and $j=0.25$ as data set close to the critical point to see how our network works.<jupyter_code>para = {'type':'ent','len':'12', 'run_1':'8', 'run_2':'8', 'j1':'0.21', 'j2':'0.25'}
p1, p2 = load_data(**para)
# train test split
X_train, X_test, Y_train, Y_test = data_pipeline(p1, p2, 0.99)
print(X_train.shape, X_test.shape)
# evaluate model
score = model_DNN.evaluate(X_test, Y_test, verbose=0)
score_hp = model_hp.evaluate(X_test, Y_test, verbose=0)
# print performance
print('Test loss:', score[0], score_hp[0])
print('Test accuracy:', score[1], score_hp[1])<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13560145808019441, 'fraction_numerical': 0.02656946132037262, 'mean_word_length': 3.8415686274509806, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '475788', 'license': 'MIT', 'n_tokens_mistral': 4377, 'n_tokens_neox': 4152, 'n_words': 1019, 'path': 'MBL_ETH_Transition.ipynb', 'repo_name': 'empyriumz/Supervised-learning-of-Many-Body-Localization', 'size': '12345'}