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 : //# $Id: SpectralElement.cc 21024 2011-03-01 11:46:18Z gervandiepen $
27 :
28 : #include <components/SpectralComponents/SpectralElement.h>
29 :
30 : #include <casacore/casa/BasicSL/Constants.h>
31 : #include <casacore/casa/BasicSL/String.h>
32 : #include <casacore/casa/Arrays/ArrayLogical.h>
33 : #include <casacore/casa/Exceptions/Error.h>
34 : #include <casacore/casa/Utilities/MUString.h>
35 : #include <casacore/scimath/Mathematics/AutoDiffMath.h>
36 : #include <casacore/scimath/Functionals/Function.h>
37 :
38 : //debug only
39 : #include <casacore/scimath/Functionals/CompiledFunction.h>
40 : #include <casacore/casa/IO/ArrayIO.h>
41 :
42 : #include <iostream>
43 :
44 : using namespace casacore;
45 : namespace casa { //# NAMESPACE CASA - BEGIN
46 :
47 16262 : SpectralElement::SpectralElement(SpectralElement::Types type, const Vector<Double>& parms)
48 16262 : : _type(type), _params(parms), _errors(parms.size(), 0),
49 32524 : _fixed(parms.size(), false) {}
50 :
51 :
52 2574283 : SpectralElement::SpectralElement(const SpectralElement &other)
53 2574283 : : _type(other._type), _params(other._params.copy()), _errors(other._errors.copy()),
54 2574283 : _fixed(other._fixed.copy()),
55 2574283 : _function(
56 : std::shared_ptr<Function<Double, Double> >(
57 2574283 : other._function->clone()
58 : )
59 5148566 : ) {}
60 :
61 2590545 : SpectralElement::~SpectralElement() {}
62 :
63 6 : SpectralElement &SpectralElement::operator=(
64 : const SpectralElement &other
65 : ) {
66 6 : if (this != &other) {
67 6 : _type = other._type;
68 6 : uInt n = other._params.size();
69 6 : _params.resize(n);
70 6 : _params = other._params.copy();
71 6 : _errors.resize(n);
72 6 : _errors = other._errors.copy();
73 6 : _fixed.resize(n);
74 6 : _fixed = other._fixed.copy();
75 12 : _function = std::shared_ptr<Function<Double, Double> >(
76 6 : other._function->clone()
77 6 : );
78 : }
79 6 : return *this;
80 : }
81 :
82 0 : Double SpectralElement::operator[](const uInt n) const {
83 0 : if (n >= _params.size()) {
84 0 : throw(AipsError("SpectralElement: Illegal index for parameter"));
85 : }
86 0 : return _params[n];
87 : }
88 :
89 0 : Bool SpectralElement::operator==(
90 : const SpectralElement& other
91 : ) const {
92 0 : if (this == &other) {
93 0 : return true;
94 : }
95 : return (
96 0 : _type == other._type && allNear(_params, other._params, 1e-8)
97 0 : && allNear(_errors, other._errors, 1e-8)
98 0 : && allTrue(_fixed == other._fixed)
99 0 : );
100 : }
101 :
102 47063 : const String* SpectralElement::allTypes(
103 : Int &nall, const SpectralElement::Types *&typ
104 : ) {
105 : static const String tname[SpectralElement::N_Types] = {
106 : String("GAUSSIAN"),
107 : String("POLYNOMIAL"),
108 : String("COMPILED"),
109 : String("GAUSSIAN MULTIPLET"),
110 : String("LORENTZIAN"),
111 : String("POWER LOGARITHMIC POLYNOMIAL"),
112 : String("LOGARITHMIC TRANSFORMED POLYNOMIAL")
113 :
114 47063 : };
115 :
116 : static const SpectralElement::Types oname[SpectralElement::N_Types] = {
117 : SpectralElement::GAUSSIAN,
118 : SpectralElement::POLYNOMIAL,
119 : SpectralElement::COMPILED,
120 : SpectralElement::GMULTIPLET,
121 : SpectralElement::LORENTZIAN,
122 : SpectralElement::POWERLOGPOLY,
123 : SpectralElement::LOGTRANSPOLY
124 :
125 : };
126 :
127 47063 : nall = SpectralElement::N_Types;
128 47063 : typ = oname;
129 47063 : return tname;
130 : }
131 :
132 807428 : void SpectralElement::_set(const Vector<Double>& params) {
133 807428 : _params = params.copy();
134 2434115 : for (uInt i=0; i<params.size(); i++) {
135 1626687 : (*_function)[i] = params[i];
136 : }
137 807428 : }
138 :
139 5 : void SpectralElement::_setType(const SpectralElement::Types type) {
140 5 : _type = type;
141 5 : }
142 :
143 47063 : const String &SpectralElement::fromType(SpectralElement::Types tp) {
144 : Int nall;
145 : const SpectralElement::Types *typ;
146 47063 : const String *const tname = SpectralElement::allTypes(nall, typ);
147 47063 : return tname[tp];
148 : }
149 :
150 0 : Bool SpectralElement::toType(
151 : SpectralElement::Types &tp, const String &typname
152 : ) {
153 : Int nall;
154 : const SpectralElement::Types *typ;
155 0 : const String *const tname = SpectralElement::allTypes(nall, typ);
156 :
157 : // Make sure a value returned
158 0 : tp = typ[0];
159 0 : Int i = MUString::minimaxNC(typname, SpectralElement::N_Types, tname);
160 0 : if (i >= nall) {
161 0 : return false;
162 : }
163 0 : tp = typ[i];
164 0 : return true;
165 : }
166 :
167 32498 : void SpectralElement::_setFunction(
168 : const std::shared_ptr<Function<Double, Double> >& f
169 : ) {
170 32498 : _function = f;
171 32498 : }
172 :
173 42315760 : Double SpectralElement::operator()(const Double x) const {
174 42315760 : return (*_function)(x);
175 : }
176 :
177 31921 : void SpectralElement::get(Vector<Double> ¶m) const {
178 31921 : param = _params.copy();
179 31921 : }
180 :
181 2920370 : Vector<Double> SpectralElement::get() const {
182 5840740 : return _params.copy();
183 : }
184 :
185 31921 : void SpectralElement::getError(Vector<Double> &err) const {
186 31921 : err = _errors.copy();
187 31921 : }
188 :
189 1002833 : Vector<Double> SpectralElement::getError() const {
190 2005666 : return _errors.copy();
191 : }
192 :
193 807606 : void SpectralElement::setError(const Vector<Double> &err) {
194 807606 : if (err.nelements() != _params.nelements()) {
195 0 : throw(
196 : AipsError(
197 : "SpectralElement: setting incorrect "
198 : "number of errors in the element"
199 : )
200 0 : );
201 : }
202 807606 : _errors = err.copy();
203 807606 : }
204 :
205 46 : void SpectralElement::fix(const Vector<Bool> &fix) {
206 46 : if (fix.nelements() != _params.nelements()) {
207 0 : throw(
208 : AipsError(
209 : "SpectralElement: setting incorrect number of fixed "
210 : "in the element"
211 : )
212 0 : );
213 : }
214 46 : _fixed = fix.copy();
215 167 : for (uInt i=0; i<_fixed.size(); i++) {
216 121 : _function->mask(i) = fix[i];
217 : }
218 46 : }
219 :
220 583057 : const Vector<Bool>& SpectralElement::fixed() const {
221 583057 : return _fixed;
222 : }
223 :
224 0 : Bool SpectralElement::toRecord(RecordInterface &out) const {
225 0 : out.define(RecordFieldId("type"), fromType(_type));
226 0 : Vector<Double> ptmp(_params);
227 0 : Vector<Double> etmp(_params);
228 0 : out.define(RecordFieldId("parameters"), _params);
229 0 : out.define(RecordFieldId("errors"), _errors);
230 0 : out.define(RecordFieldId("fixed"), _fixed);
231 0 : return true;
232 0 : }
233 :
234 529425 : void SpectralElement::set(const Vector<Double>& params) {
235 529425 : if (params.size() != get().size()) {
236 0 : throw AipsError(
237 : "Input number of parameters does not match "
238 : "the current number of parameters"
239 0 : );
240 : }
241 529425 : _set(params);
242 529425 : }
243 :
244 : } //# NAMESPACE CASA - END
245 :
246 :
|
