diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..8f13e81
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+**/__pycache__/
+**/.idea/
\ No newline at end of file
diff --git a/README.md b/README.md
index eabf9d3..b0f6213 100644
--- a/README.md
+++ b/README.md
@@ -91,4 +91,6 @@ test.negative
 - Each line corresponds to the line of test.rating, containing 99 negative samples.  
 - Each line is in the format: (userID,itemID)\t negativeItemID1\t negativeItemID2 ...
 
-Last Update Date: December 23, 2018
+
+
+Last Update Date: October 04, 2020
diff --git a/TF_KERAS_GMF.py b/TF_KERAS_GMF.py
new file mode 100644
index 0000000..941920d
--- /dev/null
+++ b/TF_KERAS_GMF.py
@@ -0,0 +1,180 @@
+# coding=utf-8
+
+'''
+Created on Aug 9, 2016
+
+Keras Implementation of Generalized Matrix Factorization (GMF) recommender model in:
+He Xiangnan et al. Neural Collaborative Filtering. In WWW 2017.  
+
+@author: Xiangnan He (xiangnanhe@gmail.com)
+'''
+
+# This TF_KERAS version is tested with Keras 2.3.1 and tensorflow-gpu 1.15.2
+
+# use GPU 0
+# import os
+# os.environ["CUDA_VISIBLE_DEVICES"] = "0"
+
+import numpy as np
+from keras import initializers
+from keras.models import Model
+from keras.layers import Embedding, Input, Dense, merge, Flatten
+from keras.optimizers import Adagrad, Adam, SGD, RMSprop
+from keras.regularizers import l2
+from Dataset import Dataset
+from evaluate import evaluate_model
+from time import time
+import argparse
+
+
+#################### Arguments ####################
+def parse_args():
+    parser = argparse.ArgumentParser(description="Run GMF.")
+    parser.add_argument('--path', nargs='?', default='Data/',
+                        help='Input data path.')
+    parser.add_argument('--dataset', nargs='?', default='ml-1m',
+                        help='Choose a dataset.')
+    parser.add_argument('--epochs', type=int, default=100,
+                        help='Number of epochs.')
+    parser.add_argument('--batch_size', type=int, default=256,
+                        help='Batch size.')
+    parser.add_argument('--num_factors', type=int, default=8,
+                        help='Embedding size.')
+    parser.add_argument('--regs', nargs='?', default='[0,0]',
+                        help="Regularization for user and item embeddings.")
+    parser.add_argument('--num_neg', type=int, default=4,
+                        help='Number of negative instances to pair with a positive instance.')
+    parser.add_argument('--lr', type=float, default=0.001,
+                        help='Learning rate.')
+    parser.add_argument('--learner', nargs='?', default='adam',
+                        help='Specify an optimizer: adagrad, adam, rmsprop, sgd')
+    parser.add_argument('--verbose', type=int, default=1,
+                        help='Show performance per X iterations')
+    parser.add_argument('--out', type=int, default=1,
+                        help='Whether to save the trained model.')
+    return parser.parse_args()
+
+
+def init_normal():
+    return initializers.normal(mean=0.0, stddev=0.01)
+
+
+def get_model(num_users, num_items, latent_dim, regs=[0, 0]):
+    # Input variables
+    user_input = Input(shape=(1,), dtype='int32', name='user_input')
+    item_input = Input(shape=(1,), dtype='int32', name='item_input')
+
+    MF_Embedding_User = Embedding(input_dim=num_users, output_dim=latent_dim, name='user_embedding',
+                                  init=init_normal(), W_regularizer=l2(regs[0]), input_length=1)
+    MF_Embedding_Item = Embedding(input_dim=num_items, output_dim=latent_dim, name='item_embedding',
+                                  init=init_normal(), W_regularizer=l2(regs[1]), input_length=1)
+
+    # Crucial to flatten an embedding vector!
+    user_latent = Flatten()(MF_Embedding_User(user_input))
+    item_latent = Flatten()(MF_Embedding_Item(item_input))
+
+    # Element-wise product of user and item embeddings 
+    # predict_vector = merge([user_latent, item_latent], mode = 'mul')
+    predict_vector = merge.multiply([user_latent, item_latent])
+
+    # Final prediction layer
+    # prediction = Lambda(lambda x: K.sigmoid(K.sum(x)), output_shape=(1,))(predict_vector)
+    prediction = Dense(1, activation='sigmoid', init='lecun_uniform', name='prediction')(predict_vector)
+
+    model = Model(input=[user_input, item_input],
+                  output=prediction)
+
+    return model
+
+
+def get_train_instances(train, num_negatives):
+    user_input, item_input, labels = [], [], []
+    num_users = train.shape[0]
+    for (u, i) in train.keys():
+        # positive instance
+        user_input.append(u)
+        item_input.append(i)
+        labels.append(1)
+        # negative instances
+        for t in range(num_negatives):
+            j = np.random.randint(num_items)
+            # while train.has_key((u, j)):
+            while (u, j) in train:
+                j = np.random.randint(num_items)
+            user_input.append(u)
+            item_input.append(j)
+            labels.append(0)
+    return user_input, item_input, labels
+
+
+if __name__ == '__main__':
+    args = parse_args()
+    num_factors = args.num_factors
+    regs = eval(args.regs)
+    num_negatives = args.num_neg
+    learner = args.learner
+    learning_rate = args.lr
+    epochs = args.epochs
+    batch_size = args.batch_size
+    verbose = args.verbose
+
+    topK = 10
+    evaluation_threads = 1  # mp.cpu_count()
+    print("GMF arguments: %s" % (args))
+    model_out_file = 'Pretrain/%s_GMF_%d_%d.h5' % (args.dataset, num_factors, time())
+
+    # Loading data
+    t1 = time()
+    dataset = Dataset(args.path + args.dataset)
+    train, testRatings, testNegatives = dataset.trainMatrix, dataset.testRatings, dataset.testNegatives
+    num_users, num_items = train.shape
+    print("Load data done [%.1f s]. #user=%d, #item=%d, #train=%d, #test=%d"
+          % (time() - t1, num_users, num_items, train.nnz, len(testRatings)))
+
+    # Build model
+    model = get_model(num_users, num_items, num_factors, regs)
+    if learner.lower() == "adagrad":
+        model.compile(optimizer=Adagrad(lr=learning_rate), loss='binary_crossentropy')
+    elif learner.lower() == "rmsprop":
+        model.compile(optimizer=RMSprop(lr=learning_rate), loss='binary_crossentropy')
+    elif learner.lower() == "adam":
+        model.compile(optimizer=Adam(lr=learning_rate), loss='binary_crossentropy')
+    else:
+        model.compile(optimizer=SGD(lr=learning_rate), loss='binary_crossentropy')
+    # print(model.summary())
+
+    # Init performance
+    t1 = time()
+    (hits, ndcgs) = evaluate_model(model, testRatings, testNegatives, topK, evaluation_threads)
+    hr, ndcg = np.array(hits).mean(), np.array(ndcgs).mean()
+    # mf_embedding_norm = np.linalg.norm(model.get_layer('user_embedding').get_weights())+np.linalg.norm(model.get_layer('item_embedding').get_weights())
+    # p_norm = np.linalg.norm(model.get_layer('prediction').get_weights()[0])
+    print('Init: HR = %.4f, NDCG = %.4f\t [%.1f s]' % (hr, ndcg, time() - t1))
+
+    # Train model
+    best_hr, best_ndcg, best_iter = hr, ndcg, -1
+    for epoch in range(epochs):
+        t1 = time()
+        # Generate training instances
+        user_input, item_input, labels = get_train_instances(train, num_negatives)
+
+        # Training
+        hist = model.fit([np.array(user_input), np.array(item_input)],  # input
+                         np.array(labels),  # labels
+                         batch_size=batch_size, nb_epoch=1, verbose=0, shuffle=True)
+        t2 = time()
+
+        # Evaluation
+        if epoch % verbose == 0:
+            (hits, ndcgs) = evaluate_model(model, testRatings, testNegatives, topK, evaluation_threads)
+            hr, ndcg, loss = np.array(hits).mean(), np.array(ndcgs).mean(), hist.history['loss'][0]
+            print('Iteration %d [%.1f s]: HR = %.4f, NDCG = %.4f, loss = %.4f [%.1f s]'
+                  % (epoch, t2 - t1, hr, ndcg, loss, time() - t2))
+            if hr > best_hr:
+                best_hr, best_ndcg, best_iter = hr, ndcg, epoch
+                if args.out > 0:
+                    model.save_weights(model_out_file, overwrite=True)
+
+    print("End. Best Iteration %d:  HR = %.4f, NDCG = %.4f. " % (best_iter, best_hr, best_ndcg))
+    if args.out > 0:
+        print("The best GMF model is saved to %s" % (model_out_file))
diff --git a/TF_KERAS_MLP.py b/TF_KERAS_MLP.py
new file mode 100644
index 0000000..60dbd36
--- /dev/null
+++ b/TF_KERAS_MLP.py
@@ -0,0 +1,186 @@
+# coding=utf-8
+
+'''
+Created on Aug 9, 2016
+Keras Implementation of Multi-Layer Perceptron (GMF) recommender model in:
+He Xiangnan et al. Neural Collaborative Filtering. In WWW 2017.  
+
+@author: Xiangnan He (xiangnanhe@gmail.com)
+'''
+
+# This TF_KERAS version is tested with Keras 2.3.1 and tensorflow-gpu 1.15.2
+
+# use GPU 0
+# import os
+# os.environ["CUDA_VISIBLE_DEVICES"] = "0"
+
+import numpy as np
+
+from keras import initializers
+from keras.regularizers import l2
+from keras.models import Model
+from keras.layers import Embedding, Input, Dense, merge, Flatten
+from keras.optimizers import Adagrad, Adam, SGD, RMSprop
+from evaluate import evaluate_model
+from Dataset import Dataset
+from time import time
+
+import argparse
+
+
+#################### Arguments ####################
+def parse_args():
+    parser = argparse.ArgumentParser(description="Run MLP.")
+    parser.add_argument('--path', nargs='?', default='Data/',
+                        help='Input data path.')
+    parser.add_argument('--dataset', nargs='?', default='ml-1m',
+                        help='Choose a dataset.')
+    parser.add_argument('--epochs', type=int, default=100,
+                        help='Number of epochs.')
+    parser.add_argument('--batch_size', type=int, default=256,
+                        help='Batch size.')
+    parser.add_argument('--layers', nargs='?', default='[64,32,16,8]',
+                        help="Size of each layer. Note that the first layer is the concatenation of user and item embeddings. So layers[0]/2 is the embedding size.")
+    parser.add_argument('--reg_layers', nargs='?', default='[0,0,0,0]',
+                        help="Regularization for each layer")
+    parser.add_argument('--num_neg', type=int, default=4,
+                        help='Number of negative instances to pair with a positive instance.')
+    parser.add_argument('--lr', type=float, default=0.001,
+                        help='Learning rate.')
+    parser.add_argument('--learner', nargs='?', default='adam',
+                        help='Specify an optimizer: adagrad, adam, rmsprop, sgd')
+    parser.add_argument('--verbose', type=int, default=1,
+                        help='Show performance per X iterations')
+    parser.add_argument('--out', type=int, default=1,
+                        help='Whether to save the trained model.')
+    return parser.parse_args()
+
+
+def init_normal():
+    return initializers.normal(mean=0.0, stddev=0.01)
+
+
+def get_model(num_users, num_items, layers=[20, 10], reg_layers=[0, 0]):
+    assert len(layers) == len(reg_layers)
+    num_layer = len(layers)  # Number of layers in the MLP
+    # Input variables
+    user_input = Input(shape=(1,), dtype='int32', name='user_input')
+    item_input = Input(shape=(1,), dtype='int32', name='item_input')
+
+    MLP_Embedding_User = Embedding(input_dim=num_users, output_dim=layers[0] // 2, name='user_embedding',
+                                   init=init_normal(), W_regularizer=l2(reg_layers[0]), input_length=1)
+    MLP_Embedding_Item = Embedding(input_dim=num_items, output_dim=layers[0] // 2, name='item_embedding',
+                                   init=init_normal(), W_regularizer=l2(reg_layers[0]), input_length=1)
+
+    # Crucial to flatten an embedding vector!
+    user_latent = Flatten()(MLP_Embedding_User(user_input))
+    item_latent = Flatten()(MLP_Embedding_Item(item_input))
+
+    # The 0-th layer is the concatenation of embedding layers
+    # vector = merge([user_latent, item_latent], mode = 'concat')
+    vector = merge.concatenate([user_latent, item_latent])
+
+    # MLP layers
+    for idx in range(1, num_layer):
+        layer = Dense(layers[idx], W_regularizer=l2(reg_layers[idx]), activation='relu', name='layer%d' % idx)
+        vector = layer(vector)
+
+    # Final prediction layer
+    prediction = Dense(1, activation='sigmoid', init='lecun_uniform', name='prediction')(vector)
+
+    model = Model(input=[user_input, item_input],
+                  output=prediction)
+
+    return model
+
+
+def get_train_instances(train, num_negatives):
+    user_input, item_input, labels = [], [], []
+    num_users = train.shape[0]
+    for (u, i) in train.keys():
+        # positive instance
+        user_input.append(u)
+        item_input.append(i)
+        labels.append(1)
+        # negative instances
+        for t in range(num_negatives):
+            j = np.random.randint(num_items)
+            # while train.has_key((u, j)):
+            while (u, j) in train:
+                j = np.random.randint(num_items)
+            user_input.append(u)
+            item_input.append(j)
+            labels.append(0)
+    return user_input, item_input, labels
+
+
+if __name__ == '__main__':
+    args = parse_args()
+    path = args.path
+    dataset = args.dataset
+    layers = eval(args.layers)
+    reg_layers = eval(args.reg_layers)
+    num_negatives = args.num_neg
+    learner = args.learner
+    learning_rate = args.lr
+    batch_size = args.batch_size
+    epochs = args.epochs
+    verbose = args.verbose
+
+    topK = 10
+    evaluation_threads = 1  # mp.cpu_count()
+    print("MLP arguments: %s " % (args))
+    model_out_file = 'Pretrain/%s_MLP_%s_%d.h5' % (args.dataset, args.layers, time())
+
+    # Loading data
+    t1 = time()
+    dataset = Dataset(args.path + args.dataset)
+    train, testRatings, testNegatives = dataset.trainMatrix, dataset.testRatings, dataset.testNegatives
+    num_users, num_items = train.shape
+    print("Load data done [%.1f s]. #user=%d, #item=%d, #train=%d, #test=%d"
+          % (time() - t1, num_users, num_items, train.nnz, len(testRatings)))
+
+    # Build model
+    model = get_model(num_users, num_items, layers, reg_layers)
+    if learner.lower() == "adagrad":
+        model.compile(optimizer=Adagrad(lr=learning_rate), loss='binary_crossentropy')
+    elif learner.lower() == "rmsprop":
+        model.compile(optimizer=RMSprop(lr=learning_rate), loss='binary_crossentropy')
+    elif learner.lower() == "adam":
+        model.compile(optimizer=Adam(lr=learning_rate), loss='binary_crossentropy')
+    else:
+        model.compile(optimizer=SGD(lr=learning_rate), loss='binary_crossentropy')
+
+        # Check Init performance
+    t1 = time()
+    (hits, ndcgs) = evaluate_model(model, testRatings, testNegatives, topK, evaluation_threads)
+    hr, ndcg = np.array(hits).mean(), np.array(ndcgs).mean()
+    print('Init: HR = %.4f, NDCG = %.4f [%.1f]' % (hr, ndcg, time() - t1))
+
+    # Train model
+    best_hr, best_ndcg, best_iter = hr, ndcg, -1
+    for epoch in range(epochs):
+        t1 = time()
+        # Generate training instances
+        user_input, item_input, labels = get_train_instances(train, num_negatives)
+
+        # Training        
+        hist = model.fit([np.array(user_input), np.array(item_input)],  # input
+                         np.array(labels),  # labels
+                         batch_size=batch_size, nb_epoch=1, verbose=0, shuffle=True)
+        t2 = time()
+
+        # Evaluation
+        if epoch % verbose == 0:
+            (hits, ndcgs) = evaluate_model(model, testRatings, testNegatives, topK, evaluation_threads)
+            hr, ndcg, loss = np.array(hits).mean(), np.array(ndcgs).mean(), hist.history['loss'][0]
+            print('Iteration %d [%.1f s]: HR = %.4f, NDCG = %.4f, loss = %.4f [%.1f s]'
+                  % (epoch, t2 - t1, hr, ndcg, loss, time() - t2))
+            if hr > best_hr:
+                best_hr, best_ndcg, best_iter = hr, ndcg, epoch
+                if args.out > 0:
+                    model.save_weights(model_out_file, overwrite=True)
+
+    print("End. Best Iteration %d:  HR = %.4f, NDCG = %.4f. " % (best_iter, best_hr, best_ndcg))
+    if args.out > 0:
+        print("The best MLP model is saved to %s" % (model_out_file))
diff --git a/TF_KERAS_NeuMF.py b/TF_KERAS_NeuMF.py
new file mode 100644
index 0000000..101c689
--- /dev/null
+++ b/TF_KERAS_NeuMF.py
@@ -0,0 +1,249 @@
+# coding=utf-8
+
+'''
+Created on Aug 9, 2016
+Keras Implementation of Neural Matrix Factorization (NeuMF) recommender model in:
+He Xiangnan et al. Neural Collaborative Filtering. In WWW 2017.  
+
+@author: Xiangnan He (xiangnanhe@gmail.com)
+'''
+
+# This TF_KERAS version is tested with Keras 2.3.1 and tensorflow-gpu 1.15.2
+
+# use GPU 0
+# import os
+# os.environ["CUDA_VISIBLE_DEVICES"] = "0"
+
+import numpy as np
+from keras import initializers
+from keras.regularizers import l2
+from keras.models import Model
+from keras.layers import Embedding, Input, Dense, merge, Flatten
+from keras.optimizers import Adagrad, Adam, SGD, RMSprop
+from evaluate import evaluate_model
+from Dataset import Dataset
+from time import time
+import argparse
+
+import TF_KERAS_GMF
+import TF_KERAS_MLP
+
+
+#################### Arguments ####################
+def parse_args():
+    parser = argparse.ArgumentParser(description="Run NeuMF.")
+    parser.add_argument('--path', nargs='?', default='Data/',
+                        help='Input data path.')
+    parser.add_argument('--dataset', nargs='?', default='ml-1m',
+                        help='Choose a dataset.')
+    parser.add_argument('--epochs', type=int, default=100,
+                        help='Number of epochs.')
+    parser.add_argument('--batch_size', type=int, default=256,
+                        help='Batch size.')
+    parser.add_argument('--num_factors', type=int, default=8,
+                        help='Embedding size of MF model.')
+    parser.add_argument('--layers', nargs='?', default='[64,32,16,8]',
+                        help="MLP layers. Note that the first layer is the concatenation of user and item embeddings. So layers[0]/2 is the embedding size.")
+    parser.add_argument('--reg_mf', type=float, default=0,
+                        help='Regularization for MF embeddings.')
+    parser.add_argument('--reg_layers', nargs='?', default='[0,0,0,0]',
+                        help="Regularization for each MLP layer. reg_layers[0] is the regularization for embeddings.")
+    parser.add_argument('--num_neg', type=int, default=4,
+                        help='Number of negative instances to pair with a positive instance.')
+    parser.add_argument('--lr', type=float, default=0.001,
+                        help='Learning rate.')
+    parser.add_argument('--learner', nargs='?', default='adam',
+                        help='Specify an optimizer: adagrad, adam, rmsprop, sgd')
+    parser.add_argument('--verbose', type=int, default=1,
+                        help='Show performance per X iterations')
+    parser.add_argument('--out', type=int, default=1,
+                        help='Whether to save the trained model.')
+    parser.add_argument('--mf_pretrain', nargs='?', default='',
+                        help='Specify the pretrain model file for MF part. If empty, no pretrain will be used')
+    parser.add_argument('--mlp_pretrain', nargs='?', default='',
+                        help='Specify the pretrain model file for MLP part. If empty, no pretrain will be used')
+    return parser.parse_args()
+
+
+def init_normal():
+    return initializers.normal(mean=0.0, stddev=0.01)
+
+
+def get_model(num_users, num_items, mf_dim=10, layers=[10], reg_layers=[0], reg_mf=0):
+    assert len(layers) == len(reg_layers)
+    num_layer = len(layers)  # Number of layers in the MLP
+    # Input variables
+    user_input = Input(shape=(1,), dtype='int32', name='user_input')
+    item_input = Input(shape=(1,), dtype='int32', name='item_input')
+
+    # Embedding layer
+    MF_Embedding_User = Embedding(input_dim=num_users, output_dim=mf_dim, name='mf_embedding_user',
+                                  init=init_normal(), W_regularizer=l2(reg_mf), input_length=1)
+    MF_Embedding_Item = Embedding(input_dim=num_items, output_dim=mf_dim, name='mf_embedding_item',
+                                  init=init_normal(), W_regularizer=l2(reg_mf), input_length=1)
+
+    MLP_Embedding_User = Embedding(input_dim=num_users, output_dim=layers[0] // 2, name="mlp_embedding_user",
+                                   init=init_normal(), W_regularizer=l2(reg_layers[0]), input_length=1)
+    MLP_Embedding_Item = Embedding(input_dim=num_items, output_dim=layers[0] // 2, name='mlp_embedding_item',
+                                   init=init_normal(), W_regularizer=l2(reg_layers[0]), input_length=1)
+
+    # MF part
+    mf_user_latent = Flatten()(MF_Embedding_User(user_input))
+    mf_item_latent = Flatten()(MF_Embedding_Item(item_input))
+    # mf_vector = merge([mf_user_latent, mf_item_latent], mode = 'mul') # element-wise multiply
+    mf_vector = merge.multiply([mf_user_latent, mf_item_latent])  # element-wise multiply
+
+    # MLP part 
+    mlp_user_latent = Flatten()(MLP_Embedding_User(user_input))
+    mlp_item_latent = Flatten()(MLP_Embedding_Item(item_input))
+    # mlp_vector = merge([mlp_user_latent, mlp_item_latent], mode = 'concat')
+    mlp_vector = merge.concatenate([mlp_user_latent, mlp_item_latent])
+    for idx in range(1, num_layer):
+        layer = Dense(layers[idx], W_regularizer=l2(reg_layers[idx]), activation='relu', name="layer%d" % idx)
+        mlp_vector = layer(mlp_vector)
+
+    # Concatenate MF and MLP parts
+    # mf_vector = Lambda(lambda x: x * alpha)(mf_vector)
+    # mlp_vector = Lambda(lambda x : x * (1-alpha))(mlp_vector)
+    # predict_vector = merge([mf_vector, mlp_vector], mode = 'concat')
+    predict_vector = merge.concatenate([mf_vector, mlp_vector])
+
+    # Final prediction layer
+    prediction = Dense(1, activation='sigmoid', init='lecun_uniform', name="prediction")(predict_vector)
+
+    model = Model(input=[user_input, item_input],
+                  output=prediction)
+
+    return model
+
+
+def load_pretrain_model(model, gmf_model, mlp_model, num_layers):
+    # MF embeddings
+    gmf_user_embeddings = gmf_model.get_layer('user_embedding').get_weights()
+    gmf_item_embeddings = gmf_model.get_layer('item_embedding').get_weights()
+    model.get_layer('mf_embedding_user').set_weights(gmf_user_embeddings)
+    model.get_layer('mf_embedding_item').set_weights(gmf_item_embeddings)
+
+    # MLP embeddings
+    mlp_user_embeddings = mlp_model.get_layer('user_embedding').get_weights()
+    mlp_item_embeddings = mlp_model.get_layer('item_embedding').get_weights()
+    model.get_layer('mlp_embedding_user').set_weights(mlp_user_embeddings)
+    model.get_layer('mlp_embedding_item').set_weights(mlp_item_embeddings)
+
+    # MLP layers
+    for i in range(1, num_layers):
+        mlp_layer_weights = mlp_model.get_layer('layer%d' % i).get_weights()
+        model.get_layer('layer%d' % i).set_weights(mlp_layer_weights)
+
+    # Prediction weights
+    gmf_prediction = gmf_model.get_layer('prediction').get_weights()
+    mlp_prediction = mlp_model.get_layer('prediction').get_weights()
+    new_weights = np.concatenate((gmf_prediction[0], mlp_prediction[0]), axis=0)
+    new_b = gmf_prediction[1] + mlp_prediction[1]
+    model.get_layer('prediction').set_weights([0.5 * new_weights, 0.5 * new_b])
+    return model
+
+
+def get_train_instances(train, num_negatives):
+    user_input, item_input, labels = [], [], []
+    num_users = train.shape[0]
+    for (u, i) in train.keys():
+        # positive instance
+        user_input.append(u)
+        item_input.append(i)
+        labels.append(1)
+        # negative instances
+        for t in range(num_negatives):
+            j = np.random.randint(num_items)
+            # while train.has_key((u, j)):
+
+            while (u, j) in train:
+                j = np.random.randint(num_items)
+            user_input.append(u)
+            item_input.append(j)
+            labels.append(0)
+    return user_input, item_input, labels
+
+
+if __name__ == '__main__':
+    args = parse_args()
+    num_epochs = args.epochs
+    batch_size = args.batch_size
+    mf_dim = args.num_factors
+    layers = eval(args.layers)
+    reg_mf = args.reg_mf
+    reg_layers = eval(args.reg_layers)
+    num_negatives = args.num_neg
+    learning_rate = args.lr
+    learner = args.learner
+    verbose = args.verbose
+    mf_pretrain = args.mf_pretrain
+    mlp_pretrain = args.mlp_pretrain
+
+    topK = 10
+    evaluation_threads = 1  # mp.cpu_count()
+    print("NeuMF arguments: %s " % (args))
+    model_out_file = 'Pretrain/%s_NeuMF_%d_%s_%d.h5' % (args.dataset, mf_dim, args.layers, time())
+
+    # Loading data
+    t1 = time()
+    dataset = Dataset(args.path + args.dataset)
+    train, testRatings, testNegatives = dataset.trainMatrix, dataset.testRatings, dataset.testNegatives
+    num_users, num_items = train.shape
+    print("Load data done [%.1f s]. #user=%d, #item=%d, #train=%d, #test=%d"
+          % (time() - t1, num_users, num_items, train.nnz, len(testRatings)))
+
+    # Build model
+    model = get_model(num_users, num_items, mf_dim, layers, reg_layers, reg_mf)
+    if learner.lower() == "adagrad":
+        model.compile(optimizer=Adagrad(lr=learning_rate), loss='binary_crossentropy')
+    elif learner.lower() == "rmsprop":
+        model.compile(optimizer=RMSprop(lr=learning_rate), loss='binary_crossentropy')
+    elif learner.lower() == "adam":
+        model.compile(optimizer=Adam(lr=learning_rate), loss='binary_crossentropy')
+    else:
+        model.compile(optimizer=SGD(lr=learning_rate), loss='binary_crossentropy')
+
+    # Load pretrain model
+    if mf_pretrain != '' and mlp_pretrain != '':
+        gmf_model = TF_KERAS_GMF.get_model(num_users, num_items, mf_dim)
+        gmf_model.load_weights(mf_pretrain)
+        mlp_model = TF_KERAS_MLP.get_model(num_users, num_items, layers, reg_layers)
+        mlp_model.load_weights(mlp_pretrain)
+        model = load_pretrain_model(model, gmf_model, mlp_model, len(layers))
+        print("Load pretrained GMF (%s) and MLP (%s) models done. " % (mf_pretrain, mlp_pretrain))
+
+    # Init performance
+    (hits, ndcgs) = evaluate_model(model, testRatings, testNegatives, topK, evaluation_threads)
+    hr, ndcg = np.array(hits).mean(), np.array(ndcgs).mean()
+    print('Init: HR = %.4f, NDCG = %.4f' % (hr, ndcg))
+    best_hr, best_ndcg, best_iter = hr, ndcg, -1
+    if args.out > 0:
+        model.save_weights(model_out_file, overwrite=True)
+
+        # Training model
+    for epoch in range(num_epochs):
+        t1 = time()
+        # Generate training instances
+        user_input, item_input, labels = get_train_instances(train, num_negatives)
+
+        # Training
+        hist = model.fit([np.array(user_input), np.array(item_input)],  # input
+                         np.array(labels),  # labels
+                         batch_size=batch_size, nb_epoch=1, verbose=0, shuffle=True)
+        t2 = time()
+
+        # Evaluation
+        if epoch % verbose == 0:
+            (hits, ndcgs) = evaluate_model(model, testRatings, testNegatives, topK, evaluation_threads)
+            hr, ndcg, loss = np.array(hits).mean(), np.array(ndcgs).mean(), hist.history['loss'][0]
+            print('Iteration %d [%.1f s]: HR = %.4f, NDCG = %.4f, loss = %.4f [%.1f s]'
+                  % (epoch, t2 - t1, hr, ndcg, loss, time() - t2))
+            if hr > best_hr:
+                best_hr, best_ndcg, best_iter = hr, ndcg, epoch
+                if args.out > 0:
+                    model.save_weights(model_out_file, overwrite=True)
+
+    print("End. Best Iteration %d:  HR = %.4f, NDCG = %.4f. " % (best_iter, best_hr, best_ndcg))
+    if args.out > 0:
+        print("The best NeuMF model is saved to %s" % (model_out_file))
diff --git a/evaluate.py b/evaluate.py
index 729f07a..432846e 100644
--- a/evaluate.py
+++ b/evaluate.py
@@ -44,7 +44,7 @@ def evaluate_model(model, testRatings, testNegatives, K, num_thread):
         ndcgs = [r[1] for r in res]
         return (hits, ndcgs)
     # Single thread
-    for idx in xrange(len(_testRatings)):
+    for idx in range(len(_testRatings)):
         (hr,ndcg) = eval_one_rating(idx)
         hits.append(hr)
         ndcgs.append(ndcg)      
@@ -61,7 +61,7 @@ def eval_one_rating(idx):
     users = np.full(len(items), u, dtype = 'int32')
     predictions = _model.predict([users, np.array(items)], 
                                  batch_size=100, verbose=0)
-    for i in xrange(len(items)):
+    for i in range(len(items)):
         item = items[i]
         map_item_score[item] = predictions[i]
     items.pop()
@@ -79,7 +79,7 @@ def getHitRatio(ranklist, gtItem):
     return 0
 
 def getNDCG(ranklist, gtItem):
-    for i in xrange(len(ranklist)):
+    for i in range(len(ranklist)):
         item = ranklist[i]
         if item == gtItem:
             return math.log(2) / math.log(i+2)