import sys
print("Python %d.%d.%d" % (sys.version_info.major, \
                           sys.version_info.minor, \
                           sys.version_info.micro))

import numpy as np
print("NumPy %s" % np.__version__)

import scipy
import scipy.io as sio
from scipy.optimize import fmin_cg
print("SciPy %s" % scipy.__version__)

import matplotlib
import matplotlib.pyplot as plt
print("matplotlib %s" % matplotlib.__version__)

"""
Explanation: Part of iPyMacLern project.

Copyright (C) 2015 by Eka A. Kurniawan
eka.a.kurniawan(ta)gmail(tod)com
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program.  If not, see http://www.gnu.org/licenses/.

Tested On
End of explanation
"""


import time

"""
Explanation: Import Modules
End of explanation
"""


# Display graph inline
%matplotlib inline

# Display graph in 'retina' format for Mac with retina display. Others, use PNG or SVG format.
%config InlineBackend.figure_format = 'retina'
#%config InlineBackend.figure_format = 'PNG'
#%config InlineBackend.figure_format = 'SVG'

# For displaying 3D graph
from mpl_toolkits.mplot3d import Axes3D

"""
Explanation: Display Settings
End of explanation
"""


movies = sio.loadmat('ex8_movies.mat')
# Rating: 1 to 5 stars
# R_ij: rating of movie i by user j
R = movies['R']
# Rating flag
# Y_ij = 0: move i has not been rated by user j
# Y_ij = 1: move i has been rated by user j
Y = movies['Y']

"""
Explanation: Collaborative Filtering[1]
We are going to implement collaborative filtering algorithm for making movie recommendation based on previously rated movies by other users.
Get and Plot Data
End of explanation
"""


np.mean(Y[0, R[0] == 1])

fig = plt.figure()
ax = fig.add_subplot(111)

imshow = ax.imshow(Y, aspect='auto')

plt.title('Movie Rating by Users (1 to 5 stars)\n', linespacing=1)
plt.xlabel('Users')
plt.ylabel('Movies')
plt.colorbar(imshow)

fig.show()

"""
Explanation: Based on movie_ids.txt file, Toy Story (1995) movie is on the list no 1 or index 0.
$$index = list - 1$$
So the rating mean is as follow.
End of explanation
"""


def load_movie_list(filename):
    movie_list = []
    with open(filename, 'r') as input_fh:
        for line in input_fh.readlines():
            space_idx = line.find(' ')
            if space_idx:
                movie_list.append(line[space_idx + 1:].strip())
    return movie_list

movie_list = load_movie_list('movie_ids.txt')

"""
Explanation: Movie List
End of explanation
"""


def normalize_rating(Y, R):
    [m, n] = Y.shape
    Y_mean = np.zeros([m, 1])
    Y_norm = np.zeros(Y.shape)
    for i in range(m):
        Y_mean[i] = np.mean(Y[i, R[i] == 1])
    Y_norm = Y - Y_mean
    return Y_mean, Y_norm

Y_mean, Y_norm = normalize_rating(Y, R)

Y_mean

Y_norm

"""
Explanation: Normalize Function
Normalize rating mean for each movie to 0.
End of explanation
"""


def compute_cost(X_Theta, args):
    Y, R, Lambda, n_m, n_u, n = args
    X_Theta = X_Theta.reshape([n_m + n_u, n])
    X = X_Theta[:n_m,:]
    Theta = X_Theta[n_m:,:]
    
    error = (np.dot(X, Theta.transpose()) - Y) * R
    J = (np.sum(error ** 2) / 2) + \
        ((Lambda * np.sum(Theta ** 2)) / 2) + \
        ((Lambda * np.sum(X ** 2)) / 2)
    return J

def compute_cost_wrapper(X_Theta):
    return compute_cost(X_Theta, args)

"""
Explanation: Cost and Gradient Functions
Collaborative filtering cost function:
$$J(x^{(1)},...,x^{(n_m)},\theta^{(1)},...,\theta^{(n_u)}) =
\frac{1}{2} \sum_{(i,j):r(i,j)=1} \left((\theta^{(j)})^T x^{(i)} - y^{(i,j)}\right)^2 + 
\frac{\lambda}{2} \sum_{i=1}^{n_m} \sum_{k=1}^{n} (x_k^{(i)})^2 +
\frac{\lambda}{2} \sum_{j=1}^{n_u} \sum_{k=1}^{n} (\theta_k^{(j)})^2$$
End of explanation
"""


def compute_grad(X_Theta, args):
    Y, R, Lambda, n_m, n_u, n = args
    X_Theta = X_Theta.reshape([n_m + n_u, n])
    X = X_Theta[:n_m,:]
    Theta = X_Theta[n_m:,:]
    
    error = (np.dot(X, Theta.transpose()) - Y) * R
    X_grad = np.dot(error, Theta) + (Lambda * X)
    Theta_grad = np.dot(error.transpose(), X) + (Lambda * Theta)
    return np.vstack([X_grad, Theta_grad]).flatten()

def compute_grad_wrapper(X_Theta):
    return compute_grad(X_Theta, args)

dummy_movies = sio.loadmat('ex8_movieParams.mat')

dummy_num_users = 4
dummy_num_movies = 5
dummy_num_features = 3

dummy_X = dummy_movies['X'][:dummy_num_movies, :dummy_num_features]
dummy_Theta = dummy_movies['Theta'][:dummy_num_users, :dummy_num_features]
dummy_Y = Y[:dummy_num_movies, :dummy_num_users]
dummy_R = R[:dummy_num_movies, :dummy_num_users]

Lambda = 1.5
X_Theta = np.vstack([dummy_X, dummy_Theta]).flatten()
args = np.asarray((dummy_Y, dummy_R, Lambda, dummy_num_movies, dummy_num_users, dummy_num_features))
J = compute_cost(X_Theta, args)

J

Lambda = 1.5
X_Theta = np.vstack([dummy_X, dummy_Theta])
args = np.asarray((dummy_Y, dummy_R, Lambda, dummy_num_movies, dummy_num_users, dummy_num_features))
grad = compute_grad(X_Theta, args)

grad

"""
Explanation: Collaborative filtering gradient function:
$$x_k^{(i)} := x_k^{(i)} - \alpha \left(\sum_{j:r(i,j)=1} \left((\theta^{(j)})^T x^{(i)} - y^{(i,j)}\right) \theta_k^{(i)} + \lambda x_k^{(i)} \right)$$
$$\theta_k^{(j)} := \theta_k^{(j)} - \alpha \left(\sum_{i:r(i,j)=1} \left((\theta^{(j)})^T x^{(i)} - y^{(i,j)}\right) x_k^{(i)} + \lambda \theta_k^{(j)} \right)$$
End of explanation
"""


# Get total movies (n_m) and total users (n_u)
n_m, n_u = Y.shape

my_ratings = np.zeros([n_m, 1])

my_ratings[0]   = 4
my_ratings[97]  = 2
my_ratings[6]   = 3
my_ratings[11]  = 5
my_ratings[53]  = 4
my_ratings[63]  = 5
my_ratings[65]  = 3
my_ratings[68]  = 5
my_ratings[182] = 4
my_ratings[225] = 5
my_ratings[354] = 5

for i, rating in enumerate(my_ratings):
    if rating:
        print('Rated %5s for %s' % ('*' * int(rating), movie_list[i]))

my_Y = np.hstack([my_ratings, Y])
my_R = np.hstack([(my_ratings != 0).astype(int), R])

# Normalize
my_Y_mean, my_Y_norm = normalize_rating(my_Y, my_R)

"""
Explanation: Rate Movies
End of explanation
"""


# Get total movies (n_m) and total users (n_u)
n_m, n_u = my_Y.shape

# Total parameters
n = 100

# Initialize random parameters X and Theta
X_init = np.random.rand(n_m, n)
Theta_init = np.random.rand(n_u, n)
X_Theta_init = np.vstack([X_init, Theta_init]).flatten()

# Regularization
Lambda = 1.5

# Arguments
args = np.asarray((my_Y_norm, my_R, Lambda, n_m, n_u, n))

tic = time.time()
X_Theta_result = fmin_cg(compute_cost_wrapper, X_Theta_init, fprime=compute_grad_wrapper)
toc = time.time()
print('Runtime: %s seconds' % int(toc - tic))

X_Theta_init

X_Theta_result = X_Theta_result.reshape([n_m + n_u, n])
X_result = X_Theta_result[:n_m,:]
Theta_result = X_Theta_result[n_m:,:]

X_result

Theta_result

# Prediction
p = np.dot(X_result, Theta_result.transpose()) + Y_mean
p = np.rint(((p - np.min(p)) / (np.abs(np.min(p)) + np.abs(np.max(p)))) * 4) + 1
my_predictions = p[:,1].reshape(p.shape[0],1)
sorted_my_predictions = np.argsort(my_predictions, axis=0)[::-1]

# Display top ten recommendation
for i in sorted_my_predictions[:10]:
    print('Rated %5s for %s' % (('*' * int(my_predictions[i])), movie_list[i]))

"""
Explanation: Learn Rating
End of explanation
"""


Y = np.asarray([[5, 5, 0, 0],
                [5, 0, 0, 0],
                [0, 4, 0, 0],
                [0, 0, 5, 4],
                [0, 0, 5, 0]])

R = np.asarray([[1, 1, 1, 1],
                [1, 0, 0, 1],
                [0, 1, 1, 0],
                [1, 1, 1, 1],
                [1, 1, 1, 0]])

movie_list = ['Love at Last',
              'Romance Forever',
              'Cute Puppies of Love',
              'Nonstop Car Chases',
              'Swords vs Karate']

def normalize_rating(Y, R):
    [m, n] = Y.shape
    Y_mean = np.zeros([m, 1])
    Y_norm = np.zeros(Y.shape)
    for i in range(m):
        Y_mean[i] = np.mean(Y[i, R[i] == 1])
    Y_norm = Y - Y_mean
    return Y_mean, Y_norm

# Normalize
Y_mean, Y_norm = normalize_rating(Y, R)

# Get total movies (n_m) and total users (n_u)
n_m, n_u = Y.shape

# Total parameters
n = 100

# Initialize random parameters X and Theta
X_init = np.random.rand(n_m, n)
Theta_init = np.random.rand(n_u, n)
X_Theta_init = np.vstack([X_init, Theta_init]).flatten()

# Regularization
Lambda = 1.5

# Arguments
args = np.asarray((Y_norm, R, Lambda, n_m, n_u, n))

X_Theta_result = fmin_cg(compute_cost_wrapper, X_Theta_init, fprime=compute_grad_wrapper)

X_Theta_result = X_Theta_result.reshape([n_m + n_u, n])
X_result = X_Theta_result[:n_m,:]
Theta_result = X_Theta_result[n_m:,:]

# Prediction
p = np.dot(X_result, Theta_result.transpose()) + Y_mean
p = np.rint(((p - np.min(p)) / (np.abs(np.min(p)) + np.abs(np.max(p)))) * 5)

p

"""
Explanation: Testing with Small Sample
End of explanation
"""


