source: branches/casa-release-4_3-test02/external-alma/components/SpectralComponents/test/tProfileFit1D.cc@ 3032

Last change on this file since 3032 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: 16.5 KB
Line 
1//# tProfileFit1D.cc: test the ProfileFit1D class
2//# Copyright (C) 1995,1996,1998,1999,2000,2001,2002,2004
3//# Associated Universities, Inc. Washington DC, USA.
4//#
5//# This program is free software; you can redistribute it and/or modify it
6//# under the terms of the GNU General Public License as published by the Free
7//# Software Foundation; either version 2 of the License, or (at your option)
8//# any later version.
9//#
10//# This program 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 General Public License for
13//# more details.
14//#
15//# You should have received a copy of the GNU General Public License along
16//# with this program; if not, write to the Free Software Foundation, Inc.,
17//# 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: tProfileFit1D.cc 21451 2014-06-10 07:48:08Z gervandiepen $
27
28#include <casa/aips.h>
29#include <components/SpectralComponents/ProfileFit1D.h>
30
31#include <casa/Containers/Record.h>
32#include <casa/Arrays/ArrayLogical.h>
33#include <casa/Arrays/ArrayMath.h>
34#include <casa/Arrays/IPosition.h>
35#include <casa/Arrays/Vector.h>
36#include <casa/Utilities/Assert.h>
37#include <components/SpectralComponents/GaussianMultipletSpectralElement.h>
38
39#include <casa/iostream.h>
40
41#include <casa/namespace.h>
42void makeData (Vector<Double>& x, Vector<Double>& y, Vector<Bool>& m,
43 Double& amp, Double& cen, Double& sigma,
44 Double& p0, Double& p1);
45void check (Double amp, Double cen, Double sigma, Double p0, Double p1,
46 const SpectralList& l);
47void checkMasks (uInt n, const ProfileFit1D<Double>& fitter, Int start,
48 Int end);
49
50GaussianMultipletSpectralElement makeMultiplet (
51 Vector<Double>& x, Vector<Double>& y,
52 const Vector<Double>& amp, const Vector<Double>& cen,
53 const Vector<Double>& sigma
54);
55
56vector<LorentzianSpectralElement> makeLorentzians (
57 Vector<Double>& x, Vector<Double>& y,
58 const Vector<Double>& amp, const Vector<Double>& cen,
59 const Vector<Double>& fwhm
60);
61
62PowerLogPolynomialSpectralElement makePowerLogPoly(
63 Vector<Double>&x, Vector<Double>& y, const Vector<Double>& coeffs
64);
65
66PolynomialSpectralElement makePoly(
67 Vector<Double>&x, Vector<Double>& y, const Vector<Double>& coeffs
68);
69
70int main() {
71
72 try {
73 {
74 // Data
75 Vector<Double> x,y;
76 Vector<Bool> m;
77 Double amp, cen, sig, p0, p1;
78 makeData(x, y, m, amp, cen, sig, p0, p1);
79 const uInt n = x.nelements();
80
81 // Make fitter, set data and fit
82 ProfileFit1D<Double> fitter;
83 fitter.setData (x,y,m);
84 fitter.setGaussianElements (1);
85 const SpectralElement *firstEl = fitter.getList(False)[0];
86 AlwaysAssert(
87 firstEl->getType() == SpectralElement::GAUSSIAN, AipsError
88 );
89 PolynomialSpectralElement p(1);
90 fitter.addElement(p);
91 AlwaysAssert(fitter.fit(), AipsError);
92 // Check ok
93
94 AlwaysAssert(fitter.getDataMask().nelements()==n, AipsError);
95 AlwaysAssert(allEQ(fitter.getDataMask(), True), AipsError);
96 AlwaysAssert(fitter.getRangeMask().nelements()==0, AipsError);
97 AlwaysAssert(fitter.getTotalMask().nelements()==n, AipsError);
98 AlwaysAssert(allEQ(fitter.getTotalMask(), True), AipsError);
99 {
100 const SpectralList& fitList = fitter.getList(True);
101 check (amp, cen, sig, p0, p1, fitList);
102 }
103 {
104 ProfileFit1D<Double> fitter2(fitter);
105 const SpectralList& fitList = fitter2.getList(True);
106 check (amp, cen, sig, p0, p1, fitList);
107 }
108 {
109 ProfileFit1D<Double> fitter2;
110 fitter2 = fitter;
111 const SpectralList& fitList = fitter2.getList(True);
112 check (amp, cen, sig, p0, p1, fitList);
113 }
114 // Set a range mask via indices
115
116 {
117 Vector<uInt> start(1), end(1);
118 start(0) = n/2; end(0) = start(0) + n/10;
119 fitter.setRangeMask (start, end, True);
120
121 // Check masks
122
123 checkMasks (n, fitter, start(0), end(0));
124
125 // Now set range mask via abcissa values
126
127 Vector<Double> startF(1), endF(1);
128 startF(0) = x(start(0)); endF(0) = x(end(0));
129 fitter.setRangeMask (startF, endF, True);
130
131 // Check masks
132
133 checkMasks (n, fitter, start(0), end(0));
134 }
135 }
136
137 {
138 Vector<Double> x, y, amp(2), cen(2), sigma(2);
139 amp[0] = 1;
140 amp[1] = 1;
141 cen[0] = 10;
142 cen[1] = 60;
143 sigma[0] = 6;
144 sigma[1] = 6;
145 GaussianMultipletSpectralElement gm0 = makeMultiplet (x, y, amp, cen, sigma);
146 ProfileFit1D<Double> fitter;
147 Vector<Bool> m (x.size(), True);
148 fitter.setData (x,y,m);
149 fitter.addElement(gm0);
150 AlwaysAssert(
151 fitter.getList(False).nelements() == 1, AipsError
152 );
153 const SpectralElement *firstEl = fitter.getList(False)[0];
154 AlwaysAssert(
155 firstEl->getType() == SpectralElement::GMULTIPLET, AipsError
156 );
157 AlwaysAssert(fitter.fit(), AipsError);
158 // Check ok
159 AlwaysAssert(fitter.getDataMask().nelements() == x.size(), AipsError);
160 AlwaysAssert(allEQ(fitter.getDataMask(), True), AipsError);
161 AlwaysAssert(fitter.getRangeMask().nelements() == 0, AipsError);
162 AlwaysAssert(fitter.getTotalMask().nelements() == x.size(), AipsError);
163 AlwaysAssert(allEQ(fitter.getTotalMask(), True), AipsError);
164 AlwaysAssert(
165 *dynamic_cast<const GaussianMultipletSpectralElement*>(
166 fitter.getList(True)[0]
167 )
168 == *dynamic_cast<const GaussianMultipletSpectralElement*>(
169 fitter.getList(False)[0]
170 ),
171 AipsError
172 );
173 Matrix<Double> r(1, 3);
174 GaussianMultipletSpectralElement gm = gm0;
175
176 cout << "amplitude ratio constrained" << endl;
177 fitter = ProfileFit1D<Double>();
178 fitter.setData (x,y,m);
179 r = 0;
180 r(0, 0) = 0.9;
181 gm = GaussianMultipletSpectralElement(gm0.getGaussians(), r);
182 fitter.addElement(gm);
183 AlwaysAssert(fitter.fit(), AipsError);
184 cout << *dynamic_cast<const GaussianMultipletSpectralElement *>(
185 fitter.getList(True)[0]
186 )
187 << endl;
188 cout << "niter " << fitter.getNumberIterations() << endl;
189 cout << endl;
190
191 cout << "amplitude ratio constrained, sigma of reference fixed" << endl;
192 fitter = ProfileFit1D<Double>();
193 fitter.setData (x,y,m);
194 r = 0;
195 r(0, 0) = 0.9;
196 gm = GaussianMultipletSpectralElement(gm0.getGaussians(), r);
197 Vector<Bool> fixed(5, False);
198 fixed[2] = True;
199 gm.fix(fixed);
200 fitter.addElement(gm);
201 AlwaysAssert(fitter.fit(), AipsError);
202 cout << *dynamic_cast<const GaussianMultipletSpectralElement *>(
203 fitter.getList(True)[0]
204 )
205 << endl;
206 cout << "niter " << fitter.getNumberIterations() << endl;
207 cout << endl;
208
209 cout << "center offset constrained" << endl;
210 fitter = ProfileFit1D<Double>();
211 fitter.setData (x,y,m);
212 r = 0;
213 r(0, 1) = 48.5;
214 gm = GaussianMultipletSpectralElement(gm0.getGaussians(), r);
215 fitter.addElement(gm);
216 AlwaysAssert(fitter.fit(), AipsError);
217 cout << *dynamic_cast<const GaussianMultipletSpectralElement *>(
218 fitter.getList(True)[0]
219 )
220 << endl;
221 cout << "niter " << fitter.getNumberIterations() << endl;
222 cout << endl;
223
224 cout << "sigma ratio constrained" << endl;
225 fitter = ProfileFit1D<Double>();
226 fitter.setData (x,y,m);
227 r = 0;
228 r(0, 2) = 0.9;
229 gm = GaussianMultipletSpectralElement(gm0.getGaussians(), r);
230 fitter.addElement(gm);
231 AlwaysAssert(fitter.fit(), AipsError);
232 cout << *dynamic_cast<const GaussianMultipletSpectralElement *>(
233 fitter.getList(True)[0]
234 )
235 << endl;
236 cout << "niter " << fitter.getNumberIterations() << endl;
237 cout << endl;
238 }
239 {
240 Vector<Double> x, y, amp(2), cen(2), fwhm(2);
241 amp[0] = 1.5;
242 amp[1] = 4;
243 cen[0] = 10;
244 cen[1] = 60;
245 fwhm[0] = 6;
246 fwhm[1] = 6.5;
247 vector<LorentzianSpectralElement> lse = makeLorentzians (x, y, amp, cen, fwhm);
248
249 ProfileFit1D<Double> fitter;
250 Vector<Bool> m (x.size(), True);
251 fitter.setData (x,y,m);
252 for (uInt i=0; i<lse.size(); i++) {
253 // perturb the initial estimates
254 LorentzianSpectralElement z = lse[i];
255 z.setAmpl(amp[i] + (i*0.25));
256 z.setCenter(cen[i] + 3 + 0.5*i);
257 z.setFWHM(fwhm[i] + i*2);
258 fitter.addElement(z);
259 }
260 AlwaysAssert(
261 fitter.getList(False).nelements() == lse.size(), AipsError
262 );
263 const SpectralElement *firstEl = fitter.getList(False)[0];
264 AlwaysAssert(
265 firstEl->getType() == SpectralElement::LORENTZIAN, AipsError
266 );
267 AlwaysAssert(fitter.fit(), AipsError);
268 // Check ok
269 AlwaysAssert(fitter.getDataMask().nelements() == x.size(), AipsError);
270 AlwaysAssert(allEQ(fitter.getDataMask(), True), AipsError);
271 AlwaysAssert(fitter.getRangeMask().nelements() == 0, AipsError);
272 AlwaysAssert(fitter.getTotalMask().nelements() == x.size(), AipsError);
273 AlwaysAssert(allEQ(fitter.getTotalMask(), True), AipsError);
274 for (uInt i=0; i<lse.size(); i++) {
275 const LorentzianSpectralElement *got = dynamic_cast<
276 const LorentzianSpectralElement*
277 >(
278 fitter.getList(True)[i]
279 );
280 LorentzianSpectralElement exp = lse[i];
281 AlwaysAssert(nearAbs(*got, exp, 1e-15), AipsError);
282 }
283
284 cout << "niter " << fitter.getNumberIterations() << endl;
285 cout << endl;
286 }
287 cout << "*** Fit a power log polynomial" << endl;
288 {
289 ProfileFit1D<Double> fitter;
290 Vector<Double> x, y;
291 Vector<Double> estimates(2);
292 estimates[0] = 0.5;
293 estimates[1] = 2;
294 makePowerLogPoly(x, y, estimates);
295 Vector<Bool> mask(x.size(), True);
296 fitter.setData(x, y, mask);
297 SpectralList list;
298 estimates[1] = 1;
299 estimates[0] = 1;
300 list.add(PowerLogPolynomialSpectralElement(estimates));
301 fitter.setElements(list);
302 AlwaysAssert(fitter.fit(), AipsError);
303 Vector<Double> parms = fitter.getList(True)[0]->get();
304 cout << "parms " << parms << endl;
305 AlwaysAssert(near(parms[0], 0.5) && near(parms[1], 2.0), AipsError);
306 }
307 {
308 ProfileFit1D<Double> fitter;
309 Vector<Double> x, y;
310 Vector<Double> estimates(3);
311 estimates[0] = 0.5;
312 estimates[1] = 2;
313 estimates[2] = 0;
314 makePowerLogPoly(x, y, estimates);
315 Vector<Bool> mask(x.size(), True);
316 fitter.setData(x, y, mask);
317 SpectralList list;
318 estimates[0] = 0.55;
319 estimates[1] = 1.93;
320 estimates[2] = 0.4;
321 list.add(PowerLogPolynomialSpectralElement(estimates));
322 fitter.setElements(list);
323 AlwaysAssert(fitter.fit(), AipsError);
324 Vector<Double> parms = fitter.getList(True)[0]->get();
325 cout << "parms " << parms << endl;
326 AlwaysAssert(near(parms[0], 0.5) && near(parms[1], 2.0), AipsError);
327 }
328 {
329 ProfileFit1D<Double> fitter;
330 Vector<Double> x, y;
331 Vector<Bool> mask(x.size(), True);
332 Vector<Double> estimates(3);
333 estimates[0] = 0.5;
334 estimates[1] = 2;
335 estimates[2] = 1;
336 makePowerLogPoly(x, y, estimates);
337 fitter.setData(x, y, mask);
338 estimates[1] = 1.99;
339 estimates[2] = 0.999;
340 SpectralList list;
341 list.add(PowerLogPolynomialSpectralElement(estimates));
342 fitter.setElements(list);
343 AlwaysAssert(fitter.fit(), AipsError);
344 Vector<Double> parms = fitter.getList()[0]->get();
345 cout << "parms " << parms << endl;
346 AlwaysAssert(
347 near(parms[0], 0.5, 1e-4)
348 && near(parms[1], 2.0, 1e-4)
349 && nearAbs(parms[2], 1.0, 1e-4),
350 AipsError
351 );
352 }
353 {
354 cout << "*** Fit a polynomial" << endl;
355 ProfileFit1D<Double> fitter;
356 Vector<Double> x, y;
357 Vector<Bool> mask(x.size(), True);
358 Vector<Double> estimates(3);
359 estimates[0] = 0.5;
360 estimates[1] = 2;
361 estimates[2] = 1;
362 makePoly(x, y, estimates);
363 fitter.setData(x, y, mask);
364 SpectralList list;
365 list.add(PolynomialSpectralElement(estimates));
366 fitter.setElements(list);
367 AlwaysAssert(fitter.fit(), AipsError);
368 Vector<Double> parms = fitter.getList()[0]->get();
369 cout << "parms " << parms << endl;
370 cout << "niter " << fitter.getNumberIterations() << endl;
371 AlwaysAssert(allNear(parms, estimates, 1e-5), AipsError);
372 Vector<Double> bad(3);
373 bad[0] = 2;
374 bad[1] = 5;
375 bad[2] = 3;
376 list.clear();
377 list.add(PolynomialSpectralElement(bad));
378 fitter.clearList();
379 fitter.setElements(list);
380 AlwaysAssert(fitter.fit(), AipsError);
381 parms = fitter.getList()[0]->get();
382 cout << "niter " << fitter.getNumberIterations() << endl;
383 cout << "parms " << parms << endl;
384 AlwaysAssert(allNear(parms, estimates, 1e-5), AipsError);
385
386 }
387 cout << "OK" << endl;
388 return 0;
389 } catch (const AipsError& err) {
390 cerr << err.getMesg() << endl;
391 }
392 return 1;
393}
394
395void makeData (Vector<Double>& x, Vector<Double>& y, Vector<Bool>& m,
396 Double& amp, Double& cen, Double& sigma,
397 Double& p0, Double& p1)
398{
399 Int n = 256;
400 x.resize(n);
401 y.resize(n);
402 m.resize(n);
403 indgen(x);
404 x *= (2.3);
405 x += (1.0);
406 m = True;
407 //
408 amp = 10.0;
409 cen = x(n/2);
410 sigma = (x[n-1] - x[0]) / 50.0;
411 p0 = 0.15;
412 p1 = 1.2;
413 GaussianSpectralElement g(amp, cen, sigma);
414 cerr << "Gaussian: " << amp << ", " << cen << ", " << sigma << endl;
415 cerr << "Polynomial: " << p0 << ", " << p1 << endl;
416 //
417 Vector<Double> pars(2);
418 pars(0) = p0;
419 pars(1) = p1;
420 PolynomialSpectralElement p(pars);
421 for (uInt i=0; i<x.nelements(); i++) {
422 y(i) = g(x[i]) + p(x[i]);
423 }
424}
425
426void check (Double amp, Double cen, Double sig, Double p0, Double p1,
427 const SpectralList& list)
428{
429 Double tol(1e-4);
430 Vector<Double> p;
431 const SpectralElement *elG = list[0];
432 const SpectralElement *elP = list[1];
433
434 elG->get(p);
435 cout << "p " << p << " amp " << amp << " tol " << tol << endl;
436 AlwaysAssert(near(amp, p[0], tol), AipsError);
437 AlwaysAssert(near(cen, p[1], tol), AipsError);
438 AlwaysAssert(near(sig, p[2], tol), AipsError);
439 p.resize(0);
440 elP->get(p);
441 AlwaysAssert(near(p0, p[0], tol), AipsError);
442 AlwaysAssert(near(p1, p[1], tol), AipsError);
443}
444
445
446void checkMasks (uInt n, const ProfileFit1D<Double>& fitter, Int start,
447 Int end)
448{
449 Vector<Bool> rangeMask = fitter.getRangeMask();
450 Vector<Bool> totalMask = fitter.getTotalMask();
451 //
452 AlwaysAssert(rangeMask.nelements()==n, AipsError);
453 AlwaysAssert(totalMask.nelements()==n, AipsError);
454 AlwaysAssert(allEQ(rangeMask, totalMask), AipsError);
455 //
456 IPosition iStart(1), iEnd(1);
457 {
458 iStart(0) = 0;
459 iEnd(0) = start-1;
460 Vector<Bool> tmp = rangeMask(iStart, iEnd);
461 AlwaysAssert(allEQ(tmp, False), AipsError);
462 }
463 {
464 iStart(0) = start;
465 iEnd(0) = end;
466 Vector<Bool> tmp = rangeMask(iStart, iEnd);
467 AlwaysAssert(allEQ(tmp, True), AipsError);
468 }
469 {
470 iStart(0) = end+1;
471 iEnd(0) = n-1;
472 Vector<Bool> tmp = rangeMask(iStart, iEnd);
473 AlwaysAssert(allEQ(tmp, False), AipsError);
474 }
475}
476
477GaussianMultipletSpectralElement makeMultiplet (
478 Vector<Double>& x, Vector<Double>& y,
479 const Vector<Double>& amp, const Vector<Double>& cen,
480 const Vector<Double>& sigma
481) {
482 Double minx = cen[0] - 5*sigma[0];
483 Double maxx = cen[0] + 5*sigma[0];
484 for (uInt i=1; i<amp.size(); i++) {
485 minx = min(minx, cen[i] - 5*sigma[i]);
486 maxx = max(maxx, cen[i] + 5*sigma[i]);
487 }
488 minx = (int)minx;
489 maxx = (int)maxx + 1;
490 x.resize((Int)(maxx-minx+1));
491 indgen(x);
492 x += minx;
493 y.resize(x.size());
494 vector<GaussianSpectralElement> g(amp.size());
495 Matrix<Double> r(g.size() - 1, 3, 0);
496 for (uInt i=0; i<amp.size(); i++) {
497 g[i] = GaussianSpectralElement(amp[i], cen[i], sigma[i]);
498 if(i > 0) {
499 r(i-1, 0) = amp[i]/amp[0];
500 }
501 }
502 GaussianMultipletSpectralElement gm(g, r);
503 for (uInt i=0; i<x.size(); i++) {
504 y[i] = gm(x[i]);
505 }
506 return gm;
507}
508
509vector<LorentzianSpectralElement> makeLorentzians (
510 Vector<Double>& x, Vector<Double>& y,
511 const Vector<Double>& amp, const Vector<Double>& cen,
512 const Vector<Double>& fwhm
513) {
514 Double minx = cen[0] - 5*fwhm[0];
515 Double maxx = cen[0] + 5*fwhm[0];
516 for (uInt i=1; i<amp.size(); i++) {
517 minx = min(minx, cen[i] - 5*fwhm[i]);
518 maxx = max(maxx, cen[i] + 5*fwhm[i]);
519 }
520 minx = (int)minx;
521 maxx = (int)maxx + 1;
522 x.resize((Int)(maxx-minx+1));
523 indgen(x);
524 x += minx;
525 y.resize(x.size());
526 vector<LorentzianSpectralElement> lse;
527 for (uInt i=0; i<amp.size(); i++) {
528 lse.push_back(LorentzianSpectralElement(amp[i], cen[i], fwhm[i]));
529 }
530 for (uInt i=0; i<x.size(); i++) {
531 y[i] = 0;
532 for (uInt j=0; j<lse.size(); j++) {
533 y[i] += lse[j](x[i]);
534 }
535 }
536 return lse;
537}
538
539
540PowerLogPolynomialSpectralElement makePowerLogPoly(
541 Vector<Double>&x, Vector<Double>& y, const Vector<Double>& coeffs
542) {
543 x.resize(10);
544 x[0] = 1;
545 x[1] = 3;
546 x[2] = 5;
547 x[3] = 6;
548 x[4] = 8;
549 x[5] = 10;
550 x[6] = 12;
551 x[7] = 15;
552 x[8] = 20;
553 x[9] = 25;
554 y.resize(x.size());
555 PowerLogPolynomialSpectralElement plp(coeffs);
556 for (uInt i=0; i<x.size(); i++) {
557 y[i] = plp(x[i]);
558 }
559 return plp;
560
561}
562
563PolynomialSpectralElement makePoly(
564 Vector<Double>&x, Vector<Double>& y, const Vector<Double>& coeffs
565) {
566 x.resize(10);
567 x[0] = 1;
568 x[1] = 3;
569 x[2] = 5;
570 x[3] = 6;
571 x[4] = 8;
572 x[5] = 10;
573 x[6] = 12;
574 x[7] = 15;
575 x[8] = 20;
576 x[9] = 25;
577 y.resize(x.size());
578 PolynomialSpectralElement poly(coeffs);
579 for (uInt i=0; i<x.size(); i++) {
580 y[i] = poly(x[i]);
581 }
582 return poly;
583
584}
585
586
Note: See TracBrowser for help on using the repository browser.