Skip to content
Snippets Groups Projects
Commit df3d0923 authored by marwan's avatar marwan
Browse files

speed enhanced

parent 6f033ff6
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 168 additions and 157 deletions
mi.m
+ 168
− 157
View file @ df3d0923
......@@ -36,7 +36,7 @@ function varargout=mi(varargin)
% Examples: x = sin(0:.2:8*pi)' + .1*randn(126,1);
% mi(x,40,100)
%
% See also HIST2, HISTN, ENTROPY.
% See also HIST2, HISTN, ENTROPY, MIGRAM.
%
% References:
% Roulston, M. S.:
......@@ -51,6 +51,9 @@ function varargout=mi(varargin)
% $Revision$
%
% $Log$
% Revision 3.10 2006/03/29 13:07:55 marwan
% problems regarding OPRPs and embedding resolved
%
% Revision 3.9 2006/02/14 11:46:15 marwan
% *** empty log message ***
%
......@@ -73,7 +76,6 @@ function varargout=mi(varargin)
global props
init_properties
action='';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% check the input
......@@ -85,7 +87,7 @@ if nargout>2, error('Too many output arguments'), end
splash_gpl('mi');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% error control
try
%try
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% read the input
......@@ -95,11 +97,13 @@ delete(findobj('Tag','msgbox'))
nogui=0;
lag=[];
nbin=10;
x={}; sum_m=0; sm_length=[];
if nargin & isnumeric(varargin{1})
x={}; sm_length=[];
if nargin && isnumeric(varargin{1})
sum_m=0;
for i=1:nargin
if isnumeric(varargin{i}) & max(size(varargin{i}))>1
if isnumeric(varargin{i}) && max(size(varargin{i}))>1
y=varargin{i};
if size(y,1)==1, y=y'; end
my=size(y,2); sum_m=sum_m+my;
......@@ -109,7 +113,7 @@ if nargin & isnumeric(varargin{1})
x(i)={y};
end
end
m=size(x,2); if sum_m==1; lag=0; x(2)=x(1); sum_m=2; m=2; end
m=size(x,2); if sum_m==1; lag=0; x(2)=x(1); m=2; end
if isempty(x), error('Not a valid input vector.'), end
i_double=find(cellfun('isclass',varargin,'double'));
......@@ -139,15 +143,15 @@ if nargin & isnumeric(varargin{1})
nogui=-1;
end
end
if nargout & (nogui<1)
if nargout && (nogui<1)
nogui=1;
action='compute';
end
if nogui==-1; nogui=0; end
elseif nargin & ischar(varargin{1}) & ~isempty(varargin{1})
elseif nargin && ischar(varargin{1}) && ~isempty(varargin{1})
if isempty(findobj('Tag','MI_Fig')) & isempty(findobj('Tag','MI_Fig'))
if isempty(findobj('Tag','MI_Fig')) && isempty(findobj('Tag','MI_Fig'))
% h=errordlg('Start the programme with mi(vector,vector).');
% set(h,'Tag','msgbox')
action='end';
......@@ -156,17 +160,17 @@ elseif nargin & ischar(varargin{1}) & ~isempty(varargin{1})
h=findobj('Tag','MI_Fig');
x=get(h(1),'UserData');
h=findobj('Tag','nBin');
nbin=str2num(get(h(1),'String'));
nbin=str2double(get(h(1),'String'));
m=size(x,2);
h=findobj('Tag','maxLag');
lag=str2num(get(h(1),'String'));
lag=str2double(get(h(1),'String'));
end
else
action='init';
end
if ~nargin
x={sin(0:.2:8*pi)'+.1*randn(126,1)}; x(2)=x; m=2; sum_m=2;
x={sin(0:.2:8*pi)'+.1*randn(126,1)}; x(2)=x; m=2;
nbin=10;
lag=40;
nogui=0;
......@@ -197,11 +201,11 @@ case 'init'
h=findobj('Label','&Help','Type','uimenu');
if isempty(h)
h=uimenu('Label','&Help');
h2=uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi');
uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi');
else
h1=flipud(get(h(1),'Children'));
set(h1(1),'Separator','on')
h2=uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi');
uimenu('Parent',h(1),'Label','&Help Mutual Information','Callback','helpwin mi');
copyobj(h1,h(1))
delete(h1)
end
......@@ -219,59 +223,59 @@ case 'init'
h2=textwrap(h,{[char(169),' AGNLD'],'University of Potsdam','2002-2006'});
set(h,'String',h2)
h=axes(props.axes,...
axes(props.axes,...
'Tag','mi_axes',...
'Box','On',...
'Position',[9 3.5 72.8333 23.0714]);
h=uicontrol(props.frame,...
uicontrol(props.frame,...
'Tag','frame',...
'Position',[91.5000 1.3571 23.5000 20.7857]);
h=uicontrol(props.text,...
uicontrol(props.text,...
'Tag','text',...
'String','Number of bins',...
'Position',[93.8333 20 17.8333 1.5000]);
h=uicontrol(props.edit,...
uicontrol(props.edit,...
'Tag','nBin',...
'String',num2str(nbin),...
'ToolTip','Number of bins for estimation of the probability function.',...
'Position',[93.8333 18.5714 17.8333 1.5000]);
h=uicontrol(props.text,...
uicontrol(props.text,...
'Tag','text',...
'String','max. Lag',...
'Position',[93.8333 16.5 17.8333 1.5000]);
h=uicontrol(props.edit,...
uicontrol(props.edit,...
'Tag','maxLag',...
'String',num2str(lag),...
'ToolTip','Maximal lag.',...
'Position',[93.8333 15.0714 17.8333 1.5000]);
h=uicontrol(props.text,...
uicontrol(props.text,...
'Tag','text_show',...
'String','Show components',...
'Enable','off',...
'Position',[93.8333 12.5 21.8333 1.5000]);
h=uicontrol(props.popup,...
uicontrol(props.popup,...
'Tag','edit_show1',...
'String','1|2',...
'Enable','off',...
'Callback','mi plot_mi',...
'Position',[93.8333 11.0714 7 1.5000]);
h=uicontrol(props.text,...
uicontrol(props.text,...
'Tag','text_show',...
'String','vs.',...
'Enable','off',...
'HorizontalAlign','center',...
'Position',[101.24995 11.0714-.2 4 1.5000]);
h=uicontrol(props.popup,...
uicontrol(props.popup,...
'Tag','edit_show2',...
'String','1|2',...
'Callback','mi plot_mi',...
......@@ -279,7 +283,7 @@ case 'init'
'Position',[105.6666 11.0714 7 1.5000]);
h=uicontrol(props.button,...
uicontrol(props.button,...
'String','Store',...
'Tag','button_store',...
'Enable','Off',...
......@@ -287,21 +291,21 @@ case 'init'
'Callback','mi store',...
'Position',[103.8666 8.2143 8.8 2.2143]);
h=uicontrol(props.button,...
uicontrol(props.button,...
'Tag','button_print',...
'CallBack','mi print',...
'ToolTip','Prints the MI window.',...
'String','Print',...
'Position',[93.8333 8.2143 8.8 2.2143]);
h=uicontrol(props.button,...
uicontrol(props.button,...
'Tag','button_close',...
'CallBack','mi close',...
'ToolTip','Closes the MI window.',...
'String','Close',...
'Position',[93.8333 2.6429 18.8333 2.2143]);
h=uicontrol(props.button,...
uicontrol(props.button,...
'String','Apply',...
'Tag','button_apply',...
'ToolTip','Starts the computation.',...
......@@ -380,6 +384,9 @@ case 'store'
case 'compute'
warning off
MI = zeros(m,m,lag);
MI_sigma = zeros(m,m,lag);
if ~nogui
h_fig=findobj('tag','MI_Fig');
setptr(gcf,'watch'),
......@@ -396,82 +403,83 @@ case 'compute'
if nogui~=2; hw=waitbar(0,'Estimation progress'); end
for t=0:lag,
if ~nogui
set(0,'ShowHidden','on')
h=findobj('tag','button_apply','Parent',h_fig(1));
if strcmpi(get(h(1),'string'),'stopped')
MI((t:lag)+1,1)=NaN;
MI_sigma((t:lag)+1,1)=NaN;
break
end
set(0,'ShowHidden','on')
h=findobj('tag','button_apply','Parent',h_fig(1));
if strcmpi(get(h(1),'string'),'stopped')
MI((t:lag)+1,1)=NaN;
MI_sigma((t:lag)+1,1)=NaN;
break
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the distributions
for k=1:m; for l=1:m;
P2=histn(x{k},x{l},nbin,t,'silent'); % 2D distribution
mP=length(size(P2));
if sum(~P2(:))/prod(size(P2)) > 0.99
warning on
warning('Too less data points for the estimation of the joint distribution.');
MI(k,l,(t:lag)+1)=NaN;
MI_sigma(k,l,(t:lag)+1)=NaN;
if ~nogui,
h=findobj('tag','button_apply');
set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','mi compute')
for j=1:length(obj);
h=findobj('Tag',obj{j},'Parent',h_fig(1));
if ~isempty(h)
set(h,'Enable','On')
for k=1:m;
for l=1:m;
P2=hist2(x{k},x{l},nbin,t,'silent'); % 2D distribution
mP=length(size(P2));
if sum(~P2(:))/numel(P2) > 0.99
warning on
warning('Too less data points for the estimation of the joint distribution.');
MI(k,l,(t:lag)+1)=NaN;
MI_sigma(k,l,(t:lag)+1)=NaN;
if ~nogui,
h=findobj('tag','button_apply');
set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','mi 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')
end
warning on
if nogui~=2, delete(hw); end
return
end
P2=P2/sum(P2(:)); % normalization
clear P1x
if mP==2
P1x=sum(P2,2); % distribution for x
P1y=sum(P2); % distribution for y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the mutual information
I1=[-sum((P1x(P1x~=0)).*log(P1x(P1x~=0))), -sum((P1y(P1y~=0)).*log(P1y(P1y~=0)))]; % entropies of Px and Py
I2=-sum(P2(P2~=0).*log(P2(P2~=0))); % entropy of joint distribution Pxy
I2_syserr=(length(P1x(P1x~=0))+length(P1y(P1y~=0))-length(P2(P2~=0).*log(P2(P2~=0)))-1)/(2*length(x{k})); % standard error for estimation of entropy
MI(k,l,t+1)=I1(1)+I1(2)-I2+I2_syserr;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the standard errors
P2xy=P1x*P1y;
i=(P2xy~=0 & P2~=0);
MI_sigma(k,l,t+1)=sqrt( sum( (log(P2xy(i)./P2(i))+MI(t+1)).^2 .*(P2(i).*(1-P2(i)))) /length(x{k}));
else
I1=[];
for i=1:mP;
P21=permute(P2,[i,i+1:mP, 1:(i-1)]);
P1x(:,i)=sum(reshape(P21,size(P21,2),size(P21,2)^(mP-1)),2); % distribution for x_i
I1 = -sum((P1x(P1x~=0)).*log(P1x(P1x~=0))); % sum of entropies of Px_i
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the mutual information
I2=-sum(P2(P2~=0).*log(P2(P2~=0))); % entropy of joint distribution Pxy
I2_syserr=(length(P1x(P1x~=0))-length(P2(P2~=0).*log(P2(P2~=0)))-1)/(2*length(x{k})); % standard error for estimation of entropy
MI(k,l,t+1)=I1-I2+I2_syserr;
end
end
setptr(h_fig(1),'arrow')
end
warning on
if nogui~=2, delete(hw); end
return
end
P2=P2/sum(P2(:)); % normalization
clear P1x
if mP==2
P1x=sum(P2')'; % distribution for x
P1y=sum(P2); % distribution for y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the mutual information
I1=[-sum((P1x(P1x~=0)).*log(P1x(P1x~=0))), -sum((P1y(P1y~=0)).*log(P1y(P1y~=0)))]; % entropies of Px and Py
I2=-sum(P2(P2~=0).*log(P2(P2~=0))); % entropy of joint distribution Pxy
I2_syserr=(length(P1x(P1x~=0))+length(P1y(P1y~=0))-length(P2(P2~=0).*log(P2(P2~=0)))-1)/(2*length(x{k})); % standard error for estimation of entropy
MI(k,l,t+1)=I1(1)+I1(2)-I2+I2_syserr;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the standard errors
P2xy=P1x*P1y;
i=(P2xy~=0 & P2~=0);
MI_sigma(k,l,t+1)=sqrt( sum( (log(P2xy(i)./P2(i))+MI(t+1)).^2 .*(P2(i).*(1-P2(i)))) /length(x{k}));
else
I1=[];
for i=1:mP;
P21=permute(P2,[i,i+1:mP, 1:(i-1)]);
P1x(:,i)=sum(reshape(P21,size(P21,2),size(P21,2)^(mP-1))')'; % distribution for x_i
I1 = -sum((P1x(P1x~=0)).*log(P1x(P1x~=0))); % sum of entropies of Px_i
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% compute the mutual information
I2=-sum(P2(P2~=0).*log(P2(P2~=0))); % entropy of joint distribution Pxy
I2_syserr=(length(P1x(P1x~=0))-length(P2(P2~=0).*log(P2(P2~=0)))-1)/(2*length(x{k})); % standard error for estimation of entropy
MI(k,l,t+1)=I1-I2+I2_syserr;
end
end, end
end
if nogui~=2; waitbar(t/lag); end
if ~nogui & t*10/lag == ceil(t*10/lag)
h=findobj('Tag','button_store');
if m==2
out_str{1}=MI; out_str{2}=MI_sigma;
else
out_str{1}=MI;
end
set(h(1),'UserData',out_str)
mi('plot_mi')
if ~nogui && t*10/lag == ceil(t*10/lag)
h=findobj('Tag','button_store');
if m==2
out_str{1}=MI; out_str{2}=MI_sigma;
else
out_str{1}=MI;
end
set(h(1),'UserData',out_str)
mi('plot_mi')
end
end
if nogui~=2; delete(hw); end
if ~nogui
......@@ -506,13 +514,12 @@ case 'compute'
if ~nogui, setptr(h_fig(1),'arrow'), end
warning on
if isnan(MI(1,1)), warning('Mutual information is NaN. Use smaller bin size.'), end
try, set(0,props.root), end
try set(0,props.root), catch end
if nargout==0
mi('plot_mi')
end
if nargout==0 & nogui
assignin('base','ans', MI)
ans=MI
if nargout==0 && nogui
varargout{1} = MI;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% output
......@@ -536,6 +543,7 @@ case 'plot_mi'
if lag==0
set(h_axes,'Visible','off'),cla
set(h_show,'enable','off')
out_str{MI} = [];
if m==2;
for k=1:length(MI)
out_str(k)={sprintf(repmat('%8.3f %6.3f ',1,2*length(MI)),reshape([MI(k,:);MI_sigma(k,:)],1,2*length(MI)))};
......@@ -549,7 +557,7 @@ case 'plot_mi'
'Tag','mi_result',...
'HorizontalAlignment', 'left',...
'FontWeight', 'bold',...
'String',['Mutual Information: ']);
'String','Mutual Information: ');
ex=get(h,'Extent'); set(h,'Position',[11 21 ex(3:4)]);
h=uicontrol(props.listbox,...
'Tag','mi_result',...
......@@ -576,14 +584,14 @@ case 'plot_mi'
patch(s3,s1,-10,'FaceColor',[.9 .9 1],'EdgeColor',[.85 .85 1])
end
hold on
end
end
plot(0:size(MI,3)-1,permute(MI(dimx,dimy,:),[3,1,2])), grid on
xlabel('Lag'), ylabel('Mutual Information')
hold off
set(gca,'Tag','mi_axes','layer','top','xlim',[0 lag])
end
drawnow
try, set(0,props.root), end
try set(0,props.root), catch end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the end
......@@ -596,55 +604,55 @@ set(0,'ShowHidden','off')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% error handling
%if 0
catch
if 0
%catch
if ~isempty(findobj('Tag','TMWWaitbar')), delete(findobj('Tag','TMWWaitbar')), end
cmd={['mean'];['var'];['std'];['median'];['squmean'];['geomean'];['bias'];['skewness'];['kurtosis']};
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
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 ==> mi:',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 ==>');
fprintf(fid,'%s\n','Please send us the following error report. Provide a brief');
fprintf(fid,'%s\n','description of what you were doing when this problem occurred.');
fprintf(fid,'%s\n','E-mail or FAX this information to us at:');
fprintf(fid,'%s\n',' E-mail: marwan@agnld.uni-potsdam.de');
fprintf(fid,'%s\n',' Fax: ++49 +331 977 1142');
fprintf(fid,'%s\n\n\n','Thank you for your assistance.');
fprintf(fid,'%s\n',repmat('-',50,1));
fprintf(fid,'%s\n',datestr(now,0));
fprintf(fid,'%s\n',['Matlab ',char(version),' on ',computer]);
fprintf(fid,'%s\n',repmat('-',50,1));
fprintf(fid,'%s\n',x_lasterr);
fprintf(fid,'%s\n',y_lastwarn);
fprintf(fid,'%s\n',[' during ==> mi:',action]);
fprintf(fid,'%s',[' input ==> ',in.class]);
if ~isempty(in2), fprintf(fid,'\t%s\n',[' (',in2,')']); end
fprintf(fid,'%s\n',' errorcode ==> no errorcode available');
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});
fprintf(fid,'%s\n',['Name',char(9),'Size',char(9),'Bytes',char(9),'Class']);
for j=1:length(z_whos);
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])));
fprintf(fid,'\t%s',content(1:min([length(content),500])));
elseif strcmp(z_whos(j).class,'cell')
content=eval(z_whos(j).name);
fprintf(fid,'\t');
for j2=1:min([length(content),500])
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';
fprintf(fid,'\t%s',content(:)');
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);
fprintf(fid,'\n');
end
end
fclose(fid);
disp('----------------------------');
disp(' ERROR OCCURED');
disp(' during executing mi');
......@@ -663,18 +671,21 @@ catch
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','mi 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
try
if ~nogui
setptr(h_fig(1),'arrow')
h=findobj('tag','button_apply');
set(h(1),'ToolTip','Starts the computation.','String','Apply','Callback','mi 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
catch
end
if nargout, varargout={NaN}; end
try, set(0,props.root), end
try set(0,props.root),catch end
set(0,'ShowHidden','off')
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment