% function [P,p]= Cal2AFC(x,y)
% Computes the sensitivity of a system using a 2Alternative Forced
% Choice Test.  (2AFC)
% 
% x is a real valued vector
% y is a binary vector of labels. It should be of the same size as x
% f is an optional vector with the frequency of each (x,y) pair
% if f is ommitted then f = 1 for all pairs
%
% P is percent correct on the 2AFC task
% p is the 2AFC percent correct for each individual (x,y) pair
% It is useful to detect (x,y) pairs are more difficult 
%
% Example 
% x= [ 1 2 3 2 3]'
% y= [ 0 0 0 1 1]'
% f= [ 3 1 1 3 4]'
%
% [P, p] = Calc2AFC(x,y)
%
% In this case P= 0.8143
%
%  p = [ 1 0.7875 0.2857 0.7 0.9]
% 
% This tells us that the pair (x=3, y=0) was very difficult, when
% presented as a representative of category 0, the probability of
% getting it right was 0.2857, quite less than chance, which is 0.5
%
% On the other hand the pair (x=3, y=1) was easy. when presented as
% a representative of category zero, the probability of getting it
% right was 0.9.
%
% Javier R. Movellan
% April 2010, May 2011 
%
function [P,p] =Calc2AFC(x,y,f)
  if(size(x,2)> 1) x = x'; end
  if(size(y,2)> 1) y = y'; end
  if(size(x,2) ~=1) 
    error('x should have only one column');
  end
  if(size(y,2) ~=1) 
    error('y should have only one column');
  end
  if nargin ==2
    f= ones(size(x));
  end
  
  c = unique(y);   
  if ( length(c) ~=2)
    error('Cannot calculate 2AFC: There are no 2 categories');
  end
  x0i= find(y==c(1));
  x0 = x(x0i);
  f0= f(x0i);
  f0 = f0/sum(f0);
  
  x1i = find(y==c(2));
  x1 = x(x1i);
  f1 = f(x1i);
  f1=f1/sum(f1);
  
  n0 = length(x0);
  n1 = length(x1);
  for k=1: n0
    % if I get item x0(k) to represent category 0 in combination with
    % any randomly select x1 item representing category 1 what's the probability that I
    % classify x0(k) correctly
    c0(k,1) = (sum(f1(find(x1>x0(k))))  + 0.5*sum( f1(find(x1 == x0(k)))));
  end
  % if I get item x1(k) to represent category 1 in combination with
  % any randomly select x0 item representing category 0 what's the probability that I
  % classify x1(k) correctly
  
  for k=1: n1
    c1(k,1) = (sum(f0(find(x0<x1(k))))  + 0.5*sum( f0(find(x0 == x1(k)))));
  end

  p(x0i,1) = c0;
  p(x1i,1) = c1;
  
  P = sum(c0.*f0); % equivalently sum(c1.*f1)
  


     