%%%% Tutorial on the basic structure of using a planar decision boundary
%%%% to divide a collection of data-points into two classes. 
%%%% by Rajeev Raizada, Jan.2010
%%%%
%%%% Please mail any comments or suggestions to: raizada at cornell dot edu
%%%%
%%%% Probably the best way to look at this program is to read through it
%%%% line by line, and paste each line into the Matlab command window
%%%% in turn. That way, you can see what effect each individual command has.
%%%%
%%%% Alternatively, you can run the program directly by typing 
%%%%
%%%%   classification_plane_tutorial
%%%%
%%%% into your Matlab command window. 
%%%% Do not type ".m" at the end
%%%% If you run the program all at once, all the Figure windows
%%%% will get made at once and will be sitting on top of each other.
%%%% You can move them around to see the ones that are hidden beneath.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%% Let's look at a toy example: classifying people as either 
%%%% sumo wrestlers or basketball players, depending on their height and weight.
%%%% Let's call the x-axis height and the y-axis weight

sumo_wrestlers =     [ 4 8; ...
                       2 6; ...
                       2 2; ...
                       3 5; ...
                       4 7];
             
basketball_players = [ 3 2; ...
                       4 5; ...
                       5 3; ...
                       5 7; ...
                       3 3];
                   

%%% Let's plot this
figure(1);
clf;
set(gca,'FontSize',14);
plot(sumo_wrestlers(:,1),sumo_wrestlers(:,2),'ro','LineWidth',2);
hold on;
plot(basketball_players(:,1),basketball_players(:,2),'bx','LineWidth',2);
axis([0 6 0 10]);
xlabel('Height');
ylabel('Weight');
legend('Sumo wrestlers','Basketball players',2); % The 2 at the end means
                                                 % put the legend in top-left corner

%%%% In order to be able to train a classifier on the input vectors,
%%%% we need to know what the desired output categories are for each one.
%%%% Let's set this to be +1 for sumo wrestlers, and -1 for basketball players

desired_output_sumo =        [ 1; ... % sumo_wrestlers =     [ 4 8; ...
                               1; ... %                        2 6; ...
                               1; ... %                        2 2; ...
                               1; ... %                        3 5; ...
                               1 ];   %                        4 7];              
desired_output_basketball = [ -1; ... % basketball_players = [ 3 2; ...
                              -1; ... %                        4 5; ...
                              -1; ... %                        5 3; ...
                              -1; ... %                        5 7; ...
                              -1  ];  %                        3 3 ];

all_desired_output = [ desired_output_sumo; ...
                       desired_output_basketball ];

%%%%%% We want to find a linear decision boundary,
%%%%%% i.e. simply a straight line, such that all the data points
%%%%%% on one side of the line get classified as sumo wrestlers,
%%%%%% i.e. get mapped onto the desired output of +1,
%%%%%% and all the data points on the other side get classified
%%%%%% as basketball players, i.e. get mapped onto the desired output of -1.
%%%%%%
%%%%%% The equation for a straight line has this form:
%%%%%% weight_vector * data_coords  -  offset_from_origin = 0;
%%%%%% 
%%%%%% We're not so interested for now in the offset_from_origin term,
%%%%%% so we can get rid of that by subtracting the mean from our data,
%%%%%% so that it is all centered around the origin.

%%%%%% Let's stack up the sumo data on top of the bastetball players data
all_data = [ sumo_wrestlers; ...
             basketball_players ];
          
%%%%%% Now let's subtract the mean from the data, 
%%%%%% so that it is all centered around the origin.
%%%%%% Each dimension (height and weight) has its own column.
mean_column_vals = mean(all_data);

%%%%%% To subtract the mean from each column in Matlab,
%%%%%% we need to make a matrix full of column-mean values
%%%%%% that is the same size as the whole data matrix.
matrix_of_mean_vals = ones(size(all_data,1),1) * mean_column_vals;

zero_meaned_data = all_data - matrix_of_mean_vals;

%%%% Now, having gotten rid of that annoying offset_from_origin term,
%%%% we want to find a weight vector which gives us the best solution
%%%% that we can find to this equation:
%%%% zero_meaned_data * weights = all_desired_output;
%%%% But, there is no such perfect set of weights. 
%%%% We can only get a best fit, such that
%%%% zero_meaned_data * weights = all_desired_output + error
%%%% where the error term is as small as possible.
%%%%
%%%% Note that our equation 
%%%% zero_meaned_data * weights = all_desired_output
%%%% 
%%%% has exactly the same form as the equation
%%%% from the tutorial code in 
%%%% http://www.dartmouth.edu/~raj/Matlab/fMRI/design_matrix_tutorial.m
%%%% which is:
%%%% Design matrix * sensitivity vector = Voxel response 
%%%%
%%%% The way we solve the equation is exactly the same, too.
%%%% If we could find a matrix-inverse of the data matrix,
%%%% then we could pre-multiply both sides by that inverse,
%%%% and that would give us the weights:
%%%%
%%%% inv(zero_meaned_data) * zero_meaned_data * weights = inv(zero_meaned_data) * all_desired_output
%%%% The inv(zero_meaned_data) and zero_meaned_data terms on the left
%%%% would cancel each other out, and we would be left with:
%%%% weights = inv(zero_meaned_data) * all_desired_output
%%%%
%%%% However, unfortunately there will in general not exist any
%%%% matrix-inverse of the data matrix zero_meaned_data.
%%%% Only square matrices have inverses, and not even all of them do.
%%%% Luckily, however, we can use something that plays a similar role,
%%%% called a pseudo-inverse. In Matlab, this is given by the command pinv.
%%%% The pseudo-inverse won't give us a perfect solution to the equation
%%%% zero_meaned_data * weights = all_desired_output
%%%% but it will give us the best approximate solution, which is what we want.
%%%%
%%%% So, instead of 
%%%% weights = inv(zero_meaned_data) * all_desired_output
%%%% we have this equation:
weights = pinv(zero_meaned_data) * all_desired_output;

%%%% Let's have a look at how these weights carve up the input space
%%%% A useful Matlab command for making grids of points
%%%% which span a particular 2D space is called "meshgrid"
[input_space_X, input_space_Y] = meshgrid([-3:0.3:3],[-3:0.3:3]);
weighted_output_Z = input_space_X*weights(1) + input_space_Y*weights(2);

%%%% The weighted output gets turned into the category-decision +1
%%%% if it is greater than 0, and -1 if it is less than zero.
%%%% The easiest way to map positive numbers to +1
%%%% and negative numbers to -1 
%%%% is by first mapping them to 1 and 0
%%%% by the inequality-test(weighted_output_Z>0)
%%%% and then turning 1 and 0 into +1 and -1
%%%% by multipling by 2 and subtracting 1.
decision_output_Z = 2*(weighted_output_Z>0) - 1;

figure(2);
clf;
hold on;
surf(input_space_X,input_space_Y,decision_output_Z);
%%% Let's show this decision surface in gray, from a good angle
colormap gray;
caxis([-3 3]);   %%% Sets white and black values to +/-3, so +/-1 are gray
shading interp;  %%% Makes the shading look prettier
grid on;
view(-10,60);
rotate3d on;   %%% Make it so we can use mouse to rotate the 3d figure
set(gca,'FontSize',14);
title('Click and drag to rotate view');

%%%% Let's plot the zero-meaned sumo and basketball data on top of this
%%%% Each class has 5 members, in this case, so we'll subtract
%%%% a mean-column-values matrix with 5 rows, to make the matrix sizes match.
one_class_matrix_of_mean_vals = ones(5,1) * mean_column_vals;
zero_meaned_sumo_wrestlers = sumo_wrestlers - one_class_matrix_of_mean_vals;
zero_meaned_basketball_players = basketball_players - one_class_matrix_of_mean_vals;

plot3(zero_meaned_sumo_wrestlers(:,1),zero_meaned_sumo_wrestlers(:,2), ...
      desired_output_sumo,'ro','LineWidth',5);
hold on;
plot3(zero_meaned_basketball_players(:,1),zero_meaned_basketball_players(:,2), ...
      desired_output_basketball,'bx','LineWidth',5,'MarkerSize',15);


xlabel('Height');
ylabel('Weight');
zlabel('Classifier output');