clear

%% DATA
fileID = fopen('radial-profile.txt','r');
C = textscan(fileID,'%f','CommentStyle','frame');
B = reshape(C{1,1},3,[]);
[m,n] = size(B);
B = horzcat(B,[0;0;0]);
p_l = 0.030303030303030; % length(mm)/pixel
R = 3.5; % radius(mm)
Component = 4; %component number
pos_ini = 1;
frame_num = 1;
for pos = 1:n
	if B(2,pos) > B(2,pos+1)
		frame(frame_num,:) = B(3,pos_ini:pos);
		frame_num = frame_num+1;
		pos_ini = pos+1;
	end
end
frame_total = frame_num;

p = B(2,n);
zero = zeros(B(2,1)-1,frame_num-1);
frame2 = vertcat(zero,frame.');
attempt_limits = 1000;
random = 0.0100;
for frame_num = 1:frame_total-1
	x = zeros(p,1);
	for i=1:p
		x(i,1) = (i-1) * p_l;
	end
	Data1 = zeros(p,2);
	Data1(:,1) = x(:,1);
	Data1(:,2) = frame2(:,frame_num);
	Data1(1:B(2,1)-1,:) = [];
	count = 0;
	count_end = 10;
	y = Data1(:,2);
	t3 = Data1(:,1);

%%
for iter=1:attempt_limits
	fprintf('Frame%dL EM Iteration %d\n',frame_num,iter)
	lam0 = random*randsample(1:999,Component);
	options = optimset('LargeScale','off','Display','off','TolX',1e-16);
   	t2 = Data1(:,1);
   	y2 = Data1(:,2);
   	plothandle = plotdatapoints(t2,y2);
   	F = @(x) fitfunc1(x,Data1,R);
   	foutputfcn = @(x,optimvalues,state) fitfun2outputfcn0(x,optimvalues,state, Data1,plothandle);
	options = optimset(options,'OutputFcn',foutputfcn);
	options = optimset(options,'Algorithm','levenberg-marquardt');
	try
	[lam,resnorm,residual,exitflag,output] = lsqnonlin(F,lam0,[],[],options);
	axis([0 3.59 -0.0005 0.008]);
	%%
	A2 = zeros(length(t2),length(lam));
	for j = 1:length(lam)
		A2(:,j) = exp(-t2.*t2/lam(j))-exp(-((R+0)^2)/lam(j));
	end
	cc2 = A2\y2; % solve A*c = y for linear parameters c
    A3 = zeros(length(t3),1);
    for i=1:length(lam)
        coefficient = cc2(i)*(exp(-(t3+0).*(t3+0)/lam(i))-exp(-((R+0)^2)/lam(i)));
        A3=A3+coefficient;
    end
	C = (A3-y).^2;
	RM = mean(C);
	RMS = sqrt(RM);
	catch
	lam=NaN(Component,1);
	cc2=NaN(Component,1);
	end

	B1 = lam>0.15;
	B2 = all(B1);
	B3 = all(cc2);
	if (exitflag >0 && B2>0 && B3>0 && count<count_end)
		count=count+1;
	Multi_result.lam(count,:)=lam;
	Multi_result.cc2(count,:)=cc2;
	Multi_result.RM(count,1)=RM;
	elseif exitflag <= 0
	else
	end

	if count == count_end
        	break    
	end
end  %iter

n = length(y);
Maximum_log_lilehood = (-n/2)*(log(2*pi)+log(Multi_result.RM)+1);
AIC = (-2*Maximum_log_lilehood)+(2*(1+2*Component+1));

for zc = 1:count_end
	if AIC(zc,1) == -Inf
		AIC(zc,1) = 1;
	end
end

number = find(AIC==min(AIC), 1);
result.lam(frame_num,:) = Multi_result.lam(number,:);
result.cc2(frame_num,:) = Multi_result.cc2(number,:);

filename = ['Approx./',num2str(Component),'Term num./',num2str(frame_num),'frame','.bmp'];
saveas(gcf,filename,'bmp')

end
save(sprintf('result1.mat'))