source: trunk/external-alma/components/SpectralComponents/GaussianMultipletSpectralElement.cc@ 3150

Last change on this file since 3150 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.6 KB
Line 
1//# SpectralElement.cc: Describes (a set of related) spectral lines
2//# Copyright (C) 2001,2004
3//# Associated Universities, Inc. Washington DC, USA.
4//#
5//# This library is free software; you can redistribute it and/or modify it
6//# under the terms of the GNU Library General Public License as published by
7//# the Free Software Foundation; either version 2 of the License, or (at your
8//# option) any later version.
9//#
10//# This library 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 Library General Public
13//# License for more details.
14//#
15//# You should have received a copy of the GNU Library General Public License
16//# along with this library; if not, write to the Free Software Foundation,
17//# Inc., 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 <components/SpectralComponents/GaussianMultipletSpectralElement.h>
28
29#include <casa/Arrays/ArrayIO.h>
30#include <casa/Arrays/ArrayLogical.h>
31#include <casa/Arrays/ArrayMath.h>
32#include <casa/Containers/Record.h>
33#include <components/SpectralComponents/GaussianSpectralElement.h>
34
35#include <casa/iostream.h>
36
37namespace casa { //# NAMESPACE CASA - BEGIN
38
39#define _ORIGIN String("GaussianMultipletSpectralElement::") + __FUNCTION__ + ":" + String::toString(__LINE__) + ": "
40
41GaussianMultipletSpectralElement::GaussianMultipletSpectralElement(
42 const vector<GaussianSpectralElement>& estimates,
43 const Matrix<Double>& constraints
44) : CompiledSpectralElement(SpectralElement::GMULTIPLET),
45 _gaussians(estimates),_constraints(constraints),
46 _paramIndices(estimates.size(), 3, 0) {
47 if(estimates.size() != constraints.nrow()+1) {
48 throw AipsError(
49 _ORIGIN
50 + "Mismatch between size of estimates and constraints"
51 );
52 }
53 if (constraints.ncolumn() != 3) {
54 throw AipsError(_ORIGIN + "constraints does not have 3 columns");
55 }
56 Matrix<Bool> fixed(LogicalArray(constraints != 0.0));
57 for (uInt i=1; i<estimates.size(); i++) {
58 if (! estimates[0].fixedAmpl()
59 && estimates[i].fixedAmpl()
60 && fixed(i-1, 0)
61 ) {
62 throw AipsError(_ORIGIN + "You cannot fix the amplitude of a "
63 "non-reference Gaussian if the reference Gaussian's "
64 "amplitude is not fixed and there is a relationship "
65 "between the two amplitudes."
66 );
67 }
68 if (! estimates[0].fixedCenter()
69 && estimates[i].fixedCenter()
70 && fixed(i-1, 1)
71 ) {
72 throw AipsError(_ORIGIN + "You cannot fix the center of a non-reference "
73 "Gaussian if the reference Gaussian's center is not fixed and there "
74 "is a relationship between the two centers."
75 );
76 }
77 if (! estimates[0].fixedWidth()
78 && estimates[i].fixedWidth()
79 && fixed(i-1, 2)
80 ) {
81 throw AipsError(_ORIGIN + "You cannot fix the width of a non-reference "
82 "Gaussian if the reference Gaussian's width is not fixed and there "
83 "is a relationship between the two widths."
84 );
85 }
86 }
87 ostringstream myfunc;
88 myfunc << "p0*exp(-0.5 * (x-p1)*(x-p1) / p2/p2)";
89 Vector<Double> parm(3 + nfalse(fixed), 0);
90 Vector<Double> errs = parm.copy();
91 parm[0] = _gaussians[0].getAmpl();
92 parm[1] = _gaussians[0].getCenter();
93 parm[2] = _gaussians[0].getSigma();
94 errs[0] = _gaussians[0].getAmplErr();
95 errs[1] = _gaussians[0].getCenterErr();
96 errs[2] = _gaussians[0].getSigmaErr();
97 Vector<Bool> f(parm.size(), True);
98 f[0] = _gaussians[0].fixedAmpl();
99 f[1] = _gaussians[0].fixedCenter();
100 f[2] = _gaussians[0].fixedSigma();
101 _paramIndices(0, 0) = 0;
102 _paramIndices(0, 1) = 1;
103 _paramIndices(0, 2) = 2;
104 uInt p = 3;
105 for (uInt i=0; i<constraints.nrow(); i++) {
106 String amp;
107 if (constraints(i, 0) != 0) {
108 amp = String::toString(constraints(i, 0)) + "*p0";
109 }
110 else {
111 amp = "p" + String::toString(p);
112 parm[p] = _gaussians[i+1].getAmpl();
113 errs[p] = _gaussians[i+1].getAmplErr();
114 f[p] = _gaussians[i+1].fixedAmpl();
115 _paramIndices(i+1, 0) = p;
116 p++;
117 }
118 String center;
119 if (constraints(i, 1) != 0) {
120 center = "x-(p1+(" + String::toString(constraints(i, 1)) + "))";
121 }
122 else {
123 center = "x-p" + String::toString(p);
124 parm[p] = _gaussians[i+1].getCenter();
125 errs[p] = _gaussians[i+1].getCenterErr();
126 f[p] = _gaussians[i+1].fixedCenter();
127 _paramIndices(i+1, 1) = p;
128 p++;
129
130 }
131 String sigma;
132 if (constraints(i, 2) != 0) {
133 sigma = String::toString(constraints(i, 2)) + "*p2";
134 }
135 else {
136 sigma = "p" + String::toString(p);
137 parm[p] = _gaussians[i+1].getSigma();
138 errs[p] = _gaussians[i+1].getSigmaErr();
139 f[p] = _gaussians[i+1].fixedSigma();
140 _paramIndices(i+1, 2) = p;
141 p++;
142 }
143 myfunc << " + " << amp << "*exp(-0.5 * (" << center << ")*("
144 << center << ") / (" << sigma << ") / (" << sigma << "))";
145 }
146 _setFunction(myfunc.str());
147 // have to set the GaussianSpectralElement parameters
148 set(parm);
149 setError(errs);
150 fix(f);
151}
152
153GaussianMultipletSpectralElement::GaussianMultipletSpectralElement(
154 const GaussianMultipletSpectralElement& other
155) : CompiledSpectralElement(other), _gaussians(other._gaussians),
156 _constraints(other._constraints),
157 _paramIndices(other._paramIndices) {}
158
159
160GaussianMultipletSpectralElement::~GaussianMultipletSpectralElement() {}
161
162SpectralElement* GaussianMultipletSpectralElement::clone() const {
163 return new GaussianMultipletSpectralElement(*this);
164}
165
166GaussianMultipletSpectralElement& GaussianMultipletSpectralElement::operator=(
167 const GaussianMultipletSpectralElement &other
168) {
169 if (this != &other) {
170 CompiledSpectralElement::operator=(other);
171 _gaussians = other._gaussians;
172 _constraints.resize(other._constraints.shape());
173 _constraints = other._constraints.copy();
174 _paramIndices.resize(other._paramIndices.shape());
175 _paramIndices = other._paramIndices.copy();
176 }
177 return *this;
178}
179
180Bool GaussianMultipletSpectralElement::operator==(
181 const GaussianMultipletSpectralElement& other
182) const {
183 return(
184 CompiledSpectralElement::operator==(other)
185 && allTrue(Vector<GaussianSpectralElement>(_gaussians) == Vector<GaussianSpectralElement>(other._gaussians))
186 && allTrue(_constraints == other._constraints)
187 );
188}
189
190const vector<GaussianSpectralElement>&
191GaussianMultipletSpectralElement::getGaussians() const {
192 return _gaussians;
193}
194
195const Matrix<Double>& GaussianMultipletSpectralElement::getConstraints() const {
196 return _constraints;
197}
198
199Bool GaussianMultipletSpectralElement::toRecord(RecordInterface& out) const {
200 out.define(RecordFieldId("type"), fromType(getType()));
201 Record gaussians;
202 for (uInt i=0; i<_gaussians.size(); i++) {
203 Record gaussian;
204 _gaussians[i].toRecord(gaussian);
205 gaussians.defineRecord(
206 "*" + String::toString(i), gaussian
207 );
208 }
209 out.defineRecord("gaussians", gaussians);
210 out.define("fixedMatrix", _constraints);
211 return True;
212}
213
214void GaussianMultipletSpectralElement::set(const Vector<Double>& param) {
215 if (get().size() > 0 && param.size() != get().size()) {
216 ostringstream x;
217 x << _ORIGIN << "Inconsistent number of parameters. Got "
218 << param.size() << ". Must be " << get().size();
219 throw AipsError(x.str());
220 }
221 SpectralElement::_set(param);
222 Double amp0 = param[0];
223 Double center0 = param[1];
224 Double sigma0 = param[2];
225 _gaussians[0].setAmpl(amp0);
226 _gaussians[0].setCenter(center0);
227 _gaussians[0].setSigma(sigma0);
228 uInt p = 3;
229 for (uInt i=3; i<_paramIndices.size(); i++) {
230 uInt gNumber = i/3;
231 uInt pNumber = i%3;
232 uInt pIndex = _paramIndices(gNumber, pNumber);
233 Double fRel = _constraints(gNumber-1, pNumber);
234 if (pNumber == 0) {
235 Double amp = pIndex == 0 ? fRel*amp0 : param[p];
236 _gaussians[gNumber].setAmpl(amp);
237 }
238 else if (pNumber == 1) {
239 Double center = pIndex == 0 ? fRel+center0 : param[p];
240 _gaussians[gNumber].setCenter(center);
241 }
242 else if (pNumber == 2) {
243 Double sigma = pIndex == 0 ? fRel*sigma0 : param[p];
244 _gaussians[gNumber].setSigma(sigma);
245 }
246 if (pIndex > 0) {
247 p++;
248 }
249 }
250}
251
252void GaussianMultipletSpectralElement::setError(const Vector<Double> &err) {
253 SpectralElement::setError(err);
254 Double amp0 = err[0];
255 Double center0 = err[1];
256 Double sigma0 = err[2];
257 Vector<Double> errors(3);
258 errors[0] = err[0];
259 errors[1] = err[1];
260 errors[2] = err[2];
261 _gaussians[0].setError(errors);
262 uInt p = 3;
263 for (uInt i=3; i<_paramIndices.size(); i++) {
264 uInt gNumber = i/3;
265 uInt pNumber = i%3;
266 uInt pIndex = _paramIndices(gNumber, pNumber);
267 Double fRel = _constraints(gNumber-1, pNumber);
268 if (pNumber == 0) {
269 errors[0] = pIndex == 0 ? fRel*amp0 : err[p];
270 }
271 else if (pNumber == 1) {
272 errors[1] = pIndex == 0 ? center0 : err[p];
273 }
274 else if (pNumber == 2) {
275 errors[2] = pIndex == 0 ? fRel*sigma0 : err[p];
276 }
277 _gaussians[gNumber].setError(errors);
278
279 if (pIndex > 0) {
280 p++;
281 }
282 }
283}
284
285void GaussianMultipletSpectralElement::fix(const Vector<Bool>& fix) {
286 SpectralElement::fix(fix);
287 Bool amp0 = fix[0];
288 Bool center0 = fix[1];
289 Bool sigma0 = fix[2];
290 Vector<Bool> fixed(3);
291 fixed[0] = fix[0];
292 fixed[1] = fix[1];
293 fixed[2] = fix[2];
294 _gaussians[0].fix(fixed);
295 uInt p = 3;
296 for (uInt i=3; i<_paramIndices.size(); i++) {
297 uInt gNumber = i/3;
298 uInt pNumber = i%3;
299 uInt pIndex = _paramIndices(gNumber, pNumber);
300 if (pNumber == 0) {
301 fixed[0] = pIndex == 0 ? amp0 : fix[p];
302 }
303 else if (pNumber == 1) {
304 fixed[1] = pIndex == 0 ? center0 : fix[p];
305 }
306 else if (pNumber == 2) {
307 fixed[2] = pIndex == 0 ? sigma0 : fix[p];
308 }
309 _gaussians[gNumber].fix(fixed);
310 if (pIndex > 0) {
311 p++;
312 }
313 }
314}
315
316ostream &operator<<(ostream& os, const GaussianMultipletSpectralElement& elem) {
317 os << SpectralElement::fromType((elem.getType())) << " element: " << endl;
318 os << " Function: " << elem.getFunction() << endl;
319 os << " Gaussians:" << endl;
320 Vector<GaussianSpectralElement> gaussians = elem.getGaussians();
321 for (uInt i=0; i<gaussians.size(); i++) {
322 os << "Gaussian " << i << ": " << gaussians[i] << endl;
323 }
324 Matrix<Double> r = elem.getConstraints();
325 os << "Constraints: " << endl;
326 for (uInt i=1; i<gaussians.size(); i++) {
327 for (uInt j=0; j<3; j++) {
328 Double v = r(i-1, j);
329 if (v != 0) {
330 switch (j) {
331 case 0:
332 os << " Amplitude ratio of component " << i
333 << " to the reference component is fixed at " << v;
334 break;
335 case 1:
336 os << " Center offset of component " << i
337 << " to the reference component is fixed at " << v;
338 break;
339 case 2:
340 os << " FWHM ratio of component " << i
341 << " to the reference component is fixed at " << v;
342 break;
343 default:
344 throw AipsError(_ORIGIN + "Logic Error");
345 }
346 }
347 }
348 }
349 return os;
350}
351
352} //# NAMESPACE CASA - END
353
354
Note: See TracBrowser for help on using the repository browser.