source: trunk/external-alma/components/SpectralComponents/test/tSpectralFit.cc@ 3006

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