Line data Source code
1 : //# ProfileFit1D.cc: Class to fit Spectral components to vectors
2 : //# Copyright (C) 2004
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: ProfileFit1D.tcc 20652 2009-07-06 05:04:32Z Malte.Marquarding $
27 :
28 : #include <components/SpectralComponents/ProfileFit1D.h>
29 :
30 : #include <casacore/casa/Arrays/ArrayLogical.h>
31 : #include <casacore/casa/Exceptions/Error.h>
32 : #include <components/SpectralComponents/SpectralEstimate.h>
33 : #include <components/SpectralComponents/SpectralElement.h>
34 : #include <casacore/casa/Utilities/DataType.h>
35 : #include <casacore/casa/Utilities/Assert.h>
36 :
37 : namespace casa {
38 :
39 : template <class T>
40 545856 : ProfileFit1D<T>::ProfileFit1D()
41 : {
42 545856 : checkType();
43 545856 : }
44 :
45 : template <class T>
46 : ProfileFit1D<T>::ProfileFit1D(const ProfileFit1D& other)
47 : {
48 : checkType();
49 : copy(other);
50 : }
51 :
52 : template <class T>
53 545770 : ProfileFit1D<T>& ProfileFit1D<T>::operator=(const ProfileFit1D& other)
54 : {
55 545770 : if (this != &other) {
56 545770 : copy(other);
57 : }
58 545770 : return *this;
59 : }
60 :
61 : template <class T>
62 545856 : ProfileFit1D<T>::~ProfileFit1D()
63 545856 : {;}
64 :
65 : template <class T>
66 545684 : bool ProfileFit1D<T>::setData (const casacore::Vector<casacore::Double>& x, const casacore::Vector<T>& y,
67 : const casacore::Vector<casacore::Bool>& mask, const casacore::Vector<casacore::Double>& weight)
68 : {
69 545684 : if (x.nelements()==0) {
70 0 : itsError = "The X vector must have some elements";
71 0 : return false;
72 : }
73 545684 : const casacore::uInt n = x.nelements();
74 : //
75 545684 : if (y.nelements() != n) {
76 0 : itsError = "The Y vector must have the same number of elements as the X vector";
77 0 : return false;
78 : }
79 : //
80 545684 : if (weight.nelements() != n && weight.nelements()!=0) {
81 0 : itsError = "The weights vector must have the same number of elements as the X vector";
82 0 : return false;
83 : }
84 : //
85 545684 : if (mask.nelements() != n && mask.nelements() != 0) {
86 0 : itsError = "The mask vector must have the same number of elements (or zero) as the data";
87 0 : return false;
88 : }
89 : //
90 545684 : itsX.resize(n);
91 545684 : itsY.resize(n);
92 545684 : itsDataMask.resize(n);
93 : //
94 545684 : itsX = x;
95 545684 : itsY = y;
96 : //
97 545684 : if (weight.nelements()==0) {
98 0 : itsWeight.resize(0);
99 : } else {
100 545684 : itsWeight.resize(n);
101 545684 : itsWeight = weight;
102 : }
103 : //
104 545684 : if (mask.nelements()==0) {
105 0 : itsDataMask = true;
106 : } else {
107 545684 : itsDataMask = mask;
108 : }
109 545684 : return true;
110 : }
111 :
112 :
113 : template <class T>
114 0 : bool ProfileFit1D<T>::setData (const casacore::Vector<casacore::Double>& x, const casacore::Vector<T>& y,
115 : const casacore::Vector<casacore::Bool>& mask)
116 : {
117 0 : casacore::Vector<casacore::Double> weight;
118 0 : return setData (x, y, mask, weight);
119 : }
120 :
121 : template <class T>
122 : bool ProfileFit1D<T>::setData (const casacore::Vector<casacore::Double>& x, const casacore::Vector<T>& y)
123 : {
124 : casacore::Vector<casacore::Bool> mask(x.nelements(), true);
125 : return setData (x, y, mask);
126 : }
127 :
128 :
129 : template <class T>
130 546314 : void ProfileFit1D<T>::setElements (const SpectralList& list)
131 : {
132 546314 : itsList = list;
133 546314 : }
134 :
135 : template <class T>
136 16211 : bool ProfileFit1D<T>::setGaussianElements (casacore::uInt nGauss)
137 : {
138 16211 : if (nGauss==0) {
139 0 : itsError = "You must specify some Gaussian components";
140 0 : return false;
141 : }
142 : //
143 16211 : if (itsY.nelements()==0) {
144 0 : itsError = "You must call function setData to set some data first";
145 0 : return false;
146 : }
147 :
148 : // Clear list
149 :
150 16211 : itsList.clear();
151 :
152 : // Make estimate for Gaussians.
153 :
154 32422 : SpectralEstimate estimator (nGauss);
155 16211 : estimator.setQ(5);
156 16211 : SpectralList listGauss = estimator.estimate (itsX, itsY); // Ignores masked data
157 16211 : itsList.add (listGauss);
158 16211 : return true;
159 : }
160 :
161 : template <class T>
162 : void ProfileFit1D<T>::addElements (const SpectralList& list)
163 : {
164 : itsList.add (list);
165 : }
166 :
167 : template <class T>
168 529010 : void ProfileFit1D<T>::addElement (const SpectralElement& el)
169 : {
170 529010 : itsList.add (el);
171 529010 : }
172 :
173 :
174 : template <class T>
175 545684 : void ProfileFit1D<T>::clearList ()
176 : {
177 545684 : itsList.clear();
178 545684 : }
179 :
180 :
181 : template <class T>
182 : bool ProfileFit1D<T>::setXRangeMask (const casacore::Vector<casacore::uInt>& start,
183 : const casacore::Vector<casacore::uInt>& end,
184 : casacore::Bool insideIsGood)
185 : {
186 : AlwaysAssert(start.nelements()==end.nelements(), casacore::AipsError);
187 : if (itsX.nelements()==0) {
188 : itsError = "You must call function setData to set some data first";
189 : return false;
190 : }
191 : //
192 : const casacore::uInt n = itsX.nelements();
193 : itsRangeMask.resize(n);
194 : casacore::Bool value = !insideIsGood;
195 : itsRangeMask = value;
196 : //
197 : for (casacore::uInt i=0; i<start.nelements(); i++) {
198 : if (start[i] > end[i]) {
199 : itsError = "The start index must be < the end index";
200 : return false;
201 : }
202 : if (start[i]>=n) {
203 : itsError = "The start index must be in the range 0->nElements-1";
204 : return false;
205 : }
206 : if (end[i]>=n) {
207 : itsError = "The end index must be in the range 0->nElements-1";
208 : return false;
209 : }
210 : //
211 : for (casacore::uInt j=start[i]; j<end[i]+1; j++) {
212 : itsRangeMask[j] = !value;
213 : }
214 : }
215 : return true;
216 : }
217 :
218 :
219 :
220 : template <class T>
221 : bool ProfileFit1D<T>::setXRangeMask (const casacore::Vector<T>& start,
222 : const casacore::Vector<T>& end,
223 : casacore::Bool insideIsGood)
224 : {
225 : if (start.nelements()!=end.nelements()) {
226 : itsError = "Start and end vectors must be the same length";
227 : return false;
228 : }
229 : if (itsX.nelements()==0) {
230 : itsError = "You must call function setData to set some data first";
231 : return false;
232 : }
233 : //
234 : const casacore::uInt n = itsX.nelements();
235 : itsRangeMask.resize(n);
236 : casacore::Bool value = !insideIsGood;
237 : itsRangeMask = value;
238 : //
239 : casacore::Vector<casacore::uInt> startIndex(start.nelements());
240 : casacore::Vector<casacore::uInt> endIndex(end.nelements());
241 :
242 : for (casacore::uInt i=0; i<start.nelements(); i++) {
243 : if (start[i] > end[i]) {
244 : itsError = "The start range must be < the end range";
245 : return false;
246 : }
247 : if (start[i]<itsX[0] || start[i]>itsX[n-1]) {
248 : itsError = "The start range must be in the X-range of the data";
249 : return false;
250 : }
251 : if (end[i]<itsX[0] || end[i]>itsX[n-1]) {
252 : itsError = "The end range must be in the X-range of the data";
253 : return false;
254 : }
255 :
256 : // Find the indices for this range
257 :
258 : casacore::Bool doneStart = false;
259 : casacore::Bool doneEnd = false;
260 : for (casacore::uInt j=0; j<n; j++) {
261 : if (!doneStart && itsX[j] >= start[i]) {
262 : startIndex[i] = j;
263 : doneStart = true;
264 : }
265 : if (!doneEnd && itsX[j] >= end[i]) {
266 : endIndex[i] = j;
267 : doneEnd = true;
268 : }
269 : if (!doneEnd) endIndex[i] = n-1;
270 : }
271 : }
272 : //
273 : return setXRangeMask (startIndex, endIndex);
274 : }
275 :
276 : template <class T>
277 4641 : bool ProfileFit1D<T>::setXMask(const std::set<casacore::uInt>& indices, casacore::Bool specifiedPixelsAreGood) {
278 4641 : const casacore::uInt n = itsX.nelements();
279 4641 : ThrowIf(n == 0, "Logic Error: setData() must be called prior to setRangeMask()");
280 4641 : itsRangeMask.resize(n);
281 4641 : itsRangeMask = ! specifiedPixelsAreGood;
282 4641 : if (indices.empty()) {
283 0 : return true;
284 : }
285 4641 : std::set<casacore::uInt>::const_iterator iter = indices.begin();
286 4641 : std::set<casacore::uInt>::const_iterator end = indices.end();
287 :
288 32487 : while (iter != end && *iter < n) {
289 27846 : itsRangeMask[*iter] = specifiedPixelsAreGood;
290 27846 : ++iter;
291 : }
292 4641 : return true;
293 : }
294 :
295 :
296 : template <class T>
297 545637 : bool ProfileFit1D<T>::fit ()
298 : {
299 : /*
300 : if (itsX.nelements()==0) {
301 : itsError = "You must call function setData to set some data first";
302 : return false;
303 : }
304 : if (itsList.nelements()==0) {
305 : itsError = "You must call function setElements to set some fit components first";
306 : return false;
307 : }
308 : */
309 :
310 : // Set list in fitter
311 545637 : itsFitter.clear();
312 545637 : itsFitter.addFitElement (itsList);
313 : // Do the fit with the total mask
314 :
315 545637 : casacore::Bool converged = itsWeight.empty()
316 1091368 : ? itsFitter.fit (itsY, itsX, makeTotalMask())
317 545637 : : itsFitter.fit (itsWeight, itsY, itsX, makeTotalMask());
318 545543 : return converged;
319 : }
320 :
321 : template <class T>
322 1137070 : const SpectralList& ProfileFit1D<T>::getList (casacore::Bool fit) const
323 : {
324 1137070 : if (fit) {
325 558959 : return itsFitter.list();
326 : } else {
327 578111 : return itsList;
328 : }
329 : }
330 :
331 :
332 : template <class T>
333 : casacore::Vector<T> ProfileFit1D<T>::getEstimate (casacore::Int which) const
334 : {
335 : casacore::Vector<T> tmp;
336 : if (itsX.nelements()==0) {
337 : itsError = "You must call function setData to set some data first";
338 : return tmp;
339 : }
340 : //
341 : if (which<0) {
342 : itsList.evaluate(tmp, itsX);
343 : } else {
344 : SpectralList list = getSubsetList (itsList, which);
345 : list.evaluate(tmp, itsX);
346 : }
347 : //
348 : return tmp;
349 : }
350 :
351 : template <class T>
352 529069 : casacore::Vector<T> ProfileFit1D<T>::getFit (casacore::Int which) const
353 : {
354 529069 : casacore::Vector<T> tmp;
355 529069 : if (itsX.nelements()==0) {
356 0 : itsError = "You must call function setData to set some data first";
357 0 : return tmp;
358 : }
359 : //
360 529069 : const SpectralList& fitList = itsFitter.list();
361 529069 : if (fitList.nelements()==0) {
362 0 : itsError = "You must call function fit first";
363 0 : return tmp;
364 : }
365 : //
366 :
367 529069 : if (which<0) {
368 529069 : fitList.evaluate(tmp, itsX);
369 : } else {
370 0 : SpectralList list = getSubsetList (fitList, which);
371 0 : list.evaluate(tmp, itsX);
372 : }
373 : //
374 529069 : return tmp;
375 : }
376 :
377 : template <class T>
378 529097 : casacore::Vector<T> ProfileFit1D<T>::getResidual (casacore::Int which, casacore::Bool fit) const
379 : {
380 529097 : casacore::Vector<T> tmp;
381 529097 : if (itsX.nelements()==0) {
382 0 : itsError = "You must call function setData to set some data first";
383 0 : return tmp;
384 : }
385 : //
386 1058194 : SpectralList list;
387 529097 : if (fit) {
388 529097 : list = itsFitter.list();
389 529097 : if (list.nelements()==0) {
390 0 : throw (casacore::AipsError("You must call function fit first"));
391 : }
392 : } else {
393 0 : list = itsList;
394 : }
395 : //
396 529097 : tmp = itsY;
397 529097 : if (which<0) {
398 529097 : list.residual(tmp, itsX);
399 : } else {
400 0 : SpectralList subList = getSubsetList (list, which);
401 0 : subList.residual(tmp, itsX);
402 : }
403 : //
404 529097 : return tmp;
405 : }
406 :
407 :
408 : // Private functions
409 :
410 :
411 : template <class T>
412 0 : SpectralList ProfileFit1D<T>::getSubsetList (const SpectralList& list, casacore::Int which) const
413 : {
414 0 : const casacore::Int n = list.nelements();
415 0 : if (which+1 > n) {
416 0 : throw casacore::AipsError("Illegal spectral element index");
417 : }
418 0 : SpectralList listOut;
419 0 : listOut.add(*list[which]);
420 0 : return listOut;
421 : }
422 :
423 : template <class T>
424 545637 : casacore::Vector<casacore::Bool> ProfileFit1D<T>::makeTotalMask () const
425 : {
426 545637 : casacore::Vector<casacore::Bool> mask;
427 545637 : if (itsRangeMask.nelements()==0) {
428 540996 : mask = itsDataMask;
429 : } else {
430 4641 : mask = itsDataMask && itsRangeMask;
431 : }
432 545637 : return mask;
433 : }
434 :
435 : template <class T>
436 545770 : void ProfileFit1D<T>::copy(const ProfileFit1D& other)
437 : {
438 545770 : itsX.resize(0);
439 545770 : itsX = other.itsX;
440 : //
441 545770 : itsY.resize(0);
442 545770 : itsY = other.itsY;
443 : //
444 545770 : itsWeight.resize(0);
445 545770 : itsWeight = other.itsWeight;
446 : //
447 545770 : itsDataMask.resize(0);
448 545770 : itsDataMask = other.itsDataMask;
449 : //
450 545770 : itsRangeMask.resize(0);
451 545770 : itsRangeMask = other.itsRangeMask;
452 : //
453 545770 : itsList = other.itsList;
454 : //
455 545770 : itsFitter = other.itsFitter;
456 : //
457 545770 : itsError = other.itsError;
458 545770 : }
459 :
460 : template <class T>
461 545856 : void ProfileFit1D<T>::checkType() const
462 : {
463 545856 : AlwaysAssert(casacore::whatType<T>()==casacore::TpDouble,casacore::AipsError);
464 545856 : }
465 :
466 : } //#End casa namespace
|
