diff --git a/fnn.m b/fnn.m new file mode 100644 index 0000000000000000000000000000000000000000..db26c20695861432a37b36dc4aa1e3587d3fe78c --- /dev/null +++ b/fnn.m @@ -0,0 +1,640 @@ +function out=fnn(varargin) +%FNN Dimension estimation by false nearest neighbours. +% Y=FNN(X) computes the vector Y of the amount of false +% nearest neighbours (FNN) as a function of the embedding dimension. +% +% Y=FNN(X,M), where M is a scalar, computes the FNN up to dimension M. +% The defeault is M=10. +% +% Y=FNN(X,M,R,S), where R and S are scalars, applies the neighbourhood +% criterion R and the size of the neighbourhood S. The defeault is R=2 +% and S=Inf. +% +% Y=FNN(X,M,R,S,N), where N is a scalar, uses N random samples for +% the determination of the FNNs. This speeds up the estimation, +% especially for long data series. The defeault is N=length(X) if +% the data length is smaller than 500, else N=200. +% +% FNN(...) without any output arguments opens a GUI for interactively +% changing the parameters. +% +% By using the GUI, the FNNs can be stored into the +% workspace. +% +% FNN(...,param) additional parameters according to the GUI are +% available: +% gui - Creates the GUI. +% nogui - Suppresses the GUI. +% silent - Suppresses all output. +% +% FNN without any arguments calls a demo (the same as the example +% below). +% +% Examples: x = sin(0:.2:8*pi)' + .1*randn(126,1); +% fnn(x,10,5) +% +% See also PHASESPACE, PSS, MI. +% +% References: +% Kennel, M. B., Brown, R., Abarbanel, H. D. I.: +% Determining embedding dimension for phase-space reconstruction +% using a geometrical construction, Phys. Rev. A, 45, 1992. + +% Copyright (c) 2006 +% Norbert Marwan, Potsdam University, Germany +% http://www.agnld.uni-potsdam.de +% +% $Date$ +% $Revision$ +% +% $Log$ +% +% +% This program is part of the new generation XXII series. +% +% This program 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 2 +% of the License, or any later version. + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% programme properties + +global props + +init_properties +action=''; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input + +error(nargchk(0,99,nargin)); +if nargout>2, error('Too many output arguments'), end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% splash the GPL + +splash_gpl('fnn'); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% error control +try + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% read the input + +set(0,'ShowHidden','on') + +delete(findobj('Tag','msgbox')) +nogui=0; +maxM_init = 10; +maxM=maxM_init; % maximal dimension +r_init = 2; +r=r_init; % r-criterion for neighbours distance +s_init = Inf; +s=s_init; % size of neighbourhood +maxN=[]; % number of samples +x={}; +method = []; +if nargin & isnumeric(varargin{1}) + + % check the numerical input + i_double=find(cellfun('isclass',varargin,'double')); + i_length = cellfun('prodofsize',varargin(i_double)); + i_x = i_double(i_length > 1); + x = varargin{i_x(1)}(:); + N = length(x); + if isempty(x), error('Not a valid input vector.'), end + i_x = i_double(i_length < 2); + switch length(i_x) + case 1 + maxM = varargin{i_x}; + case 2 + maxM = varargin{i_x(1)}; + r = varargin{i_x(2)}; + case 3 + maxM = varargin{i_x(1)}; + r = varargin{i_x(2)}; + s = varargin{i_x(3)}; + case 4 + maxM = varargin{i_x(1)}; + r = varargin{i_x(2)}; + s = varargin{i_x(3)}; + maxN = varargin{i_x(4)}; + end + + % check the char input + i_char=find(cellfun('isclass',varargin,'char')); + check_meth={'ma','eu','mi'}; % maxnorm, euclidean, nrmnorm, fan, distance + check_gui={'gui','nog','sil'}; % gui, nogui, silent + temp_meth=0; + temp_gui=0; + + if ~isempty(i_char) + for i=1:length(i_char), + varargin{i_char(i)}(4)='0'; + temp_meth = temp_meth+strcmpi(varargin{i_char(i)}(1:2),check_meth'); + temp_gui = temp_gui+strcmpi(varargin{i_char(i)}(1:3),check_gui'); + end + end + method = min(find(temp_meth)); + nogui = min(find(temp_gui))-1; + if isempty(nogui); nogui = 0; end + + + action = 'init'; + if nogui > 0, action = 'compute'; end + if nargout & nogui==0 + nogui=1; + action='compute'; + end + +elseif nargin & ischar(varargin{1}) & ~isempty(varargin{1}) + + if isempty(findobj('Tag','FNN_Fig')) & isempty(findobj('Tag','FNN_Fig')) + action='end'; + else + + action=varargin{1}; + h=findobj('Tag','FNN_Fig'); + x=get(h(1),'UserData'); + N=size(x,2); + h=findobj('Tag','maxM'); + maxM=str2num(get(h(1),'String')); + h=findobj('Tag','r'); + r=str2num(get(h(1),'String')); + h=findobj('Tag','s'); + s=str2num(get(h(1),'String')); + h=findobj('Tag','maxN'); + maxN=str2num(get(h(1),'String')); + h=findobj('Tag','method'); + method=get(h(1),'Value'); + end +else + action='init'; +end + +if ~nargin + x=sin(0:.2:8*pi)'+.1*randn(126,1); + maxM=10; + r=5; s = Inf; + nogui=0; +else + if ischar(varargin{1}) + action=varargin{1}; + end +end + + +% correct invalid inputs +N = length(x); +if isempty(maxM) | maxM > N/2, maxM = maxM_init; end +if maxM < 2, maxM = 2; end +if isempty(r), r = r_init; end +if isempty(s), s = s_init; end +if isempty(maxN) + maxN = N; + if N > 500, maxN = 200; end +end +if isempty(method), method = 1; end + +% main part begins +switch(action) + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% create gui + +case 'init' + + h8=figure(props.window,... % Plot Figure + 'Tag','FNN_Fig',... + 'MenuBar','Figure',... + 'Position',[69.5000 39.6429 120.0000 30.0714],... + 'PaperType','a4',... + 'PaperPosition',[0.25 0.25 7.7677 11.193],... + 'PaperOrientation','portrait',... + 'UserData',x,... + 'Name','False Nearest Neighbours'); + + set(0,'showhidden','on') + h=findobj('Label','&Help','Type','uimenu'); + if isempty(h) + h=uimenu('Label','&Help'); + h2=uimenu('Parent',h(1),'Label','&Help False Nearest Neighbours','Callback','helpwin fnn'); + else + h1=flipud(get(h(1),'Children')); + set(h1(1),'Separator','on') + h2=uimenu('Parent',h(1),'Label','&Help False Nearest Neighbours','Callback','helpwin fnn'); + copyobj(h1,h(1)) + delete(h1) + end + + h=axes(props.axes,... + 'Position',[89 24.8 6.8 3.5]); + logo=load('logo'); + h2=imagesc([logo.logo fliplr(logo.logo)]); + set(h2,'Tag','uniLogo') + set(h,props.logo,'Tag','axes_logo') + h=uicontrol(props.text,... + 'Tag','text_logo',... + 'String','Uni Potsdam',... + 'Position',[97 24.2143 22 3.5714]); + h2=textwrap(h,{[char(169),' AGNLD'],'University of Potsdam','2006'}); + set(h,'String',h2) + + h=axes(props.axes,... + 'Tag','fnn_axes',... + 'Box','On',... + 'Position',[9 3.5 72.8333 23.0714]); + + % frame around parameter settings + h=uicontrol(props.frame,... + 'Tag','frame',... + 'Position',[91.5000 1.3571 23.5000 21.5]); + + % dimension + h=uicontrol(props.text,... + 'Tag','text',... + 'String','Max. Dim.',... + 'Position',[93.8333 20.3 14 1.5000]); + + h=uicontrol(props.edit,... + 'Tag','maxM',... + 'String',num2str(maxM),... + 'ToolTip','False nearest neighbours will be computed up to this dimension.',... + 'Position',[105.1666 20.55 6.5 1.5000]); + + % falseness criterion + h=uicontrol(props.text,... + 'Tag','text',... + 'String','Falseness',... + 'Position',[93.8333 18.3 20.5 1.5000]); + + h=uicontrol(props.edit,... + 'Tag','r',... + 'String',num2str(r),... + 'ToolTip','Criterion value (suggested: >1).',... + 'Position',[105.1666 18.55 6.5 1.5000]); + + % neigbourhood range + h=uicontrol(props.text,... + 'Tag','text',... + 'String','Neigbourh.',... + 'Position',[93.8333 16.3 20.5 1.5000]); + + h=uicontrol(props.edit,... + 'Tag','s',... + 'String',num2str(s),... + 'ToolTip','Maximal distance between neighbours (in multiples of standard deviation).',... + 'Position',[105.1666 16.55 6.5 1.5000]); + + % norm + h=uicontrol(props.text,... + 'Tag','text',... + 'String','Norm',... + 'Position',[93.8333 14.1 15 1.5000]); + + h=uicontrol(props.popup,... + 'Tag','method',... + 'String','Max|Euc|Min',... + 'Value',method,... + 'ToolTip','Norm used for the computation of the distance.',... + 'Position',[101.1666 14.35 10.5 1.5000]); + + % number of samples + h=uicontrol(props.text,... + 'Tag','text',... + 'String','Number of samples',... + 'Position',[93.8333 11.75+.4 19 1.5000]); + + h=uicontrol(props.edit,... + 'Tag','maxN',... + 'String',num2str(maxN),... + 'ToolTip','Number of randomly chosen samples.',... + 'Position',[93.8333 10.5+.4 17.8333 1.5000]); + + + % buttons + offset = -.5; + h=uicontrol(props.button,... + 'String','Store',... + 'Tag','button_store',... + 'Enable','Off',... + 'ToolTip','Stores the number of false nearest neighbours into a variable in the workspace.',... + 'Callback','fnn store',... + 'Position',[103.8666 8.2143+offset 8.8 2.2143]); + + h=uicontrol(props.button,... + 'Tag','button_print',... + 'CallBack','fnn print',... + 'ToolTip','Prints the FNN window.',... + 'String','Print',... + 'Position',[93.8333 8.2143+offset 8.8 2.2143]); + + h=uicontrol(props.button,... + 'Tag','button_close',... + 'CallBack','fnn close',... + 'ToolTip','Closes the FNN window.',... + 'String','Close',... + 'Position',[93.8333 2.6429+offset 18.8333 2.2143]); + + h=uicontrol(props.button,... + 'String','Apply',... + 'Tag','button_apply',... + 'ToolTip','Starts the computation.',... + 'Callback','fnn compute',... + 'Position',[93.8333 5.4286+offset 18.8333 2.2143]); + + set(h8, 'HandleVis','CallBack') + tags={'FNN_Fig';'axes_logo';'text_logo';'frame';'text';'fnn_axes';'maxM';'r';'s';'method';'maxN';'button_store';'button_print';'button_close';'button_apply';}; + h=[]; + for i=1:length(tags); h=[h; findobj('Tag',tags{i})]; end + set(h,'Units','Norm') + fnn('compute') + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close windows + +case 'close' + set(0,props.root) + h=findobj('Tag','FNN_Fig'); + if ~isempty(h), close(h(1)), end + clear all + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% store + +case 'store' + if ~isempty(findobj('Tag','button_store')) + h=findobj('Tag','button_store'); + h1=findobj('Tag','button_close'); + if ~isempty(h1), vname_old=get(h1(1),'UserData'); else vname_old=''; end + if isempty(vname_old), vname_old=''; end + vname=char(inputdlg('Choose a variable name.','Store output',1,{vname_old})); + if isempty(vname) + return + else + FNN=get(h(1),'UserData'); + assignin('base',vname, [FNN]) + h=helpdlg(['False nearest neighbours have been assigned to the workspace variable ''',vname,'''.',... + ''],'Store output'); + set(h,'Tag','msgbox') + set(h1(1),'UserData',vname) + end + end + set(0,props.root) + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% print + case 'print' + + h=findobj('Tag','uniLogo'); + h_axes=findobj('Tag','fnn_axes','Parent',gcf); h_axes=h_axes(1); + h=[h; findobj('Tag','text','Parent',gcf)]; + h=[h; findobj('Tag','text_logo','Parent',gcf)]; + h=[h; findobj('Tag','frame','Parent',gcf)]; + h=[h; findobj('Tag','maxM','Parent',gcf)]; + h=[h; findobj('Tag','r','Parent',gcf)]; + h=[h; findobj('Tag','s','Parent',gcf)]; + h=[h; findobj('Tag','method','Parent',gcf)]; + h=[h; findobj('Tag','maxN','Parent',gcf)]; + h=[h; findobj('Tag','button_store','Parent',gcf)]; + h=[h; findobj('Tag','button_print','Parent',gcf)]; + h=[h; findobj('Tag','button_close','Parent',gcf)]; + h=[h; findobj('Tag','button_apply','Parent',gcf)]; + set(h,'Visible','Off') + + set(h_axes,'Units','Character') + old_pos=get(h_axes,'Position'); + set(h_axes,'Units','normalize','Position',[0.1300 0.1100 0.7750 0.8150]) + h_dlg=printdlg; + waitfor(h_dlg) + + set(h_axes,'Units','Character','Position',old_pos) + set(h_axes,'Units','normalize') + set(h,'Visible','On') + set(0,props.root) + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% computation + +case 'compute' + warning off + if ~nogui + h_fig=findobj('tag','FNN_Fig'); + setptr(gcf,'watch'), + obj=({'maxM','method','r','s','maxN','button_store','button_print','button_close','text'}); + for j=1:length(obj); + h=findobj('Tag',obj{j},'Parent',h_fig(1)); + if ~isempty(h), set(h,'Enable','Off'), end + end + h=findobj('tag','button_apply'); + set(h(1),'ToolTip','Stops the computation.','String','Stop','Callback','set(0,''ShowHidden'',''on'');h=findobj(''tag'',''button_apply'');set(h(1),''String'',''Stopped'');set(0,''ShowHidden'',''off'')') + end + if nogui~=2; hw=waitbar(0,'Estimation progress'); end + + FNN(1,1) = 1; + x_s = s * std(x); + for m = 2:maxM + % check if stop button was pressed + if ~nogui + set(0,'ShowHidden','on') + h=findobj('tag','button_apply','Parent',h_fig(1)); + if strcmpi(get(h(1),'string'),'stopped') + FNN((m:maxM)+1,1)=NaN; + break + end + end + + FNN(m,1) = 0; + t=1; + NX=N-t*(m-1); + + % create phase space vectors + clear jx,for mi=1:m; + jx(1+NX*(mi-1):NX+NX*(mi-1))=1+t*(mi-1):NX+t*(mi-1); + end + + x1=reshape(repmat(x,1,m),N*m,1); + x2=reshape(x1(jx),NX,m); % phasespace vector + jx = reshape(jx,NX,m); + if maxN >= NX; + idx = 1:NX; + else + idx = ceil(NX * rand(maxN,1)); %random samples + end + + % count the false neighbours + cnt = 0; + for i = 1:length(idx), if nogui~=2 & i*10/length(idx) == ceil(i*10/length(idx)); waitbar((length(idx)*(m-1)+i)/(length(idx)*maxM)), end + % distance in the m-dimensional space + switch method + case 1 + distance = max(abs(x2-repmat(x(jx(idx(i),:))',length(x2),1))'); + case 2 + distance = sqrt(sum((x2-repmat(x(jx(idx(i),:))',length(x2),1)).^2,2)); + case 3 + distance = sum(abs(x2-repmat(x(jx(idx(i),:))',length(x2),1)),2); + end + [distance is] = sort(distance); + + if is(2) < NX & idx(i) < NX & distance(2) < x_s + FNN(m,1) = FNN(m,1) + double((abs(x(idx(i)+m) - x(is(2)+m))/distance(2) > r)); + cnt = cnt + 1; + end + % check if stop button was pressed + if ~nogui + set(0,'ShowHidden','on') + h=findobj('tag','button_apply','Parent',h_fig(1)); + if strcmpi(get(h(1),'string'),'stopped') + FNN((m:maxM)+1,1)=NaN; + break + end + end + + end + FNN(m,1) = FNN(m,1) / cnt; + + if ~nogui, + if ~isempty(findobj('Tag','FNN_Fig')) + h=findobj('Tag','button_store'); + out_str=FNN; + set(h(1),... + 'UserData',out_str) + end + fnn('plot_fnn') + end + + end + + if nogui~=2 & ishandle(hw), delete(hw); end + if ~nogui, + set(0,'ShowHidden','on') + h=findobj('tag','button_apply'); + set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','fnn compute') + for j=1:length(obj); + h=findobj('Tag',obj{j},'Parent',h_fig(1)); + if ~isempty(h), set(h,'Enable','On'), end + end + setptr(h_fig(1),'arrow') + if ~isempty(findobj('Tag','FNN_Fig')) + h=findobj('Tag','button_store'); + set(h(1),'Enable','On') + end + + if ~nogui, setptr(h_fig(1),'arrow'), end + warning on + try, set(0,props.root), end + if nargout==0 + fnn('plot_fnn') + end + end + if nogui + out=FNN; + end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output + +case 'plot_fnn' + set(0,'ShowHidden','on') + h_fig=findobj('tag','FNN_Fig'); + if ~isempty(h_fig) + h=findobj('Tag','button_store','Parent',h_fig(1)); + if isempty(h), return, end + FNN=get(h(1),'UserData'); + + h_axes=findobj('Tag','fnn_axes','Parent',h_fig(1)); h_axes=h_axes(1); + h_res=findobj('Tag','fnn_result','Parent',h_fig(1)); delete(h_res); + set(0,'current',h_fig(1)) + + set(h_axes,'Visible','on'); + plot(1:length(FNN),FNN), grid on + xlabel('Dimension'), ylabel('False Nearest Neighbours') + set(gca,'Tag','fnn_axes','layer','top','xlim',[1 maxM],'ylim',[0 1],'xtick',[1:1:maxM]) + drawnow + end + try, set(0,props.root), end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the end + +case 'end' + +end +set(0,'ShowHidden','off') + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% error handling +%if 0 +catch + if ~isempty(findobj('Tag','TMWWaitbar')), delete(findobj('Tag','TMWWaitbar')), end + cmd={['mean'];['var'];['std'];['median'];['squmean'];['geomean'];['bias'];['skewness'];['kurtosis']}; + z_whos=whos;x_lasterr=lasterr;y_lastwarn=lastwarn;if nargin;in=varargin{1};else in.class='no input given';end + if ischar(in), in2=in; else, in2=[]; end + in=whos('in'); + if ~strcmpi(lasterr,'Interrupt') + fid=fopen('error.log','w'); + err=fprintf(fid,'%s\n','Please send us the following error report. Provide a brief'); + err=fprintf(fid,'%s\n','description of what you were doing when this problem occurred.'); + err=fprintf(fid,'%s\n','E-mail or FAX this information to us at:'); + err=fprintf(fid,'%s\n',' E-mail: marwan@agnld.uni-potsdam.de'); + err=fprintf(fid,'%s\n',' Fax: ++49 +331 977 1142'); + err=fprintf(fid,'%s\n\n\n','Thank you for your assistance.'); + err=fprintf(fid,'%s\n',repmat('-',50,1)); + err=fprintf(fid,'%s\n',datestr(now,0)); + err=fprintf(fid,'%s\n',['Matlab ',char(version),' on ',computer]); + err=fprintf(fid,'%s\n',repmat('-',50,1)); + err=fprintf(fid,'%s\n',x_lasterr); + err=fprintf(fid,'%s\n',y_lastwarn); + err=fprintf(fid,'%s\n',[' during ==> fnn:',action]); + err=fprintf(fid,'%s',[' input ==> ',in.class]); + if ~isempty(in2), err=fprintf(fid,'\t%s\n',[' (',in2,')']); end + err=fprintf(fid,'%s\n',[' errorcode ==> no errorcode available']); + err=fprintf(fid,'%s\n',' workspace dump ==>'); + if ~isempty(z_whos), + err=fprintf(fid,'%s\n',['Name',char(9),'Size',char(9),'Bytes',char(9),'Class']); + for j=1:length(z_whos); + err=fprintf(fid,'%s',[z_whos(j).name,char(9),num2str(z_whos(j).size),char(9),num2str(z_whos(j).bytes),char(9),z_whos(j).class]); + if ~strcmp(z_whos(j).class,'cell') & ~strcmp(z_whos(j).class,'struct') + content=eval(z_whos(j).name); + content=mat2str(content(1:min([size(content,1),500]),1:min([size(content,2),500]))); + err=fprintf(fid,'\t%s',content(1:min([length(content),500]))); + elseif strcmp(z_whos(j).class,'cell') + content=eval(z_whos(j).name); + err=fprintf(fid,'\t'); + for j2=1:min([length(content),500]) + content2=content{j2}; + err=fprintf(fid,'{%s} ',content{j2}); + end + elseif strcmp(z_whos(j).class,'struct') + content=fieldnames(eval(z_whos(j).name)); + content=char(content); content(:,end+1)=' '; content=content'; + err=fprintf(fid,'\t%s',content(:)'); + end + err=fprintf(fid,'\n'); + end, end + err=fclose(fid); + disp('----------------------------'); + disp(' ERROR OCCURED'); + disp(' during executing fnn'); + disp('----------------------------'); + disp(x_lasterr); + disp([' during ',action]); + disp('----------------------------'); + disp(' Please send us the error report. For your convenience, ') + disp(' this information has been recorded in: ') + disp([' ',fullfile(pwd,'error.log')]), disp(' ') + disp(' Provide a brief description of what you were doing when ') + disp(' this problem occurred.'), disp(' ') + disp(' E-mail or FAX this information to us at:') + disp(' E-mail: marwan@agnld.uni-potsdam.de') + disp(' Fax: ++49 +331 977 1142'), disp(' ') + disp(' Thank you for your assistance.') + warning('on') + end + try,if ~nogui + setptr(h_fig(1),'arrow') + h=findobj('tag','button_apply'); + set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','fnn compute') + for j=1:length(obj); + h=findobj('Tag',obj{j},'Parent',h_fig(1)); + if ~isempty(h) + set(h,'Enable','On') + end + end + end,end + if nargout, out=NaN; end + try, set(0,props.root), end + set(0,'ShowHidden','off') +end