# HG changeset patch
# User samer
# Date 1444000057 -3600
# Node ID ce9021da6ce24010ba704321e6605e5e526670f1
# Parent  a5b8bd686246461c78a1291930d2fe541b1f8bd6
Added some histogram utilities.

diff -r a5b8bd686246 -r ce9021da6ce2 histogram/accumhist.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/histogram/accumhist.m	Mon Oct 05 00:07:37 2015 +0100
@@ -0,0 +1,15 @@
+function H=accumhist(I,W,D)
+% accumhist - Just like accumarray but filters out rows with nans or infs
+%
+% accumhist ::
+%    [[N,E]->(natural|nan|inf|-inf)] ~'N rows of E-dim array indices',
+%    ([[N]->real] | real) ~'the weights associated with each row',
+%    D:[[E]->natural]  ~'the size of the array to return'
+% -> [[D]].
+%
+% Note that unlike with accumarray and other Matlab builtins, the size 
+% argument D can be a single element, eg [M], indicating that a vector
+% is desired. accumarry would require [M 1].
+
+H = accumarray(I(all(isfinite(I),2),:),W,tosize(D));
+
diff -r a5b8bd686246 -r ce9021da6ce2 histogram/hist1d.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/histogram/hist1d.m	Mon Oct 05 00:07:37 2015 +0100
@@ -0,0 +1,38 @@
+function [NN,map1]=hist1d(X,m1spec,varargin)
+% hist1d - 1D histogram
+%
+% hist1d :: 
+%    [[N,1]]	~'N rows of 1D data',
+%    (K1:natural | dmap(K1)) ~'num bins for dim 1, or bins'
+% -> [[K1]->natural] ~'bin counts',
+%    dmap(K1)        ~'the map'.
+
+mins=min(X); maxs=max(X)+eps;
+map1=mkdmap(mins(1),maxs(1),m1spec);
+
+N=accumhist([map1(X(:,1))],1,[cardr(map1)]);
+if nargout==0,
+	opts=prefs('log',1,'plot','stairs','plotopts',{},varargin{:}); 
+	switch opts.plot
+		case 'stairs', plothist_stairs(map1,N,opts.plotopts{:});
+		case 'bars',   plothist_bars(map1,N,opts.plotopts{:});
+		case 'line',   plot(centres(map1),N,opts.plotopts{:});
+	end
+	if opts.log, semilgy; end
+	xlim(domain(map1));
+else
+	NN=N;
+end
+
+
+function M=mkdmap(min,max,spec)
+
+	if length(spec)==0, spec=32; end
+	if isa(spec,'dmap') % is alread a dmap
+		M=spec;
+	elseif length(spec)==1, % spec=number of bins
+		M=linmap(min,max,spec);
+	else
+		error('invalid discretisation map');
+	end
+
diff -r a5b8bd686246 -r ce9021da6ce2 histogram/hist2d.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/histogram/hist2d.m	Mon Oct 05 00:07:37 2015 +0100
@@ -0,0 +1,37 @@
+function [NN,map1,map2]=hist2d(X,m1spec,m2spec,varargin)
+% hist2d - 2D histogram
+%
+% hist2d :: 
+%    [[N,2]]	~'N rows of 2D data',
+%    (K1:natural | dmap(K1)) ~'num bins for dim 1, or bins',
+%    (K2:natural | dmap(K2)) ~'num bins for dim 1, or bins'
+% -> [[K1,K2]->natural] ~'bin counts',
+%    dmap(K1),
+%    dmap(K2).
+%
+% if second map spec is empty, it means use the same map
+% for both dimensions.
+
+mins=min(X); mm=max(X); maxs=mm+eps(mm);
+map1=mkdmap(mins(1),maxs(1),m1spec,[]);
+map2=mkdmap(mins(2),maxs(2),m2spec,m1spec);
+
+N=accumhist([map1(X(:,1)) map2(X(:,2))],1,[cardr(map1) cardr(map2)]);
+if nargout==0,
+	opts=prefs('cmap',@(t)log(1+t),varargin{:});
+	imagexy(centres(map2),centres(map1),opts.cmap(N));
+else
+	NN=N;
+end
+
+
+function M=mkdmap(min,max,spec,def)
+	if isempty(spec), spec=def; end
+	if isa(spec,'dmap') % is alread a dmap
+		M=spec;
+	elseif length(spec)==1, % spec=number of bins
+		M=linmap(min,max,spec);
+	else
+		error('invalid discretisation map');
+	end
+
diff -r a5b8bd686246 -r ce9021da6ce2 histogram/plothist_bars.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/histogram/plothist_bars.m	Mon Oct 05 00:07:37 2015 +0100
@@ -0,0 +1,20 @@
+function h=plothist_bars(Map,Counts,varargin)
+% plothist_bars - Plot histogram using bars.
+%
+% plothist_bars :: 
+%    dmap(N)    ~'the bin map',
+%    [[N]]      ~'the bin counts'
+% -> handle.
+%
+% The width of the bars matches the width of the discretisation bins.
+% The *area* of the bars (not the height) is proportional to bin counts.
+% Any extra arguments are passed to FILL.
+
+if isvector(Counts), Counts=Counts(:)'; end
+
+x=bins(Map);
+y=Counts./diff(x);
+
+h=fill([x;flipud(x)],[zeros(2,size(y,2));y;y],1, ...
+	'EdgeColor',get(gca,'XColor'),varargin{:});
+caxis([0,2]);
diff -r a5b8bd686246 -r ce9021da6ce2 histogram/plothist_lin.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/histogram/plothist_lin.m	Mon Oct 05 00:07:37 2015 +0100
@@ -0,0 +1,28 @@
+function z=plothist_lin(h,m)
+% plothist_lin - Plot histogram assuming triangular kernel (linear interpolation)
+%
+% plothist_lin :: [[N]->natural], dmap(N-1) -> handle.
+
+	e=edges(m);
+	e=[2*e(1)-e(2),e,2*e(end)-e(end-1)];
+
+	z=ezplot(@pdf,[e(1),e(end)]);
+
+	function p=pdf(x)
+%		if ~isscalar(x), error('must be scalar'); end
+%		j=find(x>e(1:end-2) & x<e(3:end));
+		p=0;
+		for i=1:length(h) %j
+			p=p+triangle(e(i:i+2),h(i));
+		end
+
+		function y=triangle(a,h)
+			y=zeros(size(x));
+			I1=(x>a(1) & x<=a(2));
+			I2=(x>a(2) & x<a(3));
+			y(I1)=h*((x(I1)-a(1))/(a(2)-a(1)));
+			y(I2)=h*((x(I2)-a(3))/(a(2)-a(3)));
+		end
+	end
+end
+
diff -r a5b8bd686246 -r ce9021da6ce2 histogram/plothist_stairs.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/histogram/plothist_stairs.m	Mon Oct 05 00:07:37 2015 +0100
@@ -0,0 +1,12 @@
+function h=plothist_stairs(Map,Counts,varargin)
+% plothist_stairs - Plot histogram using stairs.
+%
+% plothist_stairs :: 
+%    dmap(N)    ~'the bin map',
+%    [[N]]      ~'the bin counts'
+% -> handle.
+%
+% Any extra arguments are passed to STAIRS.
+
+if isvector(Counts), Counts=Counts(:); end
+stairs(edges(Map),[Counts;Counts(end,:)],varargin{:});