source: trunk/external-alma/components/SpectralComponents/test/tProfileFit1D.cc @ 2980

Last change on this file since 2980 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.