1 |
gforget |
1.1 |
function [xh,yh,zh,varargout]=MITprof_pdf(x,xBinCenters,y,yBinCenters,varargin); |
2 |
|
|
%input: |
3 |
|
|
% x/y position vectors or 2D arrays |
4 |
|
|
% xBinCenters position bin centers |
5 |
|
|
% yBinCenters position bin centers |
6 |
|
|
% choiceStat [optional] type of statistic (pdf, hist, sum, mean, or var) |
7 |
|
|
% z [optional] data vector or 2D array |
8 |
|
|
% |
9 |
|
|
%note: all inputs are assumed to have be nan-masked accordingly |
10 |
|
|
% |
11 |
|
|
%output: |
12 |
|
|
% xh/yh are the position arrays for the historgam bin centers |
13 |
|
|
% zh is the 2D histogram |
14 |
|
|
|
15 |
|
|
warning('off','MATLAB:dsearch:DeprecatedFunction'); |
16 |
|
|
|
17 |
gforget |
1.3 |
gcmfaces_global; |
18 |
|
|
|
19 |
|
|
if ~isfield(myenv,'useDelaunayTri'); |
20 |
|
|
myenv.useDelaunayTri=~isempty(which('DelaunayTri')); |
21 |
|
|
end; |
22 |
|
|
|
23 |
gforget |
1.1 |
if nargin>4; choiceStat=varargin{1}; else; choiceStat='pdf'; end; |
24 |
|
|
if nargin>5; z=varargin{2}; else; z=[]; end; |
25 |
|
|
|
26 |
|
|
%take care of input dimensions: |
27 |
|
|
if size(x,1)==1; x=x'; end; if size(y,1)==1; y=y'; end; |
28 |
|
|
% |
29 |
|
|
if size(y,2)==1; %x is array => make sure y is so |
30 |
|
|
if size(x,2)==size(y,1); |
31 |
|
|
y=ones(size(x,1),1)*y'; |
32 |
|
|
else; |
33 |
|
|
y=y*ones(1,size(x,2)); |
34 |
|
|
end; |
35 |
|
|
end; |
36 |
|
|
% |
37 |
|
|
if size(x,2)==1; %y is array => make sure x is so |
38 |
|
|
if size(y,2)==size(x,1); |
39 |
|
|
x=ones(size(y,1),1)*x'; |
40 |
|
|
else; |
41 |
|
|
x=x*ones(1,size(y,2)); |
42 |
|
|
end; |
43 |
|
|
end; |
44 |
|
|
% |
45 |
|
|
if size(z,2)>1;%z is array => make sure x/y are so |
46 |
|
|
if size(x,2)==1&size(x,1)==size(z,1); x=x*ones(1,size(z,2)); end; |
47 |
|
|
if size(x,2)==1&size(x,1)==size(z,2); x=ones(size(z,1),1)*x'; end; |
48 |
|
|
% |
49 |
|
|
if size(y,2)==1&size(y,1)==size(z,1); y=y*ones(1,size(z,2)); end; |
50 |
|
|
if size(y,2)==1&size(y,1)==size(z,2); y=ones(size(z,1),1)*y'; end; |
51 |
|
|
end; |
52 |
|
|
if isempty(z); z=NaN*x; end;%to facilitate vector length tests |
53 |
|
|
%switch to vector form: |
54 |
|
|
x=x(:); y=y(:); z=z(:); |
55 |
|
|
%check vectors length: |
56 |
|
|
if length(x)~=length(y)|length(x)~=length(z)|length(y)~=length(z); |
57 |
|
|
error('inconsistent diemnsions of [x,y[,z]]'); |
58 |
|
|
end; |
59 |
|
|
|
60 |
|
|
%prepare output arrays: |
61 |
|
|
if size(xBinCenters,1)==1; xh=xBinCenters'; else; xh=xBinCenters; end; |
62 |
|
|
if size(yBinCenters,1)==1; yh=yBinCenters'; else; yh=yBinCenters; end; |
63 |
|
|
if size(xh,2)>1|size(yh,2)>1; error('inconsistent diemnsions of [xBinCenters yBinCenters]'); end; |
64 |
|
|
xh=xh*ones(1,size(yh,1)); |
65 |
|
|
yh=ones(size(xh,1),1)*yh'; |
66 |
|
|
|
67 |
|
|
%build the histogram |
68 |
gforget |
1.3 |
if myenv.useDelaunayTri; |
69 |
|
|
%the new way, using DelaunayTri&nearestNeighbor |
70 |
|
|
mytri.TRI=DelaunayTri(xh(:),yh(:)); |
71 |
|
|
else; |
72 |
|
|
TRI=delaunay(xh,yh,{'QJ'}); nxy = prod(size(xh)); |
73 |
|
|
Stri=sparse(TRI(:,[1 1 2 2 3 3]),TRI(:,[2 3 1 3 1 2]),1,nxy,nxy); |
74 |
|
|
end; |
75 |
gforget |
1.1 |
|
76 |
|
|
%compute grid point vector associated with lon/lat vectors |
77 |
|
|
if nargin>5; ii=find(~isnan(x.*y.*z)); else; ii=find(~isnan(x.*y)); end; |
78 |
|
|
x=x(ii); y=y(ii); z=z(ii); |
79 |
gforget |
1.3 |
if myenv.useDelaunayTri; |
80 |
|
|
ik = mytri.TRI.nearestNeighbor(x,y); |
81 |
|
|
else; |
82 |
|
|
ik = dsearch(xh,yh,TRI,x,y,Stri); |
83 |
|
|
end; |
84 |
gforget |
1.1 |
|
85 |
|
|
%compute bin sums: |
86 |
|
|
nh=zeros(size(xh)); if nargin>5; zh=nh; zh2=zh; end; |
87 |
|
|
for k=1:length(ik) |
88 |
|
|
nh(ik(k))=nh(ik(k))+1; |
89 |
|
|
if nargin>5; |
90 |
|
|
zh(ik(k))=zh(ik(k))+z(k); |
91 |
|
|
zh2(ik(k))=zh2(ik(k))+z(k)^2; |
92 |
|
|
end; |
93 |
|
|
end % k=1:length(ik) |
94 |
|
|
ii=find(nh(:)==0); nh(ii)=NaN; |
95 |
|
|
if nargin>5; zh(ii)=NaN; zh2(ii)=NaN; end; |
96 |
|
|
|
97 |
|
|
%compute output: |
98 |
|
|
if strcmp(choiceStat,'pdf'); %hist + normalization along second dimension |
99 |
|
|
zh=nh; zh(isnan(zh))=0; |
100 |
|
|
%normalize by total number of samples |
101 |
|
|
zh=zh./(nansum(zh,2)*ones(1,size(zh,2))); |
102 |
|
|
%normalize by delta(yh) => pdf |
103 |
|
|
tmp2=diff(yh,1,2); |
104 |
|
|
tmp2=0.5*(tmp2(:,1:end-1)+tmp2(:,2:end)); |
105 |
|
|
tmp2=[tmp2(:,1) tmp2 tmp2(:,end)]; |
106 |
|
|
zh=zh./tmp2; |
107 |
|
|
elseif strcmp(choiceStat,'hist'); |
108 |
|
|
zh=nh; |
109 |
|
|
elseif strcmp(choiceStat,'sum'); |
110 |
|
|
zh=zh; |
111 |
|
|
elseif strcmp(choiceStat,'mean'); |
112 |
|
|
zh=zh./nh; |
113 |
|
|
elseif strcmp(choiceStat,'var'); |
114 |
|
|
zh=zh2./nh-(zh./nh).^2; |
115 |
|
|
end; |
116 |
|
|
|
117 |
|
|
if nargout>3; |
118 |
|
|
varargout{1}=nh; |
119 |
|
|
else; |
120 |
|
|
varargout={}; |
121 |
|
|
end; |
122 |
|
|
|
123 |
|
|
%for debugging purposes: |
124 |
|
|
if 0; figure; pcolor(xh,yh,zh); colorbar; end; |
125 |
|
|
|
126 |
|
|
warning('on','MATLAB:dsearch:DeprecatedFunction'); |
127 |
|
|
|