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