[2980] | 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 | |
---|