function model = bayesian_fit_vm(angles, counts, prm)
%
% Input:   
% -- counts            N X 1 vector, a spike count per trial (N = number of trials OR number of directions)
% -- angles            N X 1 vector, direction per trial
%
% Parameters:
% -- prm.lambda        prior coefficient (default 0 for ML)
% -- prm.nRepeats      number of repeats of minimization process (default 20)
% -- prm.PriorCoef     coefficient of the exponenetial distribution (default 1)
%
% Output:
% -- model             4 X 1 vector, VM estimated parameters
  
% Hadas taubman, Gal Chechik 2013
% Based on the paper: 
%       A Bayesian approach for characterizing direction tuning curves in ...
%	the supplementary motor area of behaving monkeys
%       H. Taubman, E. Vaadia, R.Paz, G. Chechik
%       jn.physiology.org/content/early/2013/03/01/jn.00449.2012.abstract

%init
  nFails = 0; nOut = 0; goodness = inf;
  A=[]; bm=[]; Aeq=[]; beq=[]; lb=[0;eps;eps;-inf]; ub=[]; nonlcon=[];
  Options = optimset('GradObj','on','Hessian','off','Display','off',...
		     'MaxFunEvals',400,'MaxIter',400);
  
  % Repeat optimazation with different seds to avoid local minima 
  for iter = 1:prm.nRepeats
    model0 = initModel(angles, counts, prm.kpInit);
    % vonMisesError is the function to be minimized. It recieves vector w 
    % and returns scalar f, which is the objective function evaluated at w. 
    try
      [w, fval, exitflag] = fmincon(@(w) vonMisesError(angles, counts, ...
						       w, prm.lambda, ...
						       prm.PriorCoef), ...
				    model0,A,bm,Aeq,beq,lb,ub,nonlcon,Options);
    catch ME
      fprintf('fmincon failed: %s\n',ME.message);
      nFails = nFails+1; continue;
    end
    if(exitflag < 1 || ~isreal(w))
      nFails = nFails+1; continue
    end
    
    if fval < goodness - 0.001
      goodness = fval;
      model = w;
    end
  end
  if nFails > 2
    fprintf('fmincon failed %d times out of %d\n', nFails, prm.nRepeats);
  end
  if nFails + nOut == prm.nRepeats
    model=[];
  end
end

% =======================================================================
function model0 = initModel(angles, counts, kpInitMethod)
  if ~exist('kpInitMethod','var'), kpInitMethod=2; end
  model0 = zeros(4,1);
  
  model0(1) = max(0,min(counts)+abs(randn)*0.05); % Baseline 
  switch kpInitMethod
    case 1 
      model0(3) = abs(randn*4);     % kappa
    case 2 
      model0(3) = rand*20;
    case 3 
      model0(3) = rand*50;
  end
  [mx, mi] = max(counts);
  model0(2) = mx / exp(model0(3)) + randn*0.05;   % m
  %model0(4) = rand*pi;             % theta (range:[-pi,pi])
  model0(4) = angles(mi)+randn*0.05;
end

