% Copyright 2011-2012 Blerim Emruli, eblerim@gmail.com
% 
% RPM-NHRR v.0.1 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.
%
% MATLAB script that solves a Raven's Progressive Matrix similar to the one 
% considered in [1].
%
% Solving a Raven's Progressive Matrix of the form:	
%   a --> aa --> aaa
%	b --> bb --> bbb
%	c --> cc -->  ?
%
% This 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 code; if not, write to the Free Software Foundation,
% Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
% 
% References:
%  [1] Rasmussen, D. and Eliasmith, C.  (2011).
%      A neural model of rule generation in inductive reasoning
%      Topics in Cognitive Science 3 (1), 140-153.
%
%  [2] Plate, T. A. (1994).
%      Distributed Representations and Nested Compositional Structure
%      University of Toronto, Toronto, Canada.

clc;
clear;

% resultsTxt = fopen('results.txt','w');

% dimensionality of the Holographic Reduced Representations (HRRs)
dVectors = 1000;

nLoops = 100;

% save all the correlations in a matrix, in order to calculate the
% means and standard deviations for each alternative solution based on the
% predicted answer from the model
mCorr = zeros(nLoops,8);

% save the number of occurences for each alternative solution when it
% has the highest correlation with the predicted answer from the model
nOccurr(1,8) = 0;

for nLoopCount = 1:nLoops
    
    mCorr(nLoopCount,1);
    
    % define a vocabulary based on a normal distribution with mean 0 and
    % variance 1/d
    shape = random('norm',0,sqrt(1/dVectors),1,dVectors);
    square = random('norm',0,sqrt(1/dVectors),1,dVectors);
    triangle = random('norm',0,sqrt(1/dVectors),1,dVectors);
    circle = random('norm',0,sqrt(1/dVectors),1,dVectors);
    diamond = random('norm',0,sqrt(1/dVectors),1,dVectors);
    
    number = random('norm',0,sqrt(1/dVectors),1,dVectors);
    
    % Reference:
    %  [4] Rasmussen, D. (2011).
    %      A neural model of rule finding in Raven's Progressive Matrices.
    %      Technical Report: CTN-TR-20090601-004
    %      Permanent link: http://ctnsrv.uwaterloo.ca/cnrglab/node/87
    %      Page number: 13 (bottom)
    one = random('norm',0,sqrt(1/dVectors),1,dVectors);
    plusone = random('norm',0,sqrt(1/dVectors),1,dVectors);
    two = cconv_fft(one,plusone);
    three = cconv_fft(two,plusone);
    four = cconv_fft(three,plusone);
    five = cconv_fft(four,plusone);
    
    % Learning the first five mapping examples of the RPM
    % OS: One circle
    cell11 = cconv_fft(shape,circle) + cconv_fft(number,one);
    invcell11 = [cell11(1) fliplr(cell11(2:end))];
    
    % OS: Two circles
    cell12 = cconv_fft(shape,circle) + cconv_fft(number,two);
    invcell12 = [cell12(1) fliplr(cell12(2:end))];
    
    % OS: Three circles
    cell13 = cconv_fft(shape,circle) + cconv_fft(number,three);
    invcell13 = [cell13(1) fliplr(cell13(2:end))];
    
    % Rule from one circle to two circles
    Tcell1112 = cconv_fft(invcell11,cell12);
    % Rule from two circles to three circles
    Tcell1213 = cconv_fft(invcell12,cell13);
    
    % OS: One square
    cell21 = cconv_fft(shape,square) + cconv_fft(number,one);
    invcell21 = [cell21(1) fliplr(cell21(2:end))];
    
    % OS: Two squares
    cell22 = cconv_fft(shape,square) + cconv_fft(number,two);
    invcell22 = [cell22(1) fliplr(cell22(2:end))];
    
    % OS: Three squares
    cell23 = cconv_fft(shape,square) + cconv_fft(number,three);
    invcell23 = [cell23(1) fliplr(cell23(2:end))];
    
    % Rule from one square to two squares
    Tcell2122 = cconv_fft(invcell21,cell22);
    % Rule from two squares to three squares
    Tcell2223 = cconv_fft(invcell22,cell23);
    
    % OS: One triangle
    cell31 = cconv_fft(shape,triangle) + cconv_fft(number,one);
    invcell31 = [cell31(1) fliplr(cell31(2:end))];
    
    % OS: Two triangles
    cell32 = cconv_fft(shape,triangle) + cconv_fft(number,two);
    invcell32 = [cell32(1) fliplr(cell32(2:end))];
    
    % Rule from one triangle to two triangles
    Tcell3132 = cconv_fft(invcell31,cell32);
    
    % calculate the final Transformation vector
    % OS: Average of all shape rules above
    T = 1/5*(Tcell1112 + Tcell1213 + Tcell2122 + Tcell2223 + Tcell3132);
    
    % create eight different alternative solutions
    % OS: These are solutions we assume exist and can compare to
    firstAlternative = cconv_fft(shape,triangle) + cconv_fft(number,three);
    secondAlternative = cconv_fft(shape,triangle) + cconv_fft(number,two);
    thirdAlternative = cconv_fft(shape,square) + cconv_fft(number,three);
    fourthAlternative = cconv_fft(shape,triangle) + cconv_fft(number,five);
    fifthAlternative = cconv_fft(shape,triangle) + cconv_fft(number,one);
    sixthAlternative = cconv_fft(shape,triangle) + cconv_fft(number,four);
    seventhAlternative = cconv_fft(shape,square) + cconv_fft(number,four);
    eighthAlternative = cconv_fft(shape,square) + cconv_fft(number,five);
    
    % OS: This is a new input in the Raven's task and we need to specify
    % the rule/answer based on what we have learned. What is the rule for four squares?
    fourSquares = cconv_fft(shape,square) + cconv_fft(number,four);
    
    % based on the above calculated Transformation vector predict the correct answer
    predictedAnswer = cconv_fft(fourSquares,T);
    
    % compare the predicted answer with different alternative solutions
    % based on correlation of vector elements
    temp = corrcoef(predictedAnswer,firstAlternative);
    mCorr(nLoopCount,1) = temp(1,2);
    
    temp = corrcoef(predictedAnswer,secondAlternative);
    mCorr(nLoopCount,2) = temp(1,2);
    
    temp = corrcoef(predictedAnswer,thirdAlternative);
    mCorr(nLoopCount,3) = temp(1,2);
    
    temp = corrcoef(predictedAnswer,fourthAlternative);
    mCorr(nLoopCount,4) = temp(1,2);
    
    temp = corrcoef(predictedAnswer,fifthAlternative);
    mCorr(nLoopCount,5) = temp(1,2);
    
    temp = corrcoef(predictedAnswer,sixthAlternative);
    mCorr(nLoopCount,6) = temp(1,2);
    
    temp = corrcoef(predictedAnswer,seventhAlternative);
    mCorr(nLoopCount,7) = temp(1,2);
    
    temp = corrcoef(predictedAnswer,eighthAlternative);
    mCorr(nLoopCount,8) = temp(1,2);
    
    % select and save the alternative solution with the highest correlation;
    % in order to (naively) calculate the propability that the particular alternative
    % solution will have the highest correlation and be selected as the correct answer
    [C,I] = max(mCorr(nLoopCount,:));
    nOccurr(1,I) = nOccurr(1,I)+1;

    % OS: Chosen answer is I. Here it is 8 corresponding to 5 squares (which follows
    % the rule 4 squares plus 1. C is the correlation.
    
end

% OS Look at the mean correlations for each alternative
figure;
meanCorrs = [];
for i=1:8
  meanCorrs = [meanCorrs mean(mCorr(:,i))];
end
bar(meanCorrs) 

if 0
 fprintf(resultsTxt, ' The mean and the standard deviation for each alternative:');
 fprintf(resultsTxt, '\n 1 %1.4f %1.4f \n 2 %1.4f %1.4f \n 3 %1.4f %1.4f \n 4 %1.4f %1.4f \n 5 %1.4f %1.4f \n 6 %1.4f %1.4f \n 7 %1.4f %1.4f \n 8 %1.4f %1.4f', ...
    mean(mCorr(:,1)), std(mCorr(:,1)), mean(mCorr(:,2)), std(mCorr(:,2)), ...
    mean(mCorr(:,3)), std(mCorr(:,3)), mean(mCorr(:,4)), std(mCorr(:,4)), ...
    mean(mCorr(:,5)), std(mCorr(:,5)), mean(mCorr(:,6)), std(mCorr(:,6)), ...
    mean(mCorr(:,7)), std(mCorr(:,7)), mean(mCorr(:,8)), std(mCorr(:,8)));

 fprintf(resultsTxt, '\n\n The probability one of the alternative solutions to have the highest correlation and be selected as the correct answer:');
 fprintf(resultsTxt, '\n 1 %1.6f \n 2 %1.6f \n 3 %1.6f \n 4 %1.6f \n 5 %1.6f \n 6 %1.6f \n 7 %1.6f \n 8 %1.6f \n', ...
    nOccurr(1,1)/nLoops, nOccurr(1,2)/nLoops, nOccurr(1,3)/nLoops, ...
    nOccurr(1,4)/nLoops, nOccurr(1,5)/nLoops, nOccurr(1,6)/nLoops, ...
    nOccurr(1,7)/nLoops, nOccurr(1,8)/nLoops);

save workspace.mat

end
