source: trunk/external-alma/components/SpectralComponents/test/tGaussianMultipletSpectralElement.cc@ 3058

Last change on this file since 3058 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.2 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
27#include <casa/aips.h>
28#include <casa/Arrays/ArrayMath.h>
29#include <casa/Containers/Record.h>
30#include <casa/Utilities/PtrHolder.h>
31#include <components/SpectralComponents/PolynomialSpectralElement.h>
32
33#include <components/SpectralComponents/GaussianSpectralElement.h>
34#include <components/SpectralComponents/GaussianMultipletSpectralElement.h>
35#include <components/SpectralComponents/SpectralElementFactory.h>
36
37#include <casa/Utilities/Assert.h>
38#include <casa/Arrays/ArrayIO.h>
39
40#include <casa/Arrays/Vector.h>
41
42#include <casa/iostream.h>
43
44#include <casa/namespace.h>
45
46#include <vector>
47
48int main() {
49
50 {
51 cout << "Test that exception is thrown if there is "
52 << "a mismatch between number of Gaussians and "
53 << "number of relationships between them" << endl;
54 GaussianSpectralElement g1(4,4,4);
55 GaussianSpectralElement g2(5,5,5);
56 vector<GaussianSpectralElement> doublet(2);
57 doublet[0] = g1;
58 doublet[1] = g2;
59 Matrix<Double> relations(2, 3, 0);
60 relations(0,0) = 4;
61 Bool result = True;
62 try {
63 GaussianMultipletSpectralElement(doublet, relations);
64 result = False;
65 }
66 catch (AipsError) {}
67 AlwaysAssert(result, AipsError);
68 }
69 {
70 cout << "Test that exception is thrown if relations "
71 << "does not have three columns" << endl;
72 GaussianSpectralElement g1(4, 4, 4);
73 GaussianSpectralElement g2(5, 5, 5);
74 vector<GaussianSpectralElement> doublet(2);
75 doublet[0] = g1;
76 doublet[1] = g2;
77 Matrix<Double> relations(1, 4, 0);
78 relations(0,0) = 4;
79 Bool result = True;
80 try {
81 GaussianMultipletSpectralElement(doublet, relations);
82 result = False;
83 }
84 catch (AipsError) {}
85 AlwaysAssert(result, AipsError);
86 }
87 {
88 cout << "Test that exception is thrown if there is "
89 << "an amplitude fixing incompatibility" << endl;
90 GaussianSpectralElement g1(4, 4, 4);
91 GaussianSpectralElement g2(5, 5, 5);
92 g2.fixAmpl();
93 cout << "g1 " << g1 << endl;
94 cout << "g2 " << g2 << endl;
95 std::vector<GaussianSpectralElement> doublet(2);
96 doublet[0] = g1;
97 cout << "doub 0 " << doublet[0] << " " << &doublet[0] << endl;
98 cout << "doub 1 " << doublet[1] << " " << &doublet[1] << endl;
99
100 doublet[1] = g2;
101 cout << "doub 0 " << doublet[0] << " " << &doublet[0] << endl;
102 cout << "doub 1 " << doublet[1] << " " << &doublet[1] << endl;
103 cout << "g1 " << g1 << endl;
104 cout << "g2 " << g2 << endl;
105 Matrix<Double> relations(1, 3, 0);
106 relations(0,0) = 4;
107 Bool result = True;
108 try {
109 GaussianMultipletSpectralElement(doublet, relations);
110 result = False;
111 }
112 catch (AipsError x) {}
113 AlwaysAssert(result, AipsError);
114 }
115 {
116 cout << "Test that exception is thrown if there is "
117 << "a center fixing incompatibility" << endl;
118 GaussianSpectralElement g1(4, 4, 4);
119 GaussianSpectralElement g2(5, 5, 5);
120 g2.fixCenter();
121 vector<GaussianSpectralElement> doublet(2);
122 doublet[0] = g1;
123 doublet[1] = g2;
124 Matrix<Double> relations(1, 3, 0);
125 relations(0,1) = 4;
126 Bool result = True;
127 try {
128 GaussianMultipletSpectralElement(doublet, relations);
129 result = False;
130 }
131 catch (AipsError x) {}
132 AlwaysAssert(result, AipsError);
133 }
134 {
135 cout << "Test that exception is thrown if there is "
136 << "a width fixing incompatibility" << endl;
137 GaussianSpectralElement g1(4, 4, 4);
138 GaussianSpectralElement g2(5, 5, 5);
139 g2.fixFWHM();
140 vector<GaussianSpectralElement> doublet(2);
141 doublet[0] = g1;
142 doublet[1] = g2;
143 Matrix<Double> relations(1, 3, 0);
144 relations(0, 2) = 4;
145 Bool result = True;
146 try {
147 GaussianMultipletSpectralElement(doublet, relations);
148 result = False;
149 }
150 catch (AipsError x) {}
151 AlwaysAssert(result, AipsError);
152 }
153 {
154 cout << "Test gaussians were correctly constructed " << endl;
155 GaussianSpectralElement g1(4, 4.4, 4.6);
156 GaussianSpectralElement g2(5, 5.2, 5.8);
157 vector<GaussianSpectralElement> gaussians(2);
158 gaussians[0] = g1;
159 gaussians[1] = g2;
160 Matrix<Double> relations(1, 3, 0);
161 Double ampRatio = 4.5;
162 relations(0, 0) = ampRatio;
163 GaussianMultipletSpectralElement doublet(gaussians, relations);
164 for (uInt i=0; i<gaussians.size(); i++) {
165 Vector<Double> expected = gaussians[i].get();
166 if (i==1) {
167 expected[0] = ampRatio*gaussians[0].getAmpl();
168 }
169 Vector<Double> got = doublet.getGaussians()[i].get();
170 cout << "*** got " << got << endl;
171 cout << "*** exp " << expected << endl;
172 AlwaysAssert(allTrue(got == expected), AipsError);
173 }
174
175 relations = 0;
176 Double centerOff = -20;
177 relations(0, 1) = centerOff;
178 doublet = GaussianMultipletSpectralElement(gaussians, relations);
179 for (uInt i=0; i<gaussians.size(); i++) {
180 Vector<Double> expected = gaussians[i].get();
181 if (i==1) {
182 expected[1] = centerOff + gaussians[0].getCenter();
183 }
184 Vector<Double> got = doublet.getGaussians()[i].get();
185 AlwaysAssert(allTrue(got == expected), AipsError);
186 }
187 relations = 0;
188 Vector<Double> x(5, 2);
189 Double sigmaRatio = 0.5;
190 relations(0, 2) = sigmaRatio;
191 doublet = GaussianMultipletSpectralElement(gaussians, relations);
192 for (uInt i=0; i<gaussians.size(); i++) {
193 Vector<Double> expected = gaussians[i].get();
194 if (i==1) {
195 expected[2] = sigmaRatio * gaussians[0].getSigma();
196 }
197 Vector<Double> got = doublet.getGaussians()[i].get();
198 AlwaysAssert(allTrue(got == expected), AipsError);
199 }
200 }
201 {
202 cout << "Test toRecord()/fromRecord()" << endl;
203 GaussianSpectralElement g1(4.6, 4.5, 4.4);
204 GaussianSpectralElement g2(5.6, 5.2, 5.5);
205 vector<GaussianSpectralElement> gaussians(2);
206 gaussians[0] = g1;
207 gaussians[1] = g2;
208 Matrix<Double> relations(1, 3, 0);
209 relations(0, 2) = 4;
210 GaussianMultipletSpectralElement doublet(gaussians, relations);
211 Record myRec;
212 cout << "doublet real " << doublet << endl;
213 doublet.toRecord(myRec);
214 cout << "myrec " << myRec << endl;
215 cout << __FILE__ << " " << __LINE__ << endl;
216 PtrHolder<SpectralElement> ptr(SpectralElementFactory::fromRecord(myRec));
217 cout << __FILE__ << " " << __LINE__ << endl;
218 GaussianMultipletSpectralElement out = *dynamic_cast<GaussianMultipletSpectralElement*>(
219 ptr.ptr()
220 );
221 cout << __FILE__ << " " << __LINE__ << endl;
222 cout << "out " << out << endl;
223 cout << "doublet " << doublet << endl;
224 AlwaysAssert(out == doublet, AipsError);
225 }
226 {
227 cout << "Test setting/getting" << endl;
228 GaussianSpectralElement g1(4.6, 4.5, 4.4);
229 GaussianSpectralElement g2(5.6, 5.2, 5.5);
230 GaussianSpectralElement g3(6.6, 6.2, 6.5);
231
232 Vector<Double> errs(3, 0);
233 errs[0] = 0.1;
234 errs[1] = 0.2;
235 errs[2] = 0.3;
236 g1.setError(errs);
237 Vector<Bool> fixed(3, False);
238 fixed[1] = True;
239 g1.fix(fixed);
240
241
242 vector<GaussianSpectralElement> gaussians(3);
243 gaussians[0] = g1;
244 gaussians[1] = g2;
245 gaussians[2] = g3;
246
247 Matrix<Double> relations(2, 3, 0);
248 relations(0, 2) = 4;
249 relations(1, 1) = 10;
250 GaussianMultipletSpectralElement triplet(gaussians, relations);
251
252 Vector<Double> g(7, 0);
253 uInt j = 0;
254 for (uInt i=0; i<9; i++) {
255 if (i != 5 && i != 7) {
256 g[j] = gaussians[i/3].get()[i%3];
257 j++;
258 }
259 }
260 AlwaysAssert(allNear(triplet.get(), g, 1e-8), AipsError);
261 AlwaysAssert(
262 allNear(
263 gaussians[0].get(), triplet.getGaussians()[0].get(), 1e-8
264 ), AipsError
265 );
266 AlwaysAssert(
267 allNear(
268 gaussians[0].getError(), triplet.getGaussians()[0].getError(), 1e-8
269 ), AipsError
270 );
271 AlwaysAssert(
272 allTrue(
273 gaussians[0].fixed() == triplet.getGaussians()[0].fixed()
274 ), AipsError
275 );
276 Vector<Double> z = gaussians[1].get();
277 z[2] = relations(0, 2) * gaussians[0].getSigma();
278 Vector<Double> err = gaussians[1].getError();
279 err[2] = relations(0, 2) * gaussians[0].getSigmaErr();
280 AlwaysAssert(
281 allNear(
282 z, triplet.getGaussians()[1].get(), 1e-8
283 ), AipsError
284 );
285 AlwaysAssert(
286 allNear(
287 err, triplet.getGaussians()[1].getError(), 1e-8
288 ), AipsError
289 );
290 AlwaysAssert(
291 allTrue(
292 gaussians[1].fixed() == triplet.getGaussians()[1].fixed()
293 ), AipsError
294 );
295 z = gaussians[2].get();
296 z[1] = relations(1, 1) + gaussians[0].getCenter();
297 err = gaussians[2].getError();
298 err[1] = gaussians[0].getCenterErr();
299 AlwaysAssert(
300 allNear(
301 z, triplet.getGaussians()[2].get(), 1e-8
302 ), AipsError
303 );
304 AlwaysAssert(
305 allNear(
306 err, triplet.getGaussians()[2].getError(), 1e-8
307 ), AipsError
308 );
309 cout << gaussians[2].fixed() << endl;
310 cout << triplet.getGaussians()[2].fixed() << endl;
311 cout << fixed << endl;
312 AlwaysAssert(
313 allTrue(
314 triplet.getGaussians()[2].fixed() == fixed
315 ), AipsError
316 );
317 cout << triplet.getFunction() << endl;
318 Vector<Double> parms(7);
319 for (uInt i=0; i < parms.size(); i++) {
320 parms[i] = 10+i;
321 }
322 triplet.set(parms);
323
324 AlwaysAssert(
325 allNear(
326 triplet.get(), parms, 1e-5
327 ), AipsError
328 );
329 Vector<Double> exp(3);
330 exp[0] = parms[0];
331 exp[1] = parms[1];
332 exp[2] = parms[2];
333 AlwaysAssert(
334 allNear(
335 triplet.getGaussians()[0].get(), exp, 1e-5
336 ), AipsError
337 );
338 exp[0] = parms[3];
339 exp[1] = parms[4];
340 exp[2] = relations(0, 2) * parms[2];
341 AlwaysAssert(
342 allNear(
343 triplet.getGaussians()[1].get(), exp, 1e-5
344 ), AipsError
345 );
346 exp[0] = parms[5];
347 exp[1] = relations(1, 1) + parms[1];
348 exp[2] = parms[6];
349 AlwaysAssert(
350 allNear(
351 triplet.getGaussians()[2].get(), exp, 1e-5
352 ), AipsError
353 );
354
355
356
357
358
359 }
360 cout << "ok" << endl;
361 return 0;
362
363
364}
Note: See TracBrowser for help on using the repository browser.