
function pdlf_neural
    % This function implements so-called PDLF-Neural algorithm.
    % The purpose is to construct a model for a movie recommendation
    % system. To run it, you have to unzip the attached Matlab datasets
    % and place pdlf_neural.m in the same folder. Specify your parameters 
    % of choice (eg, number of hidden units for neural network) in pdlf_neural.m 
    % and type pdlf_neural in Matlab command line.
    %
    % Created by: Nik Tuzov, Purdue Computational Finance Program. 
    %             Two ROC functions borrowed from Dr Gavin C. Cawley.
    % Contact info: ntuzov@purdue.edu or http://www.linkedin.com/in/ntuzov

    %%%%%%%% Parameter section (specify your own values): %%%%%%%%%%%%%%%%%%

    nn_repl = 3;       % number of n/network replications (size of "forest")
    num_units = 10;    % number of hidden units in n/network
    max_iter = 10;     % total number of iterations for PDLF-Neural algorithm
    epochs = 30;       % number of epochs for n/network

    %%%%%%%%%%% End of parameter section. DO NOT modify the code below %%%%%%

    val_rocmax = 0;    % maximal validation ROC
    T = 14903; % training sample size
    V = 5700;  % validation sample size 

    K = 3;   % number of row clusters
    L = 3;   % number of column clusters
    M = 346; % the number of distinct users 
    N = 966; % the number of distinct movies
             % it applies to the  whole sample and to the train sample    

    % indices of observations (ratings) :
    i = 1; % user, i = 1,M
    j = 1; % movie, j = 1,N

    % row and column cluster indices:
    I = 1; % row, I = 1,K
    J = 1; % column, J = 1,L

    % randomly initialize row and column cluster distributions:
    rho = random('unid', K, 1, M);   % row 
    gma = random('unid', L, 1, N); % column 

    % load the data:
    load train_sample.mat; % 14903 * 57
    load valid_sample.mat; % 5700 * 57
    load train_matrix.mat; % M * N
    load valid_matrix.mat; % M * N
    load index_matrix.mat; % M * N

    % Patterns and targets for neural network:
    train_pat = train_sample(:, 6:55)'; % 50 * 14903
    train_targ = train_sample(:, 5)';   % 1 * 14903
    valid_pat = valid_sample(:, 6:55)'; % 50 * 5700
    valid_targ = valid_sample(:, 5)';   % 1 * 5700
    clear train_sample valid_sample;

    % to monitor deltas and cluster sizes:
    deltas_last = zeros(K, L); 
    deltas_curr = zeros(K, L);
    clust_count = zeros(K, L);

    % Cluster structure
    for I = 1 : K
        for J = 1 : L
            clusters(I, J).delta = 0; % offset       
            clusters(I, J).count = 0; % number of points in the cluster
            clusters(I, J).obs_fit(1, :) = [0 0 0 0];  % [i j observed fitted] values        
        end
    end

    %%%%%%%%%%%%%%  Start PDLF-Neural Algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    tic

    fprintf('\n Initialization... \n');

    % Initialize the fitted values of Y:
    [train_fit, valid_fit] = fit_network(train_pat, train_targ, valid_pat, valid_targ, num_units, nn_repl, epochs); 
    % [1 * 14903, 1 * 5700]

    for num_iter = 1 : max_iter

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Implement Step1: %%%%%%%%%%%%%%%
        fprintf('\n Iteration %g, Step 1. Clusters: %g * %g \n', num_iter, K, L);

        % Reset clusters:
        for I = 1 : K
            for J = 1 : L             
                clusters(I, J).count = 0; % number of points in the cluster
                clusters(I, J).obs_fit(1, :) = [0 0 0 0];  % [i j observed fitted] values        
            end
        end

        % Go through the training matrix and fill the cluster structure:
        for i = 1 : M
            for j = 1 : N
                if (~isnan(train_matrix(i, j)))             
                   obs_ij = train_targ(index_matrix(i, j)); % observed reponse         
                   fit_ij = train_fit(index_matrix(i, j));  % fitted n/network response          
                   I = rho(i); % these two were just updated in Step 3 and 4
                   J = gma(j);
                   clusters(I, J).count = clusters(I, J).count + 1;
                   clusters(I, J).obs_fit(clusters(I, J).count, :) = [i  j  obs_ij  fit_ij];             
                end
            end
        end

        % Save old deltas and recalculate deltas: 
        for I = 1 : K
            for J = 1 : L
                % Important: always use count because of redundant "tail"
                % observations
                count = clusters(I, J).count;
                deltas_last(I, J) = clusters(I, J).delta;
                resid = clusters(I, J).obs_fit(1 : count, 3) - clusters(I, J).obs_fit(1 : count, 4);            
                clusters(I, J).delta = mean(resid);
                deltas_curr(I, J) = clusters(I, J).delta;
                clust_count(I, J) = count;            
            end
        end

        %%%%%%%%%%%%%%%%%% Implement Step2 %%%%%%%%%%%%%%%%%%%%%%%%%%%

        fprintf('\n Iteration %g, Step 2. Number of n/n replications = %g, Num of units = %g, Num of epochs = %g \n',...
                num_iter, nn_repl, num_units, epochs);
        % We have to calculate updated targets for both training and validation
        % target vectors. train_targ_updt = train_targ - corresponding deltas

        train_targ_updt = train_targ * 0; % initialize, 1 * 14903 
        valid_targ_updt = valid_targ * 0; % initialize, 1 * 5700

        for i = 1 : M
            for j = 1 : N
                % training:            
                if (~isnan(train_matrix(i, j))) % (i, j) is present in the train d/s           
                   obs_ij = train_targ(index_matrix(i, j)); % observed
                   delta_ij = clusters(rho(i), gma(j)).delta;
                   train_targ_updt(index_matrix(i, j)) = obs_ij - delta_ij; % can be < 0 or >1
                   if train_targ_updt(index_matrix(i, j)) > 1  train_targ_updt(index_matrix(i, j)) = 1; end
                   if train_targ_updt(index_matrix(i, j)) < 0  train_targ_updt(index_matrix(i, j)) = 0; end
                end
                % validation:             
                if (~isnan(valid_matrix(i, j))) % (i, j) is present in the validation           
                   obs_ij = valid_targ(index_matrix(i, j)); % observed
                   delta_ij = clusters(rho(i), gma(j)).delta;
                   valid_targ_updt(index_matrix(i, j)) = obs_ij - delta_ij; % can be < 0 or >1
                   if valid_targ_updt(index_matrix(i, j)) > 1  valid_targ_updt(index_matrix(i, j)) = 1; end
                   if valid_targ_updt(index_matrix(i, j)) < 0  valid_targ_updt(index_matrix(i, j)) = 1; end
                end

            end
        end    

        % Train n/w with the new targets: 
        [train_fit, valid_fit] = fit_network(train_pat, train_targ_updt, valid_pat, valid_targ_updt,...
                                             num_units, nn_repl, epochs);

        %%%%%%%%% Implement Step3 %%%%%%%%%%%%%%%%%%%%%%%%%

        fprintf('\n Iteration %g, Step 3 \n', num_iter);

        % save old cluster assignments:
        rho_last = rho;
        gma_last = gma;

        % Update row cluster assignments: 
        for i = 1 : M    
            sse_min = +Inf;
            I_opt = 0;
            for I = 1 : K
                sse_I = 0; % sse for row assignment I            
                for j = 1 : N                 
                    if (~isnan(train_matrix(i, j)))                % (i, j) is present in the train d/s           
                       obs_ij = train_targ(index_matrix(i, j));    % observed value
                       delta_Ij = clusters(I, gma(j)).delta;       % offset from Step 1 for row assignment = I
                       fit_ij = train_fit(index_matrix(i, j));     % fitted n/n response from Step 2
                       totfit_ij = fit_ij + delta_Ij;              % total fitted value
                       if totfit_ij > 1 totfit_ij = 1; end
                       if totfit_ij < 0 totfit_ij = 0; end                  
                       sse_I = sse_I + (obs_ij - totfit_ij) ^ 2;                  
                    end
                end % for j         
                if sse_I < sse_min
                    sse_min = sse_I;
                    I_opt = I;
                end
            end  % for I
            rho(i) = I_opt; 
        end

         %%%%%%%%%%%%%%%%%% Implement Step 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

         fprintf('\n Iteration %g, Step 4 \n', num_iter);

         % Update column cluster assignments: 
        for j = 1 : N
            sse_min = +Inf;
            J_opt = 0;
            for J = 1 : L
                sse_J = 0; % sse for column assignment J               
                for i = 1 : M                
                    if (~isnan(train_matrix(i, j)))                % (i, j) is present in the train d/s           
                       obs_ij = train_targ(index_matrix(i, j));    % observed value
                       delta_iJ = clusters(rho(i), J).delta;       % offset from Step 1 for row assignment = I
                       fit_ij = train_fit(index_matrix(i, j));     % fitted response from Step 2 
                       totfit_ij = fit_ij + delta_iJ;              % total fitted value
                       if totfit_ij > 1 totfit_ij = 1; end
                       if totfit_ij < 0 totfit_ij = 0; end    
                       sse_J = sse_J + (obs_ij - totfit_ij) ^ 2;                     
                    end
                end % for i                           
                if sse_J < sse_min
                    sse_min = sse_J;
                    J_opt = J;
                end
            end  % for J
            gma(j) = J_opt; 
        end

        %%%%%%%%%%%%%%%%%%%% End of Step 4 %%%%%%%%%%%%%%%%%%%%%%%%%%

        % Compute the gloabal training and test errors:

        train_totfit = train_fit * 0; % Total training fitted values
        valid_totfit = valid_fit * 0; % Total validation fitted values

        for i = 1 : M
            for j = 1 : N
                % training:            
                if (~isnan(train_matrix(i, j))) % (i, j) is present in the train d/s     
                   delta_ij = clusters(rho(i), gma(j)).delta;
                   train_totfit(index_matrix(i, j)) = train_fit(index_matrix(i, j)) + delta_ij; % can be < 0 or >1
                   if train_totfit(index_matrix(i, j)) > 1  train_totfit(index_matrix(i, j)) = 1; end
                   if train_totfit(index_matrix(i, j)) < 0  train_totfit(index_matrix(i, j)) = 0; end
                end

                % validation:            
                if (~isnan(valid_matrix(i, j))) % (i, j) is present in the validation           
                   delta_ij = clusters(rho(i), gma(j)).delta;
                   valid_totfit(index_matrix(i, j)) = valid_fit(index_matrix(i, j)) + delta_ij; % can be < 0 or >1    
                   if valid_totfit(index_matrix(i, j)) > 1  valid_totfit(index_matrix(i, j)) = 1; end
                   if valid_totfit(index_matrix(i, j)) < 0  valid_totfit(index_matrix(i, j)) = 0; end
                end
            end
        end

        train_mse = mean((train_targ - train_totfit) .^ 2);
        train_mse_sterr = std((train_targ - train_totfit) .^ 2) / sqrt(T);
        valid_mse = mean((valid_targ - valid_totfit) .^ 2);
        valid_mse_sterr = std((valid_targ - valid_totfit) .^ 2) / sqrt(V);

        t = train_targ'; % true response
        y = train_totfit'; % fitted probabilities
        [tp, fp] = roc(t, y);
        area_train = auroc(tp, fp);  % area under curve
        t = valid_targ'; % true response
        y = valid_totfit'; % fitted probabilities
        [tp, fp] = roc(t, y);
        area_valid = auroc(tp, fp);  % area under curve    

        if area_valid > val_rocmax
            val_rocmax = area_valid;
            train_rocmax = area_train;
            clusters_opt = clusters;
            rho_opt = rho;
            gma_opt = gma;
            min_cluster = min(min(clust_count));
            max_cluster = max(max(clust_count));
            max_delta =  max(max(deltas_curr));
            min_delta = min(min(deltas_curr));
            num_iter_opt = num_iter;
        end        

        fprintf('\n ###################  RESULTS - Iteration %g ###############\n', num_iter);
        fprintf('\n Training mse = %g (%g) \n', train_mse, train_mse_sterr);
        fprintf('\n Training ROC area = %g \n', area_train);
        fprintf('\n Validation mse = %g (%g) \n', valid_mse, valid_mse_sterr);
        fprintf('\n Validation ROC area = %g \n\n', area_valid);
        fprintf('\n Row change = %g \n', sum(rho_last ~= rho));    % how many users changed row cluster assignments
        fprintf('\n Column change = %g\n', sum(gma_last ~= gma));  % how many movies changed column cluster assignments
        fprintf('\n Max current delta = %g \n', max(max(deltas_curr)));
        fprintf('\n Min current delta = %g \n', min(min(deltas_curr)));    
        fprintf('\n Sum of Squared delta changes = %g \n', sum(sum( (deltas_last - deltas_curr).^2 )) );  
        fprintf('\n Max current cluster size = %g \n', max(max(clust_count)));
        fprintf('\n Min current cluster size = %g \n', min(min(clust_count)));  
        fprintf('\n Max validation ROC so far: %g \n', val_rocmax);
        fprintf('\n ###########################################################\n\n');

    end % for num_iter

    fprintf('\n\n$$$$$$$$$$$$$$$$$$$$$$$  PDLF-Neural fitting is done $$$$$$$$$$$$$$$$$$$$\n');
    fprintf('\n Maximal validation ROC area = %g , iteration %g \n', val_rocmax, num_iter_opt);
    fprintf(' \n$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n\n');

    toc
    %%%%% End of PDLF-Neural Algorithm %%%%%%%%%%%%%%%%%%

end  % pdlf_neural

function [train_fit, valid_fit] = fit_network(train_pat, train_targ, valid_pat, valid_targ, num_units, nn_repl, epochs)
    
    % This function returns the training and validation fitted values from the network;
    
    % About TRAIN.
    % If you call train in sequence the weights will not be re-initialized
    % automatically. However, the step size is re-initialized - that's why in
    % the beginning mse jumps and then goes back. I also tried to train w/o
    % validation and then with validation - after the first step weights
    % don't re-initialize.
   
    % Important: when there's validation and the max_fail isn't reached,
    % the resulting weights correspond to the minimal validation error,
    % not to the last iteration.

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % validation patterns and targets:
    val.P = valid_pat;
    val.T = valid_targ;
    % create the n/w
    net40_rp = newff(minmax(train_pat),[num_units 1],{'logsig' 'logsig'},'trainrp');
    % parameters:
    net40_rp.trainParam.show = NaN;  % set to NaN if don't want to see anything
    net40_rp.trainParam.epochs = epochs;
    net40_rp.trainParam.min_grad = 1e-9;
    net40_rp.trainParam.max_fail = 20;
    net40_rp.trainParam.goal = 0.0;
    % performance function = mse by default

    % forest-based outputs:
    train_fit = train_targ * 0;
    valid_fit = valid_targ * 0;

    fprintf('\n Fitting neural network(s)... \n');

    for i = 1 : nn_repl  
        net40_rp = init(net40_rp);    
        [net40_rp, tr] = train(net40_rp, train_pat, train_targ, [], [], val); 
        train_fit = train_fit + sim(net40_rp, train_pat);
        valid_fit = valid_fit + sim(net40_rp, valid_pat);  
    end
    train_fit = train_fit / nn_repl;
    valid_fit = valid_fit / nn_repl;

end

function A = auroc(tp, fp)
%
% AUROC - area under ROC curve
%
%    An ROC (receiver operator characteristic) curve is a plot of the true
%    positive rate as a function of the false positive rate of a classifier
%    system.  The area under the ROC curve is a reasonable performance
%    statistic for classifier systems assuming no knowledge of the true ratio
%    of misclassification costs.
%
%    A = AUROC(TP, FP) computes the area under the ROC curve, where TP and FP
%    are column vectors defining the ROC or ROCCH curve of a classifier
%    system.
%
%    [1] Fawcett, T., "ROC graphs : Notes and practical
%        considerations for researchers", Technical report, HP
%        Laboratories, MS 1143, 1501 Page Mill Road, Palo Alto
%        CA 94304, USA, April 2004.
%
%    See also : ROC, ROCCH

%
% File        : auroc.m
%
% Date        : Wednesdaay 11th November 2004 
%
% Author      : Dr Gavin C. Cawley
%
% Description : Calculate the area under the ROC curve for a two-class
%               probabilistic classifier.
%
% References  : [1] Fawcett, T., "ROC graphs : Notes and practical
%                   considerations for researchers", Technical report, HP
%                   Laboratories, MS 1143, 1501 Page Mill Road, Palo Alto
%                   CA 94304, USA, April 2004.
%
% History     : 22/03/2001 - v1.00
%               10/11/2004 - v1.01 minor improvements to comments etc.
%
% Copyright   : (c) G. C. Cawley, November 2004.
%
%    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 2 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, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

n = size(tp, 1);
A = sum((fp(2:n) - fp(1:n-1)).*(tp(2:n)+tp(1:n-1)))/2;

% bye bye...

end

function [tp, fp] = roc(t, y)
%
% ROC - generate a receiver operating characteristic curve
%
%    [TP,FP] = ROC(T,Y) gives the true-positive rate (TP) and false positive
%    rate (FP), where Y is a column vector giving the score assigned to each
%    pattern and T indicates the true class (a value above zero represents
%    the positive class and anything else represents the negative class).  To
%    plot the ROC curve,
%
%       PLOT(FP,TP);
%       XLABEL('FALSE POSITIVE RATE');
%       YLABEL('TRUE POSITIVE RATE');
%       TITLE('RECEIVER OPERATING CHARACTERISTIC (ROC)');
%
%    See [1] for further information.
%
%    [1] Fawcett, T., "ROC graphs : Notes and practical
%        considerations for researchers", Technical report, HP
%        Laboratories, MS 1143, 1501 Page Mill Road, Palo Alto
%        CA 94304, USA, April 2004.
%
%    See also : ROCCH, AUROC

%
% File        : roc.m
%
% Date        : Friday 9th June 2005
%
% Author      : Dr Gavin C. Cawley
%
% Description : Generate an ROC curve for a two-class classifier.
%
% References  : [1] Fawcett, T., "ROC graphs : Notes and practical
%                   considerations for researchers", Technical report, HP
%                   Laboratories, MS 1143, 1501 Page Mill Road, Palo Alto
%                   CA 94304, USA, April 2004.
%
% History     : 10/11/2004 - v1.00
%               09/06/2005 - v1.10 - minor recoding
%
% Copyright   : (c) G. C. Cawley, June 2005.
%
%    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 2 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, write to the Free Software
%    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
%

% process targets

t = t > 0;

% sort by classifier output

[Y,idx] = sort(-y);
t       = t(idx);

% compute true positive and false positive rates

tp = cumsum(t)/sum(t);
fp = cumsum(~t)/sum(~t);

% add trivial end-points

tp = [0 ; tp ; 1];
fp = [0 ; fp ; 1];

% bye bye...

end

