source: trunk/external-alma/components/SpectralComponents/test/tSpectralFit.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: 11.9 KB
Line 
1//# tSpectralFit.cc -- test SpectralElement; SpectralEstimate and SpectralFit
2//# Copyright (C) 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: tSpectralFit.cc 21451 2014-06-10 07:48:08Z gervandiepen $
27
28//# Includes
29#include <casa/aips.h>
30#include <casa/Arrays/MaskArrLogi.h>
31#include <casa/Arrays/ArrayMath.h>
32#include <casa/Arrays/ArrayIO.h>
33#include <casa/Arrays/Vector.h>
34#include <casa/BasicMath/Random.h>
35#include <casa/Containers/Record.h>
36#include <casa/Exceptions/Error.h>
37#include <casa/Utilities/Assert.h>
38#include <components/SpectralComponents/SpectralEstimate.h>
39#include <components/SpectralComponents/SpectralFit.h>
40#include <components/SpectralComponents/SpectralList.h>
41#include <casa/iostream.h>
42
43#include <casa/namespace.h>
44
45int main() {
46        // RecordInterface
47        {
48                cout << "Test RecordInterface" << endl;
49                cout << "---------------------------------------------------" << endl;
50                SpectralList list;
51                GaussianSpectralElement gEl(1.0, 10.0, 2.0);
52                PolynomialSpectralElement pEl(1);
53                list.add (gEl);
54                list.add(pEl);
55
56                cout << "toRecord" << endl;
57                Record rec;
58                list.toRecord(rec);
59                cout << "fromRecord" << endl;
60                SpectralList list2;
61                String errMsg;
62                if (!list2.fromRecord(errMsg, rec)) {
63                        throw(AipsError(errMsg));
64                }
65                Double tol(1.0e-6);
66                AlwaysAssert(list.nelements()==list.nelements(),AipsError);
67                for (uInt i=0; i<list.nelements(); i++) {
68                        AlwaysAssert (list[i]->getType()==list2[i]->getType(), AipsError);
69                        Vector<Double> p1,p2;
70                        list[i]->get(p1);
71                        list2[i]->get(p2);
72                        AlwaysAssert(allNear(p1,p2,tol),AipsError);
73                }
74                cout << endl;
75        }
76
77        // get estimates
78        {
79                cout << "Test SpectralEstimate" << endl;
80                cout << "---------------------------------------------------" << endl;
81                const Int NPAR = 9;
82                const Int NPTS = 100;
83                uInt q = 2;
84                Int r;
85                Int n = NPTS;
86                Int np = NPAR;
87                Double rms = 0.5;
88                Double cutoff= 1.0;
89                Double minsig = 0.5;
90                Double par[NPAR];
91                Vector<Float> y(NPTS);
92                Vector<Float> xy(NPTS);
93
94                par[0] = 5.0; par[1] = 45.0; par[2] = 1.1;
95                par[3] = 4.0; par[4] = 50.0; par[5] = 1.1;
96                par[6] = 6.0; par[7] = 55.0; par[8] = 1.1;
97                for (Int i = 0; i<n; i++) {
98                        y(i) = 0.0;
99                        xy(i) = i;
100                        for (Int j = 0; j < np; j += 3) {
101                                if (par[j] != 0.0) {
102                                        Double a = par[j];
103                                        Double c = (par[j+1] - Double(i));
104                                        Double s = par[j+2];
105                                        y(i) += a * exp( -0.5 * c / s * c / s );
106                                };
107                        };
108                };
109                cout << "Made spectrum with components: " << endl;
110                for (uInt i=0; i<3; i++) {
111                        cout << GaussianSpectralElement(
112                                        par[3*i+0], par[3*i+1], par[3*i+2]) << endl;
113                };
114                {
115                        SpectralEstimate est(rms, cutoff, minsig);
116                        cout << "*** n " << n << endl;
117                        Vector<Float> deriv(n);
118                        const SpectralList &slist = est.estimate(y, &deriv);
119                        cout << "*** y " << y << endl;
120                        r = slist.nelements();
121                        cout << "Found (using rms, cutoff, minsig, q): (" <<
122                                        rms << ", " << cutoff << ", " << minsig << ", " << q << ") " <<
123                                        r << " estimates" << endl;
124                        for (Int i=0; i<r; i++) {
125                                cout << "Estimate " << i << ": " << *slist[i] << endl;
126                        };
127                        Vector<Float> res(NPTS);
128                        res = y;
129                        cout << "*** res " << res << endl;
130                        slist.residual(res);
131                        cout << "Minimum, Maximum residuals: " << min(res) << ", " <<
132                                        max(res) << endl;
133                        for (Int j=0; j<r; j++) {
134                                for (Int i = 0; i<n; i++) {
135                                        res(i) = slist[j]->operator()(i);
136                                }
137                        }
138                }
139                {
140                        SpectralEstimate est(rms, cutoff, minsig);
141                        est.setWindowing(True);
142                        const SpectralList &slist = est.estimate(xy, y);
143                        r = slist.nelements();
144                        cout << "Window (using rms, cutoff, minsig, q): (" <<
145                                        rms << ", " << cutoff << ", " << minsig << ", " << q << ") " <<
146                                        r << " estimates" << endl;
147                        for (Int i=0; i<r; i++) {
148                                cout << "Estimate " << i << ": " << *slist[i] << endl;
149                        };
150                }
151                {
152                        q = 1;
153                        SpectralEstimate est(rms, cutoff, minsig);
154                        est.setQ(q);
155                        const SpectralList &slist = est.estimate(xy, y);
156                        r = slist.nelements();
157                        cout << "Found (using rms, cutoff, minsig, q): (" <<
158                                        rms << ", " << cutoff << ", " << minsig << ", " << q << ") " <<
159                                        r << " estimates" << endl;
160                        for (Int i=0; i<r; i++) {
161                                cout << "Estimate " << i << ": " << *slist[i] << endl;
162                        };
163                }
164                {
165                        q = 1;
166                        rms = cutoff = minsig = 0.0;
167                        SpectralEstimate est(rms, cutoff, minsig);
168                        est.setQ(q);
169                        const SpectralList &slist = est.estimate(xy, y);
170                        r = slist.nelements();
171                        cout << "Found (using rms, cutoff, minsig, q): (" <<
172                                        rms << ", " << cutoff << ", " << minsig << ", " << q << ") " <<
173                                        r << " estimates" << endl;
174                        for (Int i=0; i<r; i++) {
175                                cout << "Estimate " << i << ": " << *slist[i] << endl;
176                        };
177                }
178                {
179                        q = 2;
180                        rms = cutoff = minsig = 0.0;
181                        SpectralEstimate est(rms, cutoff, minsig);
182                        est.setQ(q);
183                        const SpectralList &slist = est.estimate(xy, y);
184                        r = slist.nelements();
185                        cout << "Found (using rms, cutoff, minsig, q): (" <<
186                                        rms << ", " << cutoff << ", " << minsig << ", " << q << ") " <<
187                                        r << " estimates" << endl;
188                        for (Int i=0; i<r; i++) {
189                                cout << "Estimate " << i << ": " << *slist[i] << endl;
190                        };
191                }
192        }
193
194        // test fitter
195        {
196                // Number of components
197                const uInt ncomp = 4;
198                // Component data
199                Double ampl[ncomp] = { 1, 2.5, 0.5, 5 };
200                Double center[ncomp] = { 0.25, 0.5, 0.75, 0.3 };
201                Double sigma[ncomp] = { 2, 4, 8, 3 };
202                MLCG genit;
203                Normal noise(&genit, 0.0, 0.1);
204                Vector<Double> freq(1024);
205                Vector<Float> ffreq(1024);
206                for (uInt i=0; i<1024; i++) {
207                        freq(i) = 1400 + i*10.0/1024.0;
208                        ffreq(i) = freq(i);
209                };
210
211                try {
212                        cout << "Test SpectralFit" << endl;
213                        cout << "---------------------------------------------------" << endl;
214
215                        GaussianSpectralElement el[ncomp] = {
216                                        GaussianSpectralElement(ampl[0],
217                                                        center[0]*10.+1400., sigma[0]*10./1024.),
218                                                        GaussianSpectralElement(ampl[1],
219                                                                        center[1]*10.+1400., sigma[1]*10./1024.),
220                                                                        GaussianSpectralElement(ampl[2],
221                                                                                        center[2]*10.+1400., sigma[2]*10./1024.),
222                                                                                        GaussianSpectralElement(ampl[3],
223                                                                                                        center[3]*10.+1400., sigma[3]*10./1024.) };
224
225                        cout << "Spectral elements: " << endl;
226                        for (uInt j=0; j<ncomp; j++) cout << el[j] << endl;
227                        cout << "---------------------------------------------------" << endl;
228
229                        Vector<Double> dat(1024);
230                        Vector<Float> fdat(1024);
231                        for (uInt i=0; i<1024; i++) {
232                                dat(i) = 0;
233                                for (uInt j=0; j<ncomp; j++) {
234                                        dat(i) += el[j].getAmpl()*exp(-(freq(i)-el[j].getCenter())*
235                                                        (freq(i)-el[j].getCenter())*log(16.0)/
236                                                        el[j].getSigma()/
237                                                        el[j].getSigma());
238                                };
239                                dat(i) += noise();
240                                fdat(i) = dat(i);
241                        };
242                        cout << "Data for frequencies " << freq(300) << " - " << freq(315) <<
243                                        endl;
244                        for (uInt i=300; i<316; i++) cout << freq(i) << ": " << dat(i) << endl;
245                        cout << "---------------------------------------------------" << endl;
246
247                        cout << "Specify fitter on the 4 gaussians:" << endl;
248                        SpectralFit fitter;
249                        for (uInt i=0; i<ncomp; i++) fitter.addFitElement(el[i]);
250                        for (uInt i=0; i<fitter.list().nelements(); i++) {
251                                cout << *fitter.list()[i] << endl;
252                        };
253
254                        cout << "---------------------------------------------------" << endl;
255                        cout << "Execute the fitter, can do (1:Y): " <<
256                                        fitter.fit(dat,freq) << endl;
257                        cout << "The results: " << endl;
258                        Vector<Double> tmp;
259                        for (uInt i=0; i<fitter.list().nelements(); i++) {
260                                cout << *fitter.list()[i];
261                                fitter.list()[i]->get(tmp);
262                                cout << "Parameters: " << tmp << endl;
263                                fitter.list()[i]->getError(tmp);
264                                cout << "Errors:     " << tmp << endl << endl;
265                        };
266                        cout << ">>>" << endl;
267                        cout << "Number of iterations: " << fitter.nIterations() << endl;
268                        cout << "<<<" << endl;
269                        cout << "---------------------------------------------------" << endl;
270
271                        {
272                                el[3].fixAmpl();
273                                SpectralFit fitter;
274                                for (uInt i=0; i<ncomp; i++) fitter.addFitElement(el[i]);
275                                for (uInt i=0; i<fitter.list().nelements(); i++) {
276                                        cout << *fitter.list()[i] << endl;
277                                };
278
279                                cout << "---------------------------------------------------" << endl;
280                                cout << "Execute the fitter again with a fixed ampl, can do (1:Y): " <<
281                                                fitter.fit(dat,freq) << endl;
282                                cout << "The results: " << endl;
283                                Vector<Double> tmp;
284                                for (uInt i=0; i<fitter.list().nelements(); i++) {
285                                        cout << *fitter.list()[i];
286                                        fitter.list()[i]->get(tmp);
287                                        cout << "Parameters: " << tmp << endl;
288                                        fitter.list()[i]->getError(tmp);
289                                        cout << "Errors:     " << tmp << endl << endl;
290                                };
291                                cout << ">>>" << endl;
292                                cout << "Number of iterations: " << fitter.nIterations() << endl;
293                                cout << "<<<" << endl;
294                                el[3].fixAmpl(False);
295                        }
296                        cout << "---------------------------------------------------" << endl;
297
298                        cout << "Different start values: " << endl;
299                        for (uInt i=0; i<ncomp; i++) {
300                                el[i].setAmpl(0.9*el[i].getAmpl());
301                                el[i].setSigma(0.95*el[i].getSigma());
302                                fitter.setFitElement(i, el[i]);
303                        };
304                        for (uInt i=0; i<fitter.list().nelements(); i++) {
305                                cout << *fitter.list()[i] << endl;
306                        };
307                        cout << "Execute the fitter, can do (1:Y): " <<
308                                        fitter.fit(dat, freq) << endl;
309                        cout << "The results: " << endl;
310                        for (uInt i=0; i<fitter.list().nelements(); i++) {
311                                cout << *fitter.list()[i] << endl;
312                        };
313                        cout << ">>>" << endl;
314                        cout << "Number of iterations: " << fitter.nIterations() << endl;
315                        cout << "<<<" << endl;
316
317                        cout << "---------------------------------------------------" << endl;
318                        cout << "Differences: " << endl;
319                        Vector<Double> xdat(1024);
320                        Vector<Float> fxdat(1024);
321                        Double mx(-1e6);
322                        Double mn(1e6);
323                        Double avg(0);
324                        Double sg(0);
325                        xdat = dat;
326                        fitter.list().residual(xdat, freq);
327                        convertArray(fxdat, xdat);
328                        mx = max(xdat);
329                        mn = min(xdat);
330                        avg = sum(xdat);
331                        sg = sum(xdat*xdat);
332                        avg /= 1024.;
333                        sg = sqrt(sg/1024./1023.);
334                        cout << "Min difference: " << mn <<
335                                        ", max: " << mx << ", average: " << avg <<
336                                        ", sigma: " << sg << endl;
337                        cout << "---------------------------------------------------" << endl;
338
339                        cout << "---------------------------------------------------" << endl;
340                        cout << "Estimates: " << endl;
341                        SpectralEstimate mest(0.1, 0.5);
342                        mest.setQ(5);
343                        const SpectralList &slist = mest.estimate(fdat);
344                        Int mr = slist.nelements();
345                        cout << "Found " << mr << " estimates" << endl;
346                        for (Int i=0; i<mr; i++) {
347                                cout << "Estimate " << i << ": " << *slist[i] << endl;
348                        };
349                        fxdat= fdat;
350                        slist.residual(fxdat);
351                        for (Int j=0; j<mr; j++) {
352                                slist.evaluate(fxdat);
353                        };
354
355                        cout << "---------------------------------------------------" << endl;
356                        cout << "Limited number of estimates: " << endl;
357                        {
358                                SpectralEstimate mest(0.1, 0.5, 0.0, 4);
359                                mest.setQ(5);
360                                const SpectralList &slist = mest.estimate(fdat);
361                                Int mr = slist.nelements();
362                                cout << "Found " << mr << " estimates" << endl;
363                                for (Int i=0; i<mr; i++) {
364                                        cout << "Estimate " << i << ": " << *slist[i] << endl;
365                                };
366                                fxdat= fdat;
367                                slist.residual(fxdat);
368                                for (Int j=0; j<mr; j++) {
369                                        slist.evaluate(fxdat);
370                                };
371                        }
372                        cout << "---------------------------------------------------" << endl;
373
374                } catch (AipsError x) {
375                        cout << x.getMesg() << endl;
376                        return 1;
377                };
378        }
379
380        return 0;
381}
Note: See TracBrowser for help on using the repository browser.