Mercurial > hg > ltpda

comparison m-toolbox/classes/@ao/bin_data.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % BIN_DATA rebins aos data, on logarithmic scale, linear scale, or arbitrarly chosen.
2 % The rebinning is done taking the mean of the bins included in the range
3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
4 %
5 % DESCRIPTION: BIN_DATA rebins aos data, on logarithmic scale, linear scale, or arbitrarly chosen.
6 % The rebinning is done taking the mean of the bins included in the range
7 %
8 % CALL: bs = bin_data(a1,a2,a3,...,pl)
9 % bs = bin_data(as,pl)
10 % bs = as.bin_data(pl)
11 %
12 % INPUTS: aN - input analysis objects
13 % as - input analysis objects array
14 % pl - input parameter list
15 %
16 % OUTPUTS: bs - array of analysis objects, one for each input
17 %
18 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'bin_data')">Parameters Description</a>
19 %
20 % The code is inherited from D Nicolodi, UniTN
21 %
22 % VERSION: $Id: bin_data.m,v 1.20 2011/05/10 16:46:48 mauro Exp $
23 %
24 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
25
26
27
28 function varargout = bin_data(varargin)
29
30 % Check if this is a call for parameters
31 if utils.helper.isinfocall(varargin{:})
32 varargout{1} = getInfo(varargin{3});
33 return
34 end
35
36 import utils.const.*
37 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
38
39 % Collect input variable names
40 in_names = cell(size(varargin));
41 for ii = 1:nargin,in_names{ii} = inputname(ii);end
42
43 % Collect all AOs
44 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
45
46 % Decide on a deep copy or a modify
47 bs = copy(as, nargout);
48
49 % Apply defaults to plist
50 usepl = applyDefaults(getDefaultPlist(), varargin{:});
51
52 x_scale = find(usepl, 'x_scale', find(usepl, 'xscale'));
53 x_vals = find(usepl, 'x_vals', find(usepl, 'xvals'));
54 resolution = find(usepl, 'resolution');
55 range = find(usepl, 'range');
56 method = lower(find(usepl, 'method'));
57 inherit_dy = utils.prog.yes2true(find(usepl, 'inherit-dy', find(usepl, 'inherit_dy')));
58
59 % Loop over input AOs
60 for jj = 1:numel(bs)
61
62 % check input data
63 if isa(bs(jj).data, 'data2D')
64
65 w = find(usepl, 'weights');
66
67 if isa(w, 'ao')
68 w = w.y;
69 end
70
71 if isempty(w)
72 w = 1./(bs(jj).dy).^2;
73 end
74
75 if isempty(x_vals)
76 if isempty(x_scale) || isempty(resolution)
77 error('### Please specify a scale and density for binning, OR the list of the values to bin around');
78 else
79
80 switch lower(x_scale)
81 case {'lin', 'linear'}
82 % Case of linear binning
83 % number of bins in the rebinned data set
84 N = resolution;
85
86 % maximum and minimum x
87 if ~isempty(range) && isfinite(range(1))
88 xmin = range(1);
89 else
90 xmin = min(bs(jj).x);
91 end
92
93 if ~isempty(range) && isfinite(range(2))
94 xmax = range(2);
95 else
96 xmax = max(bs(jj).x);
97 end
98
99 dx = (xmax - xmin)/N;
100
101 x_min = bs(jj).x(1) + dx*(0:(N-1))';
102 x_max = bs(jj).x(1) + dx*(1:N)';
103
104 case {'log', 'logarithmic'}
105 % Case of log-based binning
106
107 % maximum and minimum x
108 if ~isempty(range) && isfinite(range(1))
109 xmin = range(1);
110 else
111 xmin = min(bs(jj).x(bs(jj).x > 0));
112 end
113
114 if ~isempty(range) && isfinite(range(2))
115 xmax = range(2);
116 else
117 xmax = max(bs(jj).x);
118 end
119
120 alph = 10^(1/resolution);
121
122 % number of bins in the rebinned data set
123 N = ceil(log10(xmax/xmin) * resolution);
124
125 % maximum and minimum x-value for each bin
126 x_min = xmin*alph.^(0:(N-1))';
127 x_max = xmin*alph.^(1:N)';
128 otherwise
129 error(['### Unknown scaling option ' x_scale '. Please choose between ''lin'' and ''log']);
130 end
131 end
132 else
133 % number of bins in the rebinned data set
134 % If the x-scale is an AO, then take the x values
135 if isa(x_vals, 'ao')
136 if eq(x_vals.xunits, bs(jj).xunits)
137 x_vals = x_vals.x;
138 else
139 error('x_vals AO and data AO have different x-units');
140 end
141 elseif ~isnumeric(x_vals)
142 error('Unsupported x_vals object');
143 end
144 N = length(x_vals) - 1;
145 x_min = x_vals(1:N);
146 x_max = x_vals(2:N+1);
147 end
148
149 x = bs(jj).x;
150 y = bs(jj).y;
151 dy = bs(jj).dy;
152
153 % preallocate output vectors
154 xr = zeros(N, 1);
155 yr = zeros(N, size(y, 2));
156 if strcmpi(method, 'mean') || strcmpi(method, 'wmean')
157 dyr = zeros(N, size(y, 2));
158 else
159 dyr = [];
160 end
161 nr = zeros(N, 1);
162
163 % compute the averages
164 for kk = 1:N
165 in = x >= x_min(kk) & x < x_max(kk);
166 if any(in)
167 nr(kk) = sum(in); % number of points averaged in this bin
168
169 switch method
170 case {'mean', 'median', 'max', 'min', 'rms'}
171 xr(kk) = feval(method, x(in)); % rebinned x bins;
172 yr(kk) = feval(method, y(in)); % rebinned y bins;
173 if strcmpi(method, 'mean')
174 dyr(kk) = std(y(in), 0)/sqrt(nr(kk));
175 % check for zeros in the uncertainty and replace it with the individual point uncertainty
176 if dyr(kk) == 0
177 if inherit_dy && ~isempty(dy)
178 dyr(kk) = mean(dy(in));
179 else
180 dyr(kk) = Inf;
181 end
182 end
183 end
184 case {'wmean'}
185 xr(kk) = mean(x(in)); % rebinned x bins;
186 yr(kk) = sum(y(in).*w(in))./sum(w(in)); % rebinned y bins;
187 dyr(kk) = 1./sqrt(sum(w(in))); % rebinned dy bins;
188 otherwise
189 error(['### Unsupported method ' method]);
190 end
191 end
192 end
193
194 % remove bins where we do not have nothing to average
195 in = nr ~= 0;
196 nr = nr(in);
197 xr = xr(in);
198 yr = yr(in,:);
199 if strcmpi(method, 'mean') || strcmpi(method, 'wmean')
200 dyr = dyr(in,:);
201 end
202
203 % set the new object data
204 bs(jj).setXY(xr, yr);
205 bs(jj).setDy(dyr);
206
207 % nr goes into the procinfo
208 bs(jj).procinfo = plist('navs', nr);
209
210 % set name
211 bs(jj).name = sprintf('bin_data(%s)', ao_invars{jj});
212 % Add history
213 bs(jj).addHistory(getInfo('None'), usepl, ao_invars(jj), bs(jj).hist);
214 else
215 warning('### Ignoring input AO number %d (%s); it is not a 2D data object.', jj, bs(jj).name)
216 end
217 end % loop over analysis objects
218
219 % Set output
220 varargout = utils.helper.setoutputs(nargout, bs);
221 end
222
223 %--------------------------------------------------------------------------
224 % Get Info Object
225 %--------------------------------------------------------------------------
226 function ii = getInfo(varargin)
227 if nargin == 1 && strcmpi(varargin{1}, 'None')
228 sets = {};
229 pl = [];
230 else
231 sets = {'Default'};
232 pl = getDefaultPlist();
233 end
234 % Build info object
235 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: bin_data.m,v 1.20 2011/05/10 16:46:48 mauro Exp $', sets, pl);
236 end
237
238 %--------------------------------------------------------------------------
239 % Get Default Plist
240 %--------------------------------------------------------------------------
241
242 function plout = getDefaultPlist()
243 persistent pl;
244 if ~exist('pl', 'var') || isempty(pl)
245 pl = buildplist();
246 end
247 plout = pl;
248 end
249
250 function pl = buildplist()
251
252 pl = plist();
253
254 % method
255 p = param({'method',['method for binning. Choose from:<ul>', ...
256 '<li>mean</li>', ...
257 '<li>median</li>', ...
258 '<li>max</li>', ...
259 '<li>min</li>', ...
260 '<li>rms</li>', ...
261 '<li>weighted mean (weights can be input or are taken from data dy)</li></ul>']}, ...
262 {1, {'MEAN', 'MEDIAN', 'MAX', 'MIN', 'RMS', 'WMEAN'}, paramValue.SINGLE});
263 pl.append(p);
264
265 % x-scale
266 p = param({'xscale',['scaling of binning. Choose from:<ul>', ...
267 '<li>log - logaritmic</li>', ...
268 '<li>lin - linear</li></ul>']}, {1, {'LOG', 'LIN'}, paramValue.SINGLE});
269 pl.append(p);
270
271 % resolution
272 p = param({'resolution',['When setting logaritmic x scale, it sets the number of points per decade.<br>' ...
273 'When setting linear x scale, it sets the number of points.']}, paramValue.DOUBLE_VALUE(10));
274 pl.append(p);
275
276 % x_vals
277 p = param({'xvals',['List of x values to evaluate the binning between.<br>', ...
278 'It may be a vector or an ao, in which case it will take the x field']}, paramValue.DOUBLE_VALUE([]));
279 pl.append(p);
280
281 % weights
282 p = param({'weights', ['List of weights for the case of weighted mean.<br>', ...
283 'If empty, weights will be taken from object(s) dy field as w = 1/dy^2']}, paramValue.DOUBLE_VALUE([]));
284 pl.append(p);
285
286 % range
287 p = param({'range', ['Range of x where to operate.<br>', ...
288 'If empty, the whole data set will be used']}, paramValue.DOUBLE_VALUE([]));
289 pl.append(p);
290
291 % inherit_dy
292 p = param({'inherit_dy', ['Choose what to do in the case of mean, and bins with only one point. Choose from:<ul>', ...
293 '<li>''yes'' - take the uncertainty from the original data, if defined</li>', ...
294 '<li>''no'' - set it to Inf so it weighs 0 in averaged means</li></ul>' ...
295 ]}, paramValue.YES_NO);
296 pl.append(p);
297 end
298 % END