tpeakvonmises

%% 2-peak von mises
clear all;
close all;

load('/Users/rysk/Dropbox/[TEMP]/matlab.mat')
mkdir('res')
ORvar = Stim.params.OR(1:end-1)';
x     = 0:0.5:360;%ORvar(end);

twopeakVonmises = @(v, x) ...
	v(1) + v(2)*exp(-1* ((x*pi/180) - v(3)).^2 * (1/(2*v(4)^2))) ...
	     + v(5)*exp(-1* ((x*pi/180) - v(3)+v(6)).^2 * (1/(2*v(4)^2)));		

options = optimoptions('lsqnonlin', 'MaxIter', 10000,...
	'Algorithm', 'trust-region-reflective', 'Tolx', 10^-15, 'TolFun', 10^-15,...
	'ScaleProblem', 'none', 'MaxFunEvals', length(ORvar)*200);
 clear s
 clear c
 clear ey
 
for c = 2:size(stimResp.meanVal, 1)
	figure;
	%set(gcf, 'Visible', 'off')
	tmp = ((stimResp.meanVal(c,1:end-1)-1));
	[tmpest, f{c}, flag{c}, lamb(c), coeff{c}, ey(:,c)] = ...
		TPgaussianFit(ORvar, tmp);
	r = coeff{c}(2);
	s(c,:) = tmpest;
		
	
	subplot(121); hold on;	
	%rectangle('Position', [0 v(2)-v(1) 360 2*std(tmp)]); 
	peak1 = 180*mod(s(c,3), 2*pi)/pi;
	peak2 = 180*mod(s(c,6)+s(c,3),2*pi)/pi;
	
	s(c,3) = mod(s(c,3), 2*pi);
	%s(c,6) = s(c,3)+mod(s(c,6)-sign(v(6))*s(c,3),2*pi);
	%s(c,6) = mod(s(c,6), 2*pi);
	
	line([0 360], s(c,1) +[s(c,2) s(c,2)])
	line(peak1*[1 1], [s(c,1) s(c,1)+s(c,2)], 'Color', [1 0 1])
	line(peak2*[1 1], [s(c,1) s(c,1)+s(c,5)], 'Color', [0 1 1])

	plot(ORvar, tmp, '--o', 'lineWidth', 0.5);
	plot(x, twopeakVonmises(s(c,:),x), '-r')
	plot(ORvar, twopeakVonmises(s(c,:),ORvar), 'dr')
	title(num2str(r))
	
	subplot(1,2,2)
	polar([ORvar 360]*pi/180, [tmp tmp(1)]); hold on;
	polar((0:360)*pi/180, twopeakVonmises(s(c,:),0:360), 'r--');
	%saveas(gcf, ['res/' num2str(c, '%03d'), '.png'])
end

% lsqnonlin

TPgaussianFit.m

function [params, f, flag, lamb, coeff, ey] =  TPgaussianFit(x, y, x0, options)
% function [params f, flag, lamb] =  vonMises(x, y, x0, options)
% func   = @(v, x) ... % Von mises function
%	v(1) + v(2)*exp(-1* ((x*pi/180) - v(3)).^2 * (1/(2*v(4)^2)))
% 
% x-range should be 0~360 

func = @(v, x) ...
	v(1) + v(2)*exp(-1* ((x*pi/180) - v(3)).^2 * (1/(2*v(4)^2))) ...
	     + v(5)*exp(-1* ((x*pi/180) - v(3)+v(6)).^2 * (1/(2*v(4)^2)));		


if nargin<3
  v(1) = min(y);	                 % baseline
  [v(2), maxidx] = max(y);	 % amplitude 1 ()
  v(3) = x(maxidx)*pi/180;   % peak deg1
  v(4) = pi/4;                              % width
  v(5) = v(2);                              % amplitude 2;
  v(6) = pi;                                  % opposite peak deg
  options = optimoptions('lsqnonlin', 'MaxIter', 10000,...
	'Algorithm', 'trust-region-reflective', 'Tolx', 10^-10, 'TolFun', 10^-10,...
	'ScaleProblem', 'none', 'MaxFunEvals', length(x)*500);  
else
  v = x0;
end

ub = [v(1)+std(y) -v(1)+v(2)+2*std(y)  v(3)+pi/4  pi/2   -v(1)+v(2)+2*std(y) pi];
lb = [v(1)-std(y) -v(1)+v(2)-2*std(y)  v(3)-pi/4  0  -v(1)+v(2)     pi];

efunc = @(v) abs((func(v,x) - y));

for i = 1:2
  [params_tmp{i}, f_tmp(i), flag_tmp{i}, lamb_tmp{i}] = ...
      lsqnonlin(efunc, [v(1); v(2); v(3); v(4); v(5); v(6)], lb,ub,options);	
  
  % opposite direction
  v(6)  = -v(6);
  ub(6) = -ub(6);
  lb(6) = -lb(6);
end

[~, o] = min(f_tmp);    % optimal pattern
params = params_tmp{o};
f = f_tmp(o);
flag = flag_tmp{o};
lamb = lamb_tmp{o};

ey = func(params, 1:360);
coeff = corrcoef(y, func(params,x));