source: branches/casa-release-4_3/external-alma/components/SpectralComponents/ProfileFit1D.tcc@ 3118

Last change on this file since 3118 was 2980, checked in by Malte Marquarding, 10 years ago

Add a copy of casacore/components/SpectralComponents to external-alma directory to prepare for its removal from casacore-trunk. DOn;t activate in SConscript yet.

File size: 10.7 KB
Line 
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 21229 2012-04-02 12:00:20Z gervandiepen $
27
28#include <components/SpectralComponents/ProfileFit1D.h>
29
30#include <casa/Arrays/ArrayLogical.h>
31#include <casa/Exceptions/Error.h>
32#include <components/SpectralComponents/SpectralEstimate.h>
33#include <components/SpectralComponents/SpectralElement.h>
34#include <casa/Utilities/DataType.h>
35#include <casa/Utilities/Assert.h>
36
37namespace casa {
38
39template <class T>
40ProfileFit1D<T>::ProfileFit1D()
41{
42 checkType();
43}
44
45template <class T>
46ProfileFit1D<T>::ProfileFit1D(const ProfileFit1D& other)
47{
48 checkType();
49 copy(other);
50}
51
52template <class T>
53ProfileFit1D<T>& ProfileFit1D<T>::operator=(const ProfileFit1D& other)
54{
55 if (this != &other) {
56 copy(other);
57 }
58 return *this;
59}
60
61template <class T>
62ProfileFit1D<T>::~ProfileFit1D()
63{;}
64
65template <class T>
66Bool ProfileFit1D<T>::setData (const Vector<Double>& x, const Vector<T>& y,
67 const Vector<Bool>& mask, const Vector<Double>& weight)
68{
69 if (x.nelements()==0) {
70 itsError = "The X vector must have some elements";
71 return False;
72 }
73 const uInt n = x.nelements();
74//
75 if (y.nelements() != n) {
76 itsError = "The Y vector must have the same number of elements as the X vector";
77 return False;
78 }
79//
80 if (weight.nelements() != n && weight.nelements()!=0) {
81 itsError = "The weights vector must have the same number of elements as the X vector";
82 return False;
83 }
84//
85 if (mask.nelements() != n && mask.nelements() != 0) {
86 itsError = "The mask vector must have the same number of elements (or zero) as the data";
87 return False;
88 }
89//
90 itsX.resize(n);
91 itsY.resize(n);
92 itsDataMask.resize(n);
93//
94 itsX = x;
95 itsY = y;
96//
97 if (weight.nelements()==0) {
98 itsWeight.resize(0);
99 } else {
100 itsWeight.resize(n);
101 itsWeight = weight;
102 }
103//
104 if (mask.nelements()==0) {
105 itsDataMask = True;
106 } else {
107 itsDataMask = mask;
108 }
109 return True;
110}
111
112
113template <class T>
114Bool ProfileFit1D<T>::setData (const Vector<Double>& x, const Vector<T>& y,
115 const Vector<Bool>& mask)
116{
117 Vector<Double> weight;
118 return setData (x, y, mask, weight);
119}
120
121template <class T>
122Bool ProfileFit1D<T>::setData (const Vector<Double>& x, const Vector<T>& y)
123{
124 Vector<Bool> mask(x.nelements(), True);
125 return setData (x, y, mask);
126}
127
128
129template <class T>
130void ProfileFit1D<T>::setElements (const SpectralList& list)
131{
132 itsList = list;
133}
134
135template <class T>
136Bool ProfileFit1D<T>::setGaussianElements (uInt nGauss)
137{
138 if (nGauss==0) {
139 itsError = "You must specify some Gaussian components";
140 return False;
141 }
142//
143 if (itsY.nelements()==0) {
144 itsError = "You must call function setData to set some data first";
145 return False;
146 }
147
148// Clear list
149
150 itsList.clear();
151
152// Make estimate for Gaussians.
153
154 SpectralEstimate estimator (nGauss);
155 estimator.setQ(5);
156 SpectralList listGauss = estimator.estimate (itsX, itsY); // Ignores masked data
157 itsList.add (listGauss);
158 return True;
159}
160
161template <class T>
162void ProfileFit1D<T>::addElements (const SpectralList& list)
163{
164 itsList.add (list);
165}
166
167template <class T>
168void ProfileFit1D<T>::addElement (const SpectralElement& el)
169{
170 itsList.add (el);
171}
172
173
174template <class T>
175void ProfileFit1D<T>::clearList ()
176{
177 itsList.clear();
178}
179
180
181template <class T>
182Bool ProfileFit1D<T>::setRangeMask (const Vector<uInt>& start,
183 const Vector<uInt>& end,
184 Bool insideIsGood)
185{
186 AlwaysAssert(start.nelements()==end.nelements(), 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 uInt n = itsX.nelements();
193 itsRangeMask.resize(n);
194 Bool value = !insideIsGood;
195 itsRangeMask = value;
196//
197 for (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 (uInt j=start[i]; j<end[i]+1; j++) {
212 itsRangeMask[j] = !value;
213 }
214 }
215 return True;
216}
217
218
219
220template <class T>
221Bool ProfileFit1D<T>::setRangeMask (const Vector<T>& start,
222 const Vector<T>& end,
223 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 uInt n = itsX.nelements();
235 itsRangeMask.resize(n);
236 Bool value = !insideIsGood;
237 itsRangeMask = value;
238//
239 Vector<uInt> startIndex(start.nelements());
240 Vector<uInt> endIndex(end.nelements());
241
242 for (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 Bool doneStart = False;
259 Bool doneEnd = False;
260 for (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 setRangeMask (startIndex, endIndex);
274}
275
276
277template <class T>
278Bool ProfileFit1D<T>::fit ()
279{
280 if (itsX.nelements()==0) {
281 itsError = "You must call function setData to set some data first";
282 return False;
283 }
284 if (itsList.nelements()==0) {
285 itsError = "You must call function setElements to set some fit components first";
286 return False;
287 }
288
289// Set list in fitter
290 itsFitter.clear();
291 itsFitter.addFitElement (itsList);
292// Do the fit with the total mask
293
294 Bool converged(False);
295 if (itsWeight.nelements()==0) {
296 converged = itsFitter.fit (itsY, itsX, makeTotalMask());
297 } else {
298 converged = itsFitter.fit (itsWeight, itsY, itsX, makeTotalMask());
299 }
300 return converged;
301}
302
303template <class T>
304const SpectralList& ProfileFit1D<T>::getList (Bool fit) const
305{
306 if (fit) {
307 return itsFitter.list();
308 } else {
309 return itsList;
310 }
311}
312
313
314template <class T>
315Vector<T> ProfileFit1D<T>::getEstimate (Int which) const
316{
317 Vector<T> tmp;
318 if (itsX.nelements()==0) {
319 itsError = "You must call function setData to set some data first";
320 return tmp;
321 }
322//
323 if (which<0) {
324 itsList.evaluate(tmp, itsX);
325 } else {
326 SpectralList list = getSubsetList (itsList, which);
327 list.evaluate(tmp, itsX);
328 }
329//
330 return tmp;
331}
332
333template <class T>
334Vector<T> ProfileFit1D<T>::getFit (Int which) const
335{
336 Vector<T> tmp;
337 if (itsX.nelements()==0) {
338 itsError = "You must call function setData to set some data first";
339 return tmp;
340 }
341//
342 const SpectralList& fitList = itsFitter.list();
343 if (fitList.nelements()==0) {
344 itsError = "You must call function fit first";
345 return tmp;
346 }
347//
348
349 if (which<0) {
350 fitList.evaluate(tmp, itsX);
351 } else {
352 SpectralList list = getSubsetList (fitList, which);
353 list.evaluate(tmp, itsX);
354 }
355//
356 return tmp;
357}
358
359template <class T>
360Vector<T> ProfileFit1D<T>::getResidual (Int which, Bool fit) const
361{
362 Vector<T> tmp;
363 if (itsX.nelements()==0) {
364 itsError = "You must call function setData to set some data first";
365 return tmp;
366 }
367//
368 SpectralList list;
369 if (fit) {
370 list = itsFitter.list();
371 if (list.nelements()==0) {
372 throw (AipsError("You must call function fit first"));
373 }
374 } else {
375 list = itsList;
376 }
377//
378 tmp = itsY;
379 if (which<0) {
380 list.residual(tmp, itsX);
381 } else {
382 SpectralList subList = getSubsetList (list, which);
383 subList.residual(tmp, itsX);
384 }
385//
386 return tmp;
387}
388
389
390// Private functions
391
392
393template <class T>
394SpectralList ProfileFit1D<T>::getSubsetList (const SpectralList& list, Int which) const
395{
396 const Int n = list.nelements();
397 if (which+1 > n) {
398 throw AipsError("Illegal spectral element index");
399 }
400 SpectralList listOut;
401 listOut.add(*list[which]);
402 return listOut;
403}
404
405template <class T>
406Vector<Bool> ProfileFit1D<T>::makeTotalMask () const
407{
408 Vector<Bool> mask;
409 if (itsRangeMask.nelements()==0) {
410 mask = itsDataMask;
411 } else {
412 mask = itsDataMask && itsRangeMask;
413 }
414 return mask;
415}
416
417template <class T>
418void ProfileFit1D<T>::copy(const ProfileFit1D& other)
419{
420 itsX.resize(0);
421 itsX = other.itsX;
422//
423 itsY.resize(0);
424 itsY = other.itsY;
425//
426 itsWeight.resize(0);
427 itsWeight = other.itsWeight;
428//
429 itsDataMask.resize(0);
430 itsDataMask = other.itsDataMask;
431//
432 itsRangeMask.resize(0);
433 itsRangeMask = other.itsRangeMask;
434//
435 itsList = other.itsList;
436//
437 itsFitter = other.itsFitter;
438//
439 itsError = other.itsError;
440}
441
442template <class T>
443void ProfileFit1D<T>::checkType() const
444{
445 T* p=0;
446 AlwaysAssert(whatType(p)==TpDouble,AipsError);
447}
448
449} //#End casa namespace
Note: See TracBrowser for help on using the repository browser.