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 | |
---|
37 | namespace casa { //# NAMESPACE CASA - BEGIN |
---|
38 | |
---|
39 | #define _ORIGIN String("GaussianMultipletSpectralElement::") + __FUNCTION__ + ":" + String::toString(__LINE__) + ": " |
---|
40 | |
---|
41 | GaussianMultipletSpectralElement::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 | |
---|
153 | GaussianMultipletSpectralElement::GaussianMultipletSpectralElement( |
---|
154 | const GaussianMultipletSpectralElement& other |
---|
155 | ) : CompiledSpectralElement(other), _gaussians(other._gaussians), |
---|
156 | _constraints(other._constraints), |
---|
157 | _paramIndices(other._paramIndices) {} |
---|
158 | |
---|
159 | |
---|
160 | GaussianMultipletSpectralElement::~GaussianMultipletSpectralElement() {} |
---|
161 | |
---|
162 | SpectralElement* GaussianMultipletSpectralElement::clone() const { |
---|
163 | return new GaussianMultipletSpectralElement(*this); |
---|
164 | } |
---|
165 | |
---|
166 | GaussianMultipletSpectralElement& 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 | |
---|
180 | Bool 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 | |
---|
190 | const vector<GaussianSpectralElement>& |
---|
191 | GaussianMultipletSpectralElement::getGaussians() const { |
---|
192 | return _gaussians; |
---|
193 | } |
---|
194 | |
---|
195 | const Matrix<Double>& GaussianMultipletSpectralElement::getConstraints() const { |
---|
196 | return _constraints; |
---|
197 | } |
---|
198 | |
---|
199 | Bool 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 | |
---|
214 | void 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 | |
---|
252 | void 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 | |
---|
285 | void 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 | |
---|
316 | ostream &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 | |
---|